Using Integrated Multi-Omics to Explore the Differences in the Three Developmental Stages of Thelephora ganbajun Zang

Abstract

Thelephora ganbajun Zang, a rare wild macrofungus, has significant culinary and medicinal value. However, it also has a high cost attributed to its inability to achieve artificial cultivation and its strict environmental requirements. To reveal the intricacies of its development, we conducted a comprehensive analysis of the proteome and metabolome in three pivotal developmental stages: the mycelium, the primordium, and the fruiting body. In our investigation, genes exhibiting various expression levels across multi-omics analyses were identified as potential candidates implicated in growth, development, or metabolic regulation. The aim of this study was to provide a clearer direction for understanding the fundamental metabolic activities and growth stages of this species. Label-free proteomic sequencing revealed a critical juncture in ectomycorrhiza formation, particularly during the transition from the mycelium to the primordium. Secreted proteins, signaling proteins, membrane proteins, and proteins with unidentified functions were rapidly synthesized, with certain amino acids contributing to the synthesis of proteins involved in signaling pathways or hormone precursor substances. In the metabolomics analysis, the classification of secondary metabolites revealed a noteworthy increase in lipid substances and organic acids, contributing to cell activity. The early mycelial development stage exhibited vigorous cell metabolism, contrasting with a decline in cell division activity during fruiting body formation. In our findings, the integration of metabolomic and transcriptomic data highlighted the potential key role of folate biosynthesis in regulating early ectomycorrhiza development. Notably, the expression of alkaline phosphatase and dihydrofolate synthase genes within this pathway was significantly up-regulated in the mycelium and fruiting body stages but down-regulated in the primordium stage. This regulation primarily influences dihydrofolate reductase activity and B vitamin synthesis.

1. Introduction

Thelephora ganbajun Zang is a representative ectomycorrhizal fungus, belonging to the Thelephora genus of the Basidiomycota family [1]. This sizable wild fungus predominantly thrives in Yunnan Province, China, and has garnered widespread public admiration for its medicinal value and distinctive flavor. However, its production remains limited due to challenging growing conditions, its inability to be cultivated artificially, and the ecological impact of human activities. Thelephora ganbajun Zang contains amino acids, proteins, polysaccharides, vitamins, and minerals that contribute to human health [2]. A Fourier-transform infrared (FTIR) spectrometer was employed to characterize the fruiting body of T. ganbajun, and the results revealed that the primary absorption peaks represented proteins and polysaccharides, with the polysaccharides containing α- and β-glycosidic bonds. A distinctive spectral feature of T. ganbajun is a clear band around 1763 cm⁻¹, indicating the presence of oil [3]. Based on T. ganbajun’s unique flavor and the presence of specific chemical components in its metabolites, investigations into the nonvolatile components identified certain novel m-nitrophenyl acetoxylated derivatives [4]. Modifying T. ganbajun’s deep fermentation process could potentially minimize the loss of active components in the fungus [5,6]. Single-nucleotide polymorphism (SNP) and allelic association analysis of T. ganbajun collected from diverse geographical areas in Yunnan Province revealed a common mtDNA recombination in the natural population of this fungus. This provides evidence of mitochondrial heterozygosity in T. ganbajun, a phenomenon not yet reported in populations of the Basidiomycota family [7].

Transcriptome sequencing results indicate that the initiation and activation of the MAPK pathway are crucial in the formation of the fruiting body. Additionally, transcripts with significant differences in the biological process between the mycelium and primordium stages primarily focus on cellular components, molecular functions, and virulence components related to the activation of plant immune defense responses. This is accompanied by an increase in electron carrier activity and antioxidant activity. Similar enrichment is observed from the mycelium to the primordium stage, as well as from the mycelium to the fruiting body stage, suggesting that the functions of ectomycorrhizal formation and the robust development of hyphae are primarily concentrated in the early stage. In the transition from the primordium to the fruiting body, secondary functions with significant proportional differences shift to reproductive development and signal transduction roles [8]. This shift indicates that the changes during this stage primarily involve cell division, material accumulation, and metabolic pathways. Moreover, at this point, the primordium has essentially completed its basic physiological functions and ectomycorrhizal organ development, enabling it to continuously obtain nutrients from the host.

In terms of evolutionary timeline, Thelephora fungi developed independently from saprophytes at an early stage [9,10]. However, when compared with saprophytes from closely related species, numerous functions related to lignin decomposition, organic material utilization, and carbohydrate-activated enzyme synthesis are notably deficient [11]. This deficiency is likely associated with the nutrient life of ectomycorrhizal fungi, in which corresponding gene functions have been lost due to the prolonged evolutionary process, involving a reduction in wood degradation and intracellular nutrient demand. Research indicates that the symbiotic plant Pinus contains a wealth of organic matter, encompassing fatty acids, sugars, terpenes, alcohols, and more [12,13]. The accumulation of these substances potentially contributes to the development of ectomycorrhizal roots, providing essential nutrients for root formation and mycelium growth. However, further exploration is needed to understand the activation of specific signaling pathways and the composition of small-molecule compounds.

At present, various mainstream sequencing analysis technologies, including genomics, transcriptomics, proteomics, and metabolomics, play pivotal roles in research [14]. The integration of multiple analysis techniques is a crucial aspect of omics. By leveraging genomic information as a backdrop, correlations among genes with varying expression levels in different life activities can be discerned [15,16]. Simultaneously, enrichment analysis using the GO [17], KEGG [18], and HMDB [19,20] databases helps to refine the screening of candidate genes, offering reliable data models to support accurate gene localization and function prediction for experiments related to gene function validation. If a gene exhibits a consistent expression trend in all omics contexts, it is highly likely to be a candidate as a housekeeping gene. Conversely, if the expression trend fluctuates between up- and down-regulation across different omics techniques, it is probable that the gene is associated with regulatory functions [21,22]. A multi-omics approach can be used to effectively screen target genes and metabolic pathways. If the genome is used as background information to annotate other omics, it is possible to use gene IDs to represent products at different omics levels to determine their correlation with genes whose information can be looked up in databases such as the NCBI.

In a laboratory setting, the mycelium undergoes gradual degeneration and subsequent demise during mycorrhizal passages, resulting in a shortened transgression time and diminished survival. This phenomenon could be attributed to the absence of specific nutrients and signaling molecules that typically stimulate relevant regulatory genes, or to metabolic dysfunctions arising from the simplified culture environment compared to the field environment. We categorized the life history of T. ganbajun into three pivotal periods: the mycelium, the primordium, and the fruiting body (see Figure 1). Subsequently, we conducted proteomic and metabolomic analyses of samples from these stages. We then carried out a comprehensive joint analysis with the transcriptome, exploring gene expression and the corresponding downstream metabolites. The aim of this approach was to uncover the process by which T. ganbajun produces substances and engages metabolic pathways crucial for ectomycorrhizal and fruiting body formation during its growth. At the same time, delving into the different components of T. ganbajun offers insights into the genes or products associated with nutritionally defective phenotypes, enhancing our understanding and providing a valuable reference for the harvesting and protection of strains in the field.

2. Materials and Methods

2.1. Sample Collection and Extraction

Samples of Thelephora ganbajun Zang were gathered from 3 pivotal developmental stages: the mycelium, the primordium, and the fruiting body. Primordium and fruiting body samples were specifically obtained from Yiliang County, Kunming City, Yunnan Province, China, and were extracted from the roots of a symbiotic plant, Pinus yunnanensis. The mycelium was collected from the same locale; the liquid was cultured in the laboratory and stored at a constant temperature. The mycelium was poured from the liquid culture into a strainer and left to stand until the liquid was completely filtered out, and the mycelium was wrapped with sterile filter paper to absorb the remaining water. The mycelium was collected in a 15 mL centrifuge tube and washed 3 times with 10 mL PBS; then, the mycelium (400–1000 g) was resuspended with 1 mL PBS and transferred to a 1.5 mL centrifuge tube for centrifugation (5–10 min at 4 °C). Then, the supernatant was aspirated with a lance tip and the wet weight of the mycelium was recorded; after that, the samples were snap-frozen in liquid nitrogen and stored at −80 °C.

2.2. Protein Extraction

To prepare proteomic sequencing of mycelium, primordium, and fruiting body samples (labeled GP, GPP, and GSP, respectively), we took 0.5–1 g samples from each, each of which comprised 4 biological replicates. The appropriate amount of liquid nitrogen for each sample was ground into powder, and 0.5 mL lysis buffer (lysis solution: protease inhibitor 50:1) was added, dissolved by vortexing, and mixed, followed by ultrasonic crushing for 1 s, then stopping for 1 s, and repeated for a total of 120 s. After that, the samples were put into 1.5 mL BP tubes and centrifuged at 14,000× g for 20 min; then, the supernatant was taken out and transferred to a new tube, 10 µL was kept for quantitative analysis, and the rest was frozen at −80 °C and stored. Protein concentration was determined using the Bradford method [23]. Samples were diluted to fall within the standard curve range. Diluted samples and standards of 10 µL each were reacted with 300 µL of protein quantitative analysis dye in the dark for 10 min. The absorbance of the standard and the sample was measured at 595 nm with an enzyme marker, a standard curve was plotted based on the absorbance and concentration of the standard in each tube, and then the concentration of the sample was calculated. Sample proteins were subjected to electrophoresis via SDS-PAGE with a 1% agarose gel, stained with Coomassie brilliant blue for 30 min, and decolorized until the background was clear (Table S1, Figure S1).

2.3. Mass Spectrometry

Capillary high-performance liquid chromatography (Thermo Scientific EASY-nLC 1000 System, Nano HPLC, Waltham, MA, USA) was performed using the following equipment: mobile phase A: 100% ultrapure water, 0.1% formic acid; mobile phase B: 100% acetonitrile, 0.1% formic acid; (formic acid no. 56302, methanol no. 14262); pre-column: Acclaim PepMap100 column (2 cm × 100 μm, C18, 3 μm); column: EASY-Spray column (12 cm × 150 μm, C18, 1.9 μm); injection vials (Thermo, 11190533); cap (Thermo, 11150635); spray needle (Thermo, PN: ES542). Experimental procedure: The lyophilized components were reconstituted with 20 μL of 2% methanol and 0.1% formic acid, followed by centrifugation at 12,000 rpm for 10 min, and then the supernatant was aspirated. The volume of sample was 10 μL, and then auto-sampling was performed. The separation flow rate of the loading pump was 300 nL/min for 15 min, and the separation gradient is shown in Table S2A. Mass Spectrometry System (Thermo, Model: Orbitrap Fusion): Ion source parameters: spray voltage 2.2 kV; capillary temperature 350 °C; ion source: NSI. Full MS: resolution: 120,000 FWHM; full-scan AGC target: 2 × 10⁵; full scan max IT: 50 ms; scan range: 350–1550 m/z. dd-MS2: resolution: 30,000 FWHM; AGC target: 5 × 10⁴; max IT: 50 ms; fragmentation Methods: HCD.

In this experiment, label-free Data-Dependent Acquisition (DDA) quantitative proteomic technology was employed for mass spectrometry [24]. Proteomics studies using a DDA scanning mode can obtain primary mass spectra and fragmentation information. However, this mode does not collect fragment ions from all parent ions; it collects the corresponding fragment ions by selecting specific parent ions into the collision cell. The limited DDA acquisition rate will cause a loss of quantification of low-abundance peptides [25,26]. Fundamentally, the mass spectrometry method involves matching experimental spectra with theoretical spectra in the database to obtain potential peptide sequences. Maxquant (v2.2.0.0; https://maxquant.net/, accessed on 10 September 2021) [27] was utilized to construct the transcriptome database and process the raw mass spectrometry files generated via label-free quantification (LFQ). We uploaded all raw data and set the grouping. On the Global Parameters screen, we clicked Sequences to add database files to be checked, set all parameters, and began searching the database. The LFQ algorithm first identified the peptide signals in the LCMS data, and then completed the qualitative data by performing a database search for all peptide signals in MS2 using the built-in Andromeda algorithm, with the specific parameter settings shown in Table S2B. The MS2 data were analyzed using Proteome Discoverer (v1.4; Thermo Fisher, Waltham, MA, USA) to quantify the peak intensity values of the peptide reporter ions (Table S2C, Figure S2).

2.4. Metabolite Extraction

For metabolomics sequencing, samples of the mycelium, the primordium, and the fruiting body (labeled GM, GPM, and GSM) were prepared, each consisting of 6 biological replicates. For this, 50 mg of each sample was weighed into an EP tube, mixed with 1000 µL of extraction solution (methanol–acetonitrile–water v/v/v = 2:2:1, internal standard concentration 2 mg/L) with internal standard (1000:2), then vortexed and mixed for 30 s. Ceramic beads were added, and the samples were processed in a 45 Hz grinder for 10 min, treated ultrasonically in an ice–water bath for 10 min, and left to stand for 1 h at −20 °C. Samples were centrifuged at 12,000 rpm, 4 °C, for 15 min; then, 500 µL of supernatant was carefully removed and placed into a new EP tube, and the extract was dried in a vacuum concentrator. Then, 160 µL of extraction solution (acetonitrile–water v/v = 1:1) was added to the dried metabolite to reconstitute the extract; this was vortexed for 30 s and sonicated in an ice–water bath for 10 min, followed by centrifuging again at 12,000 rpm, 4 °C, for 15 min. Finally, 120 µL of supernatant was carefully removed and added to a 2 mL injection bottle, and 10 µL of each sample was mixed into a QC sample for testing.

2.5. Experimental Instrumentation Platform and Reagents

The experimental instruments and reagents used in this study are detailed in Table S3. The liquid chromatography–mass spectrometry (LC-MS) system comprised an Acquity I-Class PLUS ultra-high-performance liquid phase (UHPLC) with an Xevo G2-XS QT high-resolution mass spectrometer from Waters, Milford, MA, USA. The chromatographic columns were Acquity UPLC HSS T3 (1.8 μm 2.1 × 100 mm) from Waters, USA. In both positive and negative ion modes, the mobile phases were as follows: A: 0.1% formic acid aqueous solution and B: 0.1% formic acid acetonitrile. The injection volume was maintained at 1 µL.

2.6. Data Processing and Analysis

Data processing was carried out using MassLynx (v4.2.0, Waters; https://www.medicalexpo.com.cn/). The peaks were extracted and aligned using MSConvertGUI (v3.0; https://proteowizard.sourceforge.io/download.html). The compounds were identified using the METLIN database (https://metlin.scripps.edu), including the identification of theoretical fragments. The deviation of mass counts was consistently within 100 ppm. After that, we used R to graph the processed data.

To discern differences in protein expression among samples from different groups, we employed a t-test (IBM SPSS, Armonk, NY, USA, https://www.ibm.com/products/spss-statistics), with screening conditions set at a significance value of false discovery rate (FDR), and the log2 logarithm of the quantitative ratio, log2FC (where FC is the fold change), between the two groups (up- or down-regulated) exceeding 2 (based on maximizing variation and data-driven approaches).

The correlations between all genes and metabolites were computed for each difference grouping using the Pearson correlation method (IBM SPSS). Prior to correlation calculation, the data underwent Z-value transformation for preprocessing. Subsequently, the results were filtered based on the correlation coefficient (CC) and p-value [28], with screening thresholds of |CC| > 0.80 and CCP < 0.05.

2.7. Conjoint Proteome and Metabolome Analysis

When exploring correlations between different omics techniques, the standard approach involves using KEGG co-enrichment analysis to identify shared metabolic pathways. In the integrated analysis of transcriptomic and metabolomic data, our focus was on the statistical examination of correlated pathways. Next, we filtered out genes and metabolites that might participate in the same pathway based on the extent of their correlation. Employing the ID conversion method, we transformed gene and metabolite IDs into the KEGG format and mapped them onto a metabolic pathway map. By integrating genes and metabolites with similar expression trends, we aimed to illustrate the connections between the two datasets more intuitively. During KEGG annotation, genes and metabolic pathways with a p-value < 0.05 were prioritized. This analytical approach identified pathways relevant to the study’s objectives, effectively excluding the influence of irrelevant variables [29].

Given the requirement for large sample data with multiple repetitive samplings in co-analysis, as well as the insufficient number of annotated differential proteins in the proteomic analysis, we utilized transcriptomic data from the mycelium, primordium, and fruiting body stages of T. ganbajun (corresponding to the developmental phases mentioned above). These datasets were merged and analyzed alongside the metabolomic data. The clean data from all samples in the transcriptomics analysis amounted to 5.76 Gb, with Q30 and GC content tests revealing that the percentage of bases with Q Phred > 30 in each sample exceeded 93.15%. Then, we used a multi-omics analysis platform (https://international.biocloud.net/zh/software/medical/detail/8a8300b08697667701869c3cbb920000) on the BMKCloud platform (www.biocloud.net) to generate interaction networks.

3. Results

3.1. Proteomic Test Results and Sample Analysis

Protein expression exhibits spatio-temporal specificity, and proteins displaying noteworthy variations in expression levels under distinct conditions are referred to as differentially expressed proteins (DEPs). In this study, the collection of proteins derived from the analysis of differential expression was labeled as “A vs. B”. Specifically, T. ganbajun’s mycelium was designated as GP, its primordium as GPP, and its fruiting body as GSP. The DEPs were further categorized as up- or down-regulated proteins (

Table 1. Statistics of differential protein quantity.

Diff Set	Total Number	Up- Regulated	Down- Regulated	Up-Regulated Proportion (%)	Down-Regulated Proportion (%)
GPP_vs._GSP	44	24	20	54.55	45.45
GP_vs._GPP	520	372	148	71.54	28.46
GP_vs._GSP	500	348	152	69.60	30.40

) based on the relative expression levels observed between the two groups of samples.

Significant distinctions between the two groups of DEPs are represented by volcano diagrams (Figure 2A) [30]. The volcano plot of differential protein expression across the three developmental phases revealed a pronounced difference in protein expression enrichment between the early mycelium stage and the mature primordium and fruiting body stages. Notably, the enrichment of DEPs between primordium and fruiting body was significantly diminished. To enhance the accuracy of evaluating DEPs, we repeated tests of the same samples and established correlation coefficients within the experimental design of differential expression [31]. Clustered heatmaps were used to represent the relationships of all sample proteins (Figure 2B). Principal component analysis (PCA) [32,33] was utilized to score each set group (Figure 2C). The resulting map illustrates that samples within each group were consistently categorized at the same level without aberrations. The observable compositional differences between distinct groups suggest that the principal component variances among biological replicates within the groups were minimal, indicating that the samples were in a normal condition. A Venn diagram was constructed to represent the sets of differentially expressed proteins for GP, GPP, and GSP (Figure 2D). Notably, 24 DEPs were common to both GP vs. GPP and GPP vs. GSP, while another 24 DEPs were shared between GPP vs. GSP and GP vs. GSP. Additionally, 329 DEPs were expressed in both GP vs. GPP and GP vs. GSP, and 11 DEPs were enriched in all three groups.

3.2. COG and KEGG Enrichment Analysis of DEPs

From the COG database annotation results (Figure 3A), we observed that during the developmental transition from GP to GPP, the highest number of differential proteins are enriched in post-translational modifications and protein translocation, followed by general functions and energy production transitions. Functions related to RNA processing modification, chromatin structure, prophage, transposons, and nuclear structure did not show significant enrichment. A similar trend is evident in the transition from GP to GSP, where we observed the highest number of differential proteins. It is hypothesized that this function plays a pivotal role in sclerotium formation through mycelium coalescence, nutrient accumulation, and functional proteins. The biosynthesis function was found to be weakened between the transition from the GPP to the GSP stage, indicating a cessation of biological processes related to fungal morphology construction and ectomycorrhizal formation. Instead, a considerable amount of nutrient conversion into energy was needed to support growth. Proteins involved in specific functions were not enriched during this stage. In addition to protein translocation, differentially expressed proteins increased in secondary metabolite biosynthesis, transport, and catabolism; energy production conversion; and amino acid and carbohydrate transport metabolism. During the reproductive growth of the fruiting body, intensified mitosis and a surge in expressed products were involved in cell division, promoting the construction of sporophores and further development.

Protein translocation and post-translational modifications occur throughout all stages of T. ganbajun’s lifecycle, from basic vegetative growth to reproductive growth. These processes involve a multitude of small-molecule secreted proteins, signaling proteins, and other proteins with unknown functions, facilitating the essential conditions for morphological construction, cell division, differentiation, signal transduction, and the induction of ectomycorrhizae production by plant hosts. In the reproductive growth stage, the focus shifts to the amino acid and carbohydrate metabolic processes. Ectomycorrhizal fungi, which are unable to acquire nutrients through frequent catabolism like saprophytic fungi, generate their own energy to fuel energy-consuming pathways, such as aerobic respiration, fatty acid metabolism, and the TCA cycle. This is crucial to complete the transformation from the primordium to the fruiting body in a relatively short period, which is achieved through cell division promotion and cell cycle regulation.

In the KEGG enrichment results (Figure 3B, Table S4), the most distinct pathway for differential proteins was found to be antibiotic biosynthesis between the mycelium and primordium stages (GP vs. GPP). Annotated DEPs included transaldolase, chorismate synthase, ornithine carbamoyltransferase, arginase, and methylsterol monooxygenase, suggesting a connection to the mycelium’s early competition with other fungi and bacteria for an ecological niche in the host plant, as well as defense against pathogen invasion. Following closely were amino acid biosynthesis and carbon metabolism, featuring DEPs including acetylornithine aminotransferase, dihydroxy-acid dehydratase, citrate synthase, malate dehydrogenase, and aminomethyl transferase. RNA transporters and ribosomes associated with amino acid synthesis exhibit more differential proteins, accompanied by cellular respiration and oxidative phosphorylation to generate substantial energy in this process. The pathway enrichment between the mycelium and fruiting body stages (GP vs. GSP) mirrored that of the mycelium to the primordium stage, while in the primordium to fruiting body stage (GPP vs. GSP), a marked reduction across all types of pathways was observed. The differences mainly center on protein and lipid synthesis, purine pyrimidine metabolism, and glycosylation, all intricately linked to the chromosome cycle and cell division.

In light of the evidence, it can be deduced that the initiation of ectomycorrhizal roots predominantly occurs in the early stage. The activation of numerous pathways involved in sugar and lipid synthesis lays the groundwork for the rapid expansion of the mycelium. Concurrently, through aerobic respiration, the mycelium extends into the intercellular matrix, establishing a symbiotic relationship with the host. The synthesis of various antibiotic proteins further suggests that the intricate environment encountered in the early stage and the pressure from specific environmental stresses could be influential factors affecting mycorrhizal development. Compounds generated by enzymatic action, including methyl butyrate, folic acid, and terpenoids, among others, function as precursors for ascorbic acid and acetate metabolism. These products are likely associated with the formation of certain signaling molecules in T. ganbajun, inducing the development of ectotrophic mycorrhiza, and facilitating nutrient acquisition from plants.

3.3. Metabolomic Test Results and Sample Analysis

In the metabolomic analysis, we designated T. ganbajun’s mycelium as GM, its primordium as GPM, and its fruiting body as GSM. While differences within groups were not pronounced, they were evident. PCA scores (Figure 4C) were then determined for each group, revealing that samples in the GM group were more discrete. However, samples in the GPM and GSM groups shared similar principal components, with each group being consistently categorized at the same level and devoid of anomalies. The statistics of the number of differentially expressed metabolites (DEMs) in each group are presented in

Table 2. Statistics of DEMs in different growth stages.

Groups	All Diff	Down-Regulated	Up- Regulated	Down-Regulated Proportion (%)	Up-Regulated Proportion (%)
GM_vs._GPM	3410	1129	2281	33.11	66.89
GM_vs._GSM	2925	1180	1745	40.34	59.66
GPM_vs._GSM	2235	1569	666	70.20	29.80

. Volcano plots were constructed to visually depict the differences in metabolite expression levels between two groups (Figure 4A), along with the statistical significance of these differences. Notably, significant differences in metabolite expression were observed in all three developmental stages. Comparative analyses, including mycelium vs. primordium and mycelium vs. fruiting body, revealed more DEMs, with the majority of metabolites being up-regulated and displaying a dispersed pattern, mirroring the findings in the proteomic enrichment. Conversely, the number of DEMs during the primordium and fruiting body stages was significantly lower, and they were more concentrated compared to the previous stage, suggesting a stabilization of metabolic levels. In the heatmap of all metabolites (Figure 4B), samples from the same group are clustered, and samples from different groups are significantly differentiated. The redder the color, the higher the metabolite expression, and the greener the color, the lower the metabolites expression. GM during the nutrient growth period and GPM and GSM during the reproductive growth stage were highly differentiated, with correlation coefficients for the biological replicates in the same stage of 0.9, meaning the different stages were clearly separated. According to the Venn diagram of the enrichment results (Figure 4D), 1695 DEMs were common to GM vs. GPM and GM vs. GSM, 1181 DEMs were common to GM vs. GPM and GPM vs. GSM, 991 DEMs were common to GM vs. GSM and GPM vs. GSM, and 805 DEMs were enriched in all three groups.

3.4. Multiplicative Analysis of Variance and Orthogonal Partial Least Squares–Discriminant Analysis (OPLS-DA)

Following qualitative and quantitative analyses of the detected metabolites, the quantitative information on metabolite changes within each subgroup could be effectively reflected through the log2-processed multiplicity of differences. This highlights the variability in metabolite quantities across sample stages (Figure 5A). Metabolomics data are characterized by high dimensionality (numerous types of metabolites) and small samples (limited scale). They encompass both differential variables related to categorical variables and a multitude of undifferentiated variables that may have correlations. Analyzing these high-dimensional variables with the PCA model often poses challenges. Inter-sample variance can be reflected in the principal components, which can fail to accurately reflect the specific numerical differences between samples. To address this, we employed the OPLS-DA statistical method to further process the PCA results.

OPLS-DA filters out orthogonal variables in metabolites not related to actual variables. It then cross-analyzes non-orthogonal variables and the remaining irrelevant variables that align with the expected results. This process omits low-scoring samples, enhancing the data prediction reliability. The evaluation model’s predictive parameters, R2X, R2Y, and Q2, play pivotal roles (Figure 5B). R2X and R2Y represent the explanation rate of the constructed model for the X and Y matrices, respectively. Here, the X matrix denotes the model input, representing the metabolite quantification matrix, while the Y matrix signifies the model output, representing the sample grouping matrix. Q2 represents the predictive ability of the model (Figure 5C). The closer R2Y and Q2 are to 1, the more stable and reliable the model is. Hence, the test model successfully distinguished between different sample subgroups based on differential metabolite expression. This model proved useful for screening effective models (Q2 > 0.5) and excellent models (Q2 > 0.9) across all sample groups.

3.5. Annotation Results of Differentially Expressed Metabolites

Drawing upon the metabolite categorization data from the HMDB database, our investigation entailed the statistical mapping of annotated differential metabolisms (Figure 6A). This analysis shed light on the categorization of the specific biological processes of T. ganbajun metabolites. Metabolite fractions exhibiting a significant number of differences (>100) included carboxylic acids and derivatives, fatty acyls, organooxygen compounds, and prenol lipids, with organic acids constituting over 19% and fatty acyls making up 17%. These lipids, crucial for protein recognition and signal transduction at the cell membrane, are likely associated with the formation of the fruiting body, which involves extensive cell division and heightened lipid metabolism.

KEGG annotation (Figure 6B) revealed metabolites exhibiting significant quantitative differences. Notably, there was a substantial enrichment in the biosynthesis of plant secondary metabolites (9.96%), followed closely by the purine metabolism and pyrimidine metabolism, accounting for 7.93% and 5.49%, respectively. Additionally, there were high percentages of certain amino acids and biosynthesized compounds among all differentially expressed metabolites (DEMs). Beyond those merely supporting growth, complex compounds were formed that are intricately involved in the biological processes of cell differentiation and reproductive development. By leveraging ClusterProfiler (https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html) for the KEGG annotation of differential metabolites (Figure 6C) and generating a network diagram (Figure 6D), it was found that the differential metabolic pathways were notably enriched in the GM vs. GPM phase, which included flavonoid and folate biosynthesis. This indicates a resemblance to the early development of this species. Flavonoids, which are mostly found in plants and are mostly bound to sugars in the form of glycosides or carbon sugar groups, synthesize phytochrome [34]. Folate is involved in vitamin synthesis, providing the necessary substances for producing enzymes and immune-enhancing effects, and offering sufficient trace elements for cells to fulfill all basic life activities.

In GM vs. GSM, pathways including folate biosynthesis, porphyrin metabolism, cyanoamino acid metabolism, and biosynthesis of type II polyketide products were found to be enriched. Porphyrins, complex nitrogen-containing compounds derived from changes in plant and animal pigments (chlorophyll or hemoglobin), can mimic the entry of signaling molecules into the cell membrane to participate in catalyzed peroxidase and oxidase reactions [35]. Similarly, there was significant enrichment of fatty acid degradation, which provides conditions for rapid cellular productivity. The significantly enriched metabolic pathways in GPM vs. GSM include 12-, 14-, and 16-member macrolide biosynthesis. This class of substances acts as a natural antibiotic, primarily enhancing the bacteriostatic function of organisms [36] and the antimicrobial properties of plants [37]. Arginine biosynthesis is correlated to glutamate metabolism; these pathways form a foundation for amino acid metabolism and cycling, acting as precursors for complex protein synthesis, with glutathione playing a crucial role as a multifunctional tripeptide involved in numerous catalytic reactions [38,39]. This suggests a decrease in the rate of cell division in the post-primordium stage, with fungal development leaning more toward the formation of functional organs. Therefore, late development primarily focuses on intra- and extracellular signaling, amino acid metabolism, biosynthesis of secondary metabolites, and diminished lipid metabolic activity.

3.6. Conjoint Analysis of Transcriptomes and Metabolomes

We opted for two models that satisfied the criteria of intergroup association analysis, i.e., maximizing differences in components: the mycelium vs. primordium (GM vs. GPM) and primordium vs. fruiting body (GPM vs. GSM) phases. Our analysis delved into the correlation and expression trend of each component. The outcomes of this screening, coupled with the differential folds, were visually represented in a nine-quadrant plot (note: the plot could not be generated when the metabolite differential folds threshold was 1) (Figure 7A,B) [40,41]. Given that the up- and down-regulated expressions of genes and metabolites did not meet the FC screening conditions individually, we focused on six quadrants to depict the fundamental expression trends of genes and metabolites. Quadrants 3 and 4 indicated that gene and metabolite expression trends were consistent, implying potentially positive regulation of the metabolite by the gene. Conversely, quadrants 1 and 6 indicated opposite expression trends, signifying a negative correlation between gene and metabolite. Quadrants 2 and 5 suggested that the gene remained unchanged while the metabolite underwent up- or down-regulation.

The horizontal coordinates indicate the transcriptome’s differential multiplicity (represented as log2), while the vertical coordinates indicate the differential multiplicity of the metabolome (also represented as log2). The analysis focused on trends in differential genes and metabolites. The dashed lines along the horizontal and vertical coordinates indicate the differential fold changes of the transcriptome and metabolome, respectively. Noteworthy differences in genes and metabolites are denoted outside the threshold lines, whereas non-significant differences are marked within these lines. Each fitted line corresponds to a gene or metabolite, with red lines indicating genes and metabolites that are significantly different (exhibiting either consistent or opposite expression trends), and green lines indicating metabolites that are significantly different in one of the omics but not the others.

The correlation analysis and quadrant labeling provided useful insights. In the GM vs. GPM stage, quadrants 1 and 2 suggest a higher expression abundance of metabolites than genes, indicating the up-regulation of metabolites. Quadrants 3 and 4 suggest that the gene expression trend aligns with metabolites, potentially signifying a positive regulatory effect of genes on metabolites. Quadrants 5 and 6 indicate that metabolites are expressed less abundantly than genes, corresponding to the down-regulation of metabolites. Notably, quadrants 1 and 3 have the highest number of differentially abundant metabolites (DAMs) and differentially expressed genes (DEGs), exhibiting a negative correlation in quadrant 1 and a positive correlation in quadrant 3. Moving to the GPM vs. GSM stage, the basic expression trends of genes and metabolites remain consistent with the previous stage, although the numbers of DAMs and DEGs in these quadrants are reduced. The correlation analysis shows significant differences in the genes and metabolic pathways regulating biosynthesis between different growth stages.

When conducting KEGG annotation (Figure 8A), multiple pathways are often annotated for both differential metabolites and differential genes. In such instances, prioritizing genes and metabolic pathways with a p-value < 0.05 for analysis can allow for the swift identification of pathways related to the study’s objectives. In the two-layer omics clustering, the common pathway with a p-value < 0.05 in the GM vs. GPM stage for both transcriptome and metabolome was found for the biosynthesis of amino acids. In the GPM vs. GSM stage, transcriptome pathways with a p-value < 0.05 included arginine and proline, ascorbate and aldarate, β-alanine, and biotin metabolism. Metabolomic pathways with a p-value < 0.05 included arginine and amino acids biosynthesis. Although no significant common metabolic pathways were enriched at this stage, folate biosynthesis and histidine metabolism showed high correlation with the transcriptome. KEGG co-enrichment results revealed common pathways in the transcriptome and metabolome at different stages of T. ganbajun development. This may offer insights into the positive/negative regulation of metabolites by genes, emphasizing that metabolites are not necessarily directly correlated with the gene expression profile.

Simultaneous KEGG annotation of the same grouped transcriptome unigenes and differential metabolites facilitated a deeper understanding of the relationship between genes and metabolites. Screening all enriched KEGG pathways revealed that, beyond major pathways related to growth and development (such as the MAPK signaling and sugar and lipid metabolism pathways), the folate synthesis pathway was found to play a crucial role in ectomycorrhizal root formation and substrate development. In this pathway, alkaline phosphatase (ALK) generates 7,8-dihydroneopterin-3′-triphosphate, which, through a series of reactions, ultimately synthesizes dihydrofolate (DHF). DHF can be further converted to folate or enter the folate–carbon pathway (one carbon pool by folate) (Figure 8C). In this pathway, alkaline phosphatase (PhoA) and dihydrofolate synthase (DHFS) genes were down-regulated in the GM vs. GPM stage. PhoA primarily regulates alkaline phosphatase synthesis, and weakened DHFS activity results in decreased folate and dihydrofolate, leading to reduced purine content. Both genes were up-regulated in the GPM vs. GSM stage (Figure 8B, Table S5).

4. Discussion

As a typical ectomycorrhizal fungus, Thelephora ganbajun Zang has a significantly lower capacity for breaking down plants to obtain nutrients compared to saprophytic fungi. This characteristic is intricately linked to its symbiotic behavior with the plant host. T. ganbajun has undergone evolutionary adaptations to coexist harmoniously with its host plants. This ensures the formation of mycelium in the ectomycorrhiza with the plant root system without causing damage to the cellular structure, thereby enabling the fungus to extract nutrients from the host over an extended period. In the laboratory, our observations revealed a limitation in the direct development of mycelium into ascospores. To explore this phenomenon, we delved into the specific regulatory processes at various developmental stages of T. ganbajun. Employing several omics approaches, we focused our investigation on the role of the folate metabolism pathway in regulating different growth stages. This investigation is vital, as it may have implications for the viability of mycelium. These findings illuminate the intricate dance of co-evolution between Thelephora ganbajun Zang and its host plants. The mycelium’s inability to directly transition into ascospores in non-field environments prompted a thorough exploration of the nuanced regulatory processes at play during various developmental stages. The spotlight on the folate metabolism pathway reveals its pivotal role in shaping the various growth stages of T. ganbajun, influencing the mycelium’s viability and, consequently, its overall ecological impact.

In mass spectrometry analysis, the area normalization method in mass spectrometry quantitative detection is usually suitable for samples in which all components flow out of the column and all peaks appear on the chromatogram, and the content is calculated with the normalization method. The percentage of each component is calculated by normalizing the total number of peaks based on the ratio of the peak area of each component to the total area. This method is convenient and quick, does not require a precise injection volume, and is a more common method of mass spectrometry calculation; however, the error is larger when measuring low content, especially trace impurities. This experiment was only a preliminary exploration of the content of components in T. ganbajun with higher percentages, so the chosen method was used to improve the detection efficiency. The normalization method is essential in mass spectrometry detection, as it reduces the influence of different models of instruments on the signal intensity, avoids high concentration variation, which would result in nonlinear intensities, and removes experimental anomaly [42].

In the proteomics findings, the most prominent pathways with differential protein expression, enriched from mycelium to primordium, were associated with antibiotic biosynthesis, primarily linked to defense against pathogen invasion. Subsequently, amino acid biosynthesis, carbon metabolism, and lipid synthesis pathways were evident: these play pivotal roles in cellular energy production, cell division, and the initiation of signaling pathways. Consequently, the establishment of ectotrophic mycorrhiza primarily unfolds in the early stages. Here, the activation of numerous pathways, including sugar and lipid synthesis, lays a foundation for cellular energy production. Simultaneously, through aerobic respiration, it propels the mycelium into the interstitial spaces of cells, establishing a symbiotic relationship with the host. Notable proteins such as endoplasmic reticulum proteins, small secretory molecules, membrane proteins, and enzymes involved in carbohydrate and cellulose degradation contribute to nutrient accumulation and root morphology. These proteins actively participate in nutrient accumulation and the construction of mycorrhizae.

In the metabolomics investigation, the KEGG annotation results showed that differential metabolic pathways during the mycelium vs. primordium period were primarily dominated by amino acid metabolic pathways and the biosynthesis of sugar and lipid metabolism precursors. These pathways in amino-acid-based metabolism supply precursor substances for complex macromolecules and organic compounds, emphasizing the substantial energy demand crucial for early developmental stages. Metabolic pathways enriched in mycelium vs. fruiting body, with a p-value < 0.05, included folate biosynthesis and porphyrin metabolism, both of which are involved in DNA synthesis and repair and have signaling roles. In the primordium vs. fruiting body phase, the key metabolic pathways with a p-value < 0.05 were 12-, 14-, and 16-member macrolides biosynthesis; this phase enhances fungal immunity as an endogenous antimicrobial agent, reducing the dependence of fungi on plant hosts. The remaining pathways included arginine biosynthesis. Notably, during the late development of the fruiting body, there was an enhancement in signal transduction and hormone analog synthesis pathways, likely involved in interactive processes with the host plant, while there was a decline in energy metabolism activity.

We performed a combined analysis of the metabolome and transcriptome to uncover specific substances and unique processes influencing the ectomycorrhizal development of T. ganbajun throughout its life cycle. Co-analysis results suggest that functions such as amino acid metabolism and vitamin biosynthesis may constitute key factors influencing the morphological construction of T. ganbajun during the transition from the mycelium to the primordium. The folate synthesis pathway plays a pivotal role in the synthesis of B vitamins and serves as a precursor for purines and pyrimidines, participating in both the ascorbate pathway and the methylation process [43]. Folic acid’s fundamental physiological functions are achieved through its reduced product, tetrahydrofolate (THF), and methylated derivatives of THF [44]. These derivatives primarily act as methyl donors for intracellular methylation reactions and contribute to de novo deoxyribonucleic acid synthesis [45]. In organisms, folic acid undergoes reduction by ascorbic acid, NADPH, and dihydrofolate reductase, resulting in THF production [46]. Subsequently, methylenetetrahydrofolate reductase (MTHFR) generates N-methyltetrahydrofolate, with 5-methyltetrahydrofolate (5-MTHF) being the most prevalent active form. This active form is transported to peripheral tissues for cellular metabolism [47], or it may undergo methylation to produce homocysteine by methionine synthetic reductase (MTRR) [48], ultimately lowering homocysteine levels [49]. In the biosynthetic pathway, nearly all plants, along with the majority of fungi and bacteria, synthesize folic acid [50,51]. This water-soluble B vitamin is crucial for the synthesis of thymine nucleotides and purines, as well as for methionine metabolism in fungi [50]. It serves as a vital pathway for the synthesis of micronutrients. A comprehensive analysis of both transcriptomic and metabolomic KEGG enrichment revealed a weak presence of the folate biosynthesis pathway in the GM vs. GPM phase. However, this pathway was significantly enriched in the GPM vs. GSM phase. This suggests that adequate access to this vitamin might be constrained under laboratory culture conditions. Nonetheless, the nutrient could potentially be acquired from the plant host or through self-synthesis in a field environment.

5. Conclusions

In our investigation into the development of T. ganbajun, a representative ectomycorrhizal fungus, and its potential connections with the external environment, we employed proteomic and metabolomic approaches. Our findings indicate significant up-regulation of PhoA and DHFS genes in the mycelium and fruiting body stages but down-regulation in the primordium stage. Furthermore, tetrahydrofolate and its derivatives exhibited an increase in the mycelium and primordium stages. These increases are significant for T. ganbajun’s development, contributing to the methyl cycle and participating in carbon cycling by acting as donors and acceptors of one-carbon units. Additionally, it exhibits certain antioxidant stress effects. This pathway is instrumental in promoting the formation of ectotrophic mycorrhiza, facilitating the construction of a mycorrhizal network, evading host defense mechanisms, and compensating, to some extent, for the biodegradability shortcomings of ectomycorrhizal fungi. The absence of specific nutrients in the laboratory culture environments has led to a diminished role for this pathway, consequently impeding mycelial clustering and restricting overall growth.