Mycobiome detection from a single subterranean gametophyte using metabarcoding techniques.

Premise: Several ferns and lycophytes produce subterranean gametophytes, including the Ophioglossaceae, Psilotaceae, and some members of the Schizaeaceae, Gleicheniaceae, and Lycopodiaceae. Despite the surge in plant-microbiome research, which has been particularly boosted by high-throughput sequencing techniques, the microbiomes of these inconspicuous fern gametophytes have rarely been examined. The subterranean gametophytes are peculiar due to their achlorophyllous nature, which makes them rely on fungi to obtain nutrients. Furthermore, the factors that shape the fungal communities (mycobiomes) of fern gametophytes have not been examined in depth. Methods and
Results: We present a workflow to study the mycobiome of the achlorophyllous gametophytes of Ophioderma pendulum using a high-throughput metabarcoding approach. Simultaneously, each gametophyte was investigated microscopically to detect fungal structures. Two primer sets of the nuclear ITS sequence targeting general fungi were applied, in addition to a primer set that specifically targets the nuclear small subunit ribosomal rDNA region of arbuscular mycorrhizal fungi. Both the microscopic and metabarcoding approaches revealed many diverse fungi inhabiting a single gametophyte of O. pendulum. Discussion: This study provides methodological details and discusses precautions for the mycobiome investigation of achlorophyllous gametophytes. This research is also the first to uncover the mycobiome assembly of an achlorophyllous gametophyte of an epiphytic fern.

The plant microbiome is now considered an intrinsic component of all plant lineages and plays critical roles in maintaining plant health (Vandenkoornhuyse et al., 2015; Hassani et al., 2018). Among the diverse microorganisms associated with plant roots, fungi are vital members (Anal et al., 2020; Pozo et al., 2021), with over 80% of terrestrial plants forming mycorrhizal associations (Wang and Qiu, 2006). Of these, 68% of ferns and lycophytes are estimated to be associated with mycorrhizae or root endophytes (Lehnert et al., 2017). Arbuscular mycorrhizal fungi (AMF) belonging to the Glomeromycota are the most prevalent fungi associated with the roots of ferns and lycophytes (reviewed by Lehnert et al., 2017). Recently, fine‐root endophytes (FREs) (Mucoromycota) have been shown to establish mutualistic relationships with lycophytes (Hoysted et al., 2019), while dark septate endophytes (DSEs), which are often Basidiomycota or Ascomycota fungi, have also been reported in 13% of ferns and lycophytes (Lehnert et al., 2017).

The gametophyte represents a critical life stage for ferns and lycophytes, and is the key stage for their establishment (Haufler et al., 2016); however, studies of plant–fungal associations have largely focused on the sporophytic stage. Of all types of gametophytes, subterranean gametophytes are especially peculiar due to their achlorophyllous nature. Because these subterranean gametophytes are incapable of photosynthesis, they utilize mycoheterotrophic relationships, relying on fungi for carbon translocation (Pressel et al., 2016). Their subterranean nature, however, makes them especially difficult to find in the wild, resulting in a lack of studies based on naturally occurring populations. Furthermore, subterranean gametophytes are difficult to grow axenically in culture, constraining the functional investigation of plant–fungal associations in experimental systems (Whittier, 1981). Because of the limited morphological structure of fungi associated with plant roots, fungal identification based on the barcoding region of ribosomal DNA (rDNA) is essential (Schoch et al., 2012). While morphological characterization has demonstrated fungal occurrence with the subterranean gametophytes, molecular DNA analysis has rarely been applied (Lang, 1902; Duckett and Ligrone, 2005). By Sanger sequencing the fungal nuclear ribosomal small subunit (SSU) region, Winther and Friedman (2007) detected Glomus AMF in the subterranean gametophyte of the lycophyte genus Huperzia Bernh. In addition to the SSU region, Horn et al. (2013) targeted the nuclear ribosomal large subunit and the nuclear ITS sequences of Diphasiastrum alpinum (L.) Holub gametophytes, which in combination revealed the presence of AMF, Sebacinales (Basidiomycota) fungi, and common ericoid mycorrhizal fungi of the Leotiomycetes (Ascomycota). A single gametophyte can be cut into two pieces for study, one for microscopic examination and the other for molecular identification, allowing the direct comparison between anatomical structure and molecular identification. This technique has been employed in green gametophytes using Sanger sequencing (Ogura‐Tsujita et al., 2013, 2016).

Swatzell et al. (1996) highlighted the positive effect of an aseptate fungus on the gametophyte (Schizaea pusilla Pursh) in a co‐culture experiment. In addition to typical mycorrhizal fungi, many endophytic fungi that do not produce specialized structures on the infected roots, such as DSEs, could play a critical role in improving the stress tolerance of host plants (Li et al., 2019). While DSEs have been reported from multiple case studies on fern sporophytes (Fernández et al., 2013; Muthukumar et al., 2014), their presence on subterranean gametophytes has not been demonstrated. In addition, the ferns and lycophytes include a higher percentage of epiphytes than other groups of vascular plants (Zotz, 2013). The scarcity of soil and rapidly changing moisture levels pose additional challenges to the epiphytic way of life, and can potentially influence the fungal members associated with these plants (Zotz, 2013). Compared with terrestrial ferns, Lehnert et al. (2017) hypothesized that a higher proportion of epiphytic ferns would be associated with DSEs, although this remains untested.

Recently, the advancement of high‐throughput sequencing, particularly metabarcoding (or amplicon sequencing), has enabled more comprehensive explorations of plant‐associated fungal communities (U'Ren et al., 2019). Such an approach allows the sequencing of millions of copies of the desired DNA region from targeted organisms. While the metabarcoding technique has boosted a surge of research into the fungal communities of other tracheophytes and several ferns (Benucci et al., 2020; Suetsugu et al., 2020; Oh et al., 2021), to our knowledge, it has not yet been applied to subterranean gametophytes. Two possible reasons for the lack of research on this topic are: (1) small subterranean gametophytes are often difficult to find or identify because they are embedded in the substrate (e.g., soil), and (2) there is not a sufficient amount of plant tissue for both standard DNA extraction and microscopic examination.

Missing information hinders our understanding of the abiotic and biotic interactions of fern–fungal relationships and co‐evolution. Here, we used the humus‐embedded gametophyte of Ophioderma pendulum C. Presl, an epiphyte in the Ophioglossaceae (Figure 1), as an example to showcase a workflow for simultaneous microscopic examination and fungal metabarcoding in fern gametophytes (Figure 2). The flexible barcoding of each sample provides an economic approach for studying multiple samples simultaneously. We aimed to investigate the following questions for the subterranean gametophyte mycobiome: (1) Is a single subterranean gametophyte adequate for simultaneous microscopic and metabarcoding investigations? (2) Are the common fungal primers targeting general fungi and AMF suitable for such material? To answer (1), we used one gametophyte to perform the mycobiome investigation. For (2), we tested three nuclear rDNA primer sets that targeted different fungal groups (AMF and other root‐associated fungi, including possible DSE and ericoid mycorrhizae) of interest.

Figure 1

Images of the study organism, the fern Ophioderma pendulum. (A) Ophioderma pendulum growing on a tree trunk, as indicated by a white arrow. (B) The achlorophyllous subterranean gametophyte of O. pendulum

Figure 2

Workflow of a mycobiome examination in subterranean gametophytes using a combination of DNA metabarcoding and microscopy

METHODS

An overall workflow is presented in Figure 2. The material was collected in October 2020 in a cloud forest in northern Taiwan (altitude: 500 m a.s.l., GPS coordinates: 24.84°N, 121.55°E).

Step 1: Harvest and storage of subterranean gametophytes

To search for subterranean gametophytes of O. pendulum, we targeted a few fallen tree ferns with many young O. pendulum sporophytes attached (<15 cm tall). The loose substrate (i.e., tree fern root mantle) was brought back to the lab. Using a previously published morphological description (Lang, 1902; Whittier and Braggins, 2011) and a hand lens, we were able to identify a potential Ophioderma gametophyte (<5 mm in diameter) (Figure 1B). The gametophyte was cleaned with a watercolor brush in running tap water to remove the surface debris. Following the cleaning step, the gametophyte was patted dry and cut in half with a single‐edged blade. The cut should result in two halves as symmetric and homogeneous as possible to avoid fungal detection bias. The half for microscopic examination and the other half for fungal metabarcoding were stored in 70% ethanol and in a −20°C freezer, respectively.

Step 2: Microscopic examination

We cut one gametophyte half into ca. 3 × 3‐mm2 pieces with a microtome knife and placed them into a 2‐mL Eppendorf microcentrifuge tube. To remove impurities and pigments, 1 mL 5% KOH was added to the tube, which was heated at 95°C in a dry bath for 50 min. We then removed the KOH and rinsed the sample twice with ddH2O. To balance the pH, 1 mL 5% HCl was added to the tube, and the sample was soaked at room temperature for 40 min. After removing the HCl, the sample was directly dipped into 1 mL 0.05% trypan blue (dissolved in lactoglycerol [ddH2O : lactic acid : glycerol in a 1 : 1 : 1 ratio by volume]) in a 2‐mL Eppendorf microcentrifuge tube at room temperature for 20 min. Finally, we used lactoglycerol to remove the excess dye; the first 0.5 mL of lactoglycerol was used to quickly wash trypan blue from the surface of the sample and the tube, with an additional 1 mL of lactoglycerol subsequently used to soak the sample for 1 min, which gradually removed the excess dye from inside the plant tissue. The stained gametophytes were examined under a Leica DM500 compound microscope (Leica Microsystems, Wetzlar, Germany) to visually detect fungal structures.

Step 3: DNA extraction

We then divided the gametophyte half used for DNA extraction into three pieces, with each piece representing a replicate. Due to the small size of the gametophyte (Figure 1) and the generally high concentration of phenolic compounds in underground plant tissues, a cetyltrimethylammonium bromide (CTAB) extraction method followed by a bead clean‐up (AMPure XP; Beckman Coulter, Brea, California, USA) step was applied to ensure high DNA yield and purity. The overall DNA extraction procedure was performed as described previously (Liao et al., 2014), but the volumes of the reagents were adjusted for the small size of the gametophytes. Each gametophyte sample was placed in a 2‐mL Eppendorf microcentrifuge tube containing three silica/zirconia beads (2.3 mm diameter; BioSpec Products, Bartlesville, Oklahoma, USA). The tissue was ground for 30 s at 70 Hz two times using a tissue grinder (STEP tissue grinder; ACTTR, New Taipei City, Taiwan), after which 200 μL CTAB buffer was added to each sample before they were heated in a dry bath at 65°C for 10–15 min. The tube was inverted at least three times during the incubation. Following the incubation, 200 μL of a chloroform : isoamyl alcohol (24 : 1) solution was added to the tube. The tube was vortexed immediately to allow sufficient mixing and centrifuged at 9000 × g for 10 min. After centrifugation, the supernatant was carefully transferred into a new 2‐mL Eppendorf microcentrifuge tube, and an equal amount of isopropanol was added. The tube was gently inverted a few times to mix and was then placed in a −20°C freezer for 5 min before being centrifuged at 13,000 × g for 10 min. The supernatant was carefully removed and discarded, and 800 μL of 80% ethanol was added to the tube, which was centrifuged at 13,000 × g for another 10 min before the supernatant was removed. The tube was inverted and placed on a clean paper towel at room temperature for 20 min to allow the ethanol to evaporate, after which 20 μL of nuclease‐free water was added to resuspend the pellet. The reagent recipe for the DNA extraction is provided in Appendix S1. To remove any contamination that might inhibit PCR efficiency, the extracted DNA sample was further purified with a paramagnetic bead–based system (AMPure XP). A 1 : 1 (v/v) ratio of DNA extraction : AMPure beads was used. The clean‐up procedure was performed following Chen et al. (2021).

Step 4: Host identification

Due to the limited morphological characters of fern gametophytes, confirmation of the gametophyte identity to the species level requires molecular sequencing. We amplified the trnL‐L‐F sequence (including the trnL gene and the trnL‐F intergenic spacer) of the plastid genome using primers C and F described by Taberlet et al. (1991), and confirmed that the sequence of the gametophyte was identical to what we generated from the conspecific sporophyte. The sequences were submitted to the National Center for Biotechnology Information's (NCBI) GenBank database with the accession numbers MZ723323 (gametophyte) and MZ723322 (sporophyte).

Step 5: Metabarcoding library preparation

A two‐step PCR protocol enabling flexible primer choice was applied, followed by Illumina MiSeq 300 PE sequencing (Illumina, San Diego, California, USA) (Chen et al., 2021). To begin, the targeted genomic region was amplified with ordinary primers linked to a frame‐shift region and Illumina adapters (first PCR) (Table 1, Appendices [Link], [Link], [Link], [Link]). The first PCR product was then cleaned using AMPure beads (AMPure XP). A 1 : 1 sample : AMPure beads (v/v) ratio was used. The clean first PCR product was used as the input for a second PCR, which amplified each PCR product and attached a 10‐bp barcode to each amplicon (Appendices [Link], [Link], [Link]) (Chen et al., 2021). The samples were divided into equal molar pools and were sequenced by the NGS High‐Throughput Genomics Core at Academia Sinica, Taipei, Taiwan. Because previous reports suggested the likely importance of AMF as well as general fungi of the Dikarya for the ferns and lycophytes (Horn et al., 2013; Lehnert et al., 2017; Benucci et al., 2020), three primer sets were tested here to allow the detection of a wide range of fungal diversity. We tested two general fungal primer pairs: ITS1F–ITS2, targeting the ITS1 region, and ITS3ngs1‐5–ITS4, targeting the ITS2 region (Figure 2). The other primer set, NS31–AML2, specifically targeted the SSU rDNA of AMF (Table 1). The data sets generated using ITS1F–ITS2, ITS3ngs1_5–ITS4, and NS31–AML2 will be referred to as “ITS1 data set,” “ITS2 data set,” and “SSU data set,” respectively. To amplify the ITS2 region of broader phylogenetic groups of fungi, we incorporated the primer‐cocktail approach designed by Tedersoo et al. (2015) with the two‐step PCR metabarcoding approach described by Chen et al. (2021) (Table 1, Appendix S2). Instead of using one forward primer in the first PCR reaction, five forward primers complementing each other were mixed to maximize the fungal taxa that could be targeted. The complete primer sequences designed for the two‐step PCR metabarcoding are provided in Appendix S2. All the raw reads were deposited in the Sequence Read Archive of NCBI with the BioProject accession PRJNA752349.

Table 1

Primer information used in this study

Genomic region	Direction	Primer name	Targeted fungal group	Reference	Average amplicon size (including adapter and barcode)	No. of reads received after initial demultiplexing	Read quality (Phred score) (including forward and reverse reads)	Taxonomy assignment database	Taxonomy assignment tool used with resulting sequences
ITS1	Forward	ITS1F	General fungi	White et al., 1990	447	160,925 ± 31,344	32.73 ± 0.61	UNITE	DADA2
	Reverse	ITS2	General fungi	White et al., 1990
ITS2	Forward	ITS3ngs1_5	General fungi	Tedersoo et al., 2015	512	119,189 ± 12,698	32.73 ± 0.16	UNITE	DADA2
	Reverse	ITS4	General fungi	White et al., 1990
SSU	Forward	NS31	General fungi	Kohout et al., 2014	723	73,242 ± 6,981	30.50 ± 0.73	NCBI GenBank	MEGAN
	Reverse	AML2	Arbuscular mycorrhizal fungi	Kohout et al., 2014

Step 6: Bioinformatics

Key steps of the mycobiome analyses are illustrated in Figure 2. Demultiplexing was performed by the sequencing provider using Illumina bcl2fastq Conversion software. Cutadapt (Martin, 2011) was used to remove the primers and adapters. The scripts and specific sequences required for each data set are provided in Appendix S8. For the ITS1 and ITS2 data sets, the forward and reverse reads were merged before the downstream analysis. For the SSU data set, due to the long amplicon size (Table 1), only the forward reads were subjected to further analysis. DADA2 (Callahan et al., 2016) was used for quality trimming, chimera removal, and determining the amplicon sequence variant (ASV). The same filtering criteria were applied to all three data sets and for both forward and reverse reads. After truncation (truncQ = 2), reads with the following features were discarded using the “filterAndTrim” command: (1) include Ns, (2) expected error > 2, and (3) length < 50 bp.

For samples targeting the ITS1 and ITS2 regions, the taxonomic assignment was performed using the RDP Naïve Bayesian Classifier with the UNITE fungal training set (accessed May 2021) in DADA2 (Wang et al., 2007; Callahan et al., 2016; Nilsson et al., 2019b). For samples targeting the SSU region of AMF, we first compared ASV sequences against NCBI GenBank (accessed June 2021) using BLASTN. The program optimized for searching against highly similar sequences (MegaBLAST) (Morgulis et al., 2008) was selected. Up to 100 hits with an expected value (E‐value) < 0.01 were kept. The outcomes were imported into MEGAN for the taxonomy assignment with the default lowest common ancestor algorithm (Morgulis et al., 2008; Huson et al., 2016). Alternatively, a curated AMF database is available for AMF identification (Öpik et al., 2010). The ASV data and taxa were then imported into the Phyloseq R package (McMurdie and Holmes, 2013) for downstream analysis. ASVs not belonging to the fungal kingdom were removed. The composition of the fungal community was examined at the “phylum” and “class” rank.

RESULTS

We investigated the fungal community of O. pendulum gametophytes using microscopy and NGS DNA metabarcoding to reveal the fungal structure and identity.

To confirm the gametophyte's identity, we compared its trnL‐L‐F sequence with that of a nearby O. pendulum sporophyte. The trnL‐L‐F sequences generated from the gametophyte and sporophyte shared 100% similarity, suggesting that the gametophyte was indeed O. pendulum. Trypan blue–stained structures were observed across the gametophyte, suggesting that a critical proportion of the gametophyte was colonized by fungi (Figure 3A). These fungal structures were concentrated in the center of the gametophyte and not in the protruding structures (Figure 3A). Because the gametophyte was not symmetrical, the spatial distribution of the fungi might be slightly different in the two halves of the bisected gametophyte; thus, the fungal communities detected via molecular sequencing might not completely reflect the microscopic structures in the other half. Microscopic structures were observed, including arbuscules, vesicles, and hyphal coils (Figure 3). The presence of arbuscules (Figure 3B) suggested an AMF association with the gametophyte of O. pendulum. The septate hyphae of the hyphal coils (Figure 3D) suggested that these structures were produced by non‐AMF fungi.

Figure 3

Fungal structures associated with the Ophioderma pendulum gametophyte. (A) The branches of the gametophyte were fungus‐free, but the central part was densely colonized by fungal hyphae. (B) Spherical arbuscules in the cells. (C) A vesicle produced from the middle of a hypha. (D) Hyphal coils inside the plant cells. White arrow, septa; ab, arbuscules; hc, hyphal coils; v, vesicle

All three primer sets were successfully used to amplify fungal sequences. The read quality of the SSU data set was lower than that of the ITS1 and ITS2 data sets (Table 1). The relative abundance of phyla and classes detected using the three primer sets differed (Figure 4). Both the ITS1 and ITS2 data sets indicated that >97% of reads belonged to the phylum Ascomycota. As expected, the SSU data set, which was generated using an AMF‐specific primer set, resulted in >93% of reads coming from Glomeromycota species. At the class level, both the ITS1 and ITS2 data sets suggested that Leotiomycetes dominated the gametophytic mycobiomes, while the ITS2 data set also revealed abundant Eurotiomycetes (10–33% of reads), suggesting the potential importance of eurotiomycetous fungi in the gametophytic mycobiome (Figure 4). Six fungal classes (Dothideomycetes, Eurotiomycetes, Glomeromycetes, Leotiomycetes, Agaricomycetes, and Sordariomycetes) were detected using all three primer sets (Table 2), while Tremellomycetes were detected using both the ITS1F and ITS2 data sets (Table 2). Eight classes were detected using either the ITS1 or ITS2 data sets but not the SSU data set (Table 2). All of these eight classes (Tremellomycetes, Mortierellomycetes, Saccharomycetes, Lecanoromycetes, Pucciniomycetes, Umbelopsidomycetes, Chytridiomycota‐unclassified class, Rozellomycota‐unclassified class) had a low abundance (<0.5% of reads) in all data sets, except that the Mortierellomycetes represented 1% of the reads in one sample in the ITS1 data set (Table 2, Figure 4).

Figure 4

Stack barplot showing the relative abundance of fungal taxa. (A) Phylum level. (B) Class level. ITS1, ITS1 data set; ITS2, ITS2 data set; SSU, SSU data set

Table 2

Fungal classes detected using each data set

Class	ITS1	ITS2	SSU	Phylum	Taxa of particular ecological relevance
Dothideomycetes	+	+	+	Ascomycota	Dark septate endophytes
Eurotiomycetes	+	+	+	Ascomycota	Dark septate endophytes
Glomeromycetes	+	+	+	Glomeromycota	Arbuscular mycorrhizal fungi
Leotiomycetes	+	+	+	Ascomycota	Ericoid mycorrhizal fungi
Agaricomycetes	+	+	+	Basidiomycota
Sordariomycetes	+	+	+	Ascomycota
Tremellomycetes	+	+	–	Basidiomycota
Mortierellomycetes	+	–	–	Mortierellomycota
Saccharomycetes	+	–	–	Ascomycota
Lecanoromycetes	–	+	–	Ascomycota
Pucciniomycetes	–	+	–	Basidiomycota
Umbelopsidomycetes	–	+	–	Mucoromycota
Chytridiomycota‐unclassified class	–	+	–	Chytridiomycota
Rozellomycota‐unclassified class	–	+	–	Rozellomycota

Abbreviations: ITS1, ITS1 data set; ITS2, ITS2 data set; SSU, SSU data set.

DISCUSSION

Both the microscopic and DNA metabarcoding results suggested the presence of fungi associated with the gametophyte. The obvious fungal structures and the diverse fungi detected highlight the importance of fungi during the gametophytic stage of the fern life cycle.

Metabarcoding studies can be strongly impacted by primer choices. Several studies have compared the suitability of the ITS1 and ITS2 regions for metabarcoding (Tedersoo et al., 2015; Taylor et al., 2016; Mbareche et al., 2020), suggesting that ITS1 and ITS2 each have their own advantages. The results generated using the two general fungal primer sets (ITS1F–ITS2, ITS3ngs1_5–ITS4) and the set targeting AMF (NS31–AML2) revealed distinct taxonomic abundances (Figure 4); however, fungal groups of particular interest, namely the AMF and root‐associated fungi (DSE and Leotiomycetes–ericoid fungi), were detectable across the three primer sets (Table 2). Surprisingly, Eurotiomycetes sequences were more abundant in the ITS2 libraries, suggesting that the primer pairs used to generate the ITS1 data set might have not amplified the Eurotiomycetes fungi efficiently, or that the primer pairs for the ITS2 data set preferentially amplified the Eurotiomycetes. Many fungi of the class Eurotiomycetes have been reported as DSEs, and some are beneficial to plants (Jumpponen and Trappe, 1998; Narisawa et al., 2007; Xu et al., 2020). This finding echoes the hypothesis of Lehnert et al. (2017) that DSEs might be abundant in epiphytic ferns. To further test this hypothesis, comparisons among fern mycobiomes in terrestrial and epiphytic habitats should be made using finer taxonomic investigations (e.g., delineate fungi of the class Eurotiomycetes to the level of order or genus). The SSU data set was successfully enriched with the AMF community (Figure 4, Table 2); however, using the same AMF‐specific primer set, Benucci et al. (2020) did not detect an enriched AMF‐dominated community, with only 9.7% of reads in the roots of lycophytes belonging to AMF taxa. This unexpected result might suggest that AMF‐specific primers can be non‐specific, or highlight differences in AMF prevalence within the gametophyte versus the sporophyte fungal community. These discrepancies across studies also suggest that precautions are required when using taxon‐specific primers. The differences in taxonomic abundance detected using the three primer sets indicated that the relative abundances in metabarcoding samples need to be interpreted carefully. When comparing results generated with multiple primer pairs, the presence/absence of data could be considered instead of the quantitative results (Petrolli et al., 2021). In addition to the fungi detected here, the Endogonales of the Mucoromycota (i.e., the FREs) represent another fungal group of interest. Members of the Endogonales have been reported to be inhabitants of lycophytic roots and chlorophyllous gametophytes (Hoysted et al., 2019; Ogura‐Tsujita et al., 2019), but their sequences are difficult to amplify using common universal primers. We attempted amplification using one primer set (forward: NS6, reverse: Endo18S‐1F) to detect FRE in O. pendulum (Benucci et al., 2020). While PCR amplification was successful, the PCR products sequenced using Sanger sequencing appeared to be plant host DNA (data not shown). Ophioderma pendulum may lack Endogonales, and instead the Endogonales‐specific primers may have randomly annealed to the conserved region shared by the eukaryotic lineage and amplified host DNA. More FRE primer sets need to be tested on other subterranean gametophytes to determine the associations of these fungi (Desirò et al., 2017; Ogura‐Tsujita et al., 2019).

The most commonly used metabarcoding platform today is the Illumina MiSeq, which is a highly accurate and affordable system; however, the short‐read length (300 paired ends) is not suitable for other applications beyond community assessment, such as reliable phylogenetic relationship reconstructions (Nilsson et al., 2019a). Long‐read platforms such as PacBio or Oxford Nanopore that enable longer amplicon sequencing are promising tools for uncovering novel fungal lineages and deciphering the evolutionary history of plant–fungal associations in gametophytes (Tedersoo et al., 2018; Loit et al., 2019). As an alternative to metabarcoding, shotgun sequencing methods such as metagenomics or metatranscriptomics can assess the mycobiome assembly with less primer bias (Nilsson et al., 2019a). Here, we detected AMF and other fungi using microscopy and metabarcoding; however, these techniques only allow identification of fungal structures to a shallow taxonomic level (e.g., AMF vs. non‐AMF). Techniques using species‐specific probes (e.g., fluorescence in situ hybridization [FISH]) to bind fungal structures in vivo have the potential to confirm their taxonomic identity (Vági et al., 2014). Although sequencing‐based techniques are powerful, culture‐based methods can sometimes recover fungi missed by sequencing approaches (Oita et al., 2021). Furthermore, the culturable fungi can be used to design co‐culture experiments for studying plant–fungal interactions. Because we did not perform a surface sterilization, it is possible that some surface inhabitants could have been detected. For studies that strictly focus on fungal inhabitants inside the gametophyte, surface sterilization might be possible (U'Ren et al., 2014). To verify the importance of fungal members with gametophytes, a co‐culture experiment can be set up by culturing axenic gametophytes and performing a fungal inoculation (Swatzell et al., 1996; Whittier and Braggins, 2011). Isotope tracing can then be conducted to confirm the presence of nutrient exchange between the gametophytes and the specific fungal partners (Suetsugu et al., 2020). Finally, we can try to determine the function of fungi associated with ferns and lycophytes throughout their life cycle. The comparison of the mycobiome of the gametophytes at different developmental stages with different sporophyte tissues will greatly enhance our understanding of fern and lycophyte biology, and will shed light on vascular plant and fungal ecology as a whole.

CONCLUDING REMARKS

We demonstrated the feasibility of simultaneously applying a DNA metabarcoding technique and microscopy to a single subterranean fern gametophyte. We also discussed the potential concerns in interpreting the outcomes, and pointed out several future directions for method development and biological investigation. This workflow can be applied to non‐subterranean gametophytes as well. We hope such an application can bridge the knowledge of plant–fungal associations between the two key stages—the sporophyte and gametophyte—of fern and lycophyte life cycles.

AUTHOR CONTRIBUTIONS

K.‐H.C. and L.‐Y.K. designed the project. C.‐C.C. conducted molecular work. Q.‐Y.X. conducted microscopic examination. K.‐H.C. analyzed the data and wrote the manuscript. All authors approved the final version of the manuscript.

Appendix S1. CTAB recipe for DNA extraction.

Click here for additional data file.

Appendix S2. Primer sequences for the first PCR of the ITS1, ITS2, and SSU data sets.

Click here for additional data file.

Appendix S3. The PCR program used for the first PCR.

Click here for additional data file.

Appendix S4. Reagents and volumes used for the first PCR.

Click here for additional data file.

Appendix S5. Primer sequences for the second PCR.

Click here for additional data file.

Appendix S6. Reagents and volumes used for the second PCR.

Click here for additional data file.

Appendix S7. The PCR program used for the second PCR.

Click here for additional data file.

AppendiX S8. Cutadapt scripts for the different data sets.

Click here for additional data file.