Telomere-to-Telomere Assembly of the Cordyceps militaris CH1 Genome and Integrated Transcriptomic and Metabolomic Analyses Provide New Insights into Cordycepin Biosynthesis Under Light Stress

Abstract

Cordyceps militaris, a model species in the genus Cordyceps, is widely distributed globally and is known for its significant medicinal value. It has been traditionally used in Chinese medicine to enhance immunity, alleviate fatigue, and treat tumors, among other therapeutic purposes. Here, we successfully assembled a telomere-to-telomere (T2T) level genome of C. militaris CH1 using PacBio HiFi and Hi-C technologies. The assembled genome is 32.67 Mb in size, with an N50 of 4.70 Mb. Gene prediction revealed a total of 10,749 predicted genes in the C. militaris CH1 genome, with a gene completeness of 99.20%. Phylogenetic analysis showed the evolutionary relationship between C. militaris CH1 and other Cordyceps species, suggesting that the divergence between this strain and C. militaris ATCC 34164 occurred approximately 1.36 Mya. Combined transcriptomic and metabolomic analyses identified 842 differentially expressed genes and 2052 metabolites that were significantly altered under light stress, primarily involving key pathways related to amino acid metabolism, purine metabolism, and secondary metabolite biosynthesis. Joint analysis of genes and metabolites revealed 79 genes coding for enzymes associated with the synthesis of adenine and adenosine, with the expression of 52 genes being upregulated, consistent with the accumulation trends of adenine and adenosine. Four gene clusters related to the synthesis of cordycepin were identified, with a significant upregulation of cns3 (FUN_003263), suggesting that light stress may promote cordycepin biosynthesis. This comprehensive analysis not only provides new insights into the genomics, metabolomics, and functional gene research of C. militaris CH1 but also offers a potential biological foundation for understanding the synthesis mechanisms of cordycepin and its efficient production.

Keywords: Cordyceps militaris CH1; comparative genomics; cordycepin biosynthesis; phylogenetic evolution; telomere-to-telomere; transcriptome and metabolome.

Figure 1

The telomere-to-telomere (T2T) genome assembly of C. militaris CH1. (A) The genomic survey map of C. militaris CH1 shows the distribution of different k-mers in the density plot. The black line represents the overall model fitting curve, while the red line indicates the erroneous K-mer distribution curve. (B) The Hi-C interaction heatmap for C. militaris CH1 indicates that the stronger the interaction, the redder the color, and conversely, the weaker the interaction, the more subdued the color. (C) The circos plot representing the genomic features of C. militaris CH1, with the outermost layer representing the chromosomes and the colored areas indicating coding regions. A 1 kb sliding window was used for segmentation, and the heatmap and bar chart from outer to inner layers represent GC content, gene density, density of tandem repeats, density of transposable elements, density of LTR_Gypsy, and density of LTR_Copia, respectively. The innermost lines illustrate the synteny of C. militaris CH1 itself. (D) The T2T assembly results for C. militaris CH1, including the positions of the telomeres and non-coding RNAs.

Figure 2

A bar chart for gene annotation of C. militaris CH1. (A) GO annotation statistics chart displays 12 significantly enriched terms categorized by biological process (BP), cellular component (CC), and molecular function (MF) using color coding. (B) KEGG pathway annotation chart lists the names of metabolic pathways on the vertical axis, while the horizontal axis quantifies the number of annotated genes in each pathway. (C) COG functional classification shows the functional categories of genes on the horizontal axis, with the vertical axis representing the number of genes in each functional category.

Figure 3

C. militaris CH1 gene family clustering and phylogenetic analysis with closely related species. (A) The Venn diagram depicts the gene families of C. militaris CH1 in comparison with C. pruinosa, C. fumosorosea, C. tenuipes, C. militaris ATCC, and C. javanica, highlighting the number of unique and shared gene families. (B) KEGG pathway enrichment of shared genes across six species. (C) The phylogenetic tree of C. militaris CH1 and 11 other fungi, with C. cyperi as the outgroup; the red numbers indicate divergence times, and the bar chart shows the number of different types of gene families in each fungus, including single-copy orthologs, multi-copy orthologs, unique paralogs, and other orthologs. The blue portion of the pie chart represents the number of expanded gene families, while the green portion indicates the number of contracted gene families. (D) The bar chart shows the GO functional enrichment of the expanded gene families of the C. militaris CH1 genome. (E) The dot plot illustrates the KEGG enrichment of the expanded gene families of the C. militaris CH1 genome.

Figure 4

Comparative genomics analysis of intra-species and inter-species. (A) SNP density plot with SNP counts calculated in 0.1 Mb sliding windows. (B) Comparative genomic analysis based on chromosomes. (C) Clustering analysis of all gene families. (D) Genomic similarity comparison. (E) GO enrichment analysis of genes unique. (F) KEGG enrichment analysis of genes unique. All analyses were completed between C. militaris CH1 (Cmi_CH1) and C. militaris ATCC (Cmi_ATCC), except that GCA_003332165.1 and GCF_000225605.1 were included in (D).

Figure 5

Transcriptome and metabolome analysis of mycelia growth of C. militaris CH1 under light and dark conditions. (A) The number of differentially expressed genes (upregulated and downregulated) under dark and light stress. (B) The top 20 GO enrichment items of differential genes. (C) Classification and proportion of all metabolites in C. militaris CH1 after light stress. (D) PCA analysis of control and treatment groups for metabolome analysis. (E) Volcanic maps of differential metabolites. (F) Ranked among the top 20 most significant variations in differential metabolites. (G) KEGG enrichment analysis of differential metabolites.

Figure 6

A combination of transcriptomic and metabolomic analysis of C. militaris CH1 under light stress. (A) KEGG pathways enriched by both genes and metabolites. (B) Heatmap showing the co-changes in genes and metabolites. (C) Interaction network of genes and metabolites in the purine synthesis pathway (ko00230). The red dots are target genes, the green dots are metabolites, and the connecting lines indicate the target genes required for the metabolites. (D) Interaction network of genes and metabolites involved in the synthesis of adenine and adenosine in the metabolic pathway (ko01100). (E) Heatmap of the expression levels of genes involved in the synthesis of adenine and adenosine. (F) Cordycepin biosynthesis pathway. Abbreviations: AMP: adenosine monophosphate; 3′-AMP: 3′-adenosine monophosphate; 2′,3′-cAMP: 2′,3′-cyclic adenosine monophosphate; ADK: adenosine kinase; ADEK, adenylate kinase; 3′-dAMP: 3′-deoxyadenosine monophosphate; ADP: adenosine diphosphate; 3′-dADP: 3′-deoxyadenosine diphosphates; 2′-C-3′-dA: 2 ‘-carbonyl -3′ -deoxyadenosine; RNR: ribonucleotide reductases; APRT, adenine phosphoribosyltransferase; NT5E, 5′-nucleotidase.

Figure 7

Verification of the cordycepin synthesis inferred from the transcriptome was conducted using qRT-PCR. Four genes related to cordycepin synthesis (FUN_003262, FUN_003263, FUN_003264, FUN_003266) and four genes involved in adenine synthesis (FUN_000548, FUN_000593, FUN_005768, FUN_009691) were selected for qRT-PCR analysis to validate the reliability of the RNA-seq results. The bar charts illustrate the FPKM values of these genes under light and dark conditions, as well as the relative expression levels calculated from qRT-PCR, with statistical significance between treatments denoted by asterisks: * indicates p < 0.05, ** indicates p < 0.01, and *** indicates p < 0.001.