Comparative analysis of the complete mitochondrial genomes of four cordyceps fungi.

Cordyceps is a large group of entomogenous, medicinally important fungi. In this study, we sequenced, assembled, and annotated the entire mitochondrial genome of Ophiocordyceps xuefengensis, in addition to comparing it against other three complete cordyceps mitogenomes that were previously published. Comparative analysis indicated that the four complete mitogenomes are all composed of circular DNA molecules, although their sizes significantly differ due to high variability in intron and intergenic region sizes in the Ophiocordyceps sinensis and O. xuefengensis mitogenomes. All mitogenomes contain 14 conserved genes and two ribosomal RNA genes, but varying numbers of tRNA introns. The Ka/Ks ratios for all 14 PCGs and rps3 were all less than 1, indicating that these genes have been subject to purifying selection. Phylogenetic analysis was conducted using concatenated amino acid and nucleotide sequences of the 14 PCGs and rps3 using two different methods (Maximum Likelihood and Bayesian analysis), revealing highly supported relationships between O. xuefengensis and other Ophiocordyceps species, in addition to a close relationship with O. sinensis. Further, the analyses indicated that cox1 and rps3 play important roles in population differentiation. These mitogenomes will allow further study of the population genetics, taxonomy, and evolutionary biology of medicinally important cordyceps species.

INTRODUCTION

Mitochondrial genomes are powerful tools that can be used in phylogenetic and evolutionary studies (Cameron et al., 2007; Liu et al., 2018; Saccone et al., 2000; Song et al., 2019), owing to their small size, high copy numbers, conserved orthologous genes, rare instances of recombination, and high evolutionary rates (Curole & Kocher, 1999). The rapid development of new genome sequencing technologies in recent years has led to the production of more than 180 fungal mitogenomes. Fungal mitochondria range in length from approximately 11 kbp (Hanseniaspora uvarum) to 272.2 kbp (Morchella importuna), depending on the species, but usually contain 14 conserved genes (Pramateftaki et al., 2010). In addition, mitogenome characteristics, including genome size, genome structure, gene content, gene arrangement, tRNA number, intron number, codon usage, and repeat content, can provide useful information to understand the origin, evolution, and systematics of eukaryotes (Poliseno et al., 2017; Qiang et al., 2018; Sankoff et al., 1992).

‘Cordyceps’ comprises four genera including Cordyceps, Ophiocordyceps, Metacordyceps, and Elaphocordyceps. The fungi have been used in traditional Chinese medicine, as dietary supplements, or as tonic edibles with broad pharmacological properties including antibacterial, antiviral, and antitumor activities, in addition to the ability to enhance human immunity (Yue et al., 2013; Zhong et al., 2019). Although Cordyceps have considerable economic and medicinal significance, wild cordyceps resources are on the verge of exhaustion, leading to many investigators seeking alternatives. Among these, O. xuefengensis is a newly identified cordyceps resource that we previously discovered on Xuefeng mountain in Hunan province of China (Wen et al., 2013). The novel taxa has been widely used as an ethnopharmacological invigorant by the Yao group for a long time (Jin et al., 2019). Further, we have successfully cultivated this cordyceps fungi and also characterized some of its chemical metabolites (Zhong et al., 2019). Cordyceps and Ophiocordyceps are the two most widely used groups of cordyceps fungi. In particular, O. sinensis is popularly used to treat many diseases in China. Consequently, the cultivation, novel metabolite compounds, and pharmacology of cordyceps fungi have been extensively researched, although further molecular phylogenetic and species identification studies are critically needed.

In this study, the complete mitogenome from O. xuefengensis was sequenced, assembled, and compared to other three complete cordyceps mitogenomes that were previously reported. This study thereby represents the first comparative analysis of the genomic structure, base composition, substitution, and evolutionary rates among four cordyceps species, in addition to the comprehensive molecular phylogenetic analysis of cordyceps in the Hypocreales order. This study aimed to expand our knowledge of the diversity of mitochondrial genomes and phylogenetic relationships of cordyceps.

MATERIALS AND METHODS

Fungal materials and DNA extraction

Ophiocordyceps xuefengensis strain HCMA001 was used in this study for mitogenomic analysis. In vitro mycelial fermentation in liquid culture was conducted in Erlenmeyer flasks with shaking at 150 rpm in an incubator at 25°C for 10 days. Mycelia were collected by centrifugation and quickly frozen in liquid nitrogen. Total genomic DNA was then extracted using an improved extraction method (Chen et al., 2011). Genomic DNA was quantified using an TBS‐380 fluorometer (Turner BioSystems Inc., Sunnyvale, CA). High quality DNA samples (OD260/280 = 1.8–2.0, >6 μg) were used to construct fragment libraries for genomic sequencing.

Mitochondrial DNA sequencing and genome assembly

About 8 μg of purified DNA was sequenced using a combination of third‐generation sequencing technologies including PacBio RS and Illumina sequencing platforms. The Illumina data were used to evaluate the complexity of the genome and correct the PacBio long reads. The raw paired‐end reads were trimmed, and quality controlled using Trimmomatic with parameters including SLIDINGWINDOW: 4:15 MINLEN: 75 (version 0.36). Clean data obtained from the above quality control procedures were then used for further analyses.

Annotation and comparative analysis of the Ophiocordyceps xuefengensis mitogenome

Ab initio prediction was used to generate gene models for the O. xuefengensis mitogenome. Gene models were identified using MFannot (https://megasun.bch.umontreal.ca/cgi‐bin/dev_mfa/mfannotInterface.pl). Gene models were then compared against the non‐redundant (NR) NCBI database using BLASTp, in addition to comparison against the SwissProt, KEGG, and COG databases to establish functional annotation, also using BLASTp. tRNAs were identified using the tRNAscan‐SE program (v1.23, http://lowelab.ucsc.edu/tRNAscan‐SE) and rRNAs were identified using the RNAmmer program (v1.2, http://www.cbs.dtu.dk/services/RNAmmer/). Circular maps of the four mitochondrial genomes were drawn using OGDRAW (https://chlorobox.mpimp‐golm.mpg.de/OGDraw.htmlonline) (Conant & Wolfe, 2008).

Sequence analysis

Nucleotide compositional skew was calculated according to the formulae: AT‐skew = (A − T)/(A + T) and GC‐skew = (G − C)/(G + C) (Perna & Kocher, 1995). The Sequence Manipulation Suite program package (http://www.bioinformatics.org/sms2/codon_usage.html) was used to analyze codon usage based on genetic code ‘4’ (Stothard, 2000). Genetic distances between each pair of the 14 core PCGs (atp6, atp8, atp9, cob, cox1, cox2, cox3, nad1, nad2, nad3, nad4, nad4L, nad5, and nad6) were calculated with MEGA X (Joseph, 2016), using the Kimura−2‐parameter (K2P) substitution model. DnaSP v6 (Rozas et al., 2017) was used to calculate the nonsynonymous (Ka) and synonymous (Ks) substitution rates for the 14 core PCGs among the four mitogenomes. Genomic synteny among the four mitogenomes was analyzed using the Mauve v2.4.0 program (Darling et al., 2004).

Repetitive element analysis

Tandem repeats (>10 bp in length) in the four mitogenomes were detected using the Tandem Repeats Finder program (Benson, 1999), with default parameters. In addition, repeat sequences were identified using the REPuter program (Kurtz et al., 2001) to identify forward (direct), reverse, complemented, and palindromic (reverse complemented) repeats across the four mitogenomes, using a minimum repeat size of 30 and a hamming distance set to 3.

Phylogenetic analysis of Hypocrealean species

To determine the phylogenetic relationships of the cordyceps species among the Hypocreales order, 19 mitogenomes were downloaded from the NCBI database and subjected to Maximum Likelihood (ML) and Bayesian phylogenetic analysis using the aforementioned 14 PCGs and rps3. Rhizopogon salebrosus and Rhizopogon vinicolor were used as outgroups for the analyses (Table 1). Single mitochondrial genes were first aligned using MAFFT v7.037 (Katoh et al., 2019) and these alignments were then concatenated using SequenceMatrix v1.7.8 (http://www.softpedia.com/get/Science‐CAD/Sequence‐Matrix.shtml) (Vaidya et al., 2011). The ML tree was then calculated for the combined gene set using RAxML v8.0.0 (Stamatakis, 2014), while Bayesian inference (BI) analysis was performed with MrBayes v3.2.6 (Ronquist et al., 2012). Other analyses and evaluation methods were conducted based on previously described protocols (Qiang et al., 2018).

TABLE 1

The GenBank accessions for genomes used in the phylogenetic analysis

Species	Order	Family	Genus	Genome length/bp	GenBank accession number
Rhizopogon salebrosus	Agaricomycetes	Agaricomycetidae	Boletales	66,704	MH 794152
Rhizopogon vinicolor	Agaricomycetes	Agaricomycetidae	Boletales	77,109	MH 794153
Pochonia chlamydosporia	Hypocreales	Clavicipitaceae	Pochonia	25,615	NC_022835.1
Metarhizium anisopliae	Hypocreales	Clavicipitaceae	Metarhizium	24,673	NC_008068.1
Cordyceps militaris	Hypocreales	Cordycipitaceae	Cordyceps	33,277	NC_022834.1
Cordyceps brongniartii	Hypocreales	Cordycipitaceae	Cordyceps	33,926	NC_011194.1
Beauveria bassiana	Hypocreales	Cordycipitaceae	Beauveria	29,961	NC_010652.2
Hypomyces aurantius	Hypocreales	Hypocreaceae	Hypomyces	71,638	NC_030206.1
Fusarium solani	Hypocreales	Nectriaceae	Fusarium	62,978	NC_016680.1
Fusarium oxysporum	Hypocreales	Nectriaceae	Fusarium	34,477	NC_017930.1
Fusarium graminearum	Hypocreales	Nectriaceae	Fusarium	95,676	DQ364632.1
Hirsutella thompsonii	Hypocreales	Ophiocordycipitaceae	Hirsutella	62,509	NC_040165.1
Hirsutella minnesotensis	Hypocreales	Ophiocordycipitaceae	Hirsutella	52,245	NC_027660.1
Hirsutella vermicola	Hypocreales	Ophiocordycipitaceae	Hirsutella	53,793	NC_036610.1
Hirsutella rhossiliensis	Hypocreales	Ophiocordycipitaceae	Hirsutella	62,949	MG979071.1
Ophiocordyceps sinensis	Hypocreales	Ophiocordycipitaceae	Ophiocordyceps	157,539	NC_034659.1
Tolypocladium ophioglossoides	Hypocreales	Ophiocordycipitaceae	Tolypocladium	35,159	NC_031384.1
Tolypocladium inflatum	Hypocreales	Ophiocordycipitaceae	Tolypocladium	25,328	KY924880.1
Tolypocladium cylindrosporum	Hypocreales	Ophiocordycipitaceae	Tolypocladium	34,698	MN842262.1
Ophiocordyceps xuefengensis	Hypocreales	Ophiocordycipitaceae	Ophiocordyceps	78,744	This study

Data availability statement

The newly sequenced mitogenome of O. xuefengensis strain HCMA001 was submitted to GenBank under the accession number SAMN16236787.

RESULTS

Characterization of the four cordyceps mitogenomes

Comparative analysis indicated that the complete mitogenomes of the four cordyceps all comprised circular DNA molecules. The mitogenome length of O. xuefengensis (78,744 bp) was much smaller than that of O. sinensis (157,559 bp), but longer than those of Cordyceps militaris (33,277 bp) and Cordyceps brongniartii (33,926 bp) (Figure 1). All four mitogenomes encoded 14 core PCGs involved in energy metabolism, rps3 involved in translation, and 2 rRNAs. The PCG lengths varied, but the protein coding regions themselves were conserved. The tRNA numbers of the four mitogenomes ranged from 25 (C. brongniartii) to 27 (O. sinensis). The tRNAs comprised those for all 20 amino acids and ranged in size from 71 to 86 nucleotides (Table S1). The nucleotide composition of the O. xuefengensis mitogenome was identical to that of O. sinensis, but not to Cordyceps militaris and C. brongniartii. The O. sinensis mitogenome had the highest GC content (30.20%), followed by those of O. xuefengensis (29.94%), C. brongniartii (27.34%), and C. militaris (26.79%). The high frequency of A and T usage in codons contributes to the high AT content of the four mitochondrial genomes, with AT content reaching 70.06% in O. xuefengensis, 69.80% in O. sinensis, 73.21% in C. militaris, and 72.66% in C. brongniartii. Both GC skew and AT skew were positive for all four mitogenomes (Table 2).

FIGURE 1

Circular maps of the mitochondrial genomes of four cordyceps species. Genes are represented by different colored blocks. Colored blocks outside each ring indicate that the genes are on the direct strand, while colored blocks within the ring indicate that the genes are located on the reverse strand

TABLE 2

Characteristics of cordyceps mitogenomes

Characteristic	Ophiocordyceps xuefengensis	Ophiocordyceps sinensis	Cordyceps militaris	Cordyceps brongniartii
Size (bp)	78,744	157,559	33,277	33,926
A (%)	36.08	36.70	36.98	36.51
T (%)	33.98	33.10	36.23	36.15
G (%)	16.33	15.70	15.20	15.13
C (%)	13.61	13.70	11.59	12.21
Rate of GC (%)	29.94	30.20	26.79	27.34
Rate of AT (%)	70.06	69.80	73.21	72.66
AT‐skew	0.03	0.05	0.01	0.005
GC‐skew	0.09	0.09	0.14	0.11
PCGs (n)	14	14	14	14
Introns (n)	18	54	8	6
Intronic ORFs (n)	16	45	8	5
GIY‐YIG (n)	4	10	4	4
LAGLIDADG (n)	9	27	3	1
rRNAs (n)	1	1	1	1
nRNAl (n)	1	1	1	1
tRNAs (n)	26	27	26	25
PCG regions (bp)	34,994	104,257	17,735	18,329
Intronic regions (bp) a	21,146	107,288	9,812	8,452
Exon regions (bp) b	13,848	13,287	12,978	13,098
tRNA genes (bp)	1,934	2,010	1,924	1,867
rRNA genes (bp) c	7,225	7,144	4,666	4,801
Intergenic regions (bp)	34,591	27,831	3,797	5,708

Introns in PCGs and rRNA.

Exons in PCGs.

rnl and rns including exons.

The O. sinensis mitogenome exhibited an unusually enlarged size (157,559 bp). The compositions of regions were thus analyzed for the mitogenomes, indicating that the protein coding regions among the four mitogenomes are basically identical. A significant difference was observed when comparing intronic and intergenic regions, wherein the total number of intronic nucleotides of the O. sinensis mitogenome was much higher than for those of the other cordyceps, and the total number of intergenic nucleotides of the O. xuefengensis mitogenome is the highest among the four mitogenomes (Figure 2).

FIGURE 2

The abundances of protein‐coding, intronic, intergenic, and RNA gene regions (rRNAs and tRNAs) of the entire mitochondrial genomes from four cordyceps species

The number of introns also significantly differed among the four mitogenomes. 54 introns were present in the mitogenome of O. sinensis, 18 in that of O. xuefengensis, 8 in that of C. militaris, and 6 in that of C. brongniartii (Table 2). Further, the host genes of introns also differ among the four cordyceps. Introns are distributed in the cob (no. of introns: 1), cox1 (1), cox2 (1), cox3 (1), and rnl (4) genes of the C. militaris mitogenome, and in the cob (1), cox1 (2), nad1 (1), and rnl (2) genes of the C. brongniartii mitogenome. In O. xuefengensis, the introns are distributed in cob (2), cox1 (3), cox2 (1), cox3 (3), nad1 (2), nad2 (2), nad5 (1), rnl (3), and rns (1) genes, while in the O. sinensis mitogenome, the introns are distributed in the cob (6), cox1 (14), cox2 (6), cox3 (2), nad1 (3), nad2 (2), nad4L (1), nad5 (5), nad6 (1), atp6 (2), atp9 (1), rnl (8), and rns (1) genes (Table 3). The most common intron type among the mitogenomes are group IB introns. Introns also harbor 0–2 homing endonuclease genes, including the LAGLIDADG homing endonuclease and the GIY‐YIG endonuclease.

TABLE 3

Information of introns in four cordyceps mitogenome

	O. xuefengensis	O.s sinensis	C. militaris	C. brongniartii
Intron type
IA	1	4	1	1
IB	6	20	3	2
IC1		11	2	2
IC2	4	4
ID	2	3	1	1
I(derived, A)	2	1	1
I(derived, B1)		1
I(derived, B2)		6
II	1	4
Unknown	2
Conserved domain
LAGLI‐DADG	9	25	3	1
GIY‐YIG	4	7	4	4
Hypothetical protein	2	18
RNA‐dependent DNA polymerase		1

Protein‑coding genes and codon usage among the four mitogenomes

Most PCGs and rps3 are initiated with the ATG codon, whereas cox1 and nad6 in the O. xuefengensis mitogenome are initiated with the TTG codon. Further, cox3 in the O. sinensis mitogenome and atp9 in both O. xuefengensis and O. sinensis mitogenomes are initiated with the GTG codon. In addition, most PCGs and rps3 are terminated with the TAA codon, but cox1, rps3, and cox3 in the O. sinensis mitogenome, and cox1 and cob of the O. xuefengensis mitogenome are terminated with the TAG codon (Table S2). Comparative analysis indicated that all the PCGs and rps3 encode 20 amino acids, with similar amino acid compositions among the cordyceps mitogenomes. The most overused amino acid is Leu (ranging in prevalence from 628 to 652), whereas the least used amino acid is Cys (32–48 residues). The usage of Asn, Lys, Asp, and Glu in the O. xuefengensis mitogenome is much higher than in those of the other three cordyceps (Figure 3). The relative synonymous codon usage (RSCU) of the 14 PCGs and rps3 in the four cordyceps mitogenomes was also analyzed. A total of 4,616, 4,429, 4,326 and 4,366 codons were observed among the PCGs and rps3 of the O. xuefengensis, O. sinensis, C. militaris, and C. brongniartii mitogenomes, respectively. The most used codon among all four mitogenomes is UUA, but the most frequently used codon is AGA in the O. xuefengensis and O. sinensis mitogenomes (Figure 4).

FIGURE 3

Codon usage in the mitochondrial genomes of the four cordyceps species. Codon numbers are plotted on the y‐axis for a, O. xuefengensis; b, O. sinensis; c, C. militaris; d, C. brongniartii

FIGURE 4

RSCU values for 14 PCGs and rps3 in the mitochondrial genomes of four cordyceps

Variation, genetic distance, and evolutionary rates of common genes

The lengths and/or GC contents of the 14 core PCGs and rps3 are not all consistent across the four mitogenomes, with the exception of atp8, atp9, and nad4L PCGs. In particular, the length of nad5 in O. xuefengensis is about 1,000 bp longer than in the other three cordyceps mitogenomes. Further, the GC content of all 14 PCGs and rps3 differ, indicating that the core PCGs and rps3 have been modified within different cordyceps. Across all PCGs in the mitogenomes, GC content is highest in atp9 within all four cordyceps mitogenomes, and lowest in atp8 within the O. xuefengnensis and O. sinensis mitogenomes, but lowest in rps3 within the C. militaris and C. brongniartii mitogenomes. AT skew is negative in most PCGs except rps3, while GC skew is present in most PCGs except atp8. In particular, the AT skew of atp6 in the C. militaris mitogenome and nad6 in the O. xuefengnensis mitogenome are negative, while the GC skew of only rps3 in the C. brongniartii mitogenome is negative (Figure 5).

FIGURE 5

Variation in the length and base composition of each of the 14 PCGs and rps3 among four cordyceps mitogenomes. a, PCG length variation; b, GC content across PCGs; c, AT skew; d, GC skew

Across all 14 PCGs and rps3 genes that were evaluated, rps3 has the greatest K2P genetic distance among the mitogenomes, followed by cox1. atp8 has the least genetic distance among the genes investigated, indicating a high degree of conservation. Ka analysis indicated that rps3 and cox1 have relatively high nonsynonymous substitution rates. Further, the Ks of nad1 is the highest across all four mitogenomes, while that of atp8 is the lowest. The Ka/Ks values for all 14 PCGs and rps3 are <1, suggesting that these genes are subject to purifying selection (Figure 6).

FIGURE 6

Genetic analysis of 14 PCGs and rps3 across four cordyceps mitogenomes. K2P, the overall mean Kimura‐2‐Parameter distance; Ka, the mean number of nonsynonymous substitutions per nonsynonymous site; Ks, the mean number of synonymous substitutions per synonymous site

Repetitive elements

A total of 18 repeat sequences were identified in the mitogenome of O. xuefengensis, 250 in that of O. sinensis, 6 in that of C. militaris and 23 in that of C. brongniartii (Table S3). The length of repeat sequences ranges from 38 bp to 372 bp, with pairwise nucleotide similarities ranging from 74.73% to 100%. The largest repeat region was observed in the O. sinensis mitogenome, within ORF322 and ORF308. The largest repeat region in the O. xuefengensis mitogenome is 220 bp long and located in the intronic region of cox2 and ORF301, while the largest repeat region in the C. militaris mitogenome is 259 bp long and located in the intronic region of rnl. Lastly, the largest repeat region in the C. brongniartii mitogenome is 191 bp long and located between cox1 and trnR.

Seven tandem repeats were detected in both the O. xuefengensis and C. brongniartii mitogenomes, 43 in the O. sinensis mitogenome, and 5 in the C. militaris mitogenome (Table S4). The longest tandem sequence is found in the O. sinensis mitogenome, comprises 123 bp, and is located in the in the ORF138 coding region. Of all tandem repeats in the four mitogenomes, most exist in one or two copies. REPuter identified forward (F) and palindromic (P) repeats in the mitogenome of each cordyceps (Table S5), including 69 F and 34 P repeats in the mitogenome of O. xuefengensis that account for 10.16% of the total mitogenome, along with 16 F and 1 P repeats in the C. militaris mitogenome that account for 4.95% of its total mitogenome. In addition, 38 F and 19 P repeats were identified in the C. brongniartii mitogenome that account for 14.54% of its total mitogenome, and a total of 754 F, 538 P, 31 C, and 60 R repeats were identified in the O. sinensis mitogenome that account for 65.55% of its total mitogenome.

Gene rearrangements and phylogenetic analysis

The relative positions of the mitochondrial genes (including rnl, rns, rps3, the 14 core PCGs, and the tRNAs) are highly conserved in the four cordyceps mitogenomes. The O. xuefengensis mitogenome only lacks trnI, consistent with that of O. sinensis (Figure 7). Genome synteny analysis also identified several instances of gene rearrangement in the mitogenomes of the four cordyceps. All four cordyceps mitogenomes could be divided into two homologous regions (Figure 8), with the sizes and relative positions of these homologous regions substantially differing among the four species. Based on the arrangement of homologous regions, there was a high degree of synteny. And mitogenomic rearrangements were not observed among the four species, which indicates that gene recombination did not occur in these four Cordyceps mitochondrial genomes. A complete mitogenome was produced for O. xuefengensis in this study and used with 19 other complete mitogenome sequences from six families for phylogenetic reconstruction to further investigate the phylogenetic position of O. xuefengensis within the Hypocreales order. Species and NCBI accession numbers used in these phylogenetic analyses are listed in Table 1. Bayesian inference (BI) based on mitochondrial gene datasets using the GTR +I + G nucleotide substitution model yielded an identical and well‐supported topology in which all major clades are well supported. The 20 fungal species comprised six major clusters, corresponding to the major groups of the orders Hypocreales and Agaricomycetes. The four cordyceps species were divided into two groups, wherein the close relationship between O. xuefengensis and O. sinensis was highly supported (Figure 9).

FIGURE 7

Gene order comparison among four cordyceps mitogenomes. Genes (gray), rRNA (purple), and tRNA (green) conserved across all four species are shown in the same color; and genes that varied across the four mitogenomes are indicated in yellow. Cb: Cordyceps brongniartii; Cm: Cordyceps militaris; Os: Ophiocordyceps sinensis; Ox: Ophiocordyceps xuefengensi

FIGURE 8

Collinearity analysis of four cordyceps mitogenomes, generated with Mauve 2.4.0. Two homologous regions covered by two different colored blocks were detected across the four mitochondrial genomes overall. The boundaries of colored blocks usually indicate the breakpoints of genome rearrangement, unless sequence has been gained or lost in the breakpoint region. The crossing “X” pattern of lines, which happen to occur in the vicinity of the predicted rearrangements in these organisms

FIGURE 9

Molecular phylogeny of 20 fungal species based on Bayesian inference (BI) and Maximum likelihood (ML) analysis of 14 PCGs and rps3 genes. Support values are Bayesian posterior probabilities (before slash) and bootstrap (BS) values (after slash). Species and NCBI accession numbers for genomes used in the phylogenetic analyses are provided in Table 1

DISCUSSION

Cordyceps are fungal species used in traditional medicine that are represented by the well‐known fungal resource, O. sinensis. However, natural populations of O. sinensis have been overharvested to the extent that it is an endangered species. It is consequently necessary to explore alternative cordyceps resources, because artificial cultivation of O. sinensis is difficult to achieve (Zhang, Wang, et al., 2015). Ophiocordyceps xuefengensis is a new potential cordyceps resource that was found on Xuefeng mountain in Hunan. Investigation of its artificial cultivation and biometabolite biosynthesis yielded some positive results (Zhong et al., 2019). Several previous studies have investigated the phylogenetic relationships of cordyceps, although detailed phylogenetic relationships for the species remain poorly understood, especially since increasing numbers of related species have been reported (Ghikas et al., 2010; Li et al., 2015; Zhang, Zhang, et al., 2015). It is difficult to accurately classify fungal species based on limited morphological characters and overlapping morphological features (Li, Wang, et al., 2019). Consequently, mitogenomes have been widely used in the phylogenetic analysis of eukaryotes due to many advantages, including uniparental inheritance, rapid evolutionary rates, and the presence of several molecular markers (Boore, 1999). Here, we report a new mitogenome sequence for O. xuefengensis that could provide new insights into the evolutionary histories of cordyceps fungal mitogenomes and the taxonomic relationships among cordyceps.

Mitogenome sizes significantly differ among fungi, ranging from approximately 11.09 kbp (Hanseniaspora uvarum) (Pramateftaki et al., 2010) to 272.2 kbp (Morchella importuna) (Liu et al., 2020) in length. The variation in size of mitogenomes among different fungal species primarily results from the presence and extent of intergenic regions and introns (Lang et al., 2007). The Ascomycota cordyceps fungus O. sinensis contained a mitogenome with introns accounting for 68.09% of the whole mitogenome, while introns in the mitogenomes of O. xuefengensis, C. militaris, and C. brongniartii only accounted for 26.85%, 29.49%, and 24.91% of their mitogenomes, respectively. Further, intergenic regions within the O. xuefengensis mitogenome accounted for 43.93% of the whole mitogenome, but only 11.41%–17.66% of the mitogenomes of the other three cordyceps species. Thus, considerable variation in the size of introns and intergenic regions is the major reason underlying the significant differences in the mitogenome sizes of O. sinensis and O. xuefengensis. In contrast, C. militaris and C. brongniartii exhibited similarly sized mitogenomes, due to similar numbers and sizes of intergenic regions and introns.

However, the large difference in introns did not affect the taxonomically close relationship between O. sinensis and O. xuefengensis. Maximum Likelihood and Bayesian phylogenetic analyses were constructed using 14 PCGs and rps3 from 20 fungi to reveal the taxonomic relationships of these fungi. The phylogenetic analysis highly supported the close relationship between O. xuefengensis and O. sinensis. This apparent domestication of ancestral introns may indicate an adaptation to their host genome that may be due to their decreased mobility (Novikova & Belfort, 2017), or alternatively, to their potential important role in stabilizing the gene that hosts the intron (Genetics, 2003; Korovesi et al., 2018). Other introns may be acquired late in the divergence of taxa, either through HGT events or through active transposition (Mardanov et al., 2014; Wu et al., 2015). The transposition of some introns to other genomic regions with less sequence similarity can occur more frequently under stress‐induced conditions (Coros et al., 2009; Robbins et al., 2011). A recent study demonstrated that the number and Pcl of introns are highly variable between two Rhizopogon species (Li, Ren, et al., 2019). Further, the introns in nuclear genes of S. cerevisiae play crucial roles in the survival of the organism under starvation conditions (Parenteau et al., 2019). Consequently, the abundance of introns in the mitochondrial genome of O. sinensis relative to other cordyceps species may be a consequence of the particular environments they inhabit.

The accumulation of repeated evolutionary events in fungal mitogenomes could lead to over‐dispersal of repeat sequences and the introduction of new genes through HGT, thereby contributing to dynamic changes in genome structure and gene order (Aguileta et al., 2014). Moreover, the accumulation of repeat sequences could also be highly related to gene recombination and gene loss (Zou, Jakovlić, et al., 2017). A variety of models have been proposed to investigate these mitochondrial gene rearrangements (Xia et al., 2016). Here, the mitogenome of O. sinensis harbored a larger proportion of repeat sequences than the mitogenomes of O. xuefengensis, C. militaris or C. brongniartii. Interestingly, collinearity analysis revealed that the mitogenomic gene order in these four fungi was not highly variable. In addition, positive selection of core genes was not observed by Ka/Ks analysis. This differs from observations in plants, wherein mitochondrial gene order is highly variable due to high rates of recombination (Aguileta et al., 2014). These observations collectively indicate that selection upon the genomes by environmental factors is relatively weak among the four cordyceps fungi evaluated here, resulting in highly conserved phylogenetic relationships. Thus, we could also infer that differentiation of the cordyceps genomes remains in early stages, consistent with the molecular phylogenetic analyses. One reasonable explanation for these contradictory results is that the larger proportion of repeat sequences in the genome of O. sinensis exists in intronic regions, but not in intergenic regions, leading to little influence on their evolution.

Mitochondrial genomes are obtained from a common ancestor and have been widely used in population genetics and evolutionary studies (Al‐Reedy et al., 2012; Grosemans et al., 2016). Different barcode genes, including ITS, LSU, RPB2, and EF1α, have been used for phylogenetic analysis (Buyck et al., 2014; Qiang et al., 2018). However, the large number of available molecular markers in mitogenomes and their independent evolutionary histories make them attractive tools for reconstructing phylogenetic relationships (Jiang et al., 2017). Here, we reconstructed a well‐resolved phylogeny based on the combined alignment of 14 PCGs and the rps3 gene that separated 20 fungal species into major clades. In addition, K2P and Ka/Ks analysis indicated that cox1 and rps3 divergence play important roles in population differentiation. Given those reliable genetic molecular markers are critical for further understanding the phylogenetics and classification of species, these two genes could be further evaluated as barcodes for fungal species identification.

As we know, mitogenomes have been widely used in the phylogenetic analysis of plants, fungi, and animals. Studies have shown that there is a certain correlation of genetic relationship, chemical component, and therapeutic effectiveness of medicinal plant. Generally speaking, a closer genetic relationship and higher chemical similarity in plants, and more consistent therapeutic effects expected when similarly categorized (Kim et al., 2018). For example, Rauvolfia verticillata (Lour.) Baill. was used instead of R. serpentina (L.) Benth. exKurz to develop reserpine because they have similar alkaloids and efficacy (Xiao et al., 2021). The genus Scutellaria are divided into two branches by molecular systematics, and the chemical components and pharmacological effects of the same group are similar (Shen, 2021), other similar results also occurred in Malus (Li et al., 2018) and Acer (Bi et al., 2016) genus. Through a series of comparative analyses of mitochondrial genomes, the phylogenetic analysis highly supported the close relationship between O. xuefengensis and O. sinensis, suggesting O. xuefengensis has the potential to be the alternative of O. sinensis.

CONCLUSIONS

This study extends the known mitogenomes of cordyceps fungi and identified variation in PCGs, tRNA genes, and rRNA genes among medicinal cordyceps species, including O. xuefengensis. Comparative mitochondrial genome analysis of four cordyceps fungi was conducted, revealing that their mitochondrial genomes are conserved and that differentiation of their mitogenomes is still in the early stages. Nevertheless, results indicated that the cox1 and rps3 genes likely play important roles in population differentiation among the taxa. Mitogenome analyses, like those reported here, will allow further study of the population genetics, taxonomic relationships, and evolutionary biology of medicinally important cordyceps species.

CONFLICT OF INTEREST

The authors declare no conflicts of interest.

AUTHOR CONTRIBUTIONS

Can Zhong: Project administration (lead); Writing – review & editing (lead). Jian Jin: Data curation (supporting); Methodology (supporting); Software (supporting). Rongrong Zhou: Methodology (equal). Hao Liu: Investigation (equal). Jing Xie: Software (supporting); Visualization (supporting). Dan Wan: Project administration (lead). Shengen Xiao: Conceptualization (equal); Writing – review & editing (supporting). Shuihan Zhang: Conceptualization (equal); Funding acquisition (lead); Writing – review & editing (supporting).

OPEN RESEARCH BADGES

This article has earned a Preregistered Research Designs Badge for having a preregistered research design, available at https://doi.org/10.5281/zenodo.6345492.

Table S1

Click here for additional data file.

Table S2

Click here for additional data file.

Table S3

Click here for additional data file.

Table S4

Click here for additional data file.

Table S5

Click here for additional data file.