Graphical abstract Workflows of transcriptome and proteome analyses across the differentially developmental stages of C. guangdongensis. My, Mycelial stages; P, Primordial stage; YFB, Young fruiting body stage; DFB, Developed fruiting body stage; MFB, Mature fruiting body stage. graphic file with name ga1.jpg [35]Open in a new tab Keywords: Cordyceps; Omics analysis; Correlation analysis; Bioactive ingredient biosynthesis-related enzymes, growth and development Abstract Tolypocladium guangdongense has a similar metabolite profile to Ophiocordyceps sinensis, a highly regarded fungus used for traditional Chinese medicine with high nutritional and medicinal value. Although the genome sequence of T. guangdongense has been reported, relatively little is known about the regulatory networks for fruiting body development and about the metabolite biosynthesis pathways. In order to address this, an analysis of transcriptome and proteome at differential developmental stages of T. guangdongense was performed. In total, 9076 genes were found to be expressed and 2040 proteins were identified. There were a large number of genes that were significantly differentially expressed between the mycelial stage and the stages. Interestingly, the correlation between the transcriptomic and proteomic data was low, suggesting the importance of the post-transcriptional processes in the growth and development of T. guangdongense. Among the genes/proteins that were both differentially expressed during the developmental process, there were numerous heat shock proteins and transcription factors. In addition, there were numerous proteins involved in terpenoid, ergosterol, adenosine and polysaccharide biosynthesis that also showed significant downregulation in their expression levels during the developmental process. Furthermore, both tryptophan and tryptamine were present at higher levels in the primordium stage. However, indole-3-acetic acid (IAA) levels continuously decreased as development proceeded, and the enzymes involved in IAA biosynthesis were also clearly differentially downregulated. These data could be meaningful in studying the molecular mechanisms of fungal development, and for the industrial and medicinal application of macro-fungi. 1. Introduction In China, the fruiting bodies of some Cordyceps s.l. spp. are well-known traditional medicines because of their anticancer, immunomodulation, and anti-impotence effects [36][1], [37][2], [38][3]. Up to now, only four species of Cordyceps s.l., namely C. militaris, C. cicadae, Ophiocordyceps sinensis, and Tolypocladium guangdongense previously named as C. guangdongensis, have been approved by the Chinese government as medicines or edible fungi, although hundreds of species have been found. T. guangdongense grows on the hypogeous tubers of Elaphomyces fungi, and is different from other Cordyceps s.l. species, which grow on insects [39][4]. To date, all four of these species have been successfully cultivated, and C. militaris has the largest cultivated production, and the rest have only relatively very low cultivated yields [40][5], [41][6], [42][7], [43][8], [44][9]. The fruiting body of T. guangdongense, which was approved as a novel food source by the Chinese government in 2013, contains numerous active and nutritional compounds such as high levels of cordycepic acid, adenosine, polysaccharides, micro-elements, and vitamins [45][10], [46][11]. Based on animal tests, it has been shown that its fruiting body had a clear positive effect in improving the symptoms of several diseases, such as infection with the influenza virus H9N2 and chronic renal failure [47][12], [48][13]. In addition, based on the metabolomics analysis, T. guangdongense has many similar metabolites to O. sinensis (unpublished work), but can be more easily cultivated. These data suggest that T. guangdongense could have a range of application in the functional food and medicine industries, and is worthy of further study from many aspects. In one study, it was found that the mycelia of T. guangdongense grew very slowly in PDA medium, achieving a diameter of 2.72 cm after 14 days of growth. The formation and growth of the primordia were easily affected by temperature, light and other factors [49][5]. Therefore, exploring the molecular mechanism of primordial formation and fruiting body development is very important for improving the quality and yield of T. guangdongense and accelerating its industrialization. Recently, the genome of T. guangdongense and the gene expression profiles under different light illumination times have been reported by our team [50][14], [51][15]. Nevertheless, we do not yet have a comprehensive understanding of its growth and development. The regulatory network and the function of the genes are also unknown. Advances in omics technologies such as transcriptomics and proteomics have clearly facilitated the comprehensive analysis of gene and protein expression levels under multiple experimental conditions. For example, gene and protein expression profiling have provided a viewpoint into the genes and proteins involved in nutrition transport in Morchella importuna [52][16], bioactive metabolites in Hericium erinaceus [53][17], temperature stress in Lentinula edodes [54][18], Cd^2+ stress in Pleurotus eryngii [55][19], special odor formation in Schizophyllum commune [56][20] and fruiting body development in Flammulina velutipes and Dictyophora indusiata [57][21], [58][22]. In the cases of the Cordyceps s.l. species, there have been transcriptomic and proteomic studies performed in O. sinensis, C. militaris and C. cicadae. Transcriptome analyses in O. sinensis have identified genes in the key pathways involved in fruiting body and sexual development-related genes at differential developmental stages [59][23], [60][24]. It has been proven that the levels of nutrients between artificially cultured and wild O. sinensis samples were virtually the same in terms of nucleosides, nucleotides and adenosine from the perspective of the proteome and metabolome [61][25]. Using multi-omics studies in C. militaris, researchers have uncovered gene expression differences between the mycelium and the fruiting body [62][26], the transcriptional regulation of central carbon metabolism on sucrose or glucose medium [63][27], and genes related to the biosynthesis of carotenoids [64][28]. In C. cicadae and C. kyushuensis, omics data have revealed the nature of asexual-fruiting and the genes related to cordycepin and pentostatin biosynthesis [65][29], [66][30]. These findings show that omics technologies can be effectively used for the analysis of molecular mechanisms in Cordyceps fungi. The expression patterns of genes and proteins are likely a reflection of dynamic changes in gene and protein expression during the growth and development of T. guangdongense. Nevertheless, the transcriptome or proteome of T. guangdongense have not been described. Therefore, in the present study, fungal samples from differential developmental stages were collected from T. guangdongense, and used to perform transcriptome and proteome analyses. Hub genes were detected using weighted gene co-expression network analysis (WGCNA). It is hoped that our viewpoint of the potential molecular mechanisms involved in the developmental process could be substantially improved with the help of comprehensive analyses of the transcriptome and proteome. These data should increase the knowledge of the complex mechanisms regulating primordial formation and fruiting body development, and provide theoretical support for further improvements in cultivation techniques for T. guangdongense and for its application in the health-care industry. 2. Materials and methods 2.1. Sample collection at differential developmental stages T. guangdongense strain CCTCCM206051 was used as the test strain, and was deposited in the China Center for Type Culture Collection. T. guangdongense was cultivated as described in a previous study [67][5]. After growth for 20 d on PDA medium at 23℃, mycelia from T. guangdongense were transferred into rice medium and incubated for 6 days in the dark at 23℃ and 60–70% relative humidity. For primordium formation and fruiting body growth, the cultivation room was ventilated daily for 30 min and illuminated (500 lx) for 10 h. Using a previously described method [68][30], samples of the mycelium (My), primordium (P), young fruiting body (YFB) and mature fruiting body (MFB) were collected. In addition, a sample of the stroma with length 2.5–4 cm was designated as the developed fruiting body (DFB). Three replicate samples were prepared for each developmental stage. After collection, all samples representing the five differentially developmental stages were immediately frozen in liquid nitrogen and then stored at −80 °C prior to RNA and protein isolation. 2.2. RNA isolation and sequencing Total RNA was extracted from all samples using a TRIzol Kit (Invitrogen, Dalian, China) according to the manufacturer’s instruction. After qualitative and quantitative assessment of the RNA samples using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA), the mRNAs were enriched using oligo-dT magnetic beads, and then fragmented and transcribed into cDNA using random oligonucleotides. The five cDNA libraries were sequenced on a BGISEQ-500 platform at Beijing Genomics Institute (BGI)-Shenzhen, Wuhan, China. 2.3. Transcriptome data processing After removing reads containing poly-N, sequencing adapters and low-quality reads from the raw data, the clean reads were mapped to the T. guangdongense genome using HISAT2 and Bowtie2 [69][31], [70][32]. The RSEM package was used for calculating the expression levels of genes, and the differential expression analysis was performed using DESeq2 with a Q value < 0.05 [71][33], [72][34]. Genes with a false discovery rate (FDR) of ≤ 0.001 and an absolute value of log[2]Ratio ≥ 1 were designated as differentially expressed genes (DEGs). The analysis of gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEEG) pathway enrichment was performed using Phyper based on the Hypergeometric test [73][35]. In addition, a weighted correlation network analysis (WGCNA) was performed to analyze the gene co-expression network and identify hub genes [74][36]. 2.4. Protein extraction Mycelia, primordia and developmental fruiting bodies (0.1 g) from each group (three replicates for each group) were homogenized in 0.5 mL of extraction buffer (20 mM Tris-HCl pH = 7.4, 100 mM EDTA, 2% β-mercaptoethanol, 1 mM DTT and 1% Triton X-100). Following this, 1 mL of saturated phenolic Tresturturated Phenol was added, and the mixture was shaken well for 15 min, followed by centrifugation at 25000 × g for 15 min at 4 °C to obtain the supernatant (phenol phase). Then, 5 volumes of pre-cooled 0.1 M ammonium acetate in methanol supplemented with 10 mM dithiothreitol (DTT) was added to precipitate the proteins at −20 °C for 2 h. The pelleted proteins were washed two times with pre-cooled acetone. After air-drying of the precipitate, a 1 × Cocktail (Sigma, St Louis, MO, USA) was added and the samples was incubated in an ice bath for 5 min, and then centrifuged at 25000 × g for 15 min at 4 °C. Subsequently, 10 mM DTT was added to the precipitate which was then incubated in a water bath at 56 °C for 1 h. The sample was then incubated with 55 mM iodoacetamide for 45 min in the dark. The supernatant containing proteins was quantified using the Bradford method and SDS-PAGE was used to assess the quality of the protein extract [75][18]. 2.5. Proteomic data processing The final protein solution was diluted with 0.5 M TEAB and then digested with trypsin (Promega, Madison, WI, USA) (enzyme: protein ratio of 1:20) at 37 °C for 2 h. Peptides were labelled using the IBT method as previously described [76][37]. A mass spectrometric analysis of soluble intracellular proteins was undertaken, as previously described [77][38]. Raw MS/MS data was queried against the T. guangdongense genome using the Mascot search engine MaxQuant (version 2.3.02), and was deposited in the ProteomeXchange Consortium (identifier number PXD020391). The peptides labeled with IBTs were quantitatively analyzed using the IQuant software [78][38]. Proteins with a P < 0.05 and a fold change > 1.50 or < 0.67 were defined as differentially expressed proteins (DEPs). These proteins were classified and grouped using the GO and KEGG databases, and a hypergeometric test was used to analyze the GO and KEGG pathway annotations. The interaction network and subcellular location were predicted using the STRING protein interaction database and WoLF PSORT software [79][39], [80][40]. 2.6. Quantification of tryptophan and IAA Samples (0.5 g) of mycelia, primordia and fruiting bodies were ground and then shaken for 15 min in 1 mL of buffer containing 50% methanol and 0.1% formic acid. The mixture was subjected to ultrasound for 30 min in an ice bath to ensure complete homogenization of the sample. The samples were kept at 4℃ for 30 min, and centrifuged at 10000 × g for 20 min, and the supernatant was collected. The samples were filtered using a 0.3 μm filter plate and were then analyzed by LC-MS/MS (Waters Iclass-AB Sciex 6500) at the Beijing Genomics Institute (BGI)-Shenzhen, Wuhan, China. MultiQuant was used for chromatographic peak extraction and quantitative analysis of tryptophan, tryptamine, and IAA. Standards of these three metabolites were purchased from Yuanye (Shanghai, China). 2.7. Quantitative real time polymerase chain reaction (qRT-PCR) verification of gene expression The mRNA prepared for transcriptome sequencing was also used for qRT-PCR analysis. Sixteen DEGs were selected to verify the reliability of the transcriptome. The gene-specific primers used for gene quantification were designed using Primer Primer5.0. All primers are listed in [81]Table S1. Histone H4 was used as the internal control [82][41]. The qRT-PCR reaction was performed using Applied Biosystems ABI 7500 (ABI, Foster City, CA, USA) with three biological and technical repeats. The PCR reaction was performed with an initial denaturation for 2 min at 95 °C; 40 cycles of 20 s at 95 °C and annealing at 60 °C for 20 s. The relative expression levels of target genes were calculated using the 2−^△△Ct (Ct, cycle threshold) method as described by Livak [83][42]. 3. Results 3.1. Transcriptome sequencing analysis To investigate the gene expression changes in T. guangdongense during the developmental process, triplicate samples from the five developmental stages were used to construct a total of 15 cDNA libraries after which transcriptome sequencing was performed using BGISEQ-500 sequencing, yielding a total of 334.3 million raw reads. After removing reads containing adapter sequences, poly-N, and the low-quality reads, 330.29 million clean reads were obtained, with a mean of 22.02 million clean reads for each replicate ([84]Table S2). Approximately 80% of the reads from the mycelia could be mapped to the T. guangdongense genome, whereas the mapping ratios of the reads from the other samples were>97.46% ([85]Table S2). The Q30 percentages for all the clean reads in the 15 libraries were>93%. Gene expression levels between replicates for each sample exhibited a high Pearson’s correlation coefficient value ([86]Fig. S1), indicating good repeatability between replicates.(See [87]Fig. 1) Fig. 1. [88]Fig. 1 [89]Open in a new tab Genes detected and expressed in all samples during the five developmental stages. A, Distribution diagram showing genes with different FPKM values across the five differential developmental stages. B, Venn diagram showing the overlap between genes expressed at each developmental stage. My, Mycelial stage; P, Primordial stage; YFB, Young fruiting body stage; DFB, Developmental fruiting body stage; MFB, Mature fruiting body stage. more than 99% of genes with FPKM values (>0) were detected in the 15 samples ([90]Fig. 1A and [91]Table S3). Approximately 86.31%, 98%, 98%, 98.13% and 98.37% of the genes were expressed in the My, P, YFB, DFB and MFB stages, respectively ([92]Fig. 1A). A total of 8817 genes were expressed across the five differential developmental stages ([93]Fig. 1B). Based on the gene FPKM values ([94]Table S3), the expression levels of all genes were classified into four categories. During the developmental process, the majority of genes with FPKM values (10 ≤ FPKM 〈1 0 0) were considered to be moderately expressed. Approximately 17% of genes with FPKM value > 100 were considered to be highly expressed, and approximately 25% of genes with FPKM values (0 < FPKM < 10) were considered to have low levels of expression. The number of non-expressed genes (FPKM = 0 or undetected) accounted for around 2% of the total genes at each developmental stage. 3.2. Identification and functional classification of differentially expressed genes across the five developmental stages Using a principal component analysis ([95]Fig. S2), it was noted that the samples from the five developmental stages, each with three biological replicates, could be clearly divided into two groups, with the mycelial stage separately grouped into one cluster. As shown in [96]Fig. 2A and [97]Table S4, the number of DEGs between the My and P stages was the largest (3673), accounting for 40.14% of the T. guangdongense genes, followed by the DFB vs. YFB stages (4 0 0). Fewer DEGs were found in comparisons between the MFB and DFB stages (2 5 7) and between the YFB and P stages (1 8 9). The greatest number of unique DEGs (3157) was found between the MY and P stages, whereas only 25 unique DEGs were found between the YFB and P stages ([98]Fig. 2B). In addition, only 12 DEGs were shared among the four comparisons between adjacent development stages ([99]Fig. 2B). These results suggest that the gene expression profile at the mycelial stage is unique compared to the other developmental stages, whereas the gene expression profiles across the four other developmental stages are very similar to one another with only minor differences being found. Fig. 2. [100]Fig. 2 [101]Open in a new tab Number of DEGs found between adjacent developmental stages. A, Number of DEGs found between two adjacent development stages. The number of DEGs is shown on the top of the histograms. B, Venn diagram showing the overlaps between the total number of DEGs identified in each of the four comparisons of adjacent development stages. My, Mycelial stage; P, Primordial stage; YFB, Young fruiting body stage; DFB, Developmental fruiting body stage; MFB, Mature fruiting body stage; vs, Versus. Classification of DEGs using GO annotation was performed to evaluate the potential functions of these DEGs across the five developmental stages. As shown in [102]Fig. 3A and [103]Table S5, for biological processed (BP) a total of 15 categories were found, including metabolic process (1043), cellular process (9 9 4), localization (3 1 9), biological regulation (2 2 4), cellular component organization or biogenesis, (1 8 8) and response to stimulus (1 8 3). For cellular component (CC), the largest of the 13 categories was related to membrane (1124). Other categories were related to cell (9 8 2), organelle (6 9 0) and macromolecular complex (3 3 2). For molecular function (MF), catalytic activity (1838) and binding activity (1374) were the largest categories. In comparisons between the My and P stages and between the DFB and MFB stages, the most significantly enriched GO assignment was a CC related to membrane ([104]Fig. 3B and [105]Table S5). The most significantly enriched GO terms for the P vs. YFB stages and the YFB vs. DFB stages were in the MF category involved in catalytic activity. In addition, oxidoreductase and monooxygenase activities were found in a comparison of the P vs. YFB stages. Fig. 3. [106]Fig. 3 [107]Open in a new tab GO functional classification of DEGs in each of the four comparisons of two adjacent development stages. A, Gene ontology term assignments for the DEGs. B, The most enriched GO terms. The red lines represent biological process. The green lines represent cellular component. The blue lines represent molecular function. Only the most significant GO terms with P values < 0.005 were shown. My, Mycelial stage; P, Primordial stage; YFB, Young fruiting body stage; DFB, Developmental fruiting body stage; MFB, Mature fruiting body stage; vs, Versus. (For interpretation of the references to colour in this figure legend, the