Mature and indehiscent D. officinale capsules were collected from an artificial cultivation greenhouse in Jinhua, Zhejiang Province, in November 2018. The obtained capsules were dried naturally at room temperature (25°C) for ~1 week. Next, seeds were cleaned from capsule debris, mixed, and stored in wax paper on silica gel at 4°C. Two mycorrhizal fungi, Tulasnella sp. (strain no. S6) and Serendipita sp. (strain no. 12825), were isolated from the roots of Dendrobium spp. and deposited in the Institute of Medicinal Plant Development, Chinese Academy of Medical Science. Free-living mycelium was cultivated on a PDA medium at 25°C in the dark.

The seed germination experiment was designed as shown in Figure 1. Before sowing, the seeds were immersed in sterile water and allowed to stand for 1–2 h. Seeds in the bottom of the tube were chosen for subsequent experiments. First, the seeds were surface sterilized in 1% NaClO for 3 min, rinsed three times, and diluted with sterile water into a suitable seed suspension for sowing. For symbiotic germination, 400 μL of seed suspension was dispensed with a pipette onto a square of autoclaved nylon cloth (4 × 4 cm) on oatmeal agar plates (OMA, 0.25% oatmeal, and 1% agar, 20 ml in total) in a 9 cm petri dish, and then four 5 × 5 mm inocula of Tulasnella sp. S6 (or Serendipita sp. 12825) on PDA medium were placed in each OMA for coculture. For asymbiotic germination, only 400 μL of seed suspension was dropped onto the 1/2 MS culture medium. Plates were incubated at 25°C with a 12 h/12 h light-dark cycle. Seed germination and protocorm development were evaluated under a dissecting stereomicroscope every 2 days. The seed developmental stages were defined according to the previous study by Stewart and Zettler (2002).

Symbiotic (with Tulasnella sp. or Serendipita sp.), asymbiotic germination seeds at stage 2 (germination), stage 3 (protocorm), stage 4 (seedling), and free-living mycelium of both fungi were collected either immediately frozen in liquid nitrogen or stored at −80°C for RNA extraction. There are three biological replicates for each sample. Total RNA of 11 samples including three samples for asymbiotic, three symbiotic samples with Tulasnella, three symbiotic samples with Serendipita (stages 2, 3, and 4, respectively), and two free-living mycelium fungal samples Tulasnella sp. and Serendipita sp., respectively) was extracted using the RNeasy Plant Mini Kit (QIAGEN, Hilden, Germany), and the quality, integrity, and quantification of RNA were assessed following our previous study (Chen et al., 2017). Finally, 2 μg of high-quality total RNA (RNA Quality number, RON > 7) was used for stranded RNA sequencing library preparation. The library was constructed using a KC-Digital TM Stranded mRNA-seq Library Prep Kit for Illumina® (Kangce Technology Co., LTD, Wuhan, China) following the manufacturer's instructions. Duplication bias in PCR and sequencing steps were eliminated using a unique molecular identifier (UMI) of 8 random bases to label the preamplified cDNA molecules with this kit. The library products corresponding to 200–500 bps were enriched, quantified, and finally sequenced on an Illumina HiSeq 2500 sequencer (Illumina, San Diego CA, USA) using the paired-end (PE) 150 strategy (Kangce Technology Co., LTD, Wuhan, China).

Prior to assembly and mapping, the raw sequencing data were first filtered by Trimmomatic (version.36) (Bolger et al., 2014), and then the low-quality reads and the reads containing adapter were trimmed using default parameters (PE -phred33 ILLUMINACLIP:adaptor_file:2:30:5 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 HEADCROP:0 MINLEN:36). Clean reads were further treated with in-house scripts to eliminate duplication bias introduced in library preparation and sequencing. Briefly, clean reads with the same UMI sequence were first clustered together, and then they were compared to each other using pairwise alignment. Reads with sequence identity over 95% were extracted into a new subcluster. After all subclusters had been generated, multiple sequence alignment was performed to obtain one consensus sequence for each subcluster. Using these steps, any errors and biases introduced by PCR amplification or sequencing were eliminated. The de-duplicated consensus sequences with high quality were used for downstream analyses.

Genomics mapping and alignment program analysis were performed based on the reference genome of Dendrobium catenatum (https://www.ncbi.nlm.nih.gov/assembly/GCF_001605985.2) using STAR software (version 2.5.3a) with default parameters (Dobin et al., 2013; Dobin and Gingeras, 2015). Functional annotation of predicted genes was first performed using diamond Blastx in UniProt (Universal Protein) to obtain the protein ID, and then the annotation was searched in the database of Nr (NCBI non-redundant protein sequences), Pfam (Protein family), Rfam (RNA family), eggNog (evolutionary genealogy of genes: non-supervised orthologous groups), Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) using the same ID. Gene expression was qualified using the RSEM software package (RNA-Seq by Expectation-Maximization). Genes that were differentially expressed between groups were identified using the edgeR package (version 3.12.1) (Li and Dewey, 2011). The threshold for significantly differentially expressed genes (Robinson et al., 2010) (DEGs) was set at a fold change of ≥2 and an FDR corrected p-value < 0.05 (Robinson and Oshlack, 2010).

Fungal de novo transcriptome assemblies were first reconstructed using Trinity based on the free-living mycelium of Tulasnella sp. (S6) and Serendipita sp. (12825), and the longest transcript was chosen as the unigene. Next, the unmapped reads of symbiotic samples (after mapping to the plant genome) were subjected to mapping analysis using the unigenes by Kangce Technology Co., Ltd., Wuhan, China, followed by Grabherr et al. (2011). Because we focused our analysis on the gene scale, further analysis was performed on Trinity unigenes, i.e., the longest isoform per Trinity gene. TransDecoder (version 5.5; http://transdecoder.sourceforge.net/) was used to identify putative coding regions from the assembled unigenes, and only the best-predicted protein was retained for each transcript (single_best_only option). Carbohydrate active enzyme (CAZyme) assignment was performed with dbCAN2 (Zhang et al., 2018) using HMMER, DIAMOND, and Hotpep prediction. Only domains detected by at least two tools were conserved.

Transcriptomic data were generated for the symbiotic germination of D. officinale seeds with the two mycorrhizal fungal species, asymbiotic germination seeds, and free-living mycelium of both fungi (Figure 1; Table 1). An average of 64.52% of the total 24,811 plant genes were detected to have transcriptional activity in at least one of the symbiotic germination stages (RPKM value ≥ 5) of D. officinale seeds inoculated with Tulasnella sp. A total of 1,287 genes were common differentially expressed (|fold change| ≥ 2; FDR <.05) within the three symbiotic germination stages compared to that of asymbiotic germination (Figure 2A). Similarly, an average of 66.6% of plant genes exhibited expression values (RPKM value ≥ 5) in at least one of the symbiotic samples of D. officinale seeds inoculated with Serendipita sp., and 1,943 plant genes were commonly differentially expressed in the three symbiotic stages compared to those of asymbiotic germination (Figure 2B). When we examined gene expression between both symbiotic germination groups of D. officinale inoculated with Tulasnella sp. and Serendipita sp., respectively, 1,003 common regulated plant genes were identified (Figure 2C). Among 1,003 differentially expressed genes (DEGs), 452 genes were upregulated in symbiotic groups and were clustered into nine clades based on expression profile (Figure 3; Supplementary Table S1). The genes with the highest expression level were primarily restricted to clusters I-V (Figure 3), and they displayed highly similar expression profiles in symbiotic plants inoculated with the two different mycorrhizal fungi, including the genes encoding a nodulin-like protein, Ras domain protein, sugar transporter, lysM domain protein, mannose-binding lectin, and histone-like transcription factor. These upregulated genes took part in molecular functions, biological processes, and cell components.

To further explore commonly expressed plant genes of D. officinale seeds induced by both fungal species, GO function and KEGG pathway enrichment analyses were performed (Figure 4). Compared to asymbiotic germination, all DEGs across various developmental stages of symbiotic D. officinale seeds were primarily enriched in organic or inorganic substance transport and sucrose, trehalose, and lipid synthesis and transport (Figure 4A). Several pathways, such as plant–pathogen interactions, plant hormone signal transduction, starch, and sucrose metabolism and phenylpropanoid biosynthesis, protein processing in the endoplasmic reticulum, lysosomes (for Tulasnella sp. pairs), and peroxisome (for Serendipita sp. pairs) were significantly enriched, and genes involved in these pathways exhibited a similar pattern in interaction with both Tulasnella sp. and Serendipita sp. (Figures 4B,C). Most genes in the plant-pathogen interaction pathway were encoded by calcium-binding proteins, calmodulin, LRR receptor-like, WRKY transcription factor, etc. The genes encoding auxin-responsive proteins and DELLA proteins (involved in plant hormone signal transduction) were significantly enriched during the interaction between D. officinale seeds and their mycelia. Genes encoding trehalase, beta-amylase, pectinesterase, hexokinase, glucan endo-1,3-beta-glucosidase, etc., which are implicated in starch and sucrose metabolism, were also enriched in symbiotic germination from early germination to seedling formation. Twenty-five upregulated common expression genes of the host plant (fold change ≥ 5, FDR ≤ 0.05, compared to asymbiotic germination) are listed in Supplementary Table S2, and these genes typically encoded aquaporin-like, cytochrome P450, nodulin-like, glycoside hydrolases 35 (GH35), GH18, sugar transporter, etc.

When we compared gene expression at the same germination stage of D. officinale seeds between inoculation with Tulasnella sp. and Serendipita sp., we found more DEGs in the symbiotic protocorm stage (1,852 genes, stage 3) than in the germination stage (934 genes, stage 2), implicating the protocorm stage is likely a crucial stage responding to different fungal invasions during orchid mycorrhizal development. In total, 350 plant genes were identified as significantly differentially expressed in symbiotic seeds of D. officinale across all germination stages during interaction with Tulasnella sp. compared to those symbiotic with Serendipita sp. (Figure 5). GH17, GH79, UPD-glucoronsyl transferase, UbiA prenyltransferase, CYP450, nodulin-like, hydrophobic seed protein, etc., were significantly upregulated in D. officinale seeds inoculated with Tulasnella, and upregulated plant genes in various symbiotic stages inoculated with Tulasnella sp. were primarily enriched in plant hormone signal transduction and phenylpropanoid biosynthesis, especially during the initial germination stage and photosynthesis pathway, and were significantly enriched in the protocorm stage (Figure 5).

We identified 15 common DEGs involved in plant cell wall constitution and remodeling at each symbiotic stages (Table 2) during seed germination of D. officinale inoculated with fungi. Among them, genes encoding epidermis-specific secreted glycoprotein, proline-rich receptor-like protein, and leucine-rich repeat (LRR) extensin-like protein were upregulated in the symbiotic stages and even more highly upregulated in the protocorm stage (stage 3) and seedling stage (stage 4). The gene encoding LRR extensin-like protein (LOC110096373) exhibited the highest expression of 20.2-fold in the protocorm stage of D. officinale seeds inoculated with Serendipita sp., and the gene encoding proline-rich receptor-like protein kinase displayed more highly upregulated expression with more than a 25-fold induction in the symbiotic stage compared to asymbiotic germination inoculated with either Tulasnella or Serendipita across the entire germination stage. Moreover, genes encoding enzymes that function in cell wall biosynthesis also displayed significantly differential expression between symbiotic and asymbiotic germination, such as cellulose syntheses, pectinacetylesterase, and pectinesterase. The gene encoding microtubule-associated protein RP was also significantly upregulated during symbiotic germination.

Based on our transcriptomic data, we identified 266 putative genes encoding CAZymes in the Tulasnella sp. × D. officinale symbiotic transcriptome, and among them, 66 genes were differentially expressed in at least one symbiotic stage compared to that of the FLM (Figure 6), including 3 carbohydrate-binding modules (CBM) family members (4 genes), 3 members of the carbohydrate esterase (CE) family (7 genes), 5 members of the auxiliary activity (AA) (7 genes), 13 glycoside hydrolase (GH) family members (20 genes), 14 glycosyltransferases (GT) family members (27 genes), and two pectin lyase members (PL35, PL8). Moreover, CBM43, CBM48, CBM13, CE8, AA2, AA9, GH1, GH79, GH38, and GH45 were upregulated in the entire symbiotic stage with more than 5-fold changes compared to that of the FLM.

Compared to the interaction between Tulasnella sp. × D. officinale seeds, the transcriptomics analysis of Serendipita sp. × D. officinale identified 270 putative genes encoding CAZymes. Among them, 99 genes were differentially expressed in the symbiotic stage (Figure 6). Genes encoding AA class proteins (21 genes for 7 members) targeted to pectin were significantly upregulated in symbiotic conditions by more than 5-fold compared to free living mycelium. In addition, genes in the GH family, such as GH15, GH17, GH3, GH5, GH38, GH47, and GT family members (GT2, GT22, and GT59), were upregulated in the symbiotic stage compared to FLM.

We hypothesized that the different mycorrhizal fungal species would share common gene expression patterns during interaction with seeds of the same host plant of D. officinale. To identify possible orthologous genes between the two mycorrhizal fungi, we performed a bidirectional BlastP analysis on Tulasnella sp. and Serendipita sp. at an e-value < 10⁻⁵. Additionally, we performed homology analysis using BlastP search with the published Tulasnella (AL13/4D version 1) and Serendipita (MAFF 305830 v1.0) genomes for gene annotation.

We identified 4,954 genes (25.21%) that were orthologous between Tulasnella sp. (S6) and Serendipita sp. (12825) based on the fungal transcriptomics data. Among them, 1,682 genes were expressed in symbiotic Serendipita sp. and 2,464 genes were expressed in symbiotic Tulasnella sp. in at least one symbiotic stage during the interaction with D. officinale seeds. A total of 936 orthologous genes were shared between Tulasnella sp. (S6) and Serendipita sp. (12825) in at least one symbiotic stage during colonization in seeds of D. officinale. Most of these genes exhibited similar expression profiles. According to gene expression level, 936 orthologous genes were divided into five clusters (Supplementary Figure S2; Supplementary Tables S2, S3). For example, in cluster 1, genes were primarily upregulated in both fungi during early symbiotic stages with seeds (stage 2, starting germination), and most of these genes were involved in energy production and conservation, translation and posttranslational modification, amino acid transport, and metabolism and lipid transport and metabolism. Genes encoding a fungal transcriptional regulatory protein, chitin synthase, zinc finger protein, serine/threonine-protein kinase-related, ribosomal protein, transcription factor, sugar transporter, ABC transporter-like, and acyltransferase, which participate in signal transduction, posttranslational modification, lipid transport and metabolism (androgen and estrogen metabolism), were upregulated in symbiotic Tulasnella sp. but downregulated in symbiotic Serendipita sp. during the entire process of interaction with D. officinale seeds (cluster 2), especially during the germination and protocorm stages (stages 2 and 3), indicating strikingly different expression profiles. Genes encoding ribosomal protein, protein synthesis factor (translation, ribosomal structure, and biogenesis), aminotransferase, spermine synthase, methionine synthase, glutamate-5-semialdehyde dehydrogenase, etc. (amino acid metabolism) were significantly upregulated in symbiotic Serendipita sp. compared to FLM but were not significantly differentially expressed in symbiotic Tulasnella sp. group (cluster 5). In addition, genes encoding rasGAP protein, Ca2+/calmodulin-dependent protein kinase, serine/threonine protein kinases (signal transduction mechanisms), sterol desaturase, and ergosterol biosynthesis protein (lipid metabolism) were upregulated in symbiotic Serendipita sp. compared to FLM. Notably, cluster 4 included 29 orthologous genes that were both upregulated in the two mycorrhizal fungi during interaction with D. officinale seeds in various symbiotic stages, especially in Tulasnella with the D. officinale seed group. These 29 genes primarily encode ubiquitin, short-chain dehydrogenase, ribosomal protein, histone, heat shock protein, glutathione S-transferase, carbohydrate-binding, ATPase, and amine oxidase, which are involved in carbohydrate metabolism, energy metabolism, and glutathione metabolism.

To further understand the potential function of fungal orthologous genes during interaction with orchid seeds, we analyzed the fungal genes encoding CAZymes, protease, lipase, and SSPs, which were reported to likely be involved in mycorrhizal symbiosis in previous studies. The results revealed a total of 126 of 936 commonly expressed orthologous genes encoding the three specific protein categories (61 proteases, 24 CAZymes, and 31 SPs), and the expression profiles of the two symbiotic fungi during seed symbiotic germination of D. officinale at different stages (stage 2, stage 3, and stage 4) were analyzed (Figure 7). Twenty-four CAZyme genes were differentially regulated at the transcriptional level in at least one symbiotic stage compared to FLM. Genes encoding AA9 were both highly upregulated at the early germination stage (fungal invasion and peloton formation) after inoculation with Tulasnella sp. or Serendipita sp. with more than 100-fold changes. The gene encoding AA1 was upregulated during the seedling stage (fungal peloton digestion) in D. officinale with Serendipita sp. and Tulasnella sp. Genes encoding GT66 and GT4 and AA family proteins (AA2, AA3, and AA6) were upregulated during the early stage of germination and the protocorm formation of seeds with Serendipita. Genes encoding GH family proteins (GH16, GH9, and GH5) were markedly upregulated in at least one stage of seed symbiotic germination with Tulasnella sp. Among the 31 SPs, the genes encoding peptides S8, peptides M28, and carbohydrate-binding WSC domain proteins, representing secreted proteins, were specifically upregulated in the symbiotic Tulasnella group with more than a 100-fold-change compared to FLM. GH61 was upregulated in the initial invasion stages (or early germination stage) during symbiotic germination in the two symbiotic fungi, and GH5 was highly upregulated in symbiotic Serendipita.

In orthologous and commonly expressed gene sets, genes encoding proteases, such as asparaginase, HMG-CoA, septin, ASF, and peptidase, were significantly upregulated in at least one germination stage interaction in either Tulasnella sp. or Serendipita sp. Genes encoding cytochrome b5 and major facilitator superfamily proteins were specifically upregulated in symbiotic Tulasnella compared to FLM, with more than a 100-fold-change, while ribosomal S7 and fatty acid hydroxylase were highly upregulated in symbiotic Serendipita sp.

Our comparative transcriptomics analysis revealed that an important proportion of plant genes (>60%) were expressed during the symbiotic germination of D. officinale. Most exhibited similar expression profiles at the same symbiotic germination stage with different mycorrhizal fungi. Genes encoding universal stress protein, UbiA, aquaporin-like, cytochrome P450 family, nodulin-like, plant lipid transfer protein; RAS domain family, and sugar (and other) transporter were significantly upregulated between the two mycorrhizal fungal combinations, indicating their common roles in the core mechanisms of symbiotic recognition, nutrition transport, or other aspects between the orchid host and the mycorrhizal fungi. These results are consistent with previous studies (Balestrini et al., 2014; Chen et al., 2017).

The plant cell wall is the headmost interface where interactions occur between the host plant and its mycobionts (Balestrini and Bonfante, 2014; Lionetti and Métraux, 2014). The plant cell wall is primarily composed of polysaccharides, such as cellulose, hemicelluloses, and pectins. In addition, there are some structural cell wall proteins characterized by highly repetitive sequence motifs, which are especially rich in specific amino acids, e.g., proline, glycine, and hydroxyproline glycoproteins, to modify the cell wall (Herger et al., 2019). In our present study, genes encoding epidermis-specific secreted glycoprotein, proline-rich receptor-like protein, leucine-rich repeat (LRR) extensin-like protein, and extensin-like protein were significantly upregulated during the symbiotic stage in D. officinale seeds after inoculation with Tulasnella sp. and Serendipita sp. Extensins are an abundant group of hydroxyproline-rich glycoproteins and are responsible for various processes, including embryonic development, root hair growth, cell wall assembly and structure, and biotic and abiotic stress responses (Castilleux et al., 2018). Most studies indicate that extensins exert their roles in plant defense by strengthening the cell wall, decreasing pathogen invasion, or facilitating the attachment of symbiotic organisms (Castilleux et al., 2021). In addition, new insights have been provided regarding the structure and role of extensins and their glycosylation in plant–microbe interactions using a set of anti-extensin-specific monoclonal antibodies (Castilleux et al., 2020), indicating that modeling of the cell wall architecture is likely to impact cell binding and pathogen colonization of roots. The extensins, especially arabinosylation, could even serve as markers of the plant immune response during plant-microbe interactions (Castilleux et al., 2018). In orchid symbiotic germination, previous studies revealed that the signals for the dynein monoclonal antibody JIM11 epitope were stronger in symbiotic Dendrobium protocorms than in asymbiotic ones and were primarily located in the plant walls and in special plant cells colonized by fungi rather than in the apical meristem. This indicated that the extensins were not only related to fungal attachment in a symbiotic relationship but were also likely critical for limiting fungal colonization in basal cells and preventing fungal spread inside the protocorms (Li et al., 2018). Regardless, our transcriptomics analysis indicated that the two mycorrhizal fungi both triggered plant genes encoding extensions to be significantly upregulated in the symbiotic stage during interaction with D. officinale seeds, consistent with the results of previous studies and further provided supported data.

Immunofluorescence labeling revealed that microtubules were altered and reorganized adjacent to the symbiosis interface after the fungus invaded root cells (Genre and Bonfante, 2002), although the ultrastructural changes were not illustrated in detail in the orchid mycorrhiza. In our study, the gene encoding microtubule-associated protein was upregulated in the symbiotic stage, indicating that cytoskeletal rearrangements likely occurred during the invasion of orchid seeds by the fungus. In addition, genes related to plant cell biosyntheses, such as those encoding cellulose synthase and pectinesterase, were significantly upregulated when inoculated with either Tulasnella sp. or Serendipita sp. compared to asymbiotic germination, indicating that plant cell wall biosynthesis and reinforcement was very active during the interaction between the mycorrhizal fungi and D. officinale seeds.

During interaction with the same host plant, Tulasnella sp. induced more genes involved in the symbiotic germination of D. officinale than did Serendipita sp. For example, during the initial invasion stage (correspondingly, the early germination stage), a total of 8,128 fungal genes in Tulasnella (5,268 up- and 2,960 downregulated) were differentially expressed compared to FLM, while there were 4,630 differentially expressed genes in Serendipita sp. (2,066 up to and 2,564 down). Although no genome information is available for either mycorrhizal fungus, using an overview of the reported genome of Tulasnella calospora (62.39 Mb, 19,635 genes) and Serendipita vermifera (38.09 Mb, 15,317 genes) (Kohler et al., 2015), we suspect that different taxon of mycorrhizal fungi could be contributed to the gene expression difference during interaction with the orchid plant to the extent.

In addition, we identified 4,954 orthologous genes based on the transcriptome of Serendipita sp. and Tulasnella sp. A total of 1,682 genes were expressed in symbiotic Serendipita sp. and 2,464 genes were expressed in symbiotic Tulasnella sp. compared to the FLM during interaction with D. officinale seeds (in at least a symbiotic condition). A total of 936 orthologous genes were commonly expressed, and most of the genes exhibited similar expression during mycorrhizal interactions with any mycobiont, suggesting that different mycorrhizal fungi shared a common set of genes involved in symbiotic responses.

Genome analysis revealed that orchid mycorrhizal symbionts had the largest set of CAZymes supporting their dual saprotrophic/symbiotic lifestyles (Martin et al., 2016). For example, T. calospora has 7 GH6, 27 GH7, and 33 LPMO genes for the degradation of crystalline cellulose (Kohler et al., 2015). Proteins with a cellulose-binding domain (CBM1) are also abundant in orchid mycorrhizal fungi (ORM). In our present study, 66 DEGs encoding CAZymes were identified in symbiotic Tulasnella sp., 99 genes encoding CAZymes were differentially expressed in symbiotic Serendipita sp., and AA family members, CBM, CE, and GH, were upregulated in symbiotic fungi compared with FLM. CBM 43 (β-1,3-glucanosyltransferase, former X8) and CBM48 (starch debranching enzymes) are not common in the known various fungal lifestyles, and their exact function in orchid mycorrhizal fungi is also unclear, but their significantly high expression in symbiotic Tulasnella sp. implies that they likely play a vital role in fungal cell wall construction because they are required for fungal cell wall component dual β-(1,3)-glucan elongation and branching (Aimanianda et al., 2017). In addition, chitin deacetylases of the CE4 family were upregulated in symbiotic Tulasnella sp. and chitosan, the product of Tulasnella chitin deacetylases has been identified by metabolomic analyses in symbiotic orchid protocorms (Ghirardo et al., 2020). The role of CE4 is suspected to protect the fungal cell wall from acting by plant chitin or decreasing the plant defense responses when the fungal colonization transfers chitin to chitosan in the ECM (Veneault-Fourrey et al., 2014).

Genes encoding the lignin-degrading auxiliary enzyme AA class were also upregulated in symbiotic Tulasnella sp. and Serendipita compared to FLM, especially AA9, a family of fungal origin with only lytic polysaccharide monooxygenases (Vandhana et al., 2021). AA9 is a mono-copper enzyme family and participates in the degradation of lignocellulose via the oxidative cleavage of celluloses, cello-oligosaccharides, or hemicelluloses (Zhang, 2020). Zarattini et al. (2021) employed AA9 LPMO reaction products to trigger innate immune responses in the model plant Arabidopsis thaliana during fungal pathogenesis, while results on AA17 LPMO indicated that LPMO action would help overcome plant defense barriers (Sabbadin et al., 2021). Thus, the role of AA9 in orchid symbiotic germination needs to be further investigated. In addition, the glycoside hydrolase family GH5 and GH28 were recently confirmed to be involved in cell wall remodeling for ECM symbiosis as enzymatic effectors (Zhang et al., 2021). In our study, the expression of many glycoside hydrolase family members was induced in symbiotic tissue, such as genes encoding the hemicellulose-active enzyme GH3 (β-glucosidases), which were both upregulated in the symbiotic condition of the two fungi, indicating that they likely play a key role in orchid mycorrhizal establishment and development.

Of note, several genes encoding glycosyltransferases (GTs) were differentially expressed in symbiotic Tulasnella sp. or Serendipita sp. mycelium, and there were 9 GT family proteins (GT2, GT4, GT20, GT22, GT32, GT35, GT48, GT66) with induced expression among the 24 orthologous genes encoding CAZymes. Glycosyltransferases are enzymes that catalyze the transfer of sugar moieties from activating donor molecules to specific acceptor molecules, forming glycosidic bonds and are involved in the biosynthesis of diverse carbohydrates (http://www.cazy.org/GlycosylTransferase-family). A recent result indicated that GT2 likely participates in the synthesis of extracellular or outer cell wall polysaccharides, which play a key role in facilitating many interactions between plants and fungi by enabling hyphal growth on solid matrices (King et al., 2017). However, the exact function of glycosyltransferases in orchid mycorrhizae still needs to be characterized. Adamo et al. (2020) primarily explored the roles of CAZymes of orchid mycorrhizal fungi in the development of unsuccessful and successful interactions and suspected that the expression of some key CAZymes was related to OMF switches from symbiotic to saprotrophic growth. However, the mycorrhizal relationship is particularly complex in orchids, and successful mycrrhizal interactions between orchids and their symbionts can be controlled by a balance between fungal invasion and plant defense, balancing nutritional supply and demand; therefore, the specific nutritional strategy of orchids should be considered for mycorrhizal establishment.

In addition, Tulasnella and Serendipita possessed several CAZyme, but they display different expression profiles, indicating that their function could be different in fungal species. Moreover, only a few orthologous genes encoding CAZymes were identified. Among the 24 orthologous DEGs encoding CAZymes, most exhibited distinct expression profiles between Tulasnella and Serendipita, and only genes encoding GH5, AA1, AA2, CBM1, and GT66 were upregulated in at least one symbiotic stage in both Tulasnella and Serendipita. Thus, about function of CAZymes in different mycobionts during interaction with the host plant still need to be addressed furtherly in the future.