2021 Volume 62 Issue 6 Pages 356-363
Primer bias toward Tulasnellaceae fungi during PCR is a known issue with metabarcoding analyses for the assessment of orchid mycorrhizal communities. However, this bias had not been evaluated for the fungal communities of epiphytic orchids, which account for 69% of all orchid species diversity. We compared the mycorrhizal communities detected using two primer pairs, a fungal universal primer pair (ITS86F/ITS4) and Tulasnella-specific primer pair (5.8STulngs/ITS4-Tul2), using a mock community of fungal isolates from epiphytic orchids and also environmental samples, including orchid roots and a tree bark tip from the host tree of an epiphytic orchid collected. The detected mycorrhizal communities differed widely depending on the primer pairs used. The fungal universal primer pair successfully identified Ceratobasidiaceae and Serendipitaceae fungi but did not reflect Tulasnellaceae diversity. Tulasnellaceae fungi were mainly detected using the Tulasnella-specific primer pair. These tendencies were observed in both the mock community and environmental samples. These results strongly suggest that the use of a Tulasnella-specific primer in combination with a fungal universal primer is essential for assessing the mycorrhizal communities of orchids through metabarcoding analysis, especially in epiphytic orchids. Our study contributes to further understanding of the diversity of mycorrhizal fungi in orchids.
Environmental DNA metabarcoding using next generation sequencing is a growing approach for the assessment of orchid mycorrhizal communities (Kaur, Phillips, & Sharma, 2021; Xing et al., 2020). This approach allows for the easy detection of a larger number of fungal species from a single mycorrhizal community than does the traditional Sanger sequencing-based DNA barcoding approach (Tedersoo et al., 2010). Due to the vast amounts of data available, the metabarcoding approach has rapidly replaced the Sanger sequencing-based approach. In terrestrial orchids, large numbers of studies have investigated mycorrhizal communities using the metabarcoding approach, and the findings of such studies continue to enhance our understanding of orchid mycorrhizal associations (Oja, Kohout, Tedersoo, Kull, & Kõljalg, 2015; Waud, Brys, Van Landuyt, Lievens, & Jacquemyn, 2017; Duffy, Waud, Schatz, Petanidou, & Jacquemyn, 2019; Kaur et al., 2021). However, few studies have used the metabarcoding approach to examine epiphytic orchids (Herrera et al., 2019; Xing et al., 2020), which account for 69% of all orchid species (Zotz, 2013), and thus the optimal experimental methods for accurately assessing the mycorrhizal communities of epiphytic orchids remain poorly understood.
The utility of DNA metabarcoding for orchid mycorrhizal communities using universal primers for fungal detection is limited due to primer bias, which prevents the detection of particular fungi present in a sample by primer mismatches during PCR (Waud, Busschaert, Ruyters, Jacquemyn, & Lievens, 2014). The mycorrhizal fungi of orchids mainly belong to three Basidiomycota families, the Tulasnellaceae, Ceratobasidiaceae, and Serendipitaceae, collectively described as ‘rhizoctonia-like’ (Yukawa, Ogura-Tsujita, Shefferson, & Yokoyama, 2009). Although the Tulasnellaceae are the most common partner of orchids (Dearnaley, Martos, & Selosse, 2012), its molecular detection is made difficult by mismatches with fungal universal primers (Suárez et al., 2006; Taylor & McCormick, 2008). This is because Tulasnellaceae fungi underwent accelerated evolution of the nuclear ribosomal operon (Taylor & McCormick, 2008), which increased the rate of base substitutions at the primer site compared with other Basidiomycota families. Therefore, Tulasnella-specific primers were used to prevent the underrepresentation of Tulasnellaceae fungi with the Sanger sequencing-based approach (Ogura-Tsujita, Yokoyama, Miyoshi, & Yukawa, 2012; Kartzinel, Trapnell, & Shefferson, 2013).
For metabarcoding analysis, primer bias was assessed against broad fungal taxa associated with terrestrial orchids, including Tulasnellaceae fungi in Waud et al. (2014). This study evaluated the performance of a variety of fungal universal primers and specific primers for orchid mycorrhizal fungi to characterize orchid mycorrhizal communities, and suggested several suitable primer pairs. The fungal universal primer pair ITS86F/ITS4 outperformed other primer combinations. Many studies have investigated mycorrhizal communities using ITS86F/ITS4 in both terrestrial and epiphytic orchids (Cevallos, Sánchez-Rodríguez, Decock, Declerck, & Suárez, 2017; Johnson, Gónzalez‐Chávez, Carrillo‐González, Porras‐Alfaro, & Mueller, 2020). However, Waud et al. (2014) pointed out that this primer pair excludes the Tulasnellaceae fungi including Tulasnella subgroup A defined by Girlanda et al. (2011). Furthermore, previous large-scale surveys using Sanger sequencing-based approaches have shown that the mycorrhizal communities of epiphytic orchids differ from those of terrestrial orchids (Martos et al., 2012; Xing et al., 2019). Despite these issues, no assessments of primer bias have been performed in epiphytic orchids. Appropriate primer selection is required to accurately assess the fungal communities of epiphytic orchids.
In this study, we tested whether primer choice in metabarcoding affects the diversity of the fungal communities, especially Tulasnellaceae communities, detected from epiphytic orchid roots. Two primer pairs, a fungal universal primer pair (ITS86F/ITS4) and a Tulasnella-specific primer pair, were compared by evaluating their abilities to characterize the mycorrhizal communities of bulk DNA samples, which are mock communities of fungi isolated from orchids, and environmental samples including orchid roots or a tree bark tip. Furthermore, to examine whether primer biases were related to particular Tulasnellaceae taxa, a phylogenetic analysis of Tulasnellaceae fungi was performed with additional information on the number of mismatches against used primer sequences.
Mycorrhizal fungi were isolated from nine epiphytic and two terrestrial orchid species distributed in Japan based on the methods of Rammitsu et al. (2021). Then, we identified the fungal taxa of these isolates by sequencing their internal transcribed spacer regions (ITS; spanning the ITS-1, 5.8S, and ITS-2 regions). These sequences were assigned to operational taxonomic units (OTUs) based on 97% sequence similarity across the whole ITS region (Martos et al., 2012).
We choose 20 rhizoctonia isolates (18 OTUs) including 8 Tulasnellaceae (8 OTUs), 11 Ceratobasidiaceae (9 OTUs), and 1 Serendipitaceae (1 OTU) isolates. The other Basidiomycota family, Brachybasidiaceae, isolated from epiphytic orchids, were also mixed as a non-rhizoctonia sample. To test whether we were able to identify different fungi within the same OTU, three Ceratobasidiaceae isolates (CE18A, CE18B, and CE18C) within OTU CE18 were added to the mock community. CE18B and CE18C showed 99.9% (593/597 bp) sequence similarity and had 12 or 13 base substitutions compared with CE18A, respectively. All of the isolate sequences were registered in the DNA Data Bank of Japan (LC597338-LC597358). Eighteen fungal isolates were deposited in the NITE Biological Resource Center in Japan (Table 1; NBRC114085-NBRC114087, NBRC114098, NBRC114194, NBRC114196, NBRC114198, NBRC114909-NBRC114918, NBRC114326).
| Fungal isolates | Isolation source | Information of fungal isolates | |||||
| Life form | Orchid species | Isolate ID | Fungal family | Accession No. | NBRC Accession No. | ||
| CE4 | E | Thrixspermum japonicum | F416 | Ceratobasidiaceae | LC597338 | NBRC 114914 | |
| CE5 | E | Thrixspermum japonicum | F417 | Ceratobasidiaceae | LC597339 | NBRC 114915 | |
| CE6 a | E | Vanda falcata | F313 | Ceratobasidiaceae | LC597340 | NBRC 114194 | |
| CE9 a | E | Thrixspermum japonicum | F414 | Ceratobasidiaceae | LC597341 | NBRC 114913 | |
| CE15 a | E | Vanda falcata | F210 | Ceratobasidiaceae | LC597342 | NBRC 114086 | |
| CE16 | E | Vanda falcata | F214 | Ceratobasidiaceae | LC597343 | NBRC 114087 | |
| CE18A | T | Gastrodia pubilabiata | F86 | Ceratobasidiaceae | LC597344 | − | |
| CE18B a | E | Gastrochilus japonicus | F324 | Ceratobasidiaceae | LC597345 | NBRC 114198 | |
| CE18C | E | Dendrobium officinale | F356 | Ceratobasidiaceae | LC597346 | NBRC 114326 | |
| CE21 | E | Gastrochilus japonicus | F468 | Ceratobasidiaceae | LC597347 | NBRC 114916 | |
| CE22 | E | Phalaenopsis japonica | F499 | Ceratobasidiaceae | LC597348 | NBRC 114917 | |
| TU3 a | E | Thrixspermum japonicum | F267 | Tulasnellaceae | LC597349 | NBRC 114098 | |
| TU10 a | E | Dendrobium officinale | F205 | Tulasnellaceae | LC597350 | NBRC 114085 | |
| TU11 a | E | Vanda falcata | F320 | Tulasnellaceae | LC597351 | NBRC 114196 | |
| TU12 | T | Geodorum densiflorum | F227 | Tulasnellaceae | LC597352 | NBRC 114909 | |
| TU18 a | E | Staurochilus lutchuensis | F330 | Tulasnellaceae | LC597353 | NBRC 114910 | |
| TU22 a | E | Liparis bootanensis | F372 | Tulasnellaceae | LC597354 | NBRC 114912 | |
| TU23 a | E | Luisia teres | F348 | Tulasnellaceae | LC597355 | NBRC 114911 | |
| TU36 | T | Gastrodia pubilabiata | F98 | Tulasnellaceae | LC597356 | − | |
| SE1 a | E | Dendrobium officinale | F345 | Serendipitaceae | LC597357 | − | |
| BR1a | E | Dendrobium moniliforme | F446 | Brachybasidiaceae | LC597358 | NBRC 114918 | |
To create the mock community, DNA of each isolate was adjusted to 0.1 ng/µL and pooled with 5 µL.
a DNA samples from the isolates with DNA < 0.1 ng/µL were also pooled with 5 µL. E, epiphytic species; T, terrestrial species; −, isolates without accession numbers.
Fungal DNA was extracted from the isolates using the methods of Izumitsu et al. (2012). Mycelia (0.1-1 µg) were obtained from a colony on culture medium using a sterilized toothpick and suspended in 50 µL of Tris-EDTA (TE) buffer in a 1.5-mL tube. The tubes were microwaved (600 W) for 1 min. After 30 s at room temperature, the tubes were microwaved for another minute. The tubes were cooled at −20 °C for 10 min and then centrifuged at 10,000 rpm for 5 min. The supernatant was used as the template for PCR. The DNA was quantified with a Qubit dsDNA HS Assay Kit (Invitrogen, Carlsbad, CA, USA). We diluted each DNA sample to the minimum concentration of detection (0.1 ng/µL) and mixed 5 µL of each to make a mock community of fungal isolates. The DNA concentrations in 12 samples were below the detection limit of the assay kit and failed to quantify the exact value; thus, the extracted DNA was less than 0.1 ng/µL and pooled with the 5 µL. The mixture was used as template for PCR amplification.
2.2. Creation of bulk root samples and preparation of the tree bark sampleRoots from three epiphytic orchid species (Dendrobium okinawense Hatus. et Ida, Luisia teres (Thunb.) Blume, and Thrixspermum japonicum (Miq.) Rchb.f.) and a tree bark tip were used as environmental samples (Supplementary Table S1). These samples were collected at two locations in Japan over 1,000 km apart. We collected 1-3 roots from each individual and observed hand-sliced sections of them under a microscope to assess the presence of fungal colonization. Then, the mycorrhizal roots were cut into 0.5-2 cm root fragment samples. We produced six bulk root samples, each of which mixed 2-4 fragment samples from 1-3 orchid individuals (Supplementary Table S1). We also sampled the bark surface (1 cm2) of Rhododendron indicum (L.) Sweet from near (< 1 cm) the roots of T. japonicum.
Fungal DNA was obtained from hyphal coils colonizing the root fragments using the methods of Rammitsu, Yukawa, Yamashita, Isshiki, & Ogura-Tsujita (2020) because the velamen of epiphytic orchids harbors a highly diverse group of ascomycetes (Herrera, Suárez, & Kottke, 2010), which are non-mycorrhizal endophytic fungi. The fragments were washed with TE buffer under a stereomicroscope to remove debris and crushed using tweezers to disperse the hyphal coils into a fresh volume of the same buffer. For each fragment, 100-400 hyphal coils were collected using a micro-pipette, rinsed three times in TE buffer, and triturated in 20 µL of TE buffer within a 1.5-mL tube using a BioMasher II homogenizer (Nippi, Tokyo, Japan). DNA was extracted from the suspension following the above-described method. The tree bark was homogenized with 3-mm zirconia beads at 3,000 rpm for 2 min twice and then DNA was extracted using a DNeasy Plant Mini Kit (QIAGEN, Hilden, Germany). The DNA was quantified, diluted to 0.1 ng/µL, and mixed according to the protocols used for the mock community.
2.3. PCR amplification and next-generation sequencingFungal nuclear ribosomal DNA ITS sequences were amplified using the standard fungal universal primer pair ITS86F/ITS4 (White, Bruns, Lee, & Taylor, 1990; Turenne, Sanche, Hoban, Karlowsky, & Kabani, 1999). Furthermore, we employed a Tulasnella-specific primer pair of 5.8S-Tulngs (5’-CATTCGATGAAGACCGTTGC-3’) that was slightly modified from 5.8S-Tul (Suárez et al., 2006), and ITS4-Tul2 (5’-TTCTTTTCCTCCGCTGAWTA-3’; Oja et al., 2015). To design primers to detect diverse Tulasnellaceae taxa associated with epiphytic orchids, we aligned 39 Tulasnellaceae sequences obtained from epiphytic orchids in recent studies (Martos et al., 2012; Kartzinel et al., 2013; Riofrío et al., 2013; Xing et al., 2019) or in our fungal collections. A new primer was designed for the conserved 5.8S region using Primer3 software (Koressaar & Remm, 2007; Untergasser et al., 2012) and named 5.8S-Tulngs; it has a one-base deletion from 5.8S-Tul designed by Suárez et al. (2006). Sequencing libraries were produced using a two-step protocol, as suggested by Berry, Mahfoudh, Wagner, and Loy (2011). A first PCR amplification was conducted in a total volume of 25 µL containing 1 µL of sample DNA, 12.5 µL of 2× KAPA HiFi HotStart ReadyMix (Kapa Biosystems, Wilmington, MA, USA), and 0.75 µL of each primer (0.3 µM). The PCR protocol was 95 °C for 3 min; 35 cycles of denaturing at 98 °C for 20 s, annealing at 57 °C for 15 s, and extension at 72 °C for 1 min; and a final extension at 72 °C for 1 min. For annealing, we tested the optimum temperature using both primer pairs in the range 55-60 °C in preliminary experiments, because previous metabarcoding studies annealed fungal DNA from epiphytic orchids at 55 °C (Cevallos et al., 2017) or 60 °C (Cevallos, Declerck, & Suárez, 2018; Izuddin, Srivathsan, Lee, Yam, & Webb, 2019). PCR amplifications succeeded at 55-58.4 °C, although 55 °C produced non-specific amplicons visible as a smear after gel electrophoresis. Therefore, we choose 57 °C, the lowest temperature without smears. Each PCR was run with an appropriate negative control (no template). The PCR products were extracted from 1.5% agarose gels and cleaned using a FastGene Gel/RCR Extraction Kit (Nippon Genetics, Tokyo, Japan). The purified PCR products were used as template for a second PCR to attach index and flow cell sequences to the amplicons. The PCR protocol was 98 °C for 3 min; 10 cycles of 98 °C for 20 s and 72 °C for 15 s; and a final extension at 72 °C for 5 min. After purification, the second PCR products were quantified using a Qubit fluorometer and pooled in equimolar quantities to prepare sequencing libraries. The libraries were paired-end sequenced (150 bp × 2) using the Illumina iSeq 100 sequencing platform (San Diego, CA, USA).
2.4. Data processingDue to the PCR product size ranging from 250 to 350 bp, a large number of the pair-end 150 bp reads did not overlap. Furthermore, most of the region included in the forward read was covered by the 5.8S ribosomal region, which is a highly conserved sequence. Therefore, we analyzed only the reverse read, which covered most of the ITS2 region with high sequence divergence.
The extraction of high-quality sequences was performed with the QIIME2 package (Bolyen et al., 2019). All 151 bp reads were trimmed to 81 bp by removing 64 bp and 6 bp from the 3' and 5' ends of the reads, respectively, as the inclusion of these low-quality regions caused the mean quality score to drop below Q30. The DADA2 was used for denoising and the removal of chimeras. The sequences were assigned to OTUs at 97% sequence similarity. The total number of reads for each OTU obtained with ITS86F/ITS4 and 5.8S-Tulngs/ITS4-Tul2 was calculated, and we excluded OTUs with fewer than 100 total reads to remove sequence errors and contaminants. For the remaining OTUs, the taxonomic assignment for each OTU was performed using the UNITE database ver. 8.2 (Abarenkov et al., 2020) with reference sequences of isolates included in a mock community. To ensure the identity of mycorrhizal OTUs, we also compared OTUs against the NCBI nucleotide database using BLAST. In the mock community, the sequences of assigned OTUs had 0-3 base substitutions versus the reference sequences of the isolates, probably due to sequence errors. To attribute these sequences, we performed phylogenetic analyses of the Ceratobasidiaceae and Tulasnellaceae sequences (Supplementary Fig. S1). OTUs that formed monophyletic clades with the reference sequences with high bootstrap values (81%) were regarded as the same isolates.
2.5. Phylogenetic analysisWe conducted a phylogenetic analysis of the Tulasnellaceae using the 5.8S-ITS2 regions of 76 sequences (Supplementary Table S2). Among these, 8 were from isolates within a mock community and 25 were from the main fungi detected from multiple epiphytic orchid species in previous large-scale studies (Martos et al., 2012; Herrera, Kottke, Molina, Méndez, & Suárez, 2018; Xing et al., 2019). To reveal their phylogenetic placement, sequences of representative Tulasnellaceae fungi in Girlanda et al. (2011) and Freitas et al. (2020) were added to the analysis. Two sequences of Multiclavula species were used to represent outgroup taxa, following Cruz, Suárez, Kottke, Piepenbring, and Oberwinkler (2011) and Shefferson, Kull, and Tali (2008).
The sequences were aligned using the MAFFT online server (Madeira et al., 2019). All positions with less than 90% site coverage were eliminated (i.e., < 10% alignment gaps, missing data, and ambiguous bases were allowed at any position). The final dataset comprised 273 bp positions. Phylogenetic analyses were performed using MEGA7 (Kumar, Stecher, & Tamura, 2016). The phylogenetic tree was constructed using the maximum likelihood method based on the Kimura two-parameter model (Kimura, 1980), which was the best substitution model identified by MEGA7. We estimated the relative robustness of branches using the bootstrap method (Felsenstein, Carney, Iacoviello, & Hirsch, 1985) with 1,000 replicates. Furthermore, we counted the number of mismatches against the four primers used for each sequence. The sequences lacking the primer sites were excluded from the mismatch calculation. This alignment was deposited in TreeBASE (http://www.treebase.org/) under the accession number S28425.
The number of Tulasnellaceae isolates included in the mock community was assessed more accurately using the Tulasnella-specific primer pair than using the fungal universal primer pair (Fig. 1). All eight Tulasnellaceae isolates were detected by 5.8S-Tulngs/ITS4-Tul2, while seven of them were detected by ITS86F/ITS4 (Supplementary Table S3). Tulasnellaceae isolates were relatively more abundant with 5.8S-Tulngs/ITS4-Tul2 (99.5%; 159,306 reads) than with ITS86F/ITS4 (0.05%; 124 reads). Both primer pairs detected all Ceratobasidiaceae isolates, but the relative abundances were much greater with ITS86F/ITS4 (83.9%; 207,814 reads) than with 5.8S-Tulngs/ITS4-Tul2 (0.37%; 592 reads). The Serendipitaceae isolate, SE1, was detected only with ITS86F/ITS4 amplification with 0.14% (334 reads) of all reads. Two OTUs that were not included in the mock community, one each for Tulasnellaceae and Serendipitaceae, were detected with a total of 278 reads. The sequences of CE18B and CE18C belonging to a single Ceratobasidiaceae OTU, CE18, could not be distinguished using either primer pair (Supplementary Table S3), while that of CE18A was identified by its unique substitutions located downstream of the ITS2 sequence.
The Tulasnella-specific primer pair allowed for the detection of more Tulasnellaceae OTUs than did the fungal universal primer pair (Fig. 2). Eight Tulasnellaceae OTUs were detected from seven environmental samples with 5.8S-Tulngs/ITS4-Tul2, while five of these were identified by ITS86F/ITS4 (Supplementary Table S4). Tulasnellaceae OTUs were much more abundant with 5.8S-Tulngs/ITS4-Tul2 (94.3%; 896,375 reads) than with ITS86F/ITS4 (0.09%; 1,501 reads; Supplementary Table S4). Four Ceratobasidiaceae OTUs were detected using both primer pairs, but the relative abundance of obtained reads was much larger with ITS86F/ITS4 (6.6%; 114,653 reads) than with 5.8S-Tulngs/ITS4-Tul2 (0.04%; 369 reads). Four Serendipitaceae OTUs were detected with ITS86F/ITS4 (34.0% of all reads; 585,188 reads), and two of these were also identified with 5.8S-Tulngs/ITS4-Tul2 (4.7% of all reads; 44,611 reads). For other Basidiomycota fungi, the relative abundances with ITS86F/ITS4 or 5.8STulngs/ITS4 were 22.6% (14 OTUs; 388,366 reads) and 0.28% (8 OTUs; 2,616 reads), respectively. Ascomycota fungi were detected more frequently with ITS86F/ITS4 (72 OTUs; 26.1% of all reads) than with 5.8S-Tulngs/ITS4-Tul2 (24 OTUs; 0.092% of all reads).
Based on the number of mismatches against the primers used, Tulasnellaceae sequences including the phylogenetic analysis were assigned to four groups (Fig. 3). Sequences within Group I had 3-4 and 2-3 mismatches with ITS86F and ITS4, respectively, while they perfectly matched with both 5.8S-Tulngs and ITS4-Tul2, except for one sequence (JQ240476). The sequences within Group II had 1-3 and 2-4 mismatches with ITS86F and ITS4, respectively, 0-3 mismatches with 5.8S-Tulngs, and were perfectly matched with ITS4-Tul2, except for one sequence (DQ925497). Group III sequences included 3-5 mismatches against ITS86F, ITS4, and 5.8S-Tulngs, but no mismatches for ITS4-Tul2. Group IV sequences matched completely with ITS86F, except for one sequence (DQ925537), but had 3-4 mismatches with 5.8S-Tulngs. Groups I and III corresponded to monophyletic clades supported by bootstrap values of 99% and 98%, respectively (Fig. 3). Group IV included two monophyletic clades with bootstrap values of 85% and 97%, respectively. Of the eight sequences from isolates within the mock community, seven and one fell into Groups I and II, respectively. The representative Tulasnellaceae fungi of epiphytic orchids in previous studies (Martos et al., 2012; Herrera et al., 2018; Xing et al., 2019) were positioned in Groups I (12 sequences), II (10 sequences), and III (3 sequences).
Our results reveal that the diversity of Tulasnellaceae detected using a metabarcoding analysis of the mycorrhizal communities of epiphytic orchids differed widely depending on the primer pairs used. A Tulasnella-specific primer pair, 5.8S-Tulngs/ITS4-Tul2, successfully identified all Tulasnellaceae isolates within a mock community with high relative abundances, while the fungal universal primer pair ITS86F/ITS4 did not detect one Tulasnellaceae isolates, and the abundance was much lower than that for the specific primer pair (Fig. 1; Supplementary Table S3). Conversely, all Ceratobasidiaceae isolates were detected with both primer pairs and showed high abundance with ITS86F/ITS4 in the mock community assessment. A Serendipitaceae isolate was detected only by ITS86F/ITS4 at low frequency. The identification of mycorrhizal communities from environmental samples showed similar results with that of the mock community (Fig. 2). The primer pair 5.8S-Tulngs/ITS4-Tul2 detected mainly Tulasnellaceae, while ITS86F/ITS4 detected Ceratobasidiaceae and Serendipitaceae well. These results show that using a Tulasnella-specific primer pair in combination with a fungal universal primer is essential to accurately assess the diversity of mycorrhizal fungi associated with epiphytic orchids.
Primer bias is caused by sequence mismatches between primers and their target sequences in the Tulasnellaceae taxa. Phylogenetic analysis showed that the Tulasnellaceae fungi identified from epiphytic orchids mainly fell into Group I or II, which have more than three mismatches with the fungal universal primer set ITS86F/ITS4 (Fig. 3). Primer mismatches toward Groups I and II were also observed by Waud et al. (2014) in their assessment of primer biases using a fungal mock community from terrestrial orchids. The primer pairs used by Waud et al. (2014) showed poor amplification of Tulasnella subgroup A, defined by Girlanda et al. (2011), which mainly belonged to Groups I and II in our phylogenetic analysis (Fig. 3). Because the dominant Tulasnellaceae fungi of epiphytic orchids from Réunion Island (Martos et al., 2012), Ecuador (Herrera et al., 2018), and China (Xing et al., 2019) mostly fell into Groups I and II, using 5.8S-Tulngs/ITS4-Tul2 will be effective in identifying the Tulasnellaceae fungi associated with orchids, especially in epiphytic species. For fungi within Group IV, ITS86F will be a suitable primer for detection, but the reverse primer needs to be evaluated (Fig. 3). A new forward primer will be required for the detection of fungi within Group III due to many mismatches against both ITS86F and 5.8S-Tulngs, but ITS4-Tul2 seems effective as a reverse primer.
Pair-end sequencing of 151 bp may allow us to assess the mycorrhizal diversity of orchids. In fact, the present analysis, which used only 81 out of the 151 bp, identified all OTUs, including those in a mock community (Supplementary Table S3). However, two of the three groups belonging to CE18 were indistinguishable from each other in this study; further evaluation of the sequencing conditions, e.g., the number of sequencing cycles for reverse reads, may allow read lengths longer than 81 bp.
Our results strongly suggest that it is essential to use Tulasnella-specific primer sets in combination with fungal universal primer sets to assess the mycorrhizal communities of orchids through metabarcoding analysis, especially for epiphytic orchids. A commonly used fungal universal primer pair sufficiently detected Ceratobasidiaceae and Serendipitaceae fungi but failed to detect Tulasnellaceae diversity. Our results also suggest that particular Tulasnellaceae taxa have strong primer biases toward the primers used. Some Tulasnellaceae fungi might have been overlooked in previous metabarcoding analyses of orchid mycorrhizal fungi due to such biases. Our study contributes to further understanding of the diversity of mycorrhizal fungi in epiphytic orchids.
The authors declare no conflicts of interest. All the experiments undertaken in this study comply with the current laws of the country where they were performed.
We thank S. Abe, T. Abe, M. Amano, M. Ichikawa, A. Kinoshita, N. Kotaka, M. Kudaka, N. Kudaka, N. Morita, Y. Ohshima for collecting environmental samples; Y. Nagano for DNA analyses and valuable advices; S. Inaba for depositing fungal isolates. We also thank the Handling Editor and two anonymous reviewers for their helpful comments. The DNA sequencing analyses were made using a Genetic Analyzer spectrometer at Analytical Research Center for Experimental Sciences, Saga University. This study was partially supported by TBRC Joint Usage Project Grant 2019-2020.
