Introduction

Pleurocordyceps sinensis (Q.T. Chen et al.) Y.J. Yao et al. was first isolated from the sclerotium of Ophiocordyceps sinensis collected from Kangding, Sichuan, China in June 1980, and was described as a new species named Paecilomyces sinensis Q.T. Chen et al. (Chen et al., 1984). The species gained broad scientific attention since 1980s in China. Large numbers of pharmacological studies have been carried out using the only authentic strain termed CN80-2. The species was reported to have various pharmacological activities such as anti-implantation (Lin et al., 1988), anti-inflammatory (Li et al., 1983), anti-oxidant (Liu et al., 1987, 1989, 1991), anti-tumor (Huang et al., 1988; Wu et al., 1986), fertility regulation (Li and Lin, 1991), immune modulation (Ge et al., 1989; Lin et al., 1987; Zhang et al., 1998; Zheng et al., 1983), and treating conditions of coronary arteriosclerotic heart disease (You et al., 1986) and immunological liver injury (Cheng et al., 2005; Zeng et al., 2000). Because of the similarities in the pharmacological effects and chemical components, P. sinensis has once been recognized as the possible anamorph of O. sinensis. This viewpoint was widely accepted and cited when summarizing the anamorph of O. sinensis (Fang, 1991; Liang, 1991; Liu, 1990). The idea of using Paecilomyces sinensis as a substitute of O. sinensis was thus proposed (e.g., Cheng et al., 2005; Li et al., 1983; Zeng et al., 2000). However, several independent researches based on molecular evidences rejected the anamorph-teleomorph relationship between P. sinensis and O. sinensis (Chen et al., 2001; Jiang and Yao, 2002, 2003; Li et al., 2000; Zhao et al., 1999). Wang et al. (2012) placed the species in Polycephalomyces based on morphological and molecular analyses and found Polycephalomyces species formed a new clade of clavicipitaceous fungi and stated that this new clade is distinct from the known families of Hypocreales.

The genus Polycephalomyces was identified by Kobayasi (1941) with P. formosus as the type. Only three species, that is, P. paludosus, P. cylindrosporus, and P. tomentosus, were described in the last century (Mains, 1948; Samson et al., 1981; Seifert, 1986). Until recently, after the recombination of Paecilomyces sinensis (Wang et al., 2012) and several species of Cordyceps s. l. into the genus (Kepler et al., 2013), more and more new species have been discovered and described, especially from China and Southeast Asia (Crous et al., 2017; Wang et al., 2015a, 2015b; Xiao et al., 2018; Yang et al., 2020). A total of 24 species names are currently recorded by the Index Fungorum (4 April 2022, http://www.indexfungorum.org/Names/Names.asp), among which 4 were segregated to Perennicordyceps (Matočec et al., 2014). Wang et al. (2021) recently proposed a new genus Pleurocordyceps for one of the subclades within the “Polycephalomyces clade” (Polycephalomyces sensu lato). Ten species were included in the new genus including P. sinensis. In multigene phylogenetic analyses, species of Polycephalomyces s. l. usually formed a distinct clade sister to Ophiocordycipitaceae, although this sister relationship did not receive much statistical confidence (Kepler et al., 2013; Wang et al., 2021). In other words, the family-level taxonomic position of Polycephalomyces s. l. was not fully resolved; species in this group were tentatively placed in Ophiocordycipitaceae in most researches (e.g., Kepler et al., 2013; Xiao et al., 2018), which was accepted by the Index Fungorum. Polycephalomyces s. l. may represent a new family that is different from the three existing families of clavicipitoid fungi (Wang et al., 2021), that is, Cordycipitaceae, Clavicipitaceae, and Ophiocordycipitaceae. However, more evidence from morphology and molecular phylogenetics is required to support this hypothesis.

Mitochondria play various essential roles in eukaryotic cells, including respiratory metabolism, energy production, calcium homeostasis, and are also involved in cell death and aging (Basse, 2010). Mitochondrial (mt) genomes usually have a rapid rate of evolution compared with nuclear genomes, and thus are considered as powerful tools in evolutionary biology (Berbee and Taylor, 2001; Chris et al., 1994). Previous studies revealed that the gene contents and synteny of mt genomes of hypocrealean species were largely conserved, but in the meantime, the genome sizes expanded greatly in certain species such as O. sinensis (Li et al., 2015). Complete mt genomes have been reported for a number of species of the three clavicipitaceous fungal families, that is, Cordycipitaceae (Fan et al., 2019; Kouvelis et al., 2004; Sung, 2015; Zhang et al., 2021), Ophiocordycipitaceae (Abuduaini et al., 2021; Li et al., 2015; Zhang et al., 2016; Zhang and Zhang, 2020), and Clavicipitaceae (Sun et al., 2021; Winter et al., 2018), while they have not been reported for Polycephalomyces s. l. species so far.

In this study, the complete mt genome of the type strain (CN 80-2) of the species Pleurocordyceps sinensis was sequenced, described, and compared with other hypocrealean species. Phylogenetic analyses using 14 conserved protein sequences were also performed to study the phylogenetic relationship between this species and other clavicipitaceous fungal groups.

Materials and Methods

Fungal isolation and cultivation

The ex-type strain (CN80-2) of Pleurocordyceps sinensis used in this study was isolated from a sclerotium of O. sinensis collected from Kangding, Sichuan, China in June 1980 (Chen et al., 1984). Stock strain was maintained at 4°C on potato dextrose gar (PDA) slants. Seed cultures were grown in 250-mL Erlenmeyer flasks, containing 50 mL liquid potato-dextrose medium, by shaking at 100 rpm at 25°C for 10 days. Mycelia were harvested and washed with distilled water using vacuum filtration to remove extracellular polysaccharides, frozen with liquid nitrogen, and vacuum freeze-dried using a freeze dryer. Dried mycelia were then sent to the genome sequencing company.

DNA extraction and genome sequencing

Genomic DNA was extracted by the sequencing company using TIANamp Yeast DNA Kit (TIANGEN Biotech Co., Ltd., Beijing, China) according to the manufacturer’s instruction. The amount and quality of total DNA were visualized by 1% agarose gel electrophoresis and quantified with a Qubit2.0^® Fluorometer (Life Technologies, New York, USA). A 20 K library was prepared from sheared genomic DNA (containing both mt and nuclear sequences) using a 20-Kb template library preparation workflow. Twelve single molecule real time (SMRT) sequencing cells were sequenced on PacBio RS II sequencing platform (Pacific Biosciences, Menlo Park, CA) with P6 polymerase and C4 sequencing chemistry at Tianjin Biochip Corporation (Tianjin, China).

Mitochondrial genome assembly and annotation

Mt genome of the strain CN80-2 was assembled and annotated following a procedure described in Li et al. (2015). The adapter sequences, reads with length <50 bp, or average quality <0.75 (defined as low quality) were filtered before assembling. The mt sequences were extracted from filtered reads matching each read against the fungal mt genome database (https://www.ncbi.nlm.nih.gov/genome/browse#!/organelles/), preassembled and corrected using BLASR (Chaisson and Tesler, 2012). Corrected reads were retained and then re-assembled with the Celera Assembler program (Myers et al., 2000). The assembly was further refined with Quiver (Chin et al., 2013). A circular double-stranded DNA was finally obtained and proceeded to an online annotation tool MFannot using the Mold, Protozoan, and Coelenterate Mitochondrial Code (Beck and Lang, 2010). The annotated mt genome was submitted to GenBank under the accession number OK017430. The annotated mt genetic map was generated by Circos software (Krzywinski et al., 2009) and modified with Adobe Illustrator^® CS5 (Version 15.0.0, Adobe^®, San Jose, CA).

Phylogenetic analyses and comparative genomics

All the 82 complete mt genomes available from GenBank (accessed on 28 March, 2021) within the order Hypocreales were downloaded and used for phylogenetic and/or comparative genomic analysis. Among which Ophiocordyceps camponoti-floridani EC05 (CM022976) was used only for genome comparison but not included in phylogenetic analyses, as it seems to be incorrectly assembled; an unclassified mt genome (NC_049089) was also excluded due to the large numbers of possible sequencing errors or possible assembly mistakes. A number of mt genomes were found to be incorrectly or incompletely annotated. For example, the rps3 gene was not predicted in a number of species. Those genomes were re-annotated following the same procedure used in this study with missing genes replenished and wrongly predicted genes manually corrected. A phylogenetic tree was constructed using 14 conserved protein-coding genes (cox1, cox2, cox3, cob, atp6, atp8, atp9, nad1, nad2, nad3, nad4, nad5, nad4L, and nad6). Protein sequences were aligned with BioEdit version 7.0.9.0 (Hall, 1999) and refined manually. Maximum likelihood (ML) phylogenetic analyses were performed with RAxML v.7.2.659 (Stamatakis, 2006) using the LG substitution matrix and default parameters. Bootstrap values were calculated with 1000 re-sampling iterations using an approximate likelihood ratio test. Three mt genomes from two species of the order Glomerellales, that is, Colletotrichum lindemuthianum (NC_023540) and Verticillium dahliae (NC_008248 and CM019738), were used as outgroups (Supplementary Table S1).

The gene contents and synteny of mt genomes within the order Hypocreales and related outgroup species were compared and analyzed.

Results

Genome sequencing and assembly

A total of 23,757 reads (171,039,810 bp) were identified as mt among 601,168 reads (5,005,308,071 bp) of the raw sequencing output for the whole genome of Pleurocordyceps sinensis. The lengths of the putative mt reads ranged from 276 bp to 45,245 bp with an average length of 7200 bp, reaching a coverage depth of 5371× over the mt genome of the species. The mt reads were passed through the program BLASR and assembled with Celera Assembler program and Quiver, resulting in a circular DNA of 31,841 bp (Figure 1).

Figure 1. Genetic map of the mitochondrial genome of Pleurocordyceps sinensis. Shading blocks with orange frame indicate exons of predicted coding genes, with gene names labeled on the inner side; rnl and rns genes were marked with slashes; the rps3 was in blue and nested within an intron of the rnl gene; two introns that located in rnl and cox1, respectively, were shown in blank. Thin flint lines protruding outside of the outer circle indicate tRNAs. Ten predicted ncORFs were shown as light blue blocks. The red circle represented sequence coverage with the highest at 5,200× and the lowest at 4,300× (average at 4,750×). The outer and inner diagram of curves represented GC content and GC skew, respectively.

Conserved protein genes and nonconserved open reading frames

The mt genome of P. sinensis had a low GC content of 25.46% and encoded 15 protein genes conserved within the order Hypocreales, including seven subunits of the electron transport complex I (nad1, nad2, nad3, nad4, nad4L, nad5, and nad6), cytochrome b (cob), three subunits of complex IV (cox1, cox2 and cox3), three F0 subunits of the ATP-synthase complex (atp6, atp8, and atp9), and the rps3 gene, which encodes 40S ribosomal protein S3 (Table 1, Figure 1). In addition to those genes, 10 non-conserved open reading frames (ncORFs) (7194 bp totally in length) were also predicted, among which two (ncORF3 and ncORF9) were found to encode homing endonucleases (HEs) with motif patterns GIY-YIG and LAGLIDADG, respectively (Table 1).

All conserved protein coding genes and ncORFs were found on the positive strand and oriented clockwise except for ncORF1 and ncORF2, which were on the negative strand and anticlockwise oriented. It was found that the nad2/nad3 genes were joined and nad4L/nad5 genes were fused, that is, the initial codon of the nad3 gene (ATG) followed the terminal codon of the nad2 gene (TAA), and the terminal codon of nad4L (TAA) uses the same nucleotide A with the initial codon (ATG) of nad5 (Figure 1, Table 1). Other protein-coding genes and ncORFs were separated by either long or short intergenic regions (Figure 1).

The 15 protein-coding genes and 10 predicted ncORFs employed the standard fungal mt start codon ATG, except the cox1 and ncORF10, which were initiated by ATA. In addition, 24 of those genes used TAA as the stop codon except the ncORF8, which used TAG (Table 1).

Noncoding RNAs

In addition to the 15 protein-coding genes, a large and a small ribosomal RNA (rnl and rns, respectively) and 25 tRNA genes were also identified (Table 1). The tRNA genes ranged in size from 71 to 87 bp and could correspond to 20 amino acids. A majority of amino acids were coded by only one tRNA gene; however, Serine (Ser), Arginine (Arg), Methionine (Met), and Leucine (Leu) had 2, 2, 3, and 2 tRNA genes, respectively (Table 1). All noncoding RNAs (tRNA, rRNA) were found on the positive strand and oriented clockwise.

Intronic and intergenic regions

Exons of protein-coding genes, rRNA and tRNA genes, had a total length of 20,873 bp accounting for 65.55% of the mt genome. Ten ncORFs (7,194 bp) accounted for 22.59% of the mt genome. Only two introns (group I) were predicted, including one further classified into subgroup IA (1,643 bp) in rnl and one classified into subgroup IB in cox1 (1,036 bp), respectively, making up 8.5% of the entire mt genome. The intergenic sequences had a total length of 1,070 bp covering 3.4% of the genome.

Gene component and synteny

Although different numbers of ncORFs (hypothetical proteins) would be predicted for hypocrealean fungi, the content and synteny of 15 protein-coding genes remained largely conserved, except in a few cases. For instance, location of the cox2 gene shifted in three species of Acremonium chrysogenum, A. fuci, and Clonostachys rose comparing to other Hypocreales; an additional copy of rps3 and atp9 gene was found in Beauveria malawiensis and Fusarium solani IISc-1; an extra copy of nad1 and nad4 was found from F. oxysporum UASWS AC1 (KR952337); the location of the two genes were found to be reversed in F. oxysporum f. sp. matthiolae (CM019668); and the mt genome of F. oxysporum f. sp. fragariae GL1381 (CM029251) was found to lose cox3 and nad6 genes and possess an extra reversed copy of genes of cob, cox1, nad1, nad4, atp8, and atp6; an even extreme case was found in Sarocladium implicatum in which three genes (cob, cox3, and nad6) were lost and the nad4 gene shifted its location from the nad1-atp8 junction to a position between rps3 and nad2 (Supplementary Table S1).

Phylogenetic analyses

Eighty-one complete mt genomes representing 63 distinct species from the order Hypocreales were included in phylogenetic analyses. After excluding the ambiguous aligned regions, a total of 4,345 amino acid sequences of 14 conserved proteins were retained. All species of Hypocreales formed a well-supported clade (BP = 100%) in ML analysis. Within the clade, four family-level subclades were recognized with strong supports (BP = 100%), that is, Nectriaceae, Bionectriaceae, Hypocreaceae, and Clavicipitaceae (Figure 2). Species in the family Ophiocordycipitaceae were clustered into two subclades, one subclade that consists of four Tolypocladium species showed a sister group relationship with the Clavicipitaceae clade with low bootstrap support (BP = 75%), the other highly supported (BP = 100%) subclade comprised four Hirsutella species (H. minnesotensis, H. rhossiliensis, H. thompsonii, and H. vermicola) and O. sinensis (Figure 2). It is interesting to find that P. sinensis clustered with the Clavicipitaceae clade with 100% bootstrap support rather than grouped with either two subclades of Ophiocordycipitaceae.

Figure 2. Phylogenetic relationships of Hypocreales inferred from 14 conserved protein sequences (cox1, cox2, cox3, cob, atp6, atp8, atp9, nad1, nad2, nad3, nad4, nad4L, nad5, and nad6) using the maximum likelihood method. Bootstrap values were shown above the branches. Three species in the order Glomerellares were used as outgroups.

Discussion

The complete mt genome of the ex-type strain CN 80-2 of the species P. sinensis described here is the first reported case for the newly proposed genus Pleurocordyceps (Wang et al., 2021) and Polycephalomyces s. l. It is rather compact compared with other Hypocreales species, especially O. sinensis, from which the species was isolated. The genome sizes of the two sequenced isolates of O. sinensis were 157,510 bp (KP835313) and 157,539 bp (NC_034659), respectively, almost five times larger than P. sinensis. The reason for this remarkable size variation was considered to be the presence of large numbers of repetitive regions, which mainly consisted of intronic mobile elements of HEs and reverse transcriptases (RTs) (Li et al., 2015). In the expanded mt genome of O. sinensis, 32 HEs genes (21 LAGLIDADG and 11 GIY-YIG endonuclease) and 10 RTs genes were found in group I and group II introns, respectively (Li et al., 2015), while only two HEs genes, that is, one GIY-YIG (ncORF3) and one LAGLIDADG endonuclease (ncORF9) were found in P. sinensis. After comparing mt genomes of three isolates of the same species, it was also found that larger genomes contained more introns and more intronic HEs genes in Cordyceps militaris (Zhang et al., 2015).

As mt genome sizes varied greatly in hypocrealean fungi, ranging from 22,376 bp of S. implicatum (Yao et al., 2016) to 272,497 bp of Ophiocordyceps camponoti-floridani (Will et al., 2020), it will be interesting to see whether the mt genome size variation was related to expansion of mobile elements of HEs and RTs and how those elements evolved. Megarioti and Kouvelis (2020) recently proposed an “aenaon” model for the evolution of HEs genes and their host introns; thus, free-standing introns and HEs genes were the ancestral form and could invade intron-free coding genes together; HEs genes and their host introns coevolved through recombination, transposition, and horizontal gene transfer. As observed in this case, the two HEs genes found in P. sinensis were located in intergenic regions (free-standing or sole mobile in other words) of tRNA genes (Table 1) rather than invaded into intronic regions of coding genes (intron homing). While in O. sinensis, HEs and RTs genes were all found to be intronic, either in group I or group II introns (Li et al., 2015). It may indicate that the species was earlier diverged than O. sinensis according to the “aenaon” model. Considering that the species was isolated from O. sinensis and might be a fungal parasite of the latter, and moreover, several other species in Polycephalomyces s. l. have often been found to associate with entomopathogenic Cordyceps s. l. (Kobayasi, 1941), it is reasonable to hypothesize that species of Polycephalomyces s. l. gained hyperparasitic ability to entomophagous fungi during the evolutionary process.

Despite the size variation of mt genome, the gene contents and synteny (gene order) are largely conserved within the order Hypocreales, generally encoding 15 known proteins and 2 rRNAs (rnl and rns) (Li et al., 2015). The genome size variation observed in hypocrealean fungi was probably not associated with taxonomic classification since notable variation was also observed within the same genus or even within the same species. As shown in Supplementary Table S1, mt genome sizes varied from 30,629 bp to 110,525 bp in the genus Fusarium, and from 34,477 bp to 52,424 bp within the species of F. oxysporum. This variation is largely due to the presence of various introns and the lengths of intergenic regions (Burger et al., 2003).

Although remarkable variation in terms of gene order, genome size, composition of intergenic regions, and presence of repeats, introns, and associated ncORFs have been observed between the major fungal phyla (Aguileta et al., 2014), this variance may not occur within the same fungal groups (order or below). As observed in this study, the genome size, composition of intergenic regions, and presence of repeats, introns, and associated ncORFs varied within the order Hypocreales; however, the gene content and synteny remained highly conserved even though a few exceptional cases were observed (listed in Supplementary Table S1). A part of these exceptions were probably due to the incorrect assembly. It would be interesting to know the mt genome evolutional process, that is, gene gain and loss events happened during the evolutionary history of different major fungal groups.

Most protein-coding genes and ncORFs used standard mt initial and terminal codons (ATG and TAA, respectively) in P. sinensis, except the ncORF8, cox1, and ncORF10. ncORF10 and cox1 were initiated by ATA, and ncORF8 was terminated by TAG. It is noteworthy that cox1 is usually found to use nonstandard start codons such as TCG, ACC, CGA, CTA, CCG, and AAA in insect mt genomes (Fenn et al., 2007; Wei et al., 2010), and ATA has been recorded to be used as the initial codon in organisms like Pseudocohnilembus persalinus (Gao et al., 2018), Wellcomia siamensis (Park et al., 2011), and Calanus sinicus (Wang et al., 2011). Although most hypocrealean species used ATG as the initial codon of cox1 gene, exceptional cases were also reported in Hirsutella rhossiliensis (NC_030164) and Calonectria ilicicola (Gai et al., 2020), in which TTG were used.

The mt genome released in this study provided additional evidence that P. sinensis is not the anamorph of O. sinensis but represents another fungus, and moreover, P. sinensis was found to cluster with Clavicipitaceae rather than Ophiocordycipitaceae species with very strong supports (BP = 100%) in ML phylogenetic analyses (Figure 2). It is also noteworthy that the family Ophiocordycipitaceae was paraphyletic, which contradicts previous studies applying multi-gene phylogeny (Sung et al., 2007a) although the paraphyly was not well supported (Figure 2). It should be clarified whether those contradictions were due to the incongruence of phylogenies revealed by different molecular markers since mt DNA may tell different evolutionary stories than nuclear genes (Burger et al., 2003), or were just caused by the insufficient taxon sampling or analytical difference. Sung et al. (2007b) conducted multi-gene phylogenetic analyses of clavicipitaceous fungi and compared the performance of seven loci including the nuclear ribosomal small and large subunit DNA (nrSSU and nrLSU), ß-tubulin, elongation factor 1 α (EF-1α), the largest and second largest subunits of RNA polymerase II (RPB1 and RPB2), and one mt protein-coding gene ATP Synthase subunit 6 (mtATP6), and found that seven genes gave incongruent topologies in higher-level relationships from each other and also from the combined dataset. It also showed that the only mt fragment (mtATP6) used in the study possessed localized incongruence and simultaneously provided an increased level of support for certain nodes. Phylogenetic incongruence revealed by different markers, especially those from mt and nuclear fragments, respectively, has been frequently reported and compared in different organisms (e.g. Kimball et al., 2021; Mikula et al., 2021; Zhang et al., 2021). It still remains unclear as to whether the nuclear genome sequences (fragments or whole genome data) or the mt genome sequences (fragments or complete data) could provide a better resolution of fungal phylogeny.

As in the case of clavicipitaceous fungi, almost all the later publications on taxonomy and phylogenetic studies accepted the backbone phylogeny created by Sung et al. (2007a, 2007b), and continued to use the five gene dataset (e.g. Kepler et al., 2013; Wang et al., 2021; Xiao et al., 2018). While the above two studies (Sung et al., 2007a, 2007b) failed to include species of Polycephalomyces s. l. Kepler et al. (2013) then included several species of this group and found those species represented a clade distinct from other clavicipitoid genera, and treated them as incertae sedis of Hypocreales. Further studies are needed to reconstruct a reliable phylogenetic relationship of clavicipitaceous fungi, especially the assignment of species of Polycephalomyces s. l. Since an increasing number of whole mt genomes have recently been sequenced and released for hypocrealean species, and even more are being proceeded, the plentiful phylogenetically informative sites from the conserved protein coding genes of the whole mt genome would provide valuable information for phylogenetic reconstruction of this important fungal group. While publically released data should be carefully treated since they could probably include assembly and annotation errors as observed in Ophiocordyceps camponoti-floridani EC05 (CM022976) and Ophiocordycipitaceae sp. (NC_049089), correct annotation and characterization are always necessary (Kortsinoglou et al., 2019; Megarioti and Kouvelis, 2020).