The experimental setup for preparing materials for RNA-seq is shown in Supplementary Figure S1A . Considering the challenges in collecting sufficient symbiotically germinated protocorm samples for transcriptome analyses, we cultured the symbiotically germinated protocorms and grew them. The resulting seedlings that developed roots were used as plant material for this study. First, to prepare protocorms by symbiotic germination, mature capsules of C. macranthos var. rebunense were collected from Rebun Island, Japan in September 2021, with permission from the Ministry of Environment, Japan. In vitro symbiotic germination was conducted as described previously using the OM fungus WO97, isolated from the roots of C. macranthos var. rebunense (Shimura et al., 2009; Kambara et al., 2023). We here used two subcultures of the OM fungus WO97, WO97-R and WO97-N. These two are vegetatively propagated clones derived from WO97, and their rDNA-ITS sequences are identical. These subcultures have been maintained in our laboratory and used in annual germination assays for C. macranthos var. rebunense. In preliminary experiments in symbiotic germination, we have observed that the ability of the WO97 strain to induce germination varied depending on the subcultured mycelium: WO97-R, which induced seed germination and protocorm development in vitro without differences between replicates and exhibited a stable symbiotic relationship, and WO97-N, which induced germination and protocorm growth at frequencies varying among replicates and exhibited an unstable relationship. To examine whether there were differences in the responses of the orchid when inoculated by subcultures with different colonization capabilities, WO97-R and WO97-N were selected in this experiment.

After approximately 6 months of co-cultivation with the fungus, well-grown protocorms with meristem formation were transferred to a modified T medium (based on Tsutsui and Tomita, 1990, with yeast extract excluded) supplemented with a fungicide (Tsutsui and Tomita, 1990; Shimura and Koda, 2005). Protocorms generated by inoculation with WO97-R and those by inoculation with WO97-N were transplanted separately ( Supplementary Figure S1A ). An incubation step with a fungicide was performed to promote protocorm growth and root development on a nutrient-containing medium. The step was also performed to eliminate fungal hyphae that had colonized the orchid tissues during germination, prior to the second inoculation. To prepare the inoculum plate for the second inoculation, a piece of agar containing WO97 hyphae was excised from the plate on which the protocorms had been transplanted; subsequently, the hyphae were allowed to grow again on fresh oatmeal medium and used for the second inoculation. The 20–30 seedlings cultured with fungicide for 4 months were transferred to the WO97-R or WO97-N growing oatmeal medium and further cultured at 20°C for 2 or 4 weeks.

To evaluate fungal colonization, a piece of tissue was taken from each of a few root slices prepared for RNA extraction, fixed in formalin–acetic acid–alcohol (FAA), and 10 to 11 sections from each inoculated seedling were used for observation. Fungal colonization and peloton formation were observed according to Yamamoto et al. (2017). This method has been used to observe protocorm tissue infected with fungi (Yamamoto et al., 2017), but it was also applicable to the mycorrhizal tissue of the small seedlings’ roots used in our experiment. After removing the fibrous tissues using a dissecting needle, the stained root tissues were transferred onto a microscope slide. Then, the colonization status of the cells was observed under a microscope (Nikon ECIPSE E600, Tokyo, Japan), and images were captured using a Penguin 600CL (Pixera, Osaka, Japan).

For RNA extraction, root tissue slices (approximately 0.05 g) were collected from 8–10 seedlings after 2 and 4 weeks of culture; each group of samples included three replicates. Root tissues from seedlings before fungal inoculation (pre-inoculation) were prepared in triplicates as controls. Total RNA was extracted and treated with DNase following a previously described protocol (Shimura et al., 2018), and the purified RNA was used in RNA-seq. Library preparation using the TruSeq Stranded mRNA Sample Prep Kit (Illumina) and sequencing on the Illumina NovaSeq platform were outsourced to Macrogen Inc. (Seoul, Republic of Korea). Each library was sequenced, yielding approximately 50 million 151-bp paired-end reads. Data obtained by RNA-seq were submitted to the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database (PRJNA1231755).

The raw reads were processed using fastp v0.23.4 (Chen, 2023) for base correction, adapter trimming, and quality filtering. Considering the unavailability of the C. macranthos var. rebunense genome sequence, reference sequences were constructed by de novo assembly using Trinity v2.15.1 (Grabherr et al., 2013). To expedite the assembly process, one-third of the read data (ca. 1.66 Gb) were randomly extracted from each of the three replicates in a treatment and merged using seqkit v2.1.0 (Shen et al., 2016) to create representative read data for each treatment. Five representative reads were used as input data for de novo assembly on Trinity. Redundant contigs were removed from the assembly using CD-HIT-EST v4.8.1 (parameter -c set to 0.8, Li and Godzik, 2006) and the EvidentialGene tr2aacds pipeline (Gilbert, 2013; http://arthropods.eugenes.org/EvidentialGene/about/EvidentialGene_trassembly_pipe.html). The quality of the contig set was assessed based on the contig number, mapping rate, and BUSCO (v5.6.1) completeness, following a previous report (Kambara et al., 2023).

To determine the effects of the fungus-derived reads or contigs on subsequent analyses, we prepared a contig set that excluded as many fungal sequences as possible. We used hisat2 v2.2.1 (Kim et al., 2019) to identify the reads that were significantly mapped to the WO97 draft genome sequence (our unpublished data). The mapped reads were removed from each sample’s clean reads using samtools v1.6 (Li et al., 2009), and de novo assembly was conducted using the remaining reads as described above. WO97-related contigs were further identified in the assembled contigs. Specifically, BLASTn alignment against the contig set was performed by setting the predicted CDS sequences from the WO97 draft genome as a reference and stringent parameters (e-value < 1e-50, identity ≥ 90%). Identified hit contigs were considered WO97-derived transcripts, and the corresponding contigs were removed from the contig set as described in the Results section.

The abundance of the transcripts was estimated by mapping the total clean reads to the contig set using Bowtie2 v2.5.3 (Haas et al., 2013), converted to BAM format using samtools v1.19.2 (Cock et al., 2015), and quantified by salmon v1.12.0 (Patro et al., 2017). For differentially expressed gene (DEG) analysis, contig sets without WO97-related sequences were used for mapping, because fungal transcripts in mycorrhizal tissues would certainly be detected as DEGs when comparing with non-mycorrhizal tissues. DEGs were identified using the DESeq2 package (R, v4.3.3) considering a log2 fold change (FC) threshold > |1| and an adjusted p-value (p _adj) < 0.05. The amino acid sequences of the contigs were annotated using DIAMOND BLASTP v2.1.9.163 (Buchfink et al., 2021) against the RefSeq plant database (O’Leary et al., 2016). Functional annotations were further refined using InterProScan v1.8.0 (Jones et al., 2014) for Gene Ontology (GO) term assignment and KofamKOALA v1.3.0 (Aramaki et al., 2020) for the Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology profile. For KEGG pathway identification, the most relevant pathway ID for each K number was retained based on its frequency or the highest association score. Venn analysis was performed using the online tool Jvenn (https://www.bioinformatics.com.cn/static/others/jvenn/, Bardou et al., 2014).

GO and KEGG pathway enrichment analyses were performed using R v4.3.3 with the clusterProfiler and ggplot package. GO and KEGG enrichment significance was assessed using BH multiple testing correction, with a threshold of p _adj ≤ 0.05 to define significantly enriched GO terms and KEGG pathways. The heatmap was generated based on the Python 3.9, Seaborn v0.12.2, and Matplotlib v3.7.2 libraries. PPI networks were constructed with the STRING platform (https://string-db.org/) using Dendrobium catenatum as the reference species. Proteins with sequence identities below 30% were excluded to ensure the accuracy of the analysis. When we uploaded our DEG sequences to STRING, only two orchid species (D. catenatum and Apostasia shenzhenica) returned as matched proteins. Among these, D. catenatum yielded more matched sequences, allowing for broader gene coverage in the network analysis. Therefore, we chose D. catenatum as the reference species to ensure more comprehensive and reliable interaction predictions. The SWEET genes are incorporated as anchor genes to prioritize mycorrhizal interaction-related genes, minimizing the inclusion of irrelevant genes. The confidence threshold was set at 0.2, which presents the minimum level required to include the SWEET genes in the network while maintaining a manageable scale. Unconnected nodes were excluded to enhance visualization and to focus on biologically meaningful interactions. To calculate the node weights, the refined network was further analyzed using Cytoscape v3.10.3 (Shannon et al., 2003) based on dual centrality metrics [betweenness centrality (BC) and degree centrality]. Genes ranked within the top 30% for both metrics were classified as key genes, which represented the central regulatory hubs of DEGs.

The significance of the differences in the proportion of fungal reads at different times after inoculation [2 weeks post inoculation (wpi) and 4 wpi] was evaluated using Tukey–Kramer multiple comparison. DEG detection and GO and KEGG enrichment analyses were conducted with a threshold of p _adj ≤ 0.05 to define significant genes. All statistical analyses were conducted in R v4.3.3.

At both sampling time points, the seedlings in the WO97-R and WO97-N treatments had well-developed roots and were similar in overall morphology ( Supplementary Figure S1B ). There was no visible difference in growth compared to the control (pre-inoculation) during the incubation period. To evaluate the colonization and subsequent mycorrhizal development, we observed peloton formation inside the root tissues of C. macranthos var. rebunense seedlings at 2 and 4 wpi. In the control root tissue of the seedlings before inoculation, no peloton formation was observed. At 2 wpi, the hyphae were distinctly visible and colonized in cells, and peloton formation was induced by both WO97 cultures in orchid root cells ( Figure 1 ; Supplementary Figure S2A ). At 2 wpi, there was no clear difference in the number of cells that were confirmed to be colonized with hyphae between WO97-R-inoculated and WO97-N-inoculated tissues. However, in the WO97-N-inoculated root tissues, few elongated hyphae were observed at 4 wpi, and the pelotons showed structural breakdown and light staining, indicating potential hyphal degradation. On the other hand, elongated hyphae and peloton formation were observed at 4 wpi in the root tissue inoculated with WO97-R, as observed at 2 wpi ( Figure 1 ; Supplementary Figure S2A ).

RNA-seq of 15 samples from five treatment groups (Control, WO97-R-2w, WO97-R-4w, WO97-N-2w, and WO97-N-4w) yielded 93.5%–94.2% of clean reads with a quality score of Q30 (sequencing error rate < 0.1%) from total reads. To examine the abundance of WO97-derived transcripts in each sample, we identified the short reads mapped to the WO97 draft genome sequence as WO97-derived transcripts and calculated their proportions. In the control sample, 0.0% of the reads presented WO97-derived transcripts, whereas the percentage of mapped reads was as high as 14% of the total clean reads in the inoculated treatments (WO97-R-2w, WO97-R-4w, WO97-N-2w, and WO97-N-4w) ( Figure 2 ). In WO97-R-inoculated tissues, the proportion of WO97-derived reads did not significantly differ between 2 and 4 wpi; however, in WO97-N-inoculated tissues, the proportion of WO97-derived reads decreased from 2 to 4 wpi. Since the proportion of fungal-derived reads can be considered as a reference indicator of fungal activity and colonization status, this result suggests that the transcriptional activities of WO97-R and WO97-N were similar in orchid tissues at 2 wpi, but they differed at 4 wpi. Differences in the levels of fungal transcripts (the EF1α gene) among OM fungi-inoculated tissues were also assessed by reverse transcription quantitative polymerase chain reaction (RT-qPCR), which showed similar results to the proportion of fungal reads obtained by RNA-seq ( Supplementary Figure S2B ).

To obtain gene expression data of Cypripedium from the RNA-seq reads, read mapping to the reference and count estimations were conducted. We first used a contig set derived from C. macranthos var. rebunense transcriptomes for a reference, which has been constructed in a previous study (Assembly 4, Kambara et al., 2023); however, considering the mapping rate < 90%, we decided to create a new contig set from the short reads obtained in this study. The contig set constructed by integrating the short reads obtained from the five treatments is shown in Table 1 . The number of contigs was 284,029 without removing WO97-derived short reads, while the number was 260,460 excluding the fungal sequences. The mapping rates remained at approximately 98% for both datasets. De-redundancy performed using CD-HIT-EST reduced approximately 100,000 contigs in both datasets, with improvements in complete single-copy BUSCO scores 32.0%–65.7% (including fungal sequences) or 33.8%–67.3% (excluding fungal sequences). The EvidentialGene program, used to reduce redundancy, further reduced the number of contigs and slightly decreased the mapping rate; it improved the BUSCO score by more than 87% ( Table 1 ). Among the 88,236 contigs constructed from the short reads (excluding fungal reads), 3,388 contigs (3.4%) exhibited high sequence similarity to the CDS sequences of the WO97 draft genome and, hence, were excluded as fungal genes; a final set of 84,848 orchid-derived contigs were obtained for downstream analyses. The total mapping rate of the short reads (excluding fungal sequences) to 84,848 contigs was 87.4%, and the reduction in mapping rates was minimal (0.3%), confirming the reliability of the refined contig set. Among the 84,848 contigs, 23,694 contigs (27.92%) were annotated by the RefSeq plant database, 16,549 contigs (19.50%) were assigned GO terms, and 4,437 contigs (5.23%) were assigned KEGG pathway terms.

The similarity in gene expression patterns among the samples was examined using count data via principal component analysis (PCA). As shown in Figure 3A , the five treatments (Control, WO97-R-2w, WO97-R-4w, WO97-N-2w, and WO97-N-4w) were distributed together with three replicate samples; the first principal component (PC1) explained 16.8% of the variance, and the second principal component (PC2) explained 14.8% of the variance. In the WO97-R-inoculated samples, the sample plots corresponding to 2 and 4 wpi were closely distributed, whereas they were far apart in the WO97-N-inoculated samples.

The number of DEGs derived by comparing gene expression levels between WO97-inoculated and uninoculated controls is shown in Figure 3B . In the WO97-N-2w, WO97-R-2w, and WO97-R-4w groups, the number of upregulated genes was greater than that of the downregulated genes. In contrast, in the WO97-N-4w group, the numbers of upregulated and downregulated genes were almost the same ( Figure 3B ). The finding that WO97-N-4w showed a different trend from the other three treatments was consistent with the PCA results, which is potentially related to the low abundance of fungus-derived transcripts in the WO97-N-4w group. The DEG analysis revealed that 10,139 genes expressed differently in the WO97-inoculated tissues compared with the control (before inoculation). The specificity and commonality of DEGs detected in the four comparisons are shown in Figure 3C . For further functional analyses, 2,003 and 1,926 DEGs were selected; the 2,003 DEGs shared in the four groups were predicted to be involved in the plant response to colonization by an OM fungus. The 1,926 DEGs shared only in the WO97-N-2w, WO97-R-2w, and WO97-R-4w groups were considered to be involved in adaptive mechanisms activated during continuous mycorrhizal interactions because their expression did not significantly change in the WO97-N-4w groups with reduced fungal activity. Among these DEG sets, 51.8% of the 2,003 DEGs and 30.1% of the 1,926 DEGs were annotated using the RefSeq plant database.

GO and KEGG enrichment analyses conducted for studying the functional roles of DEGs across the four comparisons revealed that genes annotated with enriched GO terms comprised 44.0% of the 2003 up- and downregulated DEGs, whereas only 10.0% of DEGs were annotated with enriched KEGG pathway terms. The top five enriched and significant GO BP terms with the lowest adjusted p-values (P _adj) were “systemic acquired resistance”, “defense response to bacterium”, “polysaccharide catabolic process”, “carbohydrate transport”, and “glucan metabolic process” ( Figure 4A ). In the KEGG pathway enrichment, “cutin, suberine and wax biosynthesis”, “platinum drug resistance”, “ABC transporters”, “synaptic vesicle cycle”, and “fructose and mannose metabolism” were detected ( Supplementary Figure S3A ), but no significance was detected due to low P _adj values. Among the 10,139 DEGs, 51 genes were annotated with the significantly enriched 10 GO-BP terms. Heatmap clustering of the 51 genes showed that the DEGs designated in Clusters 2 and 4 were strongly upregulated in OM fungi-colonized tissues when compared to noncolonized controls ( Figure 4B ); these clusters included SWEET sugar transporters, cell wall modification-related proteins, and enzymes such as expansin-like protein (EG45-like domain-containing protein), glucan 1,3-beta-glucosidase, cellulose synthase catalytic subunit, and xyloglucan endotransglucosylase/hydrolase (XTH) ( Supplementary Table S1 ). The pathogenesis-related protein (PR4) and L-type lectin receptor-like kinase (LecRLK) were also included in the upregulated gene clusters. Additionally, the 12 DEGs designated in Cluster 1 were consistently downregulated across all fungi-colonized tissues; the downregulated genes included cellulose synthase catalytic subunit and XTH genes ( Supplementary Table S1 ), suggesting a fine-tuned regulation of the physical properties of the cell wall in response to OM fungal colonization. Regarding genes related to plant hormones, auxin-responsive SAUR genes were detected in Clusters 3 and 4, which are clusters of upregulated genes.

We used the PPI network constructed using STRING to identify key genes crucial for cellular signaling events, among the detected 2,003 DEGs. The PPI network resulted in 830 nodes and 2,160 edges with significant enrichment (p < 1.0e-16) ( Supplementary Figure S4 ). Further analysis via Cytoscape using two centrality parameters (Degree and Betweenness) identified key genes ranking in the top 30% of the two centrality values, and the obtained 239 genes were used for post-PPI GO enrichment analysis. The significantly enriched GO terms (BP, CC, and MF) associated with the 239 genes are shown in Figure 4C ; in the BP terms, “defense response to bacterium”, “polysaccharide catabolic process”, “plant-type cell wall biogenesis”, and “cellulose biosynthetic process” were again detected as enriched terms. Furthermore, it included newly detected terms, such as “establishment or maintenance of cell polarity” and “actin filament bundle assembly”, and cell wall metabolism-related terms, such as “plant-type primary cell wall organization” and “plant-type secondary cell wall biogenesis”. The relative expression levels of the DEGs annotated with the enriched 11 GO-BP terms are shown in Figure 4D and Supplementary Table S2 . Genes related to the polysaccharide catabolic process or carbohydrate metabolic processes were included in both the upregulated cluster (Cluster 3, Figure 4D and Supplementary Table S2 ) and the downregulated cluster (Cluster 1, Figure 4D and Supplementary Table S2 ), indicating their involvement in the adjustment of cellular architecture to support the establishment of the mycorrhizal structure. The defense response-related genes PR-4 (Cluster 2) and LecRLK (Cluster 3) were consistently retained as upregulated genes before and after PPI processing, highlighting their potential role in regulating defense-related responses in mycorrhizal interactions. Besides BP terms, cellular component (CC) terms such as “lysosome”, “extracellular space”, and “cell wall” were retained after PPI filtering, along with MF-related terms including “cellulose synthase activity” and “dipeptidyl-peptidase activity.”

Next, we conducted GO and KEGG enrichment analyses using the 1,926 DEGs shared among the three conditions, except for the WO97-N-4w group. Genes annotated with enriched GO terms comprised 32.2% of the 1,926 DEGs, whereas only 5.3% of them were annotated with enriched KEGG pathway terms. GO-BP terms “leaf development”, “histidine biosynthetic process”, “trehalose metabolism in response to stress”, “endoplasmic reticulum to cytosol auxin transport”, and “NLS-bearing protein import into nucleus” were detected from the 1,926 DEGs ( Supplementary Figure S5 ), but no significance was obtained due to low P _adj values. In the KEGG pathway enrichment, four of the top five pathway terms with the lowest P _adj values detected in 2,003 DEGs (“Platinum drug resistance”, “ABC transporters”, “Synaptic vesicle cycle”, and “Fructose and mannose metabolism”) also appeared in 1,926 DEGs ( Supplementary Figure S3B ), but significance of the enrichment of pathway terms was not detected in the 1,926 DEGs, as in the case of the 2,003 DEGs. Using the 1,926 DEGs, a PPI network with 487 nodes and 1,386 edges was constructed (PPI enrichment, p = 5.22e-05) ( Supplementary Figure S6 ); 143 genes (top 30%) were included in the PPI network with significant centrality values and were used for GO analysis. Through the PPI filtering process, the GO terms were significantly enriched in the 143 genes ( Figure 5A ). Notably, “phosphorylation”- and “signal transduction”-related terms emerged as enriched GO-BP terms. Heatmap clustering and annotation of the 23 genes annotated with enriched 10 GO-BP terms are depicted in Figure 5B and Supplementary Table S2 . Several kinase genes, including serine/threonine-protein kinase (STY46), mitogen-activated protein kinase kinase kinase 18 (MAPKKK18), histidine kinase 1 (HK1), and calcium-dependent protein kinases were abundantly expressed in WO97-N-2w, WO97-R-2w, and WO97-R-4w, highlighting phosphorylation events in intracellular signaling cascades that coordinate host responses during continuous mycorrhizal interactions ( Figure 5B ; Supplementary Table S2 ). Although its expression level was lower than that of the kinase genes in Cluster 1, the auxin response factor ARF9 genes were included in the upregulated genes in Cluster 3 ( Supplementary Table S2 ). In CC terms, “lysosome”, “post-mRNA release spliceosomal complex”, “Prp19 complex”, and “extracellular space” were detected in enrichment after PPI filtering. In MF terms, terms related to various enzyme activities were enriched in addition to the kinases.

To investigate the potential roles of transporters in facilitating nutrient exchange under varying levels of mycorrhizal status and fungal activity, we examined the differential expression of transporter genes included in the 2,003 and 1,926 DEGs among the four conditions. A total of 58 and 31 transporter genes were identified among the sets of 2,003 and 1,926 DEGs, respectively, most of which were upregulated in both sets, highlighting their active involvement in nutrient exchange and cellular adaptation mechanisms during orchid and OM fungi interactions ( Supplementary Figure S7 ). Among the transporter genes, a notable abundance of ABC transporters across both DEG sets suggested their essential role in transferring a variety of substrates such as lipids and phytohormones in response to OM fungus colonization. The sugar transporters SWEET4 and SWEET14 were included in the set of 2,003 DEGs, but only SWEET1 was included in the set of 1,926 DEGs; hence, the functions of these SWEETs may vary depending on mycorrhizal status. Transporter genes for organic nitrogen-containing compounds such as amino acids, peptides, and nucleic acids were upregulated by OM fungal colonization, and several genes were expressed only in 1,926 DEGs. Searching for transporters of inorganic nitrogen compounds, one ammonium transporter and several nitrate transporter genes were found to be upregulated among the 2,003 DEGs. Notably, the expression levels of the transporter genes in the WO97-N-4w group were lower than those in the other three groups, with expression patterns corresponding to changes in fungal activity and nutrient availability across the different colonization stages.

To further elucidate the changes in gene expression during mycorrhizal interactions, we analyzed the DEGs across two time points (2 wpi vs. 4 wpi) ( Figure 6A ) and after inoculation with two fungal cultures (WO97-R vs. WO97-N) ( Figure 6B ). The DEG analysis revealed distinct gene expression patterns associated with temporal progression of fungal colonization and the different fungal subcultures. The temporal comparison (2 wpi vs. 4 wpi in WO97-R or WO97-N) revealed that the number of detected DEGs was significantly higher in the WO97-N-colonized tissues compared to WO97-R-colonized tissues; in the WO97-N-colonized tissues, 4,300 DEGs were upregulated and 1,326 DEGs were downregulated genes at 2 wpi compared with 4 wpi, while 447 DEGs were upregulated and 288 DEGs were downregulated at 2 wpi compared when 4 wpi in the WO97-R-colonized tissues ( Figure 6A ). The DEGs upregulated specifically in WO97-N-colonized tissues (4,061 DEGs) at 2 wpi were enriched significantly in biological process terms “defense response to bacterium”, “mitotic spindle assembly”, “polyketide biosynthetic process”, and “systemic acquired resistance” ( Figure 6A ). A total of 33 genes were associated with the top five GO terms with the lowest P _adj values, including LecRLK, PR4, kinesin-like protein, and EG45-like domain containing protein genes ( Supplementary Table S3 ). The comparative analysis of WO97-R vs. WO97-N at the same time point after the second inoculation showed that very few DEGs (222 DEGs) were detected at 2 wpi, whereas 2,461 DEGs were detected at 4 wpi; most DEGs were significantly less expressed in WO97-N-colonized tissues than in WO97-R-colonized tissues at 4 wpi ( Figure 6B ). The 2,173 genes specifically upregulated in WO97-R-colonized tissues at 4 wpi were associated with defense-related GO terms; LecRLK, PR4, and EG45-like domain containing protein genes, as well as sulfate transporter and clathrin assembly protein genes, were associated with the top five GO enriched terms with the lowest P _adj values ( Supplementary Table S3 ). On the other hand, 1,240 DEGs whose expression levels were lower at 2 wpi than at 4 wpi in WO97-N-colonized tissues (i.e., upregulated genes at 4 wpi compared to 2 wpi in WO97-N-colonized tissues) were enriched in several GO terms including “mRNA destabilization”, “lipid storage”, “photomorphogenesis”, “cellular response to heat”, and “translational elongation”, which suggests a potential reduction in metabolic demand and stress-related responses, coinciding with a decrease in fungal activity and the disappearance of peloton formation.