Genomic Sequence and Transcriptome Analysis of the Medicinal Fungus Keithomyces neogunnii.

The filamentous fungus Keithomyces neogunnii can infect the larvae of Lepidoptera (Hepialus sp.) and form an insect-fungi complex, which is utilized as an important traditional Chinese medicine. As a valuable medicinal fungus, K. neogunnii produces diverse bioactive substances (e.g., polysaccharide, vitamins, cordycepic acid, and adenosine) under cultivation conditions. Herein, we report the first high-quality genome of the K. neogunnii single-spore isolate Cg7.2a using single-molecule real-time sequencing technology in combination with Illumina sequencing. The assembled genome was 32.6 Mb in size, containing 8,641 predicted genes and having a GC content of 52.16%. RNA sequencing analysis revealed the maximum number of differentially expressed genes in the fungus during the stroma formation stage compared with those during the mycelium stage. These data are valuable to enhance our understanding of the biology, development, evolution, and physiological metabolism of K. neogunnii.

Introduction

Keithomyces neogunnii, one of the important entomopathogenic fungi, is associated with the larvae of Lepidoptera (Hepialus sp.). It is widely distributed in Southern China and has been considered a traditional Chinese medicine for more than a thousand years (Wen et al. 2017). In China, the fungus was initially denoted as Cordyceps gunnii (Liang 1983) and was later renamed as Metacordyceps neogunnii in 2017 (Wen et al. 2017). At present, it is considered a part of a new genus, namely, Keithomyces in the family Clavicipitaceae (Mongkolsamrit et al. 2020). This revised taxonomic position suggests that K. neogunnii is an essential biomaterial for studying the evolution and genetic diversity of entomopathogenic fungi. The asexual stage of the fungus (anamorph) was identified as Paecilomyces gunnii Liang sp. nov. (Liang 1985).

Keithomyces neogunnii exerts numerous pharmacological effects, including antitumor, memory-enhancing, analgesic, immunomodulatory, cardiovascular disease-preventing, and antioxidant activities, which are similar to those of the famous Chinese traditional medicine Ophiocordyceps sinensis (Zhu et al. 2016; Sun et al. 2018; Gan et al. 2019). These effects of K. neogunnii depend on its ability to synthesize various bioactive substances, including polysaccharides, nucleosides, amino acids (glutamate), vitamins (nicotinic acid), cordycepic acid, adenosine, and tyrosinase inhibitors (Zhu et al. 2012, 2013; Lu et al. 2014; Liu et al. 2019). Improving the extent of biosynthesis of these substances is of considerable important for the industrial utilization of this fungus. Till date, the lack of information regarding this fungal genome has limited further understanding of the biosynthesis of these active compounds.

Artificial cultivation of Cordyceps fungi is an important approach to solve the scarcity of wild Cordyceps resources. For some Cordyceps species, such as C. militaris, C. cicadae, and C. guangdongensis, the industrial cultivation of fruiting bodies on artificial medium has been achieved (Dong et al. 2015; Li et al. 2019). Artificial cultivation of K. neogunnii still remains a challenge, and only one successful case has been reported (Liu et al. 1990). Whole-genome sequencing provides a powerful approach to understand the molecular basis of the growth, development, and evolution of entomopathogenic fungi (Gao et al. 2011; Zheng et al. 2011; Kramer and Nodwell 2017; Lu et al. 2017; Zhang et al. 2018; Shu et al. 2020). Herein, we report for the first time the genomic sequence and the transcriptome analysis of K. neogunnii, which is immensely valuable for elucidating the genome-related biological issues of the fungus.

Results and Discussion

Genome Sequencing and General Properties

Keithomyces neogunnii Cg7 was obtained from a field fruiting body that formed a complex with a Lepidoptera larva (fig. 1). Under the artificial cultivation, the isolated fungal strain formed stroma (with a length of 1–3 cm and width of 0.2–0.4 cm) (fig. 1), which displayed a morphology that was similar to that of the field fruiting body. No ascospores were detected on the stroma in this study. The genome of the single-spore isolate Cg7.2a was sequenced using the PacBio Sequel platform, which generated 13 contigs with an N50 value of 5.08 M from a genome coverage of 131.9× (4.305 G clean data). The assembled genome was 32.6 Mb in size, containing 8,641 predicted genes and having a GC content of 52.16% (table 1). Meanwhile, the Illumina NovaSeq sequencing yielded a genome of 32.4 Mb (3.981 G clean data) with a genome coverage of 122.1× (supplementary table S1, Supplementary Material online). Analysis of the genome sequences with Benchmarking Universal Single-Copy Orthologs (BUSCO) revealed that 99.3% of fungi BUSCOs were complete (only 1.0% of these were duplicated) and 0.7% of those missed, suggesting that the assembly was relatively contiguous.

Fig. 1.

Fruiting body, MAT locus, and evolutionary analysis of Keithomyces neogunnii. (a) The field fruiting body of K. neogunnii Cg7. (b) The formed stroma of K. neogunnii Cg7 on artificial medium. (c) Comparative analysis of K. neogunnii MAT locus with C. militaris. (d) A neighbor-joining (NJ) phylogenomic tree constructed using TreeBest v1.9.2 depicting the evolutionary relationship between K. neogunnii and other fungal species. The sequence information of the fungi was collected from the NCBI genome database, and the accession numbers were GCA_012934285.1 (Ophiocordyceps sinensis), GCA_001272575.2 (O. unilateralis), GCA_001653235.2 (Pochonia chlamydosporia), GCA_000225605.1 (Cordyceps militaris), GCA_006981975.1 (C. javanica), GCA_000814975.1 (Metarhizium anisopliae), GCA_000187405.1 (M. acridum), GCA_000814955.1 (M. guizhouense), GCA_000591435.1 (M. robertsii), GCA_000814945.1 (M. majus), GCA_001653265.1 (Purpureocillium lilacinum), GCA_000804445.1 (M. album), GCA_001636745.1 (M. rileyi), GCA_000280675.1 (Beauveria bassiana), GCA_000328475.2 (Ustilago maydis), GCA_013426205.1 (M. brunneum), GCA_000687475.2 (Ustilaginoidea virens), GCA_000182925.2 (Neurospora crassa), and GCA_000011425.1 (Aspergillus nidulans).

Table 1

General Genome Properties of Keithomyces neogunnii and Other Entomopathogenic Fungi

Species	K. neogunnii	C. militaris	M. acridum	C. cicadae
Method	PacBio+Illumina	Roche 454 GS FLX system	Illumina	Illumina
Genome size (Mb)	32.61	32.2	38.05	33.9
GC content (%)	52.16	51.4	49.91	53.0
Protein-coding genes	8,641	9,684	9,849	9,701
Coverage (fold)	131.9×	147×	107×	80×
Percentage repeat rate	2.6092	3.04	1.52	3.19
Gene average length (bp)	1,613	1,742	—	—
Gene density (genes per Mb)	264.98	257	258.8	286
tRNA	124	136	122	—
Accession no.	JAIOWQ000000000	AEVU00000000	ANDI00000000	MWMN00000000
Reference	This study	Zheng et al. (2011)	Gao et al. (2011)	Lu et al. (2017)

Characteristic telomeric repeat sequences (TTAGGG/CCCTAA)n existed at either the 5ʹ or 3ʹ terminal of 11 contigs (supplementary table S2, Supplementary Material online). Among them, contig 1 contained a telomere at both ends, suggesting that it was a complete chromosome. Single-ended telomeric repeats were detected at the 3ʹ terminal of contigs 3, 5, 8, 9, and 10 or the 5ʹ terminal of contigs 2, 6, 12, and 13. The remaining two contigs lacked the telomeric repeat. A total of 1,445 (0.58% in the genome) interspersed nuclear elements and 9,092 tandem repeats (including minisatellite DNA and microsatellite DNA) were predicted in the fungal genome (supplementary tables S3 and S4, Supplementary Material online). The identified noncoding RNA in the genome included 124 transfer RNA (tRNA), 39 ribosomal RNA (rRNA), 2 small RNA (sRNA), and 17 small nuclear RNA (snRNA) (supplementary table S5, Supplementary Material online).

Mating-Type and Phylogenetic Analysis

Mating compatibility in fungi is controlled by genes located on the mating-type (MAT) loci (Ni et al. 2011). In heterothallic (self-sterile) fungal species (e.g. C. militaris), sexual propagation occurs between two strains of different MAT idiomorphs (MAT1-1 and MAT1-2) (Zheng et al. 2011). We identified only one MAT1-1 MAT locus in the genome of K. neogunnii Cg7.2a. The locus contained MAT1-1-1 (α-box), MAT1-1-2, and MAT1-1-3 (high mobility group [HMG-]-box) genes and a gene encoding DNA lyase downstream (fig. 1). Based on this result, we suggested that K. neogunnii belonged to a heterothallic filamentous ascomycete species. The taxonomic status of K. neogunnii has been revised twice based on multi-gene phylogenetic analyses and morphological characteristics (Wen et al. 2017, Mongkolsamrit et al. 2020). Phylogenomic analysis indicated that the fungal strain Cg7.2a was located on a single branch on the tree, and it was more closely related to the genera Metarhizium and Pochonia than to the genera Ophiocordyceps and Cordyceps (fig. 1).

Gene Functional Annotation

Among the predicted 8,641 genes, 96.09% (8,303), 33.11% (2,861), 33.11% (8,074), 22.34% (1,930), 65.72% (5,679), 5.22% (451), and 65.72% (5,679) genes exhibited similarity (identity ≥40%, coverage ≥40%) to the annotated proteins in Non-Redundant (NR) Protein Database, Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), Cluster of Orthologous Groups (KOG), Gene Ontology (GO), Carbohydrate-active enZYmes (CAZymes), and Pfam, respectively. Approximately, 17.96% (1,552) and 4.54% (392) of the genes encoded by the fungal genome exhibited good similarities (identity ≥40%, coverage ≥40%) with the genes in the Pathogen–Host Interactions (PHI) database and the database of fungal virulence factors (DFVF), respectively.

A total of 451 CAZyme-encoding genes were identified in the K. neogunnii genome, which contribute to the degradation of substrates (e.g., glycoside hydrolase [GH] and auxiliary activity [AA] families) and biosynthesis of oligosaccharides or polysaccharides (e.g., glycosyl tansferase [GT] family). Remarkably, the fungus contained a higher number of genes belonging to GH family (241) than other entomopathogenic fungi, such as C. militaris (134) (Zheng et al. 2011), C. cicadae (135) (Lu et al. 2017), M. aniopliae (156) (Gao et al. 2011), M. acridum (140) (Gao et al. 2011), and C. guandongensis (103) (Zhang et al. 2018) (supplementary table S6, Supplementary Material online). Several genes involved in the degradation of chitin, which is abundantly present in the insect cuticle, were detected in the fungal genome, including 4 copper-dependent lytic polysaccharide monooxygenase (AA11 family) genes and 23 chitinase (GH18 family) genes. Furthermore, the fungal genome contained 126 P450 proteins and 659 secreted proteins, which contribute to fungal physiologies and fungus–environment (or fungus–host) interactions, respectively (Staats 2014; Xu et al. 2020). A total of 36 secondary metabolism clusters (involving 386 genes), including 9 of type 1 polyketide synthases, 7 of nonribosomal peptide synthetases (NRPSs), 6 of NRPS-like, and 6 of terpene biosyntheses, were predicted to be present in the fungal genome (supplementary table S7, Supplementary Material online), which provided important information to further elucidate the biosynthesis pathways of active compounds produced by the fungus.

Gene Expression Profile of K. neogunnii during Different Growth Stages

Compared with the mycelium (Myc) growth stage, 4,230 (2,346 upregulated genes and 1,884 downregulated genes) and 4,155 (2,270 upregulated genes and 1,885 downregulated genes) genes in the K. neogunnii genome exhibited differential expression in stroma 1 (Str_1, 40 days-old culture) or stroma 2 (Str_2, 60 days-old culture) stages, respectively. Str_1 and Str_2 shared the maximum (3,423) differentially expressed genes (supplementary figs. S1 and S2, Supplementary Material online). The upregulated genes in the stroma stages are involved in sexual differentiation, secondary metabolism, and nutrition metabolism. Compared with the Myc stage, the MAT protein MAT-1-1 (A5637, α-1 HMG-box) showed a 5.7- or 5.1-fold increment in Str_1 or Str_2 stage, respectively. Meanwhile, a global secondary metabolism regulator LAE1 (methyltransferase)-encoding gene (A4589) was upregulated by 2.9-fold in Str_1 or Str_2, and a P450 protein-encoding gene (A4528) was upregulated by 3.2- or 3.6-fold during the developmental stage of Str_1 or Str_2, respectively. Compared with the Myc stage, the putative fungal hydrophobin (A8510), fatty acid desaturase (A3472), and GH18 chitinase (A2691) were downregulated in the stroma stages. The top 300 differentially expressed genes in the group of Str_1 versus Myc or Str_2 versus Myc were listed in the supplementary tables S8–S11, Supplementary Material online.

Conclusions

In this study, we presented the first high-quality genome assembly of the medicinal fungus K. neogunnii and performed the comparative transcriptome analysis of the fungus at different developmental stages. The fungus possesses abundant genes involved in nutrition metabolism, biosynthesis of active substances and fungus–host interaction. The genomic sequences and the data obtained through transcriptome analysis are immensely valuable for elucidating the biology, development, evolution, physiological metabolism, and pathogenicity of the fungus in molecular studies.

Materials and Methods

Fungal Isolate and Culture Conditions

Keithomyces neogunnii strain Cg7 was originally isolated from a wild fruiting body collected from Shibing County, Guizhou Province, China. After taxonomic identification of the strain by internal transcribed spacers sequence analysis, it was deposited in the China Centre for Type Culture Collection (CCTCC), Wuhan, China, under the number CCTCC M2021416. To obtain a pure isolate, the strain Cg7 was grown on a potato dextrose agar (PDA) plate at 25°C for sporulation, and single conidial isolates were obtained by serial dilution and spread plate method. After two rounds of single conidial isolation, the pure isolate Cg7.2a was picked and maintained on a PDA slant at 4 °C. Modified Sabouraud liquid (MSL) medium (pH 6.0) composed of 0.5% maltose, 2% glucose, 0.5% yeast extract, 1% peptone, 0.1% K2HPO4, 0.1% KCl, and 0.1% MgSO4·7H2O. The cultivation medium (per bottle) contained 20 g millet, 0.5 g silkworm pupa powder, 2.0 g sesame oil, and 40 ml nutrient solution (per liter, 30 g glucose, 10 g yeast extract, 10 g peptone, 10 g maltose, 1 g K2HPO4, 1 g MgSO4·7H2O, 0.5 g KCl, 1 g sodium citrate, 2 g sodium laurate, and 10 mg vitamin B1, pH 6.0).

Artificial cultivation of K. neogunnii Cg7 was performed as follows: A fungal block from the PDA slant was incubated in a 250-ml flask containing 50 ml MSL medium and cultured at 25 °C and 120 rpm for 3 days until the medium was filled with the fungal Myc. Next, 3 ml of seed culture was transferred to a bottle of cultivation medium, cultured in the dark at 24 °C for 3 days until the bottle was filled with the fungal Myc. To stimulate primordium growth, the cultivation room was illuminated by fluorescent light at an intensity of 800 lux for 18 h/day, for 60 days.

Genome Sequencing and Assembly

Keithomyces neogunnii isolate Cg7.2a was grown in 100 ml MSL medium at 25 °C and 120 rpm for 2 days. After collecting the Myc by filtration, the fungal genomic DNA was isolated using a standard cetyltrimethylammonium bromide method (Watanabe et al. 2010), analyzed on an agarose gel and quantified using a Qubit 2.0 Fluorometer (Thermo Scientific, USA). The fungal genome was sequenced using the PacBio Sequel platform and the Illumina NovaSeq PE150 platform at the Beijing Novogene Bioinformatics Technology Co., Ltd (Beijing, China). For single-molecule real-time (SMRT) sequencing, a DNA library with an insert size of 20 kb was generated using the SMRT bell template kit (version 1). Briefly, the DNA sample was fragmented to the target size using a g-TUBE device (Covaris, Woburn, MA, USA), repair DNA damage as well as ends, ligation of DNA fragments with hairpin connectors, purification of DNA with AMpure PB beads, size selection using the BluePippin System, repair of DNA damage, assessment of library quality on the Qubit 2.0 Fluorometer (Thermo Scientific, USA), and detection of the insert DNA size using Agilent 2100 (Agilent Technologies, USA). For the Illumina NovaSeq sequencing, the genomic DNA was fragmented to approximately 350 bp in size by sonication, and the sequence library was constructed using NEBNext Ultra DNA Library Prep Kit for Illumina (NEB, USA). Low-quality raw reads obtained from PacBio were filtered (<500 bp) to obtain clean data. Using the SMRT portal, long reads of >6,000 bp were selected as the seed sequence, the sequence accuracy was improved by alignment with other short reads, and the reads were initially assembled into contigs using SMRT Link v5.0.1 (Ardui et al. 2018). The preliminary assembly results were corrected, and the variant sites were counted using the variant Caller module of the SMRT Link software with the arrow algorithm, thus yielding the reference genome. The clean data obtained from the Illumina sequencing platform were used for the genome survey analysis. Genome completeness was evaluated with the BUSCO 3.0.2 software using the lineage data set fungi_odb9 (Seppey et al. 2019).

Analysis of Genome Component and Gene Function

Prediction of the interspersed repetitive sequences was performed using RepeatMasker (http://www.repeatmasker.org/). Tandem repeats existing in the genome were searched using the Tandem Repeats Finder (version 4.09) program (Benson 1999). tRNA genes and rRNA genes were analyzed using tRNAscan-SE (Chan and Lowe, 2019) and RNAmmer 1.2 (Lagesen et al. 2007), respectively. sRNA, snRNA, and microRNA were predicted by BLAST against the Rfam database (Kalvari et al. 2021).

A whole-genome BLAST search (e-value of <1e-5, minimal alignment length percentage of >40%) was performed against seven databases, including GO, KEGG, KOG, NR, TCDB, P450, and Swiss-Prot. The secondary metabolism gene clusters were predicted using antiSMASH version 4.02 (Medema et al., 2011). The secretory proteins were analyzed using the Signal P database. The PHI and DFVF databases were used for the analysis of related functional genes in the genome. CAZymes were predicted using the CAZymes database (Lombard et al. 2014).

Orthology and Phylogenetic Analysis

The protein sequences from the selected fungal genomes were analyzed using the BLAST software (version 2.2.26), and the gene family clustering was obtained by the comparison of similarity with Hcluster-sg (version 0.2.0) after removing redundancy using the Solar software (version 0.9.6). Single-copy orthologous genes were extracted based on the results of gene family clustering. A total of 805 orthologous proteins were obtained and subjected to multiple sequence alignment using the Multiple Sequence Comparison by Log-Expectation software. A neighbor-joining phylogenomic tree was constructed using concatenated amino acid sequences using TreeBest (Version 1.9.2).

RNA Sequencing Analysis

The undifferentiated fungal mycelia grown in MSL medium (at 25 °C and 120 rpm for 48 h) and the stroma formed after 40 or 60 days of fungal cultivation were used for RNA extraction and sequencing. The fungal samples were homogenized in liquid nitrogen, and total RNAs were extracted using the TRIzol reagent (Thermo Fisher Scientific, CA, USA), followed by DNase I (New England Biolabs, Ipswich, MA, USA) digestion. Messenger RNA was purified from the total RNAs using the poly-T oligo-attached magnetic beads method (Green and Sambrook, 2019). First-strand cDNA was synthesized using a random hexamer primer and M-MulV Reverse Transcriptase. After the synthesis of second-strand cDNA by polymerase chain reaction amplification, cDNA libraries were constructed and sequenced using the Illumina NovaSeq PE150 platform at the Beijing Novogene Bioinformatics Technology Co. Ltd (Beijing, China). Three independent biological replicates were established. Clean data were obtained by removing low-quality reads from raw data and were mapped to the fungal genome. False discovery rates (Padj ≤ 0.05) and log2-fold change of ≥1 thresholds were used to estimate the level of differential gene expression by the species under different growth stages.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.