Molecular Mechanisms of Temperature-Regulated Cordycepin Biosynthesis in Cordyceps militaris

TL;DR

This study reveals how temperature affects cordycepin production in Cordyceps militaris by analyzing genes and metabolites involved in its biosynthesis.

Contribution

The study identifies temperature-regulated molecular pathways and genes that influence cordycepin biosynthesis in Cordyceps militaris.

Findings

High temperature (25 °C) increased cordycepin content by 84% in fruiting bodies.

Transcriptomics and metabolomics showed temperature affects membrane function and key metabolic pathways.

The MEblue gene module and MFS genes are strongly linked to cordycepin biosynthesis and transport.

Abstract

Cordycepin is a key active component of Cordyceps militaris, but the molecular mechanism underlying temperature-regulated biosynthesis remains unclear. In this study, Cordyceps militaris strain KN-1 was used as experimental material, with low-temperature (15 °C), control (20 °C), and high-temperature (25 °C) treatments applied during the fruiting body stage. Transcriptomics, untargeted metabolomics, weighted gene co-expression network analysis (WGCNA), and Reverse Transcription quantitative PCR (RT-qPCR) validation were integrated to elucidate the molecular mechanism of temperature-mediated cordycepin biosynthesis. The results showed that 25 °C increased fruiting body cordycepin content by 84%, while 15 °C reduced it. Transcriptomic analysis identified differentially expressed genes (DEGs) enriched in transmembrane transport and fatty acid metabolism, and untargeted metabolomics revealed differential metabolites (DAMs) enriched in lipids and organic acids, indicating that temperature primarily affects Cordyceps militaris membrane function. WGCNA showed that the MEblue module was positively correlated with cordycepin (r = 0.93), with Major Facilitator Superfamily (MFS) members accounting for the highest proportion (47.1%) that may affect cordycepin transmembrane transport. Multi-omics analysis indicated that high temperature promotes cordycepin accumulation through the synergistic regulation of multiple pathways: upregulating genes in the pentose phosphate pathway, purine metabolism, and cordycepin biosynthetic gene cluster (Cns1–Cns3), increasing protective agent pentostatin content, downregulating cordycepin-degrading genes, and enhancing cordycepin transmembrane transport. This study clarifies the molecular mechanism of temperature-mediated cordycepin accumulation, providing a theoretical basis for improving cordycepin production via temperature regulation, optimizing Cordyceps militaris strain quality, and facilitating efficient industrial production.

Figures

Full text

1. Introduction

Cordyceps militaris, a valuable edible and medicinal fungus, is widely distributed worldwide [1,2]. Its bioactive components, including cordycepin [3], cordycepic acid [4], ergosterol [5], and cordyceps polysaccharide [6], endow it with significant application potential. Among these components, cordycepin, the first nucleoside antibiotic isolated from Cordyceps militaris, serves not only as a core quality evaluation index for Cordyceps militaris but also exhibits diverse pharmacological activities such as anti-inflammatory [7], anti-tumor [8], antioxidant [9], and blood glucose metabolism-regulating effects [10]. Cordycepin has been industrially applied in multiple fields, including pharmaceuticals [11], health products [3], functional foods, and cosmetics [12]. However, the scarcity of natural sources and extremely low cordycepin content led to high production costs. These factors severely hinder its market popularization and large-scale application. Therefore, developing efficient and low-cost strategies to enhance cordycepin yield and elucidating the molecular mechanisms underlying its biosynthesis have become research hotspots and key breakthrough directions in this field.

In recent years, the rise in multi-omics technologies has provided powerful tools for elucidating the regulatory mechanisms of cordycepin biosynthesis. Previous studies have identified a cordycepin synthesis pathway associated with 3′-AMP through the exogenous addition of xylose as the carbon source and combined with transcriptomic analysis [13]. By comparing transcriptomic differences between wild-type, high-yield, and low-yield cordycepin strains, the existence of complementary pathways for cordycepin biosynthesis in Cordyceps militaris was inferred [14]. Additionally, whole-genome sequencing has revealed the Cns1–Cns4 gene cluster related to cordycepin synthesis, and the conserved domain proteins encoded by this cluster can significantly promote cordycepin metabolism [15]. These findings have laid a foundation for an in-depth understanding of the cordycepin biosynthesis regulatory network.

Temperature is a key environmental factor influencing the growth, development, and accumulation of bioactive components in edible fungi, and it also serves as a critical process parameter in industrial production. In Ganoderma lucidum cultivation, a temperature of 27 °C significantly increases polysaccharide content in mycelia [16], while 30.1 °C promotes ganoderic acid accumulation in fermentation broth [17]. As one of the major bioactive components in Cordyceps militaris, cordycepin content in fruiting bodies is significantly enhanced by increasing the cultivation temperature and light intensity during the fruiting body culture stage [18]. However, the molecular mechanisms underlying the regulation of cordycepin biosynthesis by temperature, as well as the roles of key regulatory genes and transporter proteins, remain unclear. Therefore, clarifying the molecular mechanism of temperature-mediated cordycepin biosynthesis is of great significance for the efficient and low-cost production of cordycepin.

In this study, Cordyceps militaris strain KN-1 was used as the experimental material. During the fruiting body stage, 20 °C was set as the control, while 15 °C and 25 °C served as low- and high-temperature treatments, respectively. By integrating transcriptomics, untargeted metabolomics, and WGCNA, the molecular mechanism underlying the temperature regulation of cordycepin biosynthesis was systematically elucidated. Additionally, the core pathways and key regulatory factors involved in high-temperature-induced cordycepin accumulation were verified. This study reveals the molecular network of temperature-mediated cordycepin biosynthesis, providing new theoretical insights for improving cordycepin content via temperature optimization strategies, as well as technical support for Cordyceps militaris strain quality improvement and efficient industrial production of cordycepin.

2. Materials and Methods

2.1. Experimental Materials and Culture Media

Strain: The Cordycepin militaris strain KN-1 used in this study was provided by Jiangsu Kangneng Biotechnology Co., Ltd. (Yangzhou, Jiangsu, China). The strain was preserved in 20% glcerol at −80 °C.

Potato Dextrose Agar (PDA) medium: Used for strain activation and maintenance. The medium contained 13 g/L agar, 4 g/L potato infusion powder, 20 g/L glucose, 2 g/L peptone. All components were purchased from Sangon Biotech Co., Ltd. (Shanghai, China). The medium was sterilized at 121 °C for 30 min.

Potato Dextrose Broth (PDB) medium: Used for liquid seed culture. The medium contained 6 g/L potato infusion powder, 20 g/L glucose. All components were purchased from Sangon Biotech Co., Ltd. (Shanghai, China). The medium was sterilized at 121 °C for 30 min.

Cultivation Medium: Used for fruiting body induction and growth. Each culture vessel contained 50 goats and 60 mL deionized water, and were then sterilized at 121 °C for 1 h and cooled to 25 °C prior to inoculation.

2.2. Strain Activation and Fruiting Body Cultivation

The stock culture of Cordyceps militaris preserved in test tubes was inoculated onto the solid PDA medium under aseptic conditions for activation, followed by incubation at 20 °C in the dark for 5 d. After activation, the strain was inoculated into the liquid PDB medium at an inoculum size of 2%, and cultured at 20 °C in the dark for 5 d with shaking at 150 r/min for subsequent use.

A 6 mL aliquot of the liquid seed culture of Cordyceps militaris was inoculated into the cultivation medium and incubated at 20 °C in the dark for 10 d until the mycelium fully covered the medium surface. Subsequently, the culture was exposed to 1200 lx light to induced mycelial color transformation, which was completed after 14 d of total incubation. Primordia emerged on the surface of culture at 20 d post-incubation. Following 27 d of initial cultivation, the cultures were transferred to three temperature conditions for further fruiting body development: low temperature (15 °C), control (20 °C), and high temperature (25 °C). During the temperature treatment phase, the cultures were maintained under 1200 lx light and 80–90% ambient humidity. Fruiting bodies were harvested on 5 d, 10 d, and 15 d post-temperature treatment, corresponding to 32 d, 37 d, and 42 d post-inoculation, respectively. The fruiting bodies cultivated at 15 °C were designated as L-5, L-10, and L-15; those at 20 °C as M-5, M-10, and M-15; and those at 25 °C as H-5, H-10, and H-15.

2.3. Determination of Fruiting Body Yield and Cordycepin Content

Fruiting Body Yield Determination: The harvested fruiting bodies were placed in an oven at 55 °C and dried to a constant weight, followed by being weighed using an electronic balance.

Cordycepin Content Determination: An accurately weighed 0.5 g aliquot of the uniformly crushed sample was placed into a 100 mL volumetric flask, followed by adding approximately 80 mL of water. The flask was placed in an ultrasonic cleaner for ultrasonic extraction for 3 h. After extraction, the volume was made up to 100 mL with water and shaken well. A 1 mL aliquot of the sample solution was centrifuged, and the supernatant was passed through a 0.45 μm microporous membrane. The filtrate was used for high-performance liquid chromatography (HPLC) analysis. The filtrate was subjected to HPLC analysis under the following conditions: chromatographic column: C18 column (250 mm × 4.6 mm, 5 μm); mobile phase: acetonitrile: ultrapure water = 5:95 (v/v); flow rate: 1.0 mL/min (note: corrected from the original 10 mL/min, which is a typo consistent with standard HPLC operation); column temperature: 35 °C; detection wavelength: 260 nm; injection volume: 20 μL. A series of cordycepin standard solutions (Sigma-Aldrich, St. Louis, MO, USA) with concentrations ranging from 1.00 to 50.00 μg/mL was prepared and analyzed under the same chromatographic conditions. A standard curve was constructed with the cordycepin mass concentration as the abscissa and the corresponding peak area as the ordinate, which was then used to calculate the cordycepin content in the samples.

$$[eqn]$$

2.4. Transcriptomic Analysis

2.4.1. Transcriptome Sequencing

Fruiting body samples of Cordyceps militaris subjected to different temperature treatments for various durations were analyzed by sequencing, with three biological replicates per treatment. First, total RNA was extracted from the samples using a universal RNA extraction kit (R401; purchased from Nanjing Jisi Huiyuan Biotechnology Co., Ltd., Nanjing, Jiangsu, China). A Onedrop (OD-1000+) instrument was used to detect RNA purity (OD260/280 ratio), concentration, and the normality of nucleic acid absorption peaks, while RNA quality was assessed by agarose gel electrophoresis. After the samples passed the quality inspection, library construction was performed: eukaryotic mRNA was enriched using Oligo (dT)-coated magnetic beads, and then fragmented randomly by adding Fragmentation Buffer. Using mRNA as the template, the first strand of cDNA was synthesized with random hexamer primers; subsequently, the second strand of cDNA was synthesized by adding buffer, dNTPs, RNase H, and DNA polymerase I. Target fragment sizes were selected using AMPure XP beads, and finally, the cDNA library was obtained by PCR enrichment. After library construction, preliminary quantification was performed using Qubit2.0, and the insert size of the library was detected using Agilent 2100 (Agilent Technologies; Santa Clara, CA; USA). After confirming that the insert size met the expectations, the next step was carried out: the effective concentration of the library was accurately quantified by Q-PCR (>2 nM) to complete the library quality inspection. After passing the library quality inspection, different libraries were mixed according to the target sequencing data output, and sequencing was performed on the Illumina Novaseq X Plus platform (Illumina, Inc; San Diego, CA; USA)with a paired-end read length of 150 bp (PE150). The sequencing work was completed by Nanjing Jisi Huiyuan Biotechnology Co., Ltd.

2.4.2. Transcriptomic Data Analysis

Raw sequencing data (Raw Data) were filtered to remove adapter sequences and low-quality reads, yielding high-quality Clean Data. The Clean Data were aligned to the specified Cordyceps militarism strain KN-1 reference genome http://daehwankimlab.github.io/hisat2 (accessed on 10 September 2025) to obtain Mapped Data. Subsequently, sequencing library quality assessment was performed, including insert fragment length inspection and randomness test. Fragments Per Kilobase of transcript per Million mapped reads (FPKM) was used as the indicator to measure the expression level of transcripts or genes. Correlation coefficients between samples were calculated based on gene expression levels in the samples, and the correlation between samples was evaluated using Pearson’s correlation coefficient and Principal Component Analysis (PCA). Differentially Expressed Genes (DEGs) were screened using DESeq2 (v1.26.0) software, with the screening criteria of absolute value of Log_2_ (Fold Change) > 1 and false discovery rate (FDR) < 0.05, laying a foundation for subsequent functional annotation and functional enrichment analysis of DEGs.

2.5. Weighted Gene Co-Expression Network Analysis (WGCNA)

The weighted gene co-expression network analysis (WGCNA) method was employed for systematic transcriptomic analysis. First, the expression matrix (TPM values) was used to screen for genes with significant differences in M_VS_H and M_VS_L, as well as 5_VS_10 and 10_VS_15, under the three temperatures, and a gene co-expression matrix was constructed. Pearson correlation coefficients between different genes were calculated to establish a weighted network, and an appropriate soft threshold β of 20 was determined to ensure the network conformed to the scale-free property. Subsequently, the topological overlap matrix (TOM) was used to measure the connection relationships between genes, and hierarchical clustering was applied to classify the genes into 7 co-expression modules (gray represents genes not assigned to any co-expression module). The module classification of specific genes is shown in Supplementary Table S1: WGCNA modules. The expression pattern of each module was represented by the Module Eigengene (ME), and further correlation analysis was performed between MEs and the phenotypic traits of the samples to identify modules significantly associated with the target traits. In this study, the metabolite contents of the corresponding samples in the cordycepin metabolic pathway were used as phenotypes to calculate the correlation between module genes and metabolites.

2.6. Untargeted Metabolomic Detection and Analysis

2.6.1. Untargeted Metabolomic Detection

Fruiting body samples of Cordyceps militaris subjected to different temperature treatments for various durations were analyzed by sequencing, with three biological replicates per treatment. The sequencing work was completed by Nanjing Jisi Huiyuan Biotechnology Co., Ltd. (Nanjing, Jiangsu, China). Briefly, the samples were freeze-dried, and 25 mg ± 1 mg of uniformly crushed solid sample was accurately weighed into an EP tube under low-temperature conditions, followed by the addition of homogenization beads and 100 μL of extraction solution (methanol: acetonitrile: water = 2:2:1, v/v/v) containing isotope-labeled internal standards; after vortex mixing for 30 s, the sample was homogenized at 35 Hz for 4 min in a homogenizer and then subjected to ultrasonic extraction in an ice-water bath for 5 min, with this homogenization–ultrasonic extraction cycle repeated three times before incubating the sample at −40 °C for 1 h. A 400 μL aliquot of the supernatant was transferred to a 0.22 μm filter plate well, and the protein precipitation plate–collection plate assembly was placed in a positive pressure device to collect the filtrate under 6 psi pressure for 120 s for subsequent instrumental analysis. Quality control (QC) samples were prepared by mixing equal volumes of supernatants from all experimental samples and analyzed alongside the experimental samples to evaluate the stability and reproducibility of the experiment. Chromatographic separation was conducted using an ultra-high-performance liquid chromatography (UHPLC) system equipped with a Phenomenex Kinetex C18 column (2.1 mm × 50 mm, 2.6 μm), with the mobile phase consisting of phase A (aqueous phase containing 0.01% acetic acid) and phase B (isopropanol: acetonitrile = 1:1, v/v), the sample tray temperature maintained at 4 °C, and the injection volume set to 2 μL.

2.6.2. Untargeted Metabolomic Data Analysis

Raw mass spectrometry data were converted to mzXML format using Proteowizard software (v3.0.24054). Principal Component Analysis (PCA) of quality control (QC) samples and experimental samples was performed using the R packageropls (v3.6.2), while Pearson correlation coefficient analysis and Relative Standard Deviation (RSD) analysis of QC samples were conducted using the R package Rcorrplot (v3.6.2). Raw data were filtered using a self-developed Perl program by the company to remove data without definite substance names and without spectral matching similarity; substances with a missing rate >50% in the comparison groups were directly filtered out, and those with a missing rate <50% were imputed for missing values using the K-Nearest Neighbor (KNN) algorithm of the R package DMwR (v3.6.2). Subsequently, the Perl program was used to normalize the data based on internal standard (IS) or total ion current (TIC) of samples. After preprocessing the raw data, metabolite annotation was performed using the self-developed Perl program by the company based on the Human Metabolome Database (HMDB, v5.0) and KEGG COMPOUND database https://www.kegg.jp/kegg/compound/ (accessed on 10 September 2025). Following metabolite annotation, multivariate statistical analyses including normalized PCA and Orthogonal Partial Least Squares–Discriminant Analysis (OPLS-DA) were conducted on the data using R packages (v3.6.2); univariate statistical analyses, including Student’s t-test and Fold Change (FC) calculation, were then performed using the self-developed Perl program by the company. Finally, differential metabolites were screened by combining the results of multivariate and univariate statistical analyses, with the screening criteria of p-value < 0.05 and Variable Importance in Projection (VIP) >1. Among them, metabolites with Fold Change >1 were upregulated differential metabolites, and those with Fold Change <1 were downregulated differential metabolites. (All companies involved in this section refer to Nanjing Jisi Huiyuan Biotechnology Co., Ltd.).

2.7. Correlation Analysis of Genes and Metabolites Related to Cordycepin Biosynthesis

Genes and metabolites related to cordycepin biosynthesis were selected, and Pearson correlation coefficients were calculated based on the gene expression levels (FPKM) from transcriptomic data and the relative contents of metabolites to obtain the correlation between transcriptomic and metabolomic data.

2.8. RT-qPCR Validation

The same RNA samples as those used in Section 2.2 were employed for RT-qPCR validation. An appropriate amount of RNA (800 ng inthis experiment) was taken, and RNase-free ddH_2_O was added to a final volume of 12 μL. Then 3 μL of 5 × gDNA digester mix was added, gently mixed by pipetting, and incubated at 42 °C for 2 min. To the mixture from the previous step, 5 μL of 4 × Hifair^®^ III SuperMix plus was added and gently mixed by pipetting to obtain a 20 μL reaction system. The PCR program was set as follows: 25 °C for 5 min, 55 °C for 15 min, and 85 °C for 5 min. The product was cDNA, which could be immediately used for PCRs. The amplification conditions for two-step real-time fluorescent quantitative PCR were set as follows: amplification curve: 95 °C for 5 min (1 cycle); 95 °C for 10 s and 60 °C for 30 s (40 cycles); melting curve: 95 °C for 15 s, 60 °C for 60 s, and 95 °C for 15 s, with continuous signal detection.

2.9. Data Analysis

Each sample was subjected to three replicate experiments, and the results were expressed as mean ± standard deviation (SD). Data were analyzed using GraphPad Prism 8.0 (Northampton, MA, USA). Inter-group differences were analyzed by Duncan’s multiple range test using IBM SPSS Statistics 25.0 (SPSS Inc., Chicago, IL, USA). A p-value < 0.05 was considered statistically significant.

3. Results

3.1. Growth Characteristics and Cordycepin Content of Cordyceps militaris Under Different Temperature Treatments

Compared with the control group (20 °C), the fruiting body yield, color, and cordycepin content of Cordyceps militaris cultured at high temperature (25 °C) and low temperature (15 °C) all changed. Under low-temperature conditions, the color of Cordyceps militaris fruiting bodies was basically consistent with that of the control group, both showing a typical orange-yellow color, while under high-temperature conditions, the color of the fruiting bodies faded to light yellow (Figure 1A). Regarding yield, there was no significant difference in the fruiting body yield of Cordyceps militaris under low-temperature culture compared with the control group, whereas the fruiting body yield under high-temperature culture was significantly reduced by 10.58% compared with the control group (Figure 1B). The cordycepin content in fruiting bodies under different temperature conditions was determined by high-performance liquid chromatography (HPLC). The results showed that the cordycepin content in fruiting bodies under low-temperature culture for 5 d, 10 d, and 15 d was decreased by 46.40%, 56.33%, and 37.82% compared with the control group, respectively. In contrast, the cordycepin content in fruiting bodies under high-temperature culture for 5 d, 10 d, and 15 d was significantly increased by 42.53%, 40.90%, and 84.00% compared with the control group, respectively (Figure 1C). Although the fruiting body yield in the high-temperature group decreased by 10.58% compared with the control group, the total cordycepin yield in the high-temperature group was 6.87 ± 0.63 mg/vessel, which was 67.66% higher than that of the control group (11.52 ± 1.16 mg/vessel). This confirms that the significant increase in cordycepin concentration outweighs the 10.58% reduction in biomass, resulting in a substantial net gain in total productivity.

3.2. Transcriptomic Differential Analysis of Genes Related to Cordycepin Biosynthesis Under Different Treatments

Transcriptome sequencing was performed on Cordyceps militaris fruiting bodies under different temperature treatments for various durations, yielding a total of 197 Gb of data. After removing adapter sequences and low-quality reads, high-quality clean reads were obtained, with Q30 values > 95% and mapped reads to the reference genome exceeding 98%, ensuring the reliability of subsequent analyses. A heatmap of Pearson correlation coefficient analysis showed that the correlation coefficients between samples were all greater than 0.95, indicating good correlation (Supplementary Figure S1). Principal Component Analysis (PCA) revealed that different treatment temperatures and durations were well separated on PC1 and PC2, respectively (Figure 2A).

Based on gene-level read count data, differential expression analysis was performed using DESeq2 to calculate the Fold Change (log_2_ FC) in gene expression. DEGs were screened with the criteria of |log_2_ FC| > 1 and Q value < 0.05 after false discovery rate correction via the Benjamini–Hochberg method. Genes were compared under different temperature treatments at the same stage to reveal the expression characteristics of genes in response to temperature changes. Under low-temperature treatment, the number of DEGs increased with prolonged treatment duration. A total of 267 DEGs (196 upregulated, 71 downregulated) were identified in the M_5 VS L_5 comparison group, 114 DEGs (43 upregulated, 71 downregulated) in the M_10 VS L_10 comparison group, and 364 DEGs (214 upregulated, 150 downregulated) in the M_15 VS L_15 comparison group. Under high-temperature treatment, the number of DEGs also increased with extended treatment time. The M_5 VS H_5 comparison group yielded 131 DEGs (77 upregulated, 54 downregulated), the M_10 VS H_10 group 316 DEGs (75 upregulated, 241 downregulated), and the M_15 VS H_15 group 447 DEGs (207 upregulated, 240 downregulated). Given the focus on temperature-induced changes in gene expression, a non-redundant statistical approach was applied to analyze DEGs across different time points under various temperature treatments. A total of 1070 unique DEGs were identified across the six comparison groups. Specifically, 633 DEGs were detected in the M_ VS L group, while 664 DEGs were found in the M VS _H group, indicating that high-temperature treatment elicits a more pronounced transcriptional response. Further comparison revealed 406 and 437 DEGs that were uniquely differentially expressed under low-temperature and high-temperature stress, respectively. Additionally, 227 genes exhibited differential expression under both temperature stress conditions, suggesting the coexistence of temperature-specific regulatory mechanisms and universal stress response pathways. The statistics of DEGs among different groups are presented in Figure 2B, and detailed summaries of the differential expression analyses are provided in Supplementary Table S2: Summary of Differentially Expressed Genes.

3.3. GO and KEGG Enrichment Analysis of Differentially Expressed Genes

To comprehensively characterize the attributes of genes and their products in organisms, Gene Ontology (GO) classification and functional enrichment analysis were performed on the DEGs. A total of 131, 196, and 65 DEGs were annotated to biological process (BP), molecular function (MF), and cellular component (CC) categories, respectively. To further validate the GO annotation results, functional clustering analysis was conducted on DEGs. In the BP category, DEGs were significantly enriched in proteolysis, transmembrane transport, and rRNA processing (Figure 3A). In the MF category, DEGs were significantly enriched in transmembrane transporter activity, heme binding, and serine-type carboxypeptidase activity (Figure 3B). In the CC category, DEGs were significantly enriched in the extracellular region, membrane, and nucleus (Figure 3C). Integrating the results of GO enrichment analysis across BP, MF, and CC categories indicated that genes with transmembrane transport functions and transmembrane transporter activity localized to the membrane were significantly regulated by temperature. These fingdings suggest that temperature primarily affects membrane function, thereby influencing cordycepin content.

In organisms, different genes coordinate with each other to exert their biological functions. To better analyze the relationship between gene functions based on biological pathways, the identified DEGs were uploaded to the KEGG database for function annotation, with a significant threshold of p < 0.05 for enriched pathways. A total of 67 DEGs were enriched in 19 secondary pathways (Figure 4). In the M_VS_H comparison group, DEGs were significantly enriched in pathways including pentose and glucuronate interconversions, fatty acid metabolism, amino sugar and nucleotide sugar metabolism, and fatty acid biosynthesis. In the M_VS_L comparison group, DEGs were significantly enriched in ribosome biogenesis in eukaryotes, nitrogen metabolism, biosynthesis of secondary metabolites, and alanine, aspartate and glutamate metabolism. KEGG enrichment analysis showed that membrane-related fatty acid metabolism and fatty acid biosynthesis pathways were significantly upregulated in the M_VS_H group, whereas these two pathways were not significantly enriched in the M_VS_L group. Fatty acid metabolism and fatty acid biosynthesis can directly regulate the composition and structural stability of cell membrane lipids, and changes in membrane lipid properties are key factors affecting transmembrane substance transport [19]. Combined with the GO enrichment analysis results, these findings suggest that temperature may primarily affect intracellular transmembrane transport, thereby regulating cordycepin content.

3.4. Differential Analysis of Metabolite Responses to Temperature in Cordyceps militaris Fruiting Bodies

In this study, untargeted metabolomic analysis was employed to systematically identify changes in differential metabolites (DAMs) in Cordyceps militaris fruiting bodies under different temperature treatments for various durations. Quality control and reproducibility analysis confirmed the stability of the instrument, ensuring the reliability and reproducibility of the metabolomic data. PCA results showed that each treatment group and the control group clustered into distinct clusters with clear separation (Figure 5A). Meanwhile, Orthogonal Partial Least Squares–Discriminant Analysis (OPLS-DA) indicated significant separation between any two comparison groups (Supplementary Figures S2–S7). The above PCA and OPLS-DA results demonstrated that the data had good reproducibility and credibility, which could be used for subsequent analyses.

In this study, t-test, OPLS-DA, Variable Importance in Projection (VIP) values, and Fold Change were used to screen for DAMs (upregulated: p-value < 0.05, VIP > 1, Fold Change > 1; downregulated: p-value < 0.05, VIP > 1, Fold Change < 1). A total of 690, 706, and 758 DAMs were identified in the M_5 VS H_5, M_10 VS H_10, and M_15 VS H_15 comparison groups, respectively; 561, 671, and 604 DAMs were identified in the M_5 VS L_5, M_10 VS L_10, and M_15 VS L_15 comparison groups, respectively (Figure 5D). The number of DAMs in the high-temperature treatment group increased with the extension of treatment time, indicating that long-term high-temperature stress had a significant impact on the metabolites of Cordyceps militaris. A total of 1113 DAMs were identified in the M_VS_H comparison group (Figure 5B). Among these 1113 DAMs, 56.8% were unannotated, and the annotated metabolites with a proportion greater than 1% were lipids and lipid-like molecules (14.9%), organic acids and derivatives (14.9%), organoheterocyclic compounds, phenylpropanoids and polyketides (4.1%), organic oxygen compounds (3.5%), benzenoids (2.9%), alkaloids and derivatives (1.4%), and nucleosides, nucleotides, and analogs (1.3%). A total of 1007 DAMs were identified in the M_VS_L comparison group (Figure 5C). Among these 1007 DAMs, 58.1% were unannotated, and the annotated metabolites with a proportion greater than 1% were lipids and lipid-like molecules (15.0%), organic acids and derivatives (6.6%), organoheterocyclic compounds (5.0%), phenylpropanoids and polyketides (4.7%), benzenoids (3.4%), organic oxygen compounds (3.4%), alkaloids and derivatives (1.4%), and nucleosides, nucleotides, and analogues (1.4%). Classification and enrichment analysis of overall DAMs showed that lipids and lipid-like molecules, as well as organic acids and derivatives, accounted for the highest proportion in the M_VS_H comparison group. This is consistent with the significant enrichment of genes in fatty acid metabolism and fatty acid biosynthesis pathways in the KEGG analysis (Section 3.3), indicating that high temperature may alter cell membrane structure or energy supply by regulating lipid metabolism, thereby providing support for the transmembrane transport process related to cordycepin biosynthesis.

To explore the functions of metabolites under high-temperature treatment, functional annotation of DAMs in each high-temperature comparison group was performed using the KEGG database (Supplementary Figures S8–S13). In the M_5 VS H_5 comparison group, DAMs were significantly enriched in glucosinolate biosynthesis, Parkinson’s disease, and isoflavonoid biosynthesis. Among these, both glucosinolate biosynthesis and isoflavonoid biosynthesis are pathways related to plant secondary metabolism, while secondary metabolism in fungi often involves crosstalk among multiple metabolic pathways [20]. Therefore, the significant enrichment of glucosinolate biosynthesis and isoflavonoid biosynthesis pathways may indirectly alter the carbon metabolism flux in Cordyceps militaris, providing precursor substances for cordycepin biosynthesis. In the M_10 VS H_10 comparison group, DAMs were significantly enriched in biosynthesis of alkaloids derived from shikimate pathway and biosynthesis of various alkaloids. The shikimate pathway is a core pathway for the synthesis of aromatic amino acids and alkaloids [21], and the metabolic transformation of cordycepin synthesis precursors relies on amino acid metabolism to provide nitrogen sources [22]. This indicates that after 10 d of high-temperature stress, the alkaloid synthesis pathway may be enhanced to provide precursor substances for cordycepin biosynthesis. In the M_15 VS H_15 comparison group, DAMs were significantly enriched in glucosinolate biosynthesis and the biosynthesis of alkaloids derived from the shikimate pathway. This result suggests that under long-term, high-temperature stress, Cordyceps militaris does not randomly activate metabolic pathways but stably maintains the activity of specific secondary metabolic pathways to continuously improve the supply of precursor substances for cordycepin biosynthesis. This is consistent with the stable regulation law of secondary metabolism formed by fungi under long-term environmental stress [23].

3.5. WGCNA Related to Cordycepin Biosynthesis

To systematically clarify the synergistic regulation pattern of gene expression and its intrinsic correlation with phenotypes related to cordycepin biosynthesis, WGCNA was performed in this study. Combined with the phenotypic data of key metabolite contents such as adenosine, cordycepin, AMP, 3′-AMP, adenylosuccinate, and adenine, a co-expression regulatory network of key genes involved in cordycepin biosynthesis in Cordyceps militaris was constructed. After module hierarchical clustering and merging, a total of seven independent modules was obtained (Figure 6), providing clear targets for the subsequent screening of core regulatory modules and key genes.

Phenotype–module correlation analysis revealed that metabolites associated with cordycepin biosynthesis exhibited distinct correlations with different modules. Cordycepin was negatively correlated with adenosine, AMP, and 3′-AMP, indicating that these three metabolites, as precursor for cordycepin biosynthesis, may reduce their own contents through directional conversion to cordycepin. Among the modules, the MEblue module showed the strongest correlation with cordycepin content (p = 0.000006, r = 0.93), while the MEblack module (p = 0.002, r = −0.58) and Mered module (p = 0.01, r = −0.51) were significantly negatively correlated with cordycepin. This suggests that the MEblack and Mered modules may be involved in processes that inhibit cordycepin biosynthesis or promote cordycepin degradation. As a precursor of cordycepin biosynthesis, adenosine was only weakly positively correlated with the MEgreen module (p = 0.05, r = 0.35) and weakly negatively correlated with the MEbrown module (p = 0.07, r = −0.35), with no significant correlation with other modules. This indicates that the metabolic transformation of adenosine may depend on the regulation of specific key genes rather than the synergistic expression at the module level. As a core intermediate product of purine metabolism, AMP was significantly positively correlated with the MEred module (p = 0.02, r = 0.46) and weakly negatively correlated with the MEyellow module (p = 0.2, r = −0.28), suggesting that the MEred module may participate in the precursor supply for cordycepin biosynthesis by regulating the production or transformation of AMP. Collectively, the module correlation results indicated that the MEblue module is a key module regulating cordycepin biosynthesis, and its significant positive correlation with cordycepin content provides a critical direction for the subsequent mining of key genes related to cordycepin biosynthesis and elucidation of the molecular mechanism underlying temperature-regulated cordycepin biosynthesis.

To determine the dynamic change trends of 307 DEGs in the MEblue module, GO enrichment analysis was performed on MEblue module DEGs associated with cordycepin content (Figure 7). Among these DEGs, 69 were significantly enriched BP, along with 99 in MF and 16 CC. In the BP category, DEGs were significantly enriched in transmembrane transport (28 DEGs, accounting for 40.58%), proteolysis (10 DEGs, accounting for 14.49%), and the carbohydrate metabolic process (7 DEGs, accounting for 10.14%). In the MF category, DEGs were significantly enriched in transmembrane transporter activity (23 DEGs, accounting for 23.23%), O-methyltransferase activity (3 DEGs, accounting for 3.03%), sequence-specific DNA binding (3 DEGs, accounting for 3.03%), serine-type peptidase activity (5 DEGs, accounting for 5.05%), oxidoreductase activity (5 DEGs, accounting for 5.05%), iron ion binding (6 DEGs, accounting for 6.06%), and iron–sulfur cluster binding (2 DEGs, accounting for 2.02%). In the CC category, DEGs were only significantly enriched in the membrane (13 DEGs, accounting for 81.25%). These results indicate that the increase in cordycepin content under temperature stress is closely associated with transmembrane transport proteins.

3.6. Putative Transporters Related to Cordycepin Biosynthesis

To further analyze the molecular mechanism underlying the increase in cordycepin content under temperature stress, DEGs from the MEblue module were uploaded to the Transporter Classification Database TCDB, https://tcdb.org/browse.php (accessed on 30 September 2025) for annotation. Among the 307 DEGs in the blue module, 34 were annotated as transporters, accounting for 11.07%. Detailed information on these genes is provided in Supplementary Table S3: Transporters in the blue module. Statistical analysis of the family classification of these transporters showed that the 2.A.1 family had the highest proportion (47.1%), followed by the 3.A.1 and 8.A.47 families (each accounting for 5.9%) (Table 1). The 2.A.1 family belongs to the Major Facilitator Superfamily (MFS), which can catalyze the transport of secondary metabolites. It is speculated that these MFS transporters are associated with the increased cordycepin content under high-temperature stress.

3.7. Regulatory Relationships Between Genes and Metabolites Related to Cordycepin Biosynthesis

To further analyze the potential mechanism underlying the increase in cordycepin content under temperature stress, 22 key genes related to cordycepin biosynthesis were extracted in this study, including those involved in purine metabolism and the pentose phosphate pathway. The key cordycepin biosynthetic gene cluster Cns1~Cns4 was renamed as Cordycepin synthesis. A gene–metabolite correlation expression matrix of the cordycepin metabolic pathway was constructed using these 22 key genes and eight metabolites related to cordycepin biosynthesis (Figure 8). The result showed that in the cordycepin biosynthetic pathway, the gene with the highest correlation with cordycepin is A9K55_004233. In the pentose phosphate pathway, the genes highly correlated with cordycepin are A9K55_001711 and A9K55_005654, followed by A9K55_006418, while the gene negatively correlated with adenylosuccinate and pentostatin is A9K55_001765. Adenosine, the end product of purine metabolism, has low correlation with all enzymes in its synthetic pathway. In the purine metabolism pathway, the gene with the highest correlation with cordycepin is A9K55_006221. Adenine, an intermediate metabolite of the purine metabolism pathway related to cordycepin biosynthesis, has the highest correlation with A9K55_008988 in the purine metabolism pathway, followed by A9K55_000442 and A9K55_008038. Adenylosuccinate, another intermediate metabolite of the purine metabolism pathway related to cordycepin biosynthesis, has the highest correlation with A9K55_000049 in the purine metabolism pathway, followed by A9K55_008988. Pentostatin, a metabolite related to protecting cordycepin from degradation, has the highest correlation with A9K55_000049 and A9K55_008988 in the purine metabolism pathway related to cordycepin biosynthesis, followed by A9K55_000442.

3.8. Putative Metabolic Pathway of Cordycepin Under Temperature Stress

Given that high-temperature stress can significantly increase cordycepin content in Cordyceps militaris, this study systematically identified DEGs and DAMs associated with increased cordycepin content through KEGG enrichment analysis, WGCNA, and gene–metabolite correlation analysis. Based on the previous findings by Zheng [24] and Vongsangnak [25] that cordycepin biosynthesis mainly relies on the pentose phosphate pathway and purine metabolism pathway, as well as the research foundation of Yu [26] on the key cordycepin biosynthetic gene cluster Cns1~Cns3, this study further predicted the metabolic pathway related to cordycepin biosynthesis under high-temperature stress (Figure 9).

The pentose phosphate pathway provides core precursor substances for cordycepin biosynthesis. Key genes in this pathway including fructose-1,6-bisphosphatase (A9K55_001765), transketolase (tktA, A9K55_001711), transaldolase (TALDO1, A9K55_005654), and ribose-phosphate pyrophosphokinase (RPPS, A9K55_006418) were all upregulated. This promotion of phosphoribosyl pyrophosphate (PRPP) production supplies sufficient precursors for cordycepin biosynthesis. Concurrently, the synchronous upregulation of nucleoside diphosphate kinase (A9K55_007777) in the energy metabolism pathway effectively maintains intracellular ATP homeostasis, providing adequate energy support for the aforementioned biosynthetic processes.

In the purine metabolism pathway, critical genes such as phosphoribosylpyrophosphate amidotransferase (A9K55_008988), GAR transformylase (GART, A9K55_000442), adenylosuccinate lyase (purB, A9K55_003166), adenylosuccinate synthetase (purA, A9K55_000049 and A9K55_008138), and adenosine kinase (ADK, A9K55_003166) were upregulated, mediating the cascade catalysis of PRPP to AMP. Metabolomic data showed that adenylosuccinate content gradually decreased under temperature stress, indicating that more adenylosuccinate was converted to AMP to fuel cordycepin synthesis. Adenine phosphoribosyltransferase (A9K55_007511) was upregulated under high-temperature stress, catalyzing the conversion of adenine to AMP and providing additional substrates for cordycepin biosynthesis. Metabolomic results revealed that adenine content increased at 5 d and 10 d of high-temperature treatment but gradually decreased with prolonged exposure, suggesting that high temperature promotes adenine accumulation and most of the accumulated adenine is catalytically converted to AMP to support cordycepin production. Additionally, AMP content increased at 5 d of high-temperature stress but decreased at 10 d and 15 d, indicating enhanced conversion of AMP to adenosine. This is consistent with the elevated adenosine levels observed in metabolomic data under high-temperature conditions, thus providing more abundant precursors for cordycepin biosynthesis.

In the key cordycepin biosynthetic gene cluster, Cns3 (A9K55_004233) was upregulated under high-temperature stress, catalyzing the conversion of adenosine to 3′-AMP. Notably, Cns3 has been reported to be involved in the biosynthesis of pentostatin, a metabolite that protects cordycepin from degradation [27]. Metabolomic data showed that pentostatin content increased at 5 d and 10 d of high-temperature stress but gradually decreased with extended treatment, indicating increased utilization of pentostatin to safeguard intracellular cordycepin from degradation. Cns2 (A9K55_004234) was upregulated under high-temperature stress, catalyzing the conversion of 3′-AMP to 2′-C-3′-dA. Cns1 (A9K55_004235) was upregulated at 5 d of high-temperature stress, catalyzing the conversion of 2′-C-3′-dA to cordycepin. Adenosine deaminase (A9K55_002334) can degrade intracellular cordycepin to protect cells from the toxicity of high cordycepin concentrations [28]. However, its expression was continuously downregulated under high-temperature stress, leading to reduced intracellular cordycepin degradation and subsequent accumulation of cordycepin.

3.9. RT-qPCR Validation of Gene Expression Patterns

To verify the reliability of transcriptomic data, upregulated transporter genes A9K55_005911, A9K55_008619, and A9K55_007015 in the MEblue module and genes related to the cordycepin biosynthetic gene cluster A9K55_004233, A9K55_004234, and A9K55_004235 were randomly selected. RT-qPCR assays were performed using Cordyceps militaris fruiting bodies subjected to different temperature treatments for various durations. The results showed that the gene expression trends obtained from RT-qPCR (relative expression levels) were consistent with those from RNA-seq data (FPKM values) (Figure 10), confirming the reliability of the transcriptomic data.

4. Discussion

Cordycepin is an important active substance in Cordyceps militaris, and improving cordycepin content has long been a research focus in related fields. Currently, numerous methods have been developed to enhance cordycepin production, such as adding exogenous substances in liquid fermentation [29], optimizing medium formulations for solid-state fermentation [30], ultraviolet mutagenesis of high-yield strains [31], and improving culture conditions [32]. Among these, temperature regulation, as a key strategy for optimizing culture conditions, has been confirmed to affect cordycepin content in Cordyceps militaris. However, existing studies have only clarified the relationship between temperature and cordycepin content, without in-depth exploration of the underlying molecular regulatory mechanisms, nor have they defined the roles of key regulatory factors and transport processes. In this study, we increased cordycepin content in fruiting bodies by altering the fruiting body cultivation temperature. Combined with multi-omics technologies, we verified the effect of high temperature on cordycepin content and simultaneously revealed the regulatory network of cordycepin biosynthesis under temperature influence, thus addressing the gaps in current research.

Metabolomic data from this study showed that the contents of 3-methylxanthine (H-5: 0.72, M-5: 0.79; H-10: 0.53, M-10: 0.74; H-15: 0.68, M-15: 0.41) and 7-methylxanthine (H-5: 0.72, M-5: 0.79; H-10: 0.53, M-10: 0.74; H-15: 0.68, M-15: 0.41) significantly decreased at 10 d and 15 d of high-temperature treatment, indicating that the purine catabolic pathway may be inhibited. Adenine and adenosine, as precursor substances for cordycepin biosynthesis, were not decomposed into methylated products but were utilized for cordycepin synthesis. Colchicine can induce fungal differentiation, enhance antibacterial activity, and promote cordycepin biosynthesis, thereby increasing cordycepin content [33]. The content of colchicine significantly decreased at 5 and 10 d of high-temperature treatment, suggesting that more colchicine was utilized to promote the increase in cordycepin content after high-temperature treatment, which is consistent with the study by Yu [34]. Amino acids, the fundamental building blocks of proteins, are closely associated with the biological activities of fungi. The content of proline significantly decreased under high-temperature stress, and the content of tryptophan significantly decreasing at 10 and 15 d of high-temperature stress indicated that proline and tryptophan were used for cordycepin biosynthesis. This result is consistent with the study by Chen [35]. Rotenone is a respiratory inhibitor that can interfere with amino acid metabolism and nucleotide metabolism to promote cordycepin accumulation [36]. Its content significantly increased at 5 and 10 d of high-temperature treatment, which may lead to an increase in cordycepin content in Cordyceps militaris fruiting bodies, consistent with the study by Ma [37]. Ergosterol is another important active component in Cordyceps militaris, with anti-cancer [38] and anti-inflammatory activities [39]. Its content significantly increased at 10 d of high-temperature treatment, indicating that high-temperature treatment can also increase other active components in Cordyceps militaris, providing a new direction for the comprehensive quality improvement of Cordyceps militaris.

The typical orange-yellow color of Cordyceps militaris fruiting bodies is mainly associated with carotenoids [40]. In this study, high-temperature stress induced the fading of fruiting body color to pale yellow, implying a reduction in carotenoid content under high-temperature conditions. The conjugated double bonds in the molecular structure of carotenoids are susceptible to oxidation by reactive oxygen species (ROS), generating colorless oxidation products such as aldehydes and ketones [41] and thereby impairing their coloring capacity. Thus, it is hypothesized that high-temperature stress increases intracellular ROS levels in fruiting bodies, leading to carotenoid degradation and the subsequent color change from orange-yellow to pale yellow. Fungal secondary metabolites often act as adaptive mediators in response to environmental stress [42]. Cordycepin has been confirmed to function as a stress-protective molecule in Cordyceps militaris, and its antioxidant activity helps alleviate oxidative damage caused by elevated ROS levels during high-temperature stress [43]. ROS-induced lipid peroxidation can disrupt membrane integrity [44], and the synergistic upregulation of cordycepin biosynthesis and fatty acid metabolism under high temperature may constitute a coordinated protective mechanism. KEGG enrichment analysis revealed significant upregulation of fatty acid metabolism and biosynthesis pathways in the high-temperature group, which can modulate membrane lipid composition and stability. Meanwhile, untargeted metabolomics confirmed that lipids and lipid-like molecules were the most abundant differential metabolites, indicating active membrane remodeling by cells in response to heat stress. This membrane adaptation is crucial for maintaining cellular homeostasis under thermal stress, as high temperature may alter the ratio of saturated to unsaturated fatty acids and disrupt membrane fluidity [45]. Cordycepin may preserve membrane integrity by scavenging ROS and reducing lipid peroxidation. The dual regulation of membrane lipid remodeling and cordycepin-mediated antioxidant defense helps Cordyceps militaris cope with high-temperature stress, minimizing cellular damage and ensuring the maintenance of basic metabolic activities.

Transcription factors play a pivotal role in activating secondary metabolic pathways, mediating environmental stress responses, and regulating developmental processes in fungi [46,47]. As a class of key regulatory factors to fungi, Zn_2_Cys_6_ binuclear cluster transcription factors have been demonstrated to respond to temperature stress, modulate the biosynthesis of secondary metabolites, and participate in the regulation of fruiting body morphogenesis and developmental processes [48,49,50]. Chen [51] found that overexpressing two key Zn_2_Cys_6_ transcription factors under L-alanine stress significantly increased cordycepin content in the mycelia of Cordyceps militaris. Yang [52] revealed that Zn_2_Cys_6_ transcription factors are involved in regulating light-responsive genes in Cordyceps militaris; the light-responsive gene WC-1 can modulate the expression of adenylosuccinate synthase, thereby affecting the production of adenosine, which serves as a key precursor for cordycepin biosynthesis, and ultimately regulating cordycepin accumulation. In addition, Wang [53] conducted a genome-wide identification and expression analysis of Coriolopsis trogii under high-temperature stress, showing that Zn_2_Cys_6_ transcription factors were significantly enriched to respond to high-temperature stress. Consistent with these findings, the present study identified three Zn_2_Cys_6_ transcription factors, A9K55_007036, A9K55_007435, and A9K55_008577, that were significantly upregulated under high-temperature stress. These results suggest that Zn_2_Cys_6_ transcription factors may enhance cordycepin accumulation in Cordyceps militaris fruiting bodies under high-temperature stress by regulating specific secondary metabolic pathways and promoting the biosynthesis of cordycepin precursors.

Although existing studies have focused on the effect of temperature on cordycepin content, the transport mechanism of cordycepin synthesis and the mechanism by which temperature regulates cordycepin synthesis remain unclear. Cordycepin is synthesized in intracellular lipid droplets, and excessive intracellular cordycepin accumulation is toxic to the fungus itself. Therefore, it is necessary to excrete cordycepin extracellularly to mitigate self-toxicity. By comparing the expression of differential genes between low-yield and high-yield cordycepin strains, Zhang [14] speculated that cordycepin may be excreted extracellularly via the ATP-binding cassette (ABC) transporter encoded by Cns4 and the transporter encoded by CMM_04722, thereby maintaining a low intracellular cordycepin concentration in the cytoplasm and further promoting intracellular cordycepin biosynthesis. However, Cns4 has been reported to function as a transporter responsible for the extracellular excretion of pentostatin. As a cordycepin protectant, no transcriptional difference in Cns4 was observed among the control groups in this study. Meanwhile, pentostatin content first increased and then decreased with the extension of high-temperature treatment, indicating that more intracellular pentostatin was utilized to protect cordycepin from degradation. Chen [51] significantly increased cordycepin content by adding L-alanine to the Cordyceps militaris fermentation broth. Transcriptomic analysis of mycelia showed that differential genes were significantly enriched in the Major Facilitator Superfamily (MFS) transporters, which is consistent with the significant enrichment of MFS family transporters in the MEblue module of this study. The MFS family is an ancient superfamily of transporters, which can be divided into 76 subfamilies and has the functions of catalyzing the transport of various substrate molecules and multiple compounds [54]. It has been reported to mediate the efflux of secondary metabolites in Trichoderma species [55] and facilitate the transport of nucleosides such as adenosine and purine nucleosides in Escherichia coli [56]. As a nucleoside-derived secondary metabolite in Cordyceps militaris, the extracellular transport mechanism of cordycepin remains unclear. Combined with the correlation between MFS transporters and cordycepin content observed in this study, we speculate that MFS family transporters may be involved in the transmembrane transport of cordycepin. Future studies will perform knockout and overexpression experiments on key MFS candidate genes to clarify their functional roles in the transmembrane transport of cordycepin, and further elucidate the regulatory mechanisms underlying the extracellular efflux of cordycepin mediated by MFS transporters.

In conclusion, this study systematically clarified that high-temperature stress promotes cordycepin accumulation in Cordyceps militaris through the synergistic effects of increased supply of cordycepin synthesis precursors, upregulated expression of key genes in the synthesis pathway, inhibition of the cordycepin degradation pathway, and enhanced activity of transporters via integrated multi-omics technologies. The results of this study provide a solid theoretical basis for optimizing Cordyceps militaris cultivation techniques through temperature regulation. Subsequent research focusing on the functional verification of key transporters will further promote breakthroughs in the industrial application of cordycepin and facilitate the efficient production of Cordyceps militaris with high cordycepin content