The Cordyceps blackwelliae strain ZYJ0835 used in this study was isolated from a fresh specimen collected from Guangdong, China in 2019 and deposited at the China Center for Type Culture Collection (CCTCC) under no. M2022663 (Fan et al. 2022). The strain was grown on potato dextrose agar (PDA) overlaid with sterile cellophane at 25°C for seven days. Mycelia were harvested for genomic DNA extraction using the cetyl trimethyl ammonium bromide method (Zhang et al. 2010). To prepare samples for RNA sequencing and metabolic analysis, fungal materials were prepared at three distinct growth stages in accordance with the procedures described in our previous study (Fan et al. 2022). Briefly, the fungus was cultivated at 25°C and 160 r/min for six days using the potato dextrose broth medium to obtain the liquid culture mycelia (LM), which represents the undifferentiated, fast-growing vegetative phase. The wheat culture mycelia (WM) were obtained by inoculating the fungus to wheat media in glass bottles and cultivating the fungus under dark conditions at 20°C with humidity levels of 60–70%. Mycelia were collected when the mycelia reached full growth across the wheat substrate. This stage mimics the natural substrate colonisation phase where invasive growth initiation occurs. To induce fruiting body (FB) formation, duplicate WM cultures were exposed to 200 lux fluorescent light with alternating temperature (22°C for 14 h and 17°C for 10 h) to stimulate primordial formation. After primordia differentiated, the fungus was kept at 22°C with 80% relative humidity and received a daily photoperiod of 12 hours under a light intensity of 1,000 lux. Fruiting bodies were then ready to be harvested when they grew to 3–5 cm long. The FB stage represents the reproductive phase, crucial for studying developmental morphogenesis. All collected samples were promptly frozen in liquid nitrogen and stored at -80°C until further use. Total RNA was extracted using Trizol reagent (DP424, Tiangen Biotech Co., Ltd., Beijing, China) with quantity verification performed using an Agilent 2100 bioanalyser (Agilent Technologies, Palo Alto, CA, USA).

Genome sequencing was performed by Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China) using both Illumina NovaSeq 6000 (PE150) and Nanopore PromethION R9.4 platforms. Illumina data underwent quality control filtering by fastp v.0.23.0 (Chen et al. 2018) to remove low-quality reads. For Nanopore data, quality control was conducted using Guppy v.4.3.4 within the MinKNOW v.2.2 package for base calling and only reads with an average Q score exceeding seven were selected for subsequent analyses. Genome assembly was performed firstly from the Nanopore data using CANU assembler v.1.7 (Koren et al. 2017), followed by polishing with the Illumina data using Pilon v.1.23 (Walker et al. 2014). The quality of the genome assembly was assessed using Benchmarking Universal Single-Copy Orthologs (BUSCO) v.3.0.2, based on the fungi_odb10 database (Simão et al. 2015) and Core Eukaryotic Genes Mapping Approach (CEGMA) v.2.5 (Parra et al. 2007). Telomeres were manually annotated by identifying tandem repeat sequences (CCCTAA)n or (TTAGGG)n located at the ends of scaffolds (Lee et al. 2021). Transfer RNAs (tRNAs) were detected using tRNAscan-SE v.2.0 (Lowe et al. 1997). Genome annotation was carried out with the MAKER2 v.2.31.9 pipeline (Holt et al. 2011), employing three ab inito gene prediction methods: Augustus (Stanke et al. 2006), GeneMark-ES (Lukashin et al. 1998) and SNAP (Korf 2004). Local BLAST searches were performed against the NCBI Non-Redundant Protein Database (Nr) (Sayers et al. 2022), the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al. 2000) and the Clusters of Orthologous Groups of proteins (COG) (Tatusov et al. 2000) using DIAMOND v.0.8.35 (Buchfink et al. 2021). Additionally, HMMER v.3.1b2 (Eddy 2011) was used for local searches of protein signatures against the Pfam database (Finn et al. 2016). Gene Ontology (GO) (Ashburner et al. 2000) annotations was conducted using Blast2GO v.2.5 (Götz et al. 2008). The upset plot for multi-database annotation was generated using the online tool ChiPlot (https://www.chiplot.online/). Carbohydrate-active enzymes (CAZymes) were identified through the CAZy database (https://www.cazy.org/). Predictions of secondary metabolite gene clusters were made using the web-based fungal version of antiSMASH v.7.0 (Blin et al. 2023). The presence of signal peptides in putative proteins was predicted with SignalP v.4.1 (Petersen et al. 2011) and transmembrane domain predictions was performed using TMHMM v.2.0 (Krogh et al. 2001). Sequences harbouring a signal peptide cleavage site, but lacking a transmembrane region were categorised as putatively secreted proteins.

The aforementioned annotation methodology effectively identified APN2 (encoding a DNA lyase), one of the genes commonly located adjacent to mating-type (MAT) genes. However, the MAT genes (MAT1-1-1 and MAT1-2-1) and SLA2 (encoding a cytoskeleton assembly control protein) were not annotated. To locate these genes within the genome of C. blackwelliae, local BLAST searches were conducted with corresponding genes of C. militaris (GenBank accession nos. KP721256 for MAT1-1-1, KP721273 for MAT1-2-1 and XM_006671662 for SLA2) as queries. To assess the expression of the identified MAT genes in C. blackwelliae, RNA sequencing (RNA-seq) and real-time quantitative reverse transcription PCR (RT-qPCR) analyses were conducted, following the methodologies outlined below.

Structural comparisons at the MAT locus were performed across nine additional Cordyceps species/strains, including C. cateniannulata (FRD 24 and MBC 23), C. chanhua (ZJ1611, BA-001 and CC02), C. fumosorosea (ARSEF 2679), C. javanica (F084), C. militaris (ATCC 34164) and C. tenuipes (YF17002). All sequences were retrieved from GenBank (Suppl. material 1: table S2) and the MAT, APN2 and SLA2 genes were identified through local BLAST searches as described above. The results were visualised using IBS v.1.0 (Liu et al. 2015). Amino acid sequences of the MAT genes were aligned using the MUSCLE algorithm integrated within MEGA v.11 (Tamura et al. 2021) and the conserved domain of MAT1-2-1 was predicted via the SMART website (https://smart.embl.de/) (Schultz et al. 1998).

To investigate the evolution of MAT genes in Cordyceps, nucleotide sequences of MAT1-1-1, MAT1-2-1, as well as the concatenated sequences of APN2 and SLA2, were used for phylogenetic analyses. The nine Cordyceps species/strains mentioned earlier, along with C. blackwelliae ZYJ0835, were designated as the ingroup, with Beauveria asiatica ARSEF 4384 serving as the outgroup. Following sequence alignment by MUSCLE (Edgar 2004), poorly aligned regions were removed by trimAl (Capella-Gutiérrez et al. 2009). Phylogenetic relationships were estimated using the Maximum-Likelihood method as implemented in IQ-TREE v.1.6.12 (Nguyen et al. 2015) with 1,000 bootstrap replicates. The resulting tree files were visualised using FigTree v.1.4.4 (Rambaut 2009).

Two datasets were used for phylogenomic analyses. Dataset I comprised 3,488 single-copy nuclear genes, which were extracted using BUSCO v.5.8.0 with the hypocreales_odb10 database (Seppey et al. 2019). Dataset II consisted of 14 core mitochondrial protein-coding genes (i.e. atp6, 8-9, cob, cox1-3, nad1-6 and nad4L). For dataset I, nuclear phylogenomic matrices were generated using the BUSCO_phylogenomics pipeline (accessible at https://github.com/jamiemcg/BUSCO_phylogenomics). The analysed fungal taxa are detailed in Suppl. material 1: table S3. A Maximum-Likelihood phylogenetic tree was constructed, based on concatenated sequences of 2,039,507 amino acid positions with the method described above. For dataset II, fungal species from genera consistent with those used in dataset I were selected (Suppl. material 1: table S4). Phylogenetic relationships were estimated using a Maximum-Likelihood approach, based on the concatenated sequences of 4,274 amino acid positions using the aforementioned methodology.

Complementary DNA (cDNAs) libraries were constructed by enriching mRNA from total RNA using polyT oligo-attached magnetic beads and then reverse-transcribed into cDNA. Subsequently, the cDNA libraries were sequenced by Novogene Technology Co., Ltd. (Beijing, China) on an Illumina NovaSeq platform with PE150 reads. Raw sequencing data underwent filtration using fastp v.0.23.0 (Chen et al. 2018) and then clean reads were aligned to the C. blackwelliae genome using HISAT2 v.2.0.5 (Kim et al. 2015). The mapped reads for each sample were subsequently assembled using StringTie v.1.3.3b (Pertea et al. 2015). The number of reads mapped to each gene was quantified using FeatureCounts v.1.5.0-p3 (Liao et al. 2014), which facilitated the calculation of gene expression levels represented as fragments per kilobase exon model per million mapped reads (FPKM). Differential expression analysis amongst sample groups was performed using the R package DESeq2 v.1.20.0 (Love et al. 2014), with P-values adjusted using the Benjamini and Hochberg’s approach (Benjamini et al. 2018) for controlling the false discovery rate. Genes with an adjusted P-value of less than 0.05 and a |log2FoldChange| greater than 2 were regarded as DEGs. GO and KEGG enrichment analyses were performed using clusterProfiler v.3.8.1 (Yu et al. 2012), with an adjusted P-value threshold of less than 0.05 deemed statistically significant.

Alternative splicing (AS) events were analysed using rMATS v.4.1.0 (Shen et al. 2014). Various types of AS events, including exon skipping (SE), mutually exclusive exon (MXE), intron retention (RI), alternative 5’ splicing (A5SS) and alternative 3’ splicing (A3SS), were compared across different growth stages, with a false discovery rate (FDR) < 0.05 indicating statistical significance. The significant AS events were visualised using the Python programme rmats2sashimiplot (https://github.com/Xinglab/rmats2sashimiplot). Additionally, GO enrichment analysis was performed on all genes associated with the identified AS events using the aforementioned methodology.

To confirm the reliability of RNA-seq data, four DEGs (gene 6087, gene 0120, gene 1814 and gene 1408) and two non-DEGs (gene 5222 and gene 0234) were randomly selected and their expression levels were quantified through RT-qPCR. Complementary DNA were synthesised from 100 ng of mRNA, enriched as described above, using the TransScript^® One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech, Beijing, China) in a 20 μl reaction mixture that included 1 μl of TransScript^® RT/RI Enzyme Mix, 10 μl of 2×TS Reaction Mix, 1 μl of Anchored Oligo(dT)₁₈ Primer, 1 μl of gDNA Remover, 2 μl of RNase-free water and 5 μl of mRNA. RNA concentration and purity (i.e. OD260/OD280) were measured using NanoDrop^TM One Microvolume UV-Vis Spectrophotometers (Thermo Scientific, USA). RT-qPCR was performed on a CFX Connect^TM Real-Time PCR system (Bio-Rad, USA). The RT-qPCR reaction mixture (in 20 μl) consisted of 10 μl of 2× PerfectStart^® Green qPCR SuperMix, 0.4 μl of each primer (10 μM), 3 μl of diluted cDNA (200 ng/μl) and 6.2 μl of nuclease-free water. The RT-qPCR programme included an initial denaturation step at 95°C for 30 s, followed by 40 cycles of 10 s at 95°C, 20 s at 56°C and 15 s at 72°C. All reactions were conducted in biological triplicates. Fluorescence signals were recorded during the extension stage and a melting curve was set to detect the specificity of the amplified products. Relative gene expression was evaluated using the 2^-ΔΔCT method (Livak et al. 2001). The 18S rRNA was utilised as an internal reference, which has been previously validated as a stable internal reference gene in C. blackwelliae (Li et al. 2025b). All primers used in this study are listed in Suppl. material 1: table S5. Least significance difference (LSD) test was used to assess the statistical significance of differences observed across different stages.

The R package WGCNA was employed to conduct a weighted correlation analysis based on all genes (Langfelder et al. 2008). An adjacency matrix was calculated with the optimal soft-thresholding power (β) (Suppl. material 2: fig. S1) and transformed into a topological overlap matrix (TOM), which was further converted into a distance matrix to build a tree of gene clusters. A dynamic cutting algorithm was applied to cluster all genes into colour-coded modules. To identify modules that are strongly related to phenotypic traits, Pearson’s correlation coefficients (Benesty et al. 2009) were calculated between each module and each sample. A heat map was drawn according to these correlation coefficients, allowing for the selection of phenotype-related modules. A gene co-expression network was constructed, with a focus on edges with a weight exceeding 0.2 in the magenta and black modules and a weight greater than 0.36 in the turquoise module. The resulting network was visualised using Cytoscape 3.10.2 (Shannon et al. 2003). Hub genes were identified as the most highly connected nodes within each module.

Samples of LM, WM and FB were subjected to grinding in liquid nitrogen to produce fine powders suitable for Ultra-High-Performance Liquid Chromatography coupled with Tandem Mass Spectrometry (UHPLC-MS/MS) analysis. This analysis was conducted using a Vanquish UHPLC system (Thermo Fisher, Germany) in conjunction with an Orbitrap Q Exactive^TM HF mass spectrometer (Thermo Fisher, Germany). Samples were injected on to a Hypersil Gold column (100 × 2.1 mm, 1.9 μm) using a 12-min linear gradient at a flow rate of 0.2 ml/min. For a positive polarity mode, analyses were performed with (A) 0.1% formic acid in water and (B) methanol. For a negative polarity mode, mobile phases comprised (A) 5 mM ammonium acetate (pH = 9.0) and (B) methanol. The separation process followed the gradient: 2% B, 1.5 min; 2%–85% B, 3 min; 85%–100% B, 10 min; 100%–2% B, 10.1 min; and 2% B, 12 min. The Q Exactive^TM HF mass spectrometer was operated in positive/negative polarity mode with a spray voltage set at 3.5 KV, a capillary temperature of 320°C, a sheath gas flow rate of 35 psi, an aux gas flow rate of 10 l/min, an S-lens RF level of 60 and an aux gas heater temperature of 350°C.

The raw data files obtained from UHPLC-MS/MS were processed using Compound Discoverer 3.3 (Thermo Fisher, Germany) for data pretreatments, which included peak alignment, peak feature extraction and quantification. Subsequently, the identified peaks were matched against mzCloud (https://www.mzcloud.org/), mzVault (Thermo Fisher, Germany) and MassList (Novogene, China) databases. Statistical analyses were performed using R v.3.4.3 and Python v.2.7.6. The metabolites were annotated using HMDB (https://hmdb.ca/metabolites) and LIPIDMaps (http://www.lipidmaps.org/) databases. Partial Least Squares Discriminant Analysis (PLS-DA) was performed using metaX (Westerhuis et al. 2010; Wen et al. 2017). Univariate analysis (t-test) was applied to calculate the statistical significance (P-value). Metabolites were considered as DMs if they exhibit a variable importance in projection (VIP) > 1, a P-value < 0.05 and a fold change (FC) ≥ 2 or FC ≤ 0.5.

To investigate the correlation between DEGs and DMs, Pearson’s correlation coefficient between gene expression and metabolite abundance were calculated using the mixOmics package in R (Rohart et al. 2017). The top 50 DMs and top 50 DEGs were presented, based on their statistical significance (P-value). A correlation heatmap illustrating the relationship between genes and metabolites was built using the corrplot package in R (Wei et al. 2017). Additionally, KEGG co-enrichment diagrams were created by utilising the ggplot2 package in R (Wickham. 2009). The global metabolic pathway was visualised through iPath v.3.0. (https://pathways.embl.de/).

To establish a reference genome for C. blackwelliae, Illumina and Nanopore sequencing technologies were employed to generate short and long reads, respectively. These sequencing efforts yielded 6.55 Gb of Illumina data (21,697,168 paired reads, 211× coverage) and 5.37 Gb of Nanopore data (623,214 reads, 173× coverage). Except for a single scaffold (i.e. Scaffold10) representing the mitochondrial genome, the nuclear genome was assembled into 11 scaffolds, totally 31.06 Mb, with an N50 value of 4.14 Mb (Table 1). Notably, telomeric repeats were identified at both termini of three scaffolds (i.e. scaffolds 2, 7 and 8) and at either the 5’ or 3’ end of five scaffolds (i.e. scaffolds 1, 3, 9, 11 and 12) (Suppl. material 1: table S6). BUSCO and CEGMA analyses confirmed that the assembly captured 99.3% of fungal single-copy orthologs and 95.6% of core eukaryotic genes, respectively (Table 1).

A total of 8,136 genes were predicted, including 7,969 protein-coding genes, 129 tRNA genes and 38 rRNA genes. Functional annotations using NR, GO, KEGG, COG and Pfam databases revealed that 1,585 genes were annotated across all databases and 6,435 genes had annotation in at least two databases (Suppl. material 2: fig. S2; Suppl. material 1: table S7). The most enriched GO terms were related to “Cellular anatomical entity”, “Catalytic activity” and “Binding” (Suppl. material 2: fig. S3a). These annotated genes accounted for 54.7% of the total protein coding genes and are likely to play a significant role in developmental processes. The top three COG categories were “Carbohydrate transport and metabolism”, “General function prediction only” and “Translation, ribosomal structure and biogenesis”, which encompassed 351, 326 and 294 genes, respectively (Suppl. material 2: fig. S3b). KEGG annotation assigned 3,126 proteins to 387 distinct pathways, categorised into seven major functional groups: cellular processes (5 branches, 416 genes), environmental information processing (3 branches, 229 genes), genetic information processing (4 branches, 610 genes), human diseases (12 branches, 524 genes), metabolism (12 branches, 847 genes), organismal systems (10 branches, 363 genes) and 1,194 unclassified genes. Amongst the annotated genes, the category “Global and overview maps” contained the highest number of genes, followed by “Transport and catabolism”, “Signal transduction”, “Translation” and “Carbohydrate metabolism” (Suppl. material 2: fig. S3c).

To investigate C. blackwelliae fruiting body formation, the MAT locus structure was analysed. The genome contained a complete MAT1-2-1 gene and a truncated MAT1-1-1 gene, a pattern conserved across ten Cordyceps genomes (Fig. 1a). The two MAT genes with opposite transcriptional orientations were found to be more distantly spaced in C. blackwelliae (7.3 kb) than those in other Cordyceps species/strains (1.5–2.1 kb). The MAT1-2-1 gene within these Cordyceps genomes contained two introns (46–51 bp) and the amino acid sequence alignment revealed an 85.36% conservation rate amongst different Cordyceps species/strains (Fig. 1b). The HMG-box domain was identified through searches in the SMART database, with strong E-values (1.53e−17 to 8.76e−18).

Figure 1. 

Structural and evolutionary analysis of the MAT locus in C. blackwelliae: a schematic representation of the MAT1-2 locus across various Cordyceps species/strains. The APN2 gene, which encodes a DNA lyase and the SLA2 gene, which encodes a protein involved in cytoskeleton assemble control, flank the MAT genes. Gene orientations are indicated by arrows. Please refer to Suppl. material 1: table S2 for details of fungal strains and corresponding scaffolds harbouring MAT, APN2 and SLA2 genes; b alignment of amino acid sequences of the MAT1-2-1 protein, with the conserved HMG-box domain indicated by a red box; c, d validation of MAT gene expression at three different developmental stages through RT-qPCR analysis. Relative expression was calculated using the 2^-ΔΔCT method and normalised against the LM stage of the truncated MAT1-1-1 gene. Gene expression levels (FPKM) from RNA-seq data were presented as heatmaps e–g phylogenetic analyses based on the concatenated SLA2 and APN2 genes (e), the truncated MAT1-1-1 gene (f) and the MAT1-2-1 gene (g).

To further clarify the potential functions of MAT1-2-1 and the truncated MAT1-1-1 gene in the regulation of fruiting body formation, gene expression profiles were examined across three different developmental stages. The MAT1-2-1 gene exhibited differential expression (read counts 0–122 and FPKM values 0–8.67), with the lowest expression levels observed at the FB stage, while expression levels increased by 8.2-fold and 4.2-fold at the LM and WM stages, respectively (Suppl. material 1: table S8). In contrast, the truncated MAT1-1-1 gene demonstrated minimal expression (read counts 0–11 and FPKM values 0–1.41) throughout the three developmental stages (Suppl. material 1: table S8). RT-qPCR confirmed these expression patterns (Fig. 1c, d).

To ascertain whether the different MAT genes exhibit divergent evolutionary trajectories, phylogenetic analyses were conducted, based on nucleotide sequences of MAT1-2-1, the truncated MAT1-1-1 or concatenated sequences of SLA2 and APN2. The topologies of the three phylogenetic trees were found to be basically congruent, consistently placing C. blackwelliae in close relationship to C. chanhua and C. tenuipes (Fig. 1e–g).

To elucidate the phylogenetic placement of C. blackwelliae within the family Cordycipitaceae, phylogenetic trees were constructed, based on 3,488 nuclear single-copy orthologs and 14 mitochondrial core protein-coding genes, respectively. Both analyses robustly clustered Cordyceps species together (Fig. 2). Each of the other genera that included more than two taxa formed distinct clades, with the exception of Lecanicillium. Both phylogenetic trees consistently positioned C. blackwelliae in close phylogenetic proximity to C. chanhua and C. tenuipes, although some minor discrepancies were noted (Fig. 2). The mitochondrial phylogeny positioned C. blackwelliae closely with C. chanhua, followed by C. tenuipes, while the nuclear phylogeny grouped C. chanhua and C. tenuipes together, with C. blackwelliae as their closest relative.

Figure 2. 

Phylogenomic analyses based on both nuclear and mitochondrial genomes. The nuclear phylogeny (left) was constructed, based on 3,488 single-copy orthologs, while the mitochondrial phylogeny (right) was constructed, based on 14 core protein-coding genes. Detailed information regarding the fungal taxa and their corresponding accession numbers can be found in Suppl. material 1: tables S3, S4. The target fungus of this study is marked with a star (black star).

To elucidate the molecular mechanisms underlying fruiting body development, a total of 15 RNA-seq libraries were generated from three growth stages (i.e. LM, WM and FB), each with five biological replicates. The sequencing process yielded 631 million raw reads; following quality control and cleaning procedures, 601 million clean reads were retained, resulting in an average of 40.1 million reads per replicate (Suppl. material 1: table S9). Notably, more than 90% clean reads were successfully mapped to the C. blackwelliae genome and 553 novel genes were identified (Suppl. material 1: table S10). Principal component analysis (PCA) demonstrated clear separation amongst stages and strong replicate correlations (Suppl. material 2: fig. S4a).

To identify genes that may play a crucial role in mycelial growth and fruiting body development, three comparative analyses were conducted: WM vs. LM, FB vs. WM and FB vs. LM. The analyses revealed that 2,078 genes exhibited differential expression with a |log2FoldChange| > 2 (Suppl. material 1: table S11), with 121 DEGs shared across all three comparisons (Suppl. material 2: fig. S4b). Heatmap clustering showed that the highest number of DEGs was observed in the FB vs. LM comparison, while the lowest was noted in the FB vs. WM comparison (Fig. 3a). Transcription profiles of DEGs were visualised through volcano plots utilising -log10(padj) and log2(FoldChange) as statistical parameters (Fig. 3b–d), highlighting the most significantly regulated genes. Notably, MetB (gene 6635), encoding a putative cystathionine gamma-synthase involved in Cys/Met metabolism, was not expressed at the FB stage, but active at both LM and WM stages (FPKM values 4.6–314.85). At the LM stage, gene 4153, encoding a putative HSP20-like chaperone, was dramatically up-regulated by 385-fold compared to the WM stage and 466-fold compared to the FB stage. Conversely, gene 6087, encoding a putative RNA binding protein, was significantly down-regulated by 38,736-fold compared to the WM stage and 77,560-fold compared to the FB stage (Suppl. material 1: table S12). In addition, ribosomal genes were significantly up-regulated during the WM stage, while septin protein family genes were up-regulated during the FB stage. Furthermore, genes encoding chitin deacetylase and β-1,3-glucanase, related to cell wall remodelling, were up-regulated during the FB stage, reflecting dynamic cell wall modifications essential for fruiting body development (Suppl. material 2: fig. S5). To validate the RNA-seq results, four DEGs and two non-DEGs were randomly chosen for validation by RT-qPCR. The RT-qPCR results showed strong concordance with the RNA-seq data, confirming the reliability of our transcriptomic analysis (Suppl. material 2: fig. S6).

Figure 3. 

Gene expression profiles across three developmental stages: a clustering heat map of DEGs. Each column represents a sample, and each row corresponds to a gene. The FPKM values are indicated by a colour gradient, with red indicating high expression and green indicating low expression. The results of gene clustering are displayed as the dendrograms on the left; b–d volcano plot of DEGs between different developmental stages, with red dots indicating up-regulated genes, green dots indicating down-regulated genes and blue dots representing genes without statistically significant expression variation. The two vertical dashed lines represent |log2FoldChange| of 2, while the horizontal dashed lines represent the Padj value thresholds of 0.05; e C. blackwelliae samples at three different developmental stages that were used to perform RNA-seq and metabolomic analyses.

To identify pathways associated with phenotypes, GO and KEGG enrichment analysis was performed on the 2,078 DEGs. The findings indicated a significant overexpression of genes linked to DNA repair and recombination, including those involved in “Base excision repair”, “Mismatch repair” and “Homologous recombination”. This overexpression is crucial for maintaining genomic stability during rapid cell division and growth, particularly at the LM stage, where the mycelia exhibited rapid growth, forming dense and homogeneous mycelial pellets (Fig. 3e). Furthermore, the up-regulated expression of membrane-related GO terms, including “Integral component of membrane”, “Membrane parts” and “Transmembrane transporter activity”, also occurs at the LM stage. This up-regulation is crucial for maintaining cellular integrity, as well as for facilitating intercellular communication and signal transduction, supporting the observed rapid proliferation and biomass accumulation (Fig. 3e). At the WM stage, the up-regulated genes were significantly enriched in pathways associated with endopeptidase, oxidoreductase, peptidase and hydrolase activities. Additionally, the up-regulation of co-factor binding functions, including iron ion, heme and tetrapyrrole binding, is crucial for the proper functioning of many enzymes, as these co-factors are involved in electron transfer and catalytic reactions. This coordinated regulation of enzyme activity and co-factor binding likely supports the decomposition and utilisation of nutrients in the culture medium by the vegetative mycelia, which spread across the solid medium surface and penetrate the substrate to form an extensive hyphal network. At the FB stage, there was a significant up-regulation of genes involved in various biosynthesis and catabolism processes, including those related to amino acids, N-glycan, sphingolipids, carbohydrates, glycerophospholipids and nucleotides (Suppl. material 1: table S13). These processes are likely to be pivotal in the development of fruiting bodies and signify the transition from vegetative to reproductive growth, characterised by the differentiation of mycelia into specialised structures such as stalks and spore-producing tissues (Fig. 3e).

Alternative splicing (AS) events were examined across the three developmental stages (Suppl. material 1: table S14). A total of 707 AS events were identified in the FB vs. LM comparison, comprising 670 SE and 37 MXE. In the FB vs. WM comparison, we observed 658 AS events, consisting of 623 SE and 35 MXE. The WM vs. LM comparison revealed 704 AS events, with 668 SE and 36 MXE. Notably, we found no evidence of RI, A5SS and A3SS. Considering that several genes exhibited two or more AS events, we estimated that a total of 681 genes underwent AS with each gene being counted only once. Of particular interest, 53.23% of SE and 23.91% of MXE were present in all three comparisons. These findings indicated that AS is dynamically regulated during development. Nine significant AS events were identified according to our stringent criteria (Suppl. material 2: fig. S7, Suppl. material 1: table S14). GO annotation analysis revealed that 40 genes, associated with the “Structural molecule activity” term, exhibited alternative splicing, accounting for 27% of all genes categorised under this GO term. This was followed by “Transporter activity” and “ATP-dependent activity”, which also demonstrated a high frequency of AS events and a great proportion of affected genes (Suppl. material 2: fig. S8).

WGCNA was utilised to construct a gene co-expression network and to identify hub genes implicated in fungal development, based on 8,136 initially predicted genes and 553 novel genes. The analysis revealed 26 distinct gene modules (Fig. 4a). To identify modules strongly correlated with phenotypes, a correlation heatmap of modules and samples was generated (Fig. 4b). Based on the module-sample relationships, three modules were selected: a turquoise module (1,610 genes) that showed peak expression at the LM stage (Fig. 4c), a magenta module (333 genes) with positive correlation to the WM stage (Fig. 4d) and a black module (396 genes) that was overexpressed at the FB stage (Fig. 4e). These modules were used to construct the gene co-expression network (Suppl. material 1: table S15; Fig. 4f–h). The network revealed the following hub genes: at the LM stage—ARP3 (actin-like protein 3), HXT1 (low-affinity glucose transporter), LAC1 (L-ascorbate oxidase) and CIR2 (electron transfer flavoprotein-ubiquinone oxidoreductase) (Fig. 4f); at the WM stage—SSN2 (RNA polymerase II mediator subunit), AOX1 (alternative oxidase) and CHS7 (Chitin synthase III catalytic subunit) (Fig. 4g); and at the FB stage—gene 4088 (alpha-hydroxy acid dehydrogenase), GIT2 (glycerophosphodiester transporter) and BEA3 (ABC transporter), all likely involved in fruiting body formation (Fig. 4h). Notably, the occurrence of AS events in these hub genes (27.27%) was significantly higher than in the remaining genes (7.84%)

Figure 4. 

WGCNA investigation of modules associated with growth and fruiting body development: a display of original and merged modules within the clustering tree. The upper section of the diagram illustrates the gene clustering tree within the network. Each branch of the hierarchical tree or each vertical line in the accompanying colour bar corresponds to an individual gene. The middle section represents the original modules generated through dynamic tree clipping. The merged colours at the bottom indicate merged modules with a dissimilarity coefficient of less than 0.25. Genes not assigned to any modules are depicted in grey; b heatmap of correlations between modules and samples as indicated by a colour gradient, with red indicating a positive correlation and blue indicating a negative correlation between modules and samples. The turquoise, magenta and black modules are positively correlated with the LM, WM and FB stages, respectively; c–e gene expression patterns of the turquoise (c), magenta (d) and black (e) modules, with specific gene information available in Suppl. material 1: table S15; f–h gene co-expression networks corresponding to the turquoise (f), magenta (g) and black (h) modules. The construction of the gene co-expression network for each module involved selecting edges, based on weight thresholds of 0.36 for the turquoise module and 0.20 for both the black and magenta modules. The turquoise module comprises 58 nodes and 101 edges, the magenta module includes 56 nodes and 101 edges and the black module consists of 44 nodes and 95 edges. In this context, nodes represent individual genes, while edges represent interactions between gene pairs.

In order to elucidate the metabolic dynamics associated with morphogenesis and maturation, an analysis utilising UHPLC-MS/MS was performed. A comparative assessment of the total ion chromatograms indicated that, although the types of metabolites present across the three developmental stages were comparable, there were notable differences in the concentrations (Suppl. material 2: fig. S9). A total of 692 metabolites were identified in the positive ion mode and 469 metabolites in the negative ion mode (Suppl. material 1: table S16). These metabolites were classified into ten distinct categories, with lipids and lipid-like molecules making up the largest fraction (39.41%), followed by organic acids and their derivatives (21.19%) (Suppl. material 2: fig. S10a). Amongst the identified lipids and lipid-like molecules, the predominant subclasses included fatty acids and their conjugates, glycerophosphoserines and eicosanoids (Suppl. material 2: fig. S10b).

To characterise the global metabolic differences across the three developmental stages, PLS-DA was conducted. The first two principal components of the PLS-DA score plot accounted for over 70% of the total variance, clearly distinguishing between developmental stages. Additionally, the R²Y and Q²Y values approached one, confirming the model’s accuracy and reliability (Suppl. material 2: fig. S11a–c). A permutation test was conducted to evaluate the model’s quality, revealing that all red Q² and blue R² values on the left were lower than the original points on the right (Suppl. material 2: fig. S11d–f), which suggests that the model is not subject to overfitting. A total of 1,014 DMs were identified across the three stages (Suppl. material 2: fig. S12). Notably, amino acids were significantly more abundant at the FB stage, likely playing a critical role in the biosynthesis of structural proteins, enzymes and other nitrogen-containing compounds essential for fruiting body development (Suppl. material 2: fig. S13). Conversely, corey lactone diol (a gamma-butyrolactone) was more abundant in the mycelial stages (i.e. LM and WM), particularly at the WM stage and less so at the FB stage. Previous studies have shown that butyrolactones can increase hyphal branching and reduce hyphal extension in fungi like Aspergillus terreus (Schimmel et al. 1998). Thus, this metabolite likely promotes dense mycelial network formation during the mycelial stage, enhancing nutrient absorption and substrate colonisation (Suppl. material 2: fig. S14).

To better understand the pivotal genes and metabolites implicated in the differentiation of fruiting bodies, we analysed the similarities in gene expression and metabolite profiles across the three developmental stages. KEGG co-enrichment analysis revealed that both DMs and DEGs were significantly enriched in the pathways of “Nitrogen metabolism” and “Lysine biosynthesis” when comparing FB to LM (Fig. 5a). In the FB vs. WM comparison, the pathways enriched for both DMs and DEGs included “Starch and sucrose metabolism” and “Biosynthesis of unsaturated fatty acids” (Fig. 5b). Additionally, “Tyrosine metabolism” and “Butanoate metabolism” were enriched in the WM vs. LM comparison (Fig. 5c).

Figure 5. 

KEGG pathway enrichment analysis of DEGs and DMs across the three distinct comparisons. The horizontal axes indicate gene ratios, which reflect the proportion of DEGs/DMs relative to the total number of genes/metabolites associated with a specific KEGG pathway. The size of the data points corresponds to gene counts and the colour represents the -log10 P-value.

For a comprehensive understanding of the transitions associated with fruiting body differentiation, DEGs and DMs from the FB vs. WM comparison were mapped on to the metabolic pathways using iPath. The results were consistent with the GO and KEGG enrichment analyses of DEGs, indicating that pathways related to purine, amino sugar, nucleotide sugar, glycerophospholipid, starch, sucrose, phenylalanine and steroid metabolism and biosynthesis are essential for fruiting body development (Suppl. material 2: fig. S15a). Notably, at the FB stage, pathways for galactose, fructose, mannose, starch and sucrose metabolism were markedly up-regulated. The hydrolysis of macromolecular cellodextrin yielded sucrose and α,α-trehalose, which are readily utilised to provide energy for fruiting body formation (Suppl. material 2: fig. S15b). In contrast, metabolic pathways related to arachidonic acid, linoleic acid and α-linolenic acid metabolism were down-regulated at the FB stage (Suppl. material 2: fig. S15c).

Furthermore, a correlation plot was generated to illustrate the relationship between DEGs and DMs in the FB vs. WM comparison. This analysis identified two distinct groups of genes regulating metabolites. Group 1 consisted of 17 genes (including gene 0752 and gene 1433), exhibiting a similar regulatory pattern, while Group 2 contained the remaining 33 genes displaying an opposing regulatory pattern (Fig. 6a). Group 1 genes were significantly up-regulated at the FB stage (Fig. 6b), whereas Group 2 genes were down-regulated (Fig. 6c). Consequently, the regulation of fruiting body development was predominantly influenced by Group 1 genes. A total of 34 metabolites, positively correlated with Group 1 genes, showed significant increases at the FB stage (Table 2; Suppl. material 2: fig. S16), suggesting their crucial roles in fruiting body development.

Figure 6. 

Correlation analysis between DEGs and DMs in the FB vs. WM comparison: a the figure shows the pairwise correlations between the top 50 DMs and the top 50 DEGs, sorted by P-values from the smallest to the largest. The intensity of the colour indicates the strength of the correlation, with red hues representing strong positive correlations and blue hues indicating strong negative correlations. Correlation values are more pronounced when the ellipses are flatter and statistically significant correlations (P < 0.05) are denoted with an asterisk (*); b, c based on the results shown in panel a, metabolite-regulating genes were categorised into two groups. Group 1 comprises 17 genes that show a positive correlation with 34 metabolites (Table 2), but a negative correlation with the remaining 16 metabolites (b). In contrast, Group 2 genes exhibit an inverse regulatory pattern relative to Group 1 (c).

The authors have declared that no competing interests exist.

No ethical statement was reported.

All the fungal strains used in this study have been legally obtained, respecting the Convention on Biological Diversity (Rio Convention).

This study was supported by the National Natural Science Foundation of China (31872162, 32370014), the Fundamental Research Program of Shanxi Province (202203021221036), the Open Fund Project of Key Laboratory of Microbial Resources Collection and Preservation, Ministry of Agriculture and Rural Affairs (KLMRCP2023-01) and the ‘Sanjin Talents’ Plan of Shanxi Province.

Y.J.Z. and S.Z. designed research; J.N.L. performed research and analysed data; J.N.L., S.Z. and Y.J.Z. wrote the paper. All authors have read and agreed to the published version of the manuscript.

Shu Zhang https://orcid.org/0009-0009-3779-1700

Yong-Jie Zhang https://orcid.org/0000-0002-6331-2189

The genome assembly of Cordyceps blackwelliae ZYJ0835 has been deposited at NCBI with an accession no. JBKZWQ000000000. RNA-seq data have been deposited at NCBI with a BioProject no. PRJNA1214394.