Comparative mitochondrial genome analyses reveal conserved gene arrangement but massive expansion/contraction in two closely related Exserohilum pathogens.

Exserohilum turcicum and E. rostratum, two closely related fungal species, are both economically important pathogens but have quite different target hosts (specific to plants and cross-kingdom infection, respectively). In the present study, complete circular mitochondrial genomes of the two Exserohilum species were sequenced and de novo assembled, which mainly comprised the same set of 13 core protein-coding genes (PCGs), two rRNAs, and a certain number of tRNAs and unidentified open reading frames (ORFs). Comparative analyses indicated that these two fungi had significant mitogenomic collinearity and consistent mitochondrial gene arrangement, yet with vastly different mitogenome sizes, 264,948 bp and 64,620 bp, respectively. By contrast with the 17 introns containing 17 intronic ORFs (one-to-one) in the E. rostratum mitogenome, E. turcicum involved far more introns (70) and intronic ORFs (126), which was considered as the main contributing factors of their mitogenome expansion/contraction. Within the generally intron-rich gene cox1, a total of 18 and 10 intron position classes (Pcls) were identified separately in the two mitogenomes. Moreover, 16.16% and 10.85% ratios of intra-mitogenomic repetitive regions were detected in E. turcicum and E. rostratum, respectively. Based on the combined mitochondrial gene dataset, we established a well-supported topology of phylogeny tree of 98 ascomycetes, implying that mitogenomes may act as an effective molecular marker for fungal phylogenetic reconstruction. Our results served as the first report on mitogenomes in the genus Exserohilum, and would have significant implications in understanding the origin, evolution and pathogenic mechanisms of this fungal lineage.

Introduction

Exserohilum is an important plant-pathogenic genus belonging to the order Pleosporales (Dothideomycetes, Ascomycota), and until now, a total of 40 Exserohilum species have been recorded (https://www.indexfungorum.org). This taxonomic group was originally described as Helminthosporium, which was then subdivided into Bipolaris, Curvularia and Exserohilum based on the traditional taxonomical criteria especially by their hilum morphology [1], [2], [3]. The classification result was also supported by multigene phylogenetic analyses based on nine nuclear loci including ITS, LSU, act, tub2, cam, gapdh, his, tef1 and rpb2 [3]. Although known as Setosphaeria (teleomorphic stage) as well, the scientific name of anamorph Exserohilum was generally recommended for use according to the International Code of Nomenclature for algae, fungi and plants [3], [4].

As the type species of Exserohilum, E. turcicum, a filamentous hemibiotrophic fungus, is notorious for causing northern corn leaf blight (NCLB) [3], [5]. Of late years, this disease has become increasingly problematic worldwide with the extensive planting of susceptible varieties of maize in many countries [6], [7]. The most common symptom of NCLB is necrotic green-gray lesions on leaves with a length range from 2 to 14 cm, which could lead to 10–40% grain yield losses in severe cases [8]. Besides, the fungal pathogen is also an important blight agent of some Sorghum spp. [9]. By contrast, another noteworthy species of this genus, E. rostratum, has a more complex life style. It was first described as the pathogen of thrush grass in 1923 in the USA [10], and then reported to infect more than 30 plant species belonging to 28 genera of 11 families [11], including gramineous crops like corn and rice [12], woody plants like rubber tree [13], and some tropical cash crops like sugarcane [14]. Recently, E. rostratum seemed to be better known as an opportunistic pathogen found in humans, causing keratitis, meningitis, subdermal fungal disease, etc. [15], [16], [17]. Between the above two closely related Exserohilum species with fairly different target hosts, comparative analysis was rarely performed. Even if their nuclear genomes have been already available, it is not clear how E. turcicum and E. rostratum differ in their virulence strategies and molecular evolution mechanisms [15], [18].

Mitochondria are a kind of double membrane organelles well-known for their major role in energy supply, and therefore commonly referred to as the power factories of eukaryotes. In addition, they can also participate in the ion homeostasis, intermediate metabolism, cell senescence and apoptosis [19]. In recent years, an increasing body of evidence suggested that mitochondria were derived from the symbiotic combination of α-proteobacteria and eukaryotic cells [20], [21]. And their evolutionary processes varied among the main eukaryote groups including animals, plants and fungi [22], of which, however, fungal mitochondria shared the least spotlight. The mitochondrial genomes (mitogenomes) of different fungal species showed significant differences in genome size, gene arrangement, ratio of repetitive and intergenic regions, and number of introns and open reading frames (ORFs) [23], [24]. Among all published fungal mitogenomes, Hanseniaspora uvarum has the smallest one that is only 18.84 kb in size [25], while the mitogenome of Morchella crassipes is the largest, approximately 531.19 kb in size [26]. In spite of the varying sequence lengths, fungal mitogenomes usually contain a set of conserved protein-coding genes (PCGs) including seven NADH deaminase complex subunits (nad1-6 and nad4L), three cytochrome c oxidase subunits (cox1-3), one cytochrome b gene (cob) and three ATP synthase subunits (atp6, atp8 and atp9), two ribosomal RNA (rRNA) genes, and 22–26 transfer RNA (tRNA) genes [27], [28]. Studies have found that expansion and contraction of the mitogenomes in different fungi were mainly due to the number and length of introns, repeat rate and introduction of new genes through horizontal transfer [29], [30]. There are two types of introns in fungal mitogenomes, namely group I and group II, in which homing endonuclease genes (HEGs) were often detected and could produce LAGLIDADG/GIY-YIG motif endonucleases [31], [32]. With the development of the next-generation sequencing (NGS) technology and continuous reduction of sequencing costs, more than 12,500 complete mitogenome sequences of eukaryotes have been released to date in the NCBI database (https://www.ncbi.nlm.nih.gov/genome/browse/#!/organelles/). Nonetheless, for most fungal taxa including the genus Exserohilum, there has been no research on their mitochondrial genomes, which will limit our knowledge of the “second genome” in some significant fungal lineages.

In the current study, mitogenomes of the two Exserohilum pathogens, E. turcicum and E. rostratum, were sequenced, assembled, annotated, and compared with some other Pleosporales species. Our main objectives are as follows: 1) to depict the mitogenome contents, structures and organizations of the two Exserohilum species; 2) to perform comparative mitogenomic analysis between the different Exserohilum species and reveal their variations and conservations; 3) to reveal the intron dynamic changes of cox1 genes in 13 Pleosporales species including the two Exserohilum fungi; 4) to clarify the phylogenetic status of Exserohilum in the phylum Ascomycota based on combined mitochondrial gene sets. This study served as the first investigation on Exserohilum mitogenomes, which will lay the foundation for further understanding the genetic evolution, species differentiation, and ecological adaptation of Exserohilum species.

Materials and methods

Fungal isolates, DNA extraction and genome sequencing

With the method of tissue isolation [11], both E. turcicum and E. rostratum strains (ZM10601 and ZM170581) were obtained from corn leaves in the city of Zhengzhou (34°48′ N, 113°39′ E), Henan province, China. Their hyphal tips were transferred onto PDA plates for purification, and after sufficient growth, these culture plates were kept at 4 ℃ for short-term storage. For long-term preservation, however, the fungal colonies needed to be rinsed off twice with sterile water for collecting conidia, which were then stored in 15% glycerol at − 80 ℃ in the fungal collection of Henan Agricultural University. Species identification was based on their morphological characters and sequence analysis of the ITS, LSU, gapdh, tef1 and rpb2 genes. These sequences have been deposited in GenBank under the accession numbers listed in Supplementary Table 1. Total genomic DNA (gDNA) extraction of the two specimens collected by us was performed by using the cetyl trimethyl ammonium bromide (CTAB) method [33]. To ensure that each DNA sample was of adequate quality for PCR, UV spectrophotometry (Nanodrop 2000, Thermo, Wilmington, DE) was used to determine the concentration and purity of DNA.

High-quality gDNA was then sent to the Novogene Co., Ltd. (Tianjin, China) for library preparation and genome sequencing. Short insert libraries (350 bp) were constructed using the NEBNext® Ultra DNA Library Prep Kit for Illumina (NEB, USA). Whole-genome sequencing (WGS) was carried out on an Illumina Hiseq X Ten platform, generating 150 bp paired-end reads for each sample.

Mitogenome assembly and annotation

Approximately 5 Gb of raw data were obtained through WGS, and were then trimmed of adapters and low-quality reads (ratio of the bases less than Q20 > 30% or containing undetermined bases) by using the fastp v0.13.1 [34]. FastUniq v1.1 [35] and Musket v1.1 [36] were successively used for duplicate removing and error correction. Cleaned paired-end reads were de novo assembled using the SPAdes v3.14.1 software with k-mers of 21, 33, 55, 77, 99 and 127 [37]. Mitochondria-related contigs were identified and pooled for each sample by performing BLASTn search against the reference mitogenomes from B. cookei [38] and B. sorokiniana [39], which are closely related with the genera Exserohilum. We then filled gaps between these contigs by using the MITObim v1.9 to build closed-circular mtDNA macromolecules of the two Exserohilum samples [40]. In addition, both the assembled mitochondrial sequences were verified by NOVOPlasty [41].

The obtained complete mitogenomes of E. turcicum and E. rostratum were annotated according to the methods previously described by Wu et al. [42]. In short, we initially used Mfannot [43] and MITOS [44] to predict the PCGs, introns, unidentified open reading frames (uORFs), and non-coding RNA (ncRNA) genes among the mitogenome sequences, both based on genetic code 4. Then, the NCBI Open Reading Frame Finder (https://www.ncbi.nlm.nih.gov/orffinder) were used to modify or predict the PCGs and uORFs, which were further annotated by BLASTp searches against the NCBI non-redundant protein sequence database [45]. We verified the intron–exon borders of PCGs by using the exonerate v2.2 software [46] with the B. cookei mitogenome as reference [38]. Secondary structures of the tRNA genes were predicted by MITOS with default parameters [44], and were redrawn in Adobe Illustrator CS6. Graphical maps of the two Exserohilum mitogenomes were drawn using the online program OGDraw v1.2 [47].

In addition, for further comparative analyses, we downloaded another 11 complete mitogenomes of Pleosporales species (including B. cookei and B. sorokiniana that were previously used as references) from the NCBI GenBank database. Their detailed information (such as accession numbers) is provided in Supplementary Table 2.

Analyses of sequences and repetitive elements

Strand asymmetries of the 13 Pleosporales mitogenomes were assessed based on the following formulas: AT skew = [A - T] / [A + T], and GC skew = [G - C] / [G + C] [48]. Mitochondrial genome collinearity analysis of the two Exserohilum species and three closely related Bipolaris species was conducted by using Mauve v2.4.0 [49]. Based on genetic code 4, the sequence manipulation suite [50] was used to analyze the codon usage frequency and preference of the two Exserohilum mitogenomes. The DnaSP v6.10.01 software [51] was used to calculate synonymous (Ks) and nonsynonymous substitution rates (Ka) for 12 core PCGs (atp6, cob, cox1, cox2, cox3, nad1, nad2, nad3, nad4, nad4L, nad5, and nad6) in all 13 acquired Pleosporales mitogenomes. Based on the Kimura-2-parameter (K2P) substitution model, MEGA v6.06 [52] was used to calculate the genetic distances between each pair of the 12 core PCGs. We performed BLASTn searches [53] of the two Exserohilum mitogenomes against themselves to detect if there were interspersed repeats or intragenomic duplications of large fragments, based on an E-value of < 1e-10. Tandem Repeats Finder [54] was used to detect tandem repeats (length > 10 bp) within the two Exserohilum mitogenomes. In addition, the contribution rates of different genetic compositions to the mitogenome expansion/contraction of two Exserohilum species were calculated according to the formulas: [A - B] / [C - D] (A: size of one genetic component in the larger mitogenome; B: size of the same genetic component in the smaller mitogenome; C: size of the larger mitogenome; D: size of the smaller mitogenome).

Comparative mitogenomic and intron analyses

A comparative mitogenomic analysis was performed to evaluate the variations and conservations between different Pleosporales mitogenomes in genome size, base composition, GC content, gene number, intron number, gene arrangement, and gene content. According to the method described by Férandon et al. [55], introns within the cox1 genes of the 13 Pleosporales mitogenomes could be divided into different position classes (Pcls). Taking the cox1 gene of B. cookei as reference [38], cox1 genes from the other 12 Pleosporales species were aligned by Clustal W [56] to detect the insertion sites of introns. Pcls were named according to their insertion sites in the corresponding reference sequence. The same Pcls from different species usually have high sequence similarity and contain homologous intronic ORFs.

Mitochondrial phylogenomic analysis

In order to investigate the phylogenetic positions of the two Exserohilum species, we constructed a phylogenetic tree composed of 98 ascomycetes based on the concatenated mitochondrial gene set (nad1-6, nad4L, cox1-3, cob, rps3, atp6, atp8 and atp9). Both Taphrina deforman and T. wiesneri from Taphrinomycetes were appointed as the outgroup [57]. The software MEGA v6.06 [52] was first used to align individual mitochondrial genes, and then the SequenceMatrix v1.7.8 [58] was used to concatenated them into a combined mitochondrial gene set. Partition homogeneity test was used to detect potential phylogenetic conflicts between different mitochondrial genes. PartitionFinder v2.1.1 [59] was applied to determination of best-fit models of phylogeny and partitioning scheme of the gene set. We performed phylogenetic analysis using both Bayesian inference (BI) and maximum likelihood (ML) methods [60]. The ML analysis was implemented in the RAxML v8.2.12 with 1,000 bootstrap replications under the substitution model GTRGAMMA [61]. The BI analysis was performed using the MrBayes v3.2.6 [62]. Four simultaneous Markov chains were run starting from a random tree for 2,000,000 generations and trees were sampled every 100 generations. The first 25% of samples were discarded as burn-in, and the remaining trees were used to calculate values of Bayesian posterior probabilities (BPPs) in a 50% majority-rule consensus tree. Tree rendering was carried out using the FigTree v1.4.4 (https://tree.bio.ed.ac.uk/software/figtree/).

Data availability

The complete mitogenomes of E. turcicum and E. rostratum were deposited in the GenBank database under the accession numbers OK381869 and OK377062, respectively; and their raw sequencing data were deposited in the Sequence Read Archive (SRA) database under the accession numbers SRR16781348 and SRR16781347, respectively.

Results

Features, genetic compositions and PCGs of two Exserohilum mitogenomes

Mitochondrial genomes of the two Exserohilum species, E. turcicum and E. rostratum, were circularly assembled, with the total sizes of 264,948 bp and 64,620 bp, respectively (Fig. 1A). Their GC contents were very close, with the values of 29.76% and 29.00%, respectively (Supplementary Table 2). For the E. turcicum mitogenome, both its AT and GC skews were positive, while in E. rostratum negative AT skew and positive GC skew were presented. Intronic regions occupied the largest proportion for both the two Exserohilum mitogenomes, reaching 62.25% and 31.53%, respectively (Fig. 1B). In the mitogenome of E. turcicum, protein-coding sequence was the second largest region accounting for 19.26%, whereas it ranked the third in E. rostratum, accounting for 27.64%. Intergenic spacers accounted for 14.25% and 29.63%, respectively, of the two mitochondrial genomes. Their ncRNA genes (including tRNAs and rRNAs), took up the smallest proportions, reaching only 4.23% and 11.20%, respectively.

Fig. 1

Circular maps and genetic compositions of two Exserohilum mitogenomes. (A) Circular maps of the mitogenomes of two Exserohilum species, E. turcicum (left) and E. rostratum (right). Colored blocks along the outer rings represent the mitochondrial core protein-coding genes, rRNAs, tRNAs, and non-conserved ORFs. Inner circles show the GC contents. (B) Proportions of the protein-coding, intronic, intergenic, and ncRNA regions in the two Exserohilum mitogenomes.

Both the two newly sequenced Exserohilum mitogenomes had 12 typical core PCGs (atp6, cob, cox1-3, nad1-6, and nad4L) and a conserved core gene coding for the putative ribosomal protein S3 (rps3). Particularly, the cox1 genes in the two mitogenomes were seamlessly adjacent to the cox2, just resembling a fused gene. The mitogenome of E. turcicum contained two pseudogenes, which were respectively part of the two core PCGs, atp6 and nad2 (Supplementary Table 3). In addition, 34 free-standing ORFs were detected in the E. turcicum mitogenome, of which 21 were identified to encode homing endonucleases (HEs) and the others had unrecognized functions, whereas in E. rostratum, there were only six free-standing ORFs with unknown functions. A total of 70 introns harbouring 126 intronic ORFs were detected in the core PCGs and rRNAs of the E. turcicum mitogenomes, while in E. rostratum, only 17 introns with 17 intronic ORFs were found within the core PCGs. These intron-encoding ORFs mainly produced LAGLIDADG and GIY-YIG HEs (Supplementary Table 3). In E. turcicum, most of the ORFs in introns encoded LAGLIDADG endonucleases, which were three times as many as GIY-YIG endonucleases. Analogously, the intronic ORFs in E. rostratum encoded more LAGLIDADG endonucleases, two times as many as GIY-YIG endonucleases.

rRNA, tRNA, and codon analyses

Both the above Exserohilum mitogenomes had two rRNA genes, namely the small subunit ribosomal RNA (rns) and large subunit ribosomal RNA (rnl) (Supplementary Table 3). Nine and eight introns were respectively detected in the rnl and rns of the E. turcicum mitogenome, whereas no intron was found in both the two rRNA genes of E. rostratum. Therefore, E. turcicum possessed much longer mitochondrial rnl and rns genes (19,676 bp and 10,448 bp) than E. rostratum (3,408 bp and 1,682 bp).

A total of 32 and 29 tRNA genes were respectively identified in the mitogenomes of E. turcicum and E. rostratum, with lengths ranging from 70 bp (trnC-2) to 85 bp (trnY and trnS) (Fig. 2), which could encode 20 standard amino acids. Five (trnC, trnS, trnL, trnP and trnV) and four (trnL, trnP, trnN and trnS) tRNA genes were noticed to be duplicated in the E. turcicum and E. rostratum mitogenomes, respectively. Four copies of trnR and three copies of trnM were present in both the two mitogenomes. The copy number of trnN was three in the E. turcicum mitogenome. Overall, compared with the mitogenome of E. rostratum, the E. turcicum mitogenome contained three additional tRNA genes, i.e., trnN-3, trnC-2 and trnV-2. No introns had been found in these tRNA genes, and therefore all tRNA genes in the two mitogenomes could be folded into classical cloverleaf secondary structures. Of the 29 tRNAs shared by the two Exserohilum mitogenomes, only two variable sites were found in the trnP-1 and trnN-1, respectively.

Fig. 2

Putative secondary structures of tRNA genes in the two Exserohilum mitogenomes. The tRNA genetic structures drawn in green represent those shared by the two Exserohilum species, while the blue tRNAs are only present in E. turcicum. Residues conserved between the two mitogenomes are shown in green, while the variable sites are shown in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

ATG was the most commonly used start codon in the core PCGs of the two Exserohilum mitogenomes, except in the atp6 gene (TTG as the start codon) of E. turcicum. TAA was usually used as the stop codon in core PCGs of the two mitogenomes, followed by TAG. Codon usage analysis showed that the most frequently used codons in the two Exserohilum mitogenomes were TTA (for Leucine; Leu), AAA (for lysine; Lys), ATA (for isoleucine Ile), TCT (for serine; Ser), and AAT (for asparagine; Asn) (Fig. 3 and Supplementary Table 4). The high frequent use of A and T in codons resulted in the high AT content of the two Exserohilum mitogenomes (average 70.62%).

Fig. 3

Codon usage in the mitochondrial genomes of two Exserohilum species, E. turcicum (A) and E. rostratum (B), respectively. Frequency of codon usage is plotted on the y-axis.

Overlapping nucleotides and intergenic regions

There were six pairs of overlapping nucleotides detected in the mitogenome of E. turcicum, which were located across the neighboring genes orf139 and orf313 (−70 bp), orf70 and orf285 (−22 bp), orf130 and orf96 (−5 bp), orf329 and orf193 (−20 bp), orf555 and orf68 (−43 bp), and nad6 and orf388 (−85 bp), respectively. The E. rostratum mitogenome contained two pairs of overlapping nucleotides, of which one was located across the neighboring nad4 and orf99 genes (62 bp) and the other was located between orf99 and nad4L (1 bp) (Supplementary Table 3). A total of 37,766 bp and 19,147 bp of intergenic sequences were detected in the mitogenomes of E. turcicum and E. rostratum, respectively. Lengths of the intergenic spacers of the two Exserohilum mitogenomes ranged from 0 to 3,720 bp and from 0 to 1,446 bp, respectively. The longest intergenic segments in the E. turcicum and E. rostratum mitogenomes were located between orf168 and orf766 and between cox1 and orf138, respectively.

Analysis of repeat elements

Through BLASTn searches of the two Exserohilum mitogenomes against themselves, a total of 678 and 135 repetitive elements were separately detected in the mitogenomes of E. turcicum and E. rostratum (Supplementary Table 5). Sizes of these repeat regions ranged from 33 bp to 432 bp, with pairwise nucleotide similarities ranging from 81% to 100%. The longest repeat elements in the two mitogenomes were 432 bp and 424 bp, respectively. And they were observed to be located around the first exon of the cox1 gene and in the intergenic region between trnC and orf348, respectively. The overall length of the repetitive sequence in the E. turcicum mitogenome was 42,810 bp, accounting for 16.16% of its whole mitogenome, while the 7,012 bp of repetitive sequence in E. rostratum accounted for 10.83% of its mitogenome. In addition, a total of 50 and 26 tandem repeats were identified in the mitogenomes of E. turcicum and E. rostratum, respectively (Supplementary Table 6). Most of them in the two Exserohilum mitogenomes contained 2 to 4 copies, despite the highest copy number (65) observed in E. rostratum. The longest tandem repeat was found in E. turcicum, encompassing four repeat loci with a length of 93 bp. Tandem repeat sequences accounted for 1.43% and 3.00% of the mitogenomes of E. turcicum and E. rostratum, respectively.

Intron dynamics of the cox1 genes in Pleosporales

A total of 275 introns, most of them belonged to group Ⅰ, were detected in the mitogenomes of 13 Pleosporales species, each of which contained 1 to 70 introns (Supplementary Table 2). The great difference in intron number suggested that acquire/loss events of the introns had frequently appeared in the evolution of Pleosporales. In addition, the cox1 was estimated to be the largest host gene of introns in Pleosporales, containing 33.82% of the total introns, which was thus used for further analysis in the intron dynamic change among all the 13 Pleosporales species.

A total of 22 different Pcls were detected in the cox1 genes of 13 Pleosporales species (Fig. 4). The newly sequenced E. turcicum mtDNA was found to contain the largest number of Pcls (18) within the order Pleosporales, including two novel positions (P821 and P1307), whereas no introns were identified in the cox1 gene of Phaeosphaeria nodorum (Fig. 5). Half of the 22 Pcls were found to be widespread, among which P615, P731 and P1125, existing in nine of the 13 Pleosporales species, were the most common, followed by P1262 detected in seven Pleosporales species. P678 and P807 could be separately detected in only one of the 13 Pleosporales mitogenomes. Besides, the P678 of Stemphylium lycopersici seemed to be a species-specific Pcl that had been never detected in other ascomycetes.

Fig. 4

Insertion sites of different position classes (Pcls) in the coding regions of cox1 genes of 13 Pleosporales species. Protein sequences encoded by the cox1 genes of 12 other Pleosporales species were aligned with the cox1 of Bipolaris cookei (as reference). The Pcls were named according to their insertion sites in the reference cox1 sequence of B. cookei. The Pcls in brown, yellow, green, red, blue, purple, orange and black represent their different numbers (9, 7, 6, 5, 4, 3, 2 and 1, respectively) among the 13 Pleosporales mitogenomes. The symbols ‘+1′ and ‘+2′ refer to the different insertion positions of Pcls within triplet codons: ‘+1′ when between the 1st and 2nd nt of a codon; and ‘+2′ when between the 2nt and 3rd nt of a codon. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5

Position class (Pcl) information of the cox1 genes in 13 Pleosporales species. The Pcls were named according to their insertion sites in the reference cox1 sequence of Bipolaris cookei (MF784482). Phylogenetic positions of the 13 Pleosporales species were established by using both the Bayesian inference (BI) and Maximum Likelihood (ML) methods based on a concatenated mitochondrial gene set. All species are shown in Supplementary Table 2.

Genetic distance and evolutionary rates of 12 core PCGs

As the rps3 gene was often absent from the Pleosporales mitogenomes, here only 12 core PCGs were used to calculate the genetic distances and substitution rates between each pair of the 13 Pleosporales species. It was found that the nad3 gene had the largest K2P genetic distance (average value 0.34), followed by the nad2 (average value 0.18), which meant that they exhibited the fastest mutation rates among the 12 core PCGs in Pleosporales (Fig. 6). The cox1 gene showed the lowest genetic differentiation between the 13 Pleosporales species, with an overall mean K2P distance of 0.10, indicating that this gene was highly conserved. Across the 12 core PCGs detected, the nad3 gene had the highest non-synonymous substitution rate (Ka) (average value 0.19) between the 13 Pleosporales species, followed by nad2 (average value 0.08), while the cox1 gene exhibited the lowest Ka value (average 0.02). The nad3 gene exhibited the highest synonymous substitution (Ks) rate (average value 3.73), while the lowest Ks value (average 0.78) among the 12 PCGs was observed in the cob gene. The overall Ka/Ks values for all detected core PCGs were<1, suggesting that these genes were subject to purifying selection.

Fig. 6

Genetic and evolutionary analyses of 12 core protein-coding genes (PCGs) in the 13 Pleosporales mitogenomes. K2P, the Kimura-2-parameter distance; Ka, nonsynonymous substitution site; Ks, synonymous substitution site. The software DnaSP v6.10.01 was used to calculate Ka and Ks values. The MEGA v6.06 was used to calculate the K2P distances.

Mitochondrial gene arrangement in Pleosporales species

Here, we compared the arrangements of mitochondrial core PCGs and rRNAs genes among the 13 Pleosporales species (Fig. 7), and found that the mitochondrial gene arrangement in Pleosporales often varied greatly, but was fairly conserved within genera, such as Exserohilum and Bipolaris. There were three Pleosporales mitogenomes that had no rps3 genes, including Phaeosphaeria nodorum, Didymella pinodes and Pithomyces chartarum. As mentioned earlier in the Exserohilum mitogenomes, an uninterrupted gene pair formed by the cox1 and cox2 was also found in all other Pleosporales species. The similar situation was observed between nad2 and nad3 and between nad4L and nad5. Besides, the three genes rnl, nad6 and rns together constituted a conserved gene block in most Pleosporales species except Shiraia bambusicola. Among the three closely related species of Bipolaris, no difference in gene order was found between B. sorokiniana and B. cookei, though position exchanges of two pairs of genes (cox1&cox2 and nad6&rns) were observed in the B. oryzae mitogenome. As for the two Exserohilum species, their mitochondrial gene arrangements were completely consistent, which appeared to be exactly the same as those of B. sorokiniana and B. cookei.

Fig. 7

Mitochondrial gene arrangement analysis of the 13 Pleosporales species. A total of 13 core protein-coding genes (PCGs) and two rRNAs were involved here, which were shown in order of their occurrence in the mitochondrial genome, starting from cox1. Phylogenetic relationship of the 13 Pleosporales species were established using both the Bayesian inference (BI) and Maximum Likelihood (ML) methods based on concatenated mitochondrial genes. Species IDs and NCBI accession numbers used for this analysis are listed in Supplementary Table 2.

Mitogenome comparison

Comparative analysis showed that the mitogenome of E. turcicum was the largest one among all the 13 Pleosporales species available online (Supplementary Table 2). Its size (264,948 bp) was approximately 1.92 to 6.79 times larger than the other Pleosporales mitogenomes, indicating that E. turcicum might have experienced a huge mitogenome expansion during evolution. GC contents of the two Exserohilum mitogenomes (average 29.38%) were just slightly higher than those of other Pleosporales species (average 29.13%). Of the 13 Pleosporales species, each contained 15–49 PCGs (including core PCGs and free-standing ORFs) with the mean value of 20, which was close to the number of PCGs (19) in E. rostratum but much fewer than that (49) of E. turcicum. Besides, the E. turcicum mitogenome also had the largest numbers of introns (70) and intronic ORFs (1 2 6) among the 13 Pleosporales species, followed by the B. cookei mitogenome with 40 introns and 51 intronic ORFs. In addition, these 13 Pleosporales species generally contained 21–38 tRNA genes, but always possessed the same number of rRNA genes (2).

Whole-mitogenome collinearity analysis was carried out between five closely related Pleosporales species, i.e., two Exserohilum and three Bipolaris fungi. A total of 11 homologous regions (A to K) were detected in each of the five Pleosporales mitogenomes (Fig. 8). Compared to the E. turcicum mitogenome, most of these homologous regions in the other four species were significantly reduced in size, especially for E. rostratum that seemed to have experienced a massive mitogenome contraction during evolution. Nevertheless, between the two Exserohilum mitogenomes, ten of the homologous regions (only except B) were arranged in the identical order. In addition, it was observed that, in the genus Bipolaris, the relative locations of most homologous regions were basically consistent as well. The homologous regions in B. sorokiniana and B. cookei had almost the same order (only except E), whereas the positions of five homologous regions (D, E, G, I, and J) varied between B. oryzae and the other two Bipolaris mitogenomes. Remarkably, E. turcicum and B. cookei, belonging to two separate genera, had exactly the same order of homologous regions. In summary, the above Exserohilum and Bipolaris species, except B. oryzae, have a high degree of collinearity whether within or between genera.

Fig. 8

Mitogenomic collinearity analysis of five closely related Pleosporales species by using Mauve v2.4.0. Homologous regions between different species were represented by the same color blocks and connected by the same color lines. ET: Exserohilum turcicum, ER: E. rostratum, BC: Bipolaris cookei, BS: B. sorokiniana, BO: B. oryzae.

A phylogenetic tree involving 98 fungal species of Ascomycota was constructed based on both Bayesian inference (BI) and Maximum likelihood (ML) analyses of the 15 concatenated mitochondrial conserved PCG genes (Fig. 9). All major branches within this tree were well supported (BPP ≥ 0.95; BS ≥ 89). Both T. deforman and T. wiesneri from Taphrinomycetes were appointed as the outgroup, while the other 96 species were well assigned into five clades in this topology, including Dothideomycetes, Sordariomycetes, Eurotiomycetes, Lecanoromycetes, and Pneumocystidomycetes (Supplementary Table 7). In the Dothideomycetes, different orders could be also well separated, namely Pleosporales, Cladosporiales, Mycosphaerellales, Dothideales, and Botryosphaeriales. In the order of Pleosporales, as expected, the two Exserohilum species tested here were clustered tightly with the highest support rates (BPP = 1; BS = 100%). And the genus Exserohilum was the most contiguous to the Bipolaris cluster containing three closely related species, B. cookie, B. sorokiniana and B. oryzae.

Fig. 9

Molecular phylogeny of 98 Ascomycota species based on both Bayesian inference (BI) and Maximum Likelihood (ML) analyses of 15 protein coding genes. Bootstrap (BS) values and Bayesian posterior probabilities (BPPs) are respectively placed before and after the slashes. The asterisks indicate that the BPP and BS values are 1 and 100%, respectively. The BI and ML analyses were performed by using MrBayes v3.2.6 and RAxML v8.2.12, respectively. NCBI accession numbers of these species used in the phylogenetic analysis are provided in Supplementary Table 7.

Discussion

Previous studies have indicated that fungal mitochondrial genomes were highly variable in size from 18.84 kb to 531.19 kb [24], [25], [26], and four main factors were reported to bring about this variation, including the introns, intergenic regions, repeat sequences, and plasmid derived dynamic change regions [63], [64]. In this research, complete mtDNAs of two Exserohilum fungi were sequenced, assembled, and analyzed. Despite belonging to the same genus, these two species had vastly different mitogenome sizes, of which E. turcicum was the fourth largest (264,948 bp in length) among all reported mitogenomes of Ascomycetes, only smaller than the three macro-fungi Morchella crassipes (531,195 bp), Tuber calosporum (287,403 bp) and M. importuna (272,238 bp) [26], [65], [66]. Besides, it also became the largest one among all the 13 Pleosporales mitogenomes available online (Supplementary Table 1), which was about twice as big as the previously reported largest mitogenome of B. sorokiniana (137,775 bp) [39]. However, the size of the E. rostratum mitogenome (64,620 bp) was less than one-fourth of that of E. turcicum, even lower than the average value (∼79 kb) of the other 11 Pleosporales species.

The earlier morphological identification, later multigene phylogenetic analyses based on nuclear loci, and current phylogenetic tree based on combined mitochondrial gene set (Fig. 9) all supported that the genera Exserohilum and Bipolaris had quite a close relationship [2], [3]. Nevertheless, these two closely related taxa differed greatly in the mitogenome scale. Unlike Exserohilum, the three Bipolaris species used here, not surprisingly, possessed relatively similar mitogenome sizes, 124,887 bp, 135,790 bp and 137,775 bp, respectively, with a rangeability of only 10% ((max–min)/min). Interestingly, mitogenomes of the two Exserohilum species were either much larger or far smaller than those of the genus Bipolaris, which implied that the two Exserohilum fungi might have undergone massive mitogenome expansion and contraction, respectively, during evolution. We found that, compared with E. rostratum and the three Bipolaris species, the mitogenome of E. turcicum had expanded by 310% (vs. E. rostratum), 112% (vs. B. oryzae), 95% (vs. B. cookei), and 92% (vs. B. sorokiniana), respectively. Of all the compositions of mitogenomes (described in the Results 3.1), intronic regions seemed to contribute the most to this expansion, with the contribution rates of 72.16%, 73.23%, 63.85% and 82.30%, respectively (Supplementary Fig. 1), followed by the protein coding regions that contributed 16.56%, 14.67%, 27.78% and 18.87%, respectively. On the contrary, compared with E. turcicum and the three Bipolaris species, the mitogenome of E. rostratum had shrunk by 76% (vs. E. turcicum), 53% (vs. B. sorokiniana), 52% (vs. B. cookei), and 48% (vs. B. oryzae), respectively. Similarly, the intronic regions seemed to contribute the most to this contraction, with the contribution rates of 72.16%, 54.54%, 87.25% and 69.68%, respectively (Supplementary Fig. 1). All these results indicated that the expansion/contraction of the Exserohilum mitogenomes was primarily due to their increase/decrease of intronic regions. Besides, repetitive sequences might also contribute to the expansion/contraction of the two Exserohilum mitogenomes, though their contribution rates were far below the intronic regions. In addition, plasmids and plasmid-like elements, which had been reported in some fungi as important factors leading to the change of mitogenome size [64], [67], [68], were not detected in the two Exserohilum mitogenomes.

Despite the huge distinction of the two Exserohilum mitogenomes in size, their core PCGs were very conserved in both number and sequence, which were consistent with the gene contents of other Pleosporales mitogenomes, consisting mainly of nad1-nad6, nad4L, cob, cox1-cox3 and atp6. However, there were another two core PCGs, atp8 and atp9, that could be commonly observed in most orders of Dothideomycetes, such as Mycosphaerellales [69], but were missing in all reported mitogenomes of Pleosporales including the genus Exserohilum [38], [39], [70], [71], [72], [73], [74]. In addition, the two Exserohilum mitogenomes had a significant feature, i.e., the gene fusion between cox1 and cox2, which was also found in some other Pleosporales species, such as Corynespora cassiicola and S. lycopersici [73], [74]. Estimation of genetic distances indicated that the sequences of cox1 genes were the most conserved between the 13 Pleosporales mitogenomes, whereas the nad3 genes exhibited considerable genetic differentiation. The Ka and Ks values varied not only between different core genes, but amongst the same genes of different species, especially for nad4 and nad4L, implying that in some Pleosporales species these genes may have evolved at a faster rate. All the 12 core PCGs of the 13 Pleosporales mitogenomes were identified to have been subject to purifying selection (Ka/Ks < 1), which meant that these genes had underwent a relatively conservative evolutionary process in functionality. In addition, the numbers of tRNA genes in the two Exserohilum mitogenomes (average 30.5) were just slightly higher than the average value (28.5) of other Pleosporales species. However, we noticed that, compared with E. turcicum, three tRNA genes (trnN-3, trnC-2 and trnV-2) were lost in the E. rostratum mitogenome; and besides, two variable sites occurred in the tRNAs trnP-1 and trnN-1 between the two Exserohilumum species, of which one was located in the acceptor stem (trnN-1) and the other was in the anticodon (trnP-1). The former may be a meaningless mutation, whereas the latter, if a non-synonymous mutation, can lead to the change of specific recognition of mRNA codon, and consequently affect the correct amino acid incorporation during translation [75], [76]. The tRNA mutations that affect protein synthesis and metabolism have been found in some other eukaryotes [77], [78].

Arrangement of the mitochondrial genes could provide a lot of reference information for understanding the phylogeny and evolution of eukaryotic species [79], [80]. In animals, the mitochondrial gene rearrangement has been widely studied and a number of models associated with it have been proposed, such as the tandem duplication-random loss (TDLR) and duplication and nonrandom loss model [64], [81]. By contrast, the rearrangement of fungal mitochondrial genes became more frequent, and their order and orientation could change in a larger scale [22]. In this research, a remarkable variation in gene arrangement was observed between different members of the Pleosporales, which was consistent with previous reports [72]. Nevertheless, several conserved gene blocks were found, like cox1-cox2, nad2-nad3, nad4L-nad5 and rnl-nad6-rns. Even between some genera, namely Exserohilum and Bipolaris, the gene arrangements were fairly conserved, which accorded with the fact that these two taxa were closely related and used to be classed as the same group. In addition, between B. oryzae and the other two Bipolaris species, obvious position exchanges of two pairs of mitochondrial genes (cox1&cox2 and nad6&rns) were observed, which led to a slight gene rearrangement within the genus. This limited exchange of gene positions between related species had been also found in the mitogenomes of Boletales [82]. Additionally, according to some previous studies [22], [24], accumulation of the repetitive sequences in fungal mitogenomes correlated with the mitochondrial gene rearrangement. However, in this study, although a large number of repeat elements were detected in both the two Exserohilum mitogenomes, there was no gene rearrangement between them, which was consistent with the discovery in Boletus [82]. Overall, the mechanism of mitochondrial gene rearrangement in Pleosporales is quite important for understanding their evolution history, and deserves further investigation.

In eukaryotes, introns with various sizes and frequencies are present in nuclear, mitochondrial, and plastidial genes from different kingdoms [83]. Mitochondrial introns are classified in two main groups, namely group I and II, of which the former was mainly found in fungi and the latter was more common in plants [22]. In this study, a total of 275 introns, most of them belonged to group I, were detected in 13 Pleosporales mitogenomes, and each species had 21 introns on average. As more than a third of these introns were found in the cox1, this gene was thus selected to analyze the dynamic changes of mitochondrial introns among the 13 Pleosporales species. By looking into the 22 Pcls detected in the cox1 genes, we noticed that the 13 Pleosporales species had a great difference in the number and type of mitochondrial introns, which indicated that the gain/loss events of introns might have occurred in the evolution of Pleosporales. P615, P731 and P1125 were found to be the most conserved intron insertion sites in the order Pleosporales, and were therefore considered to be derived from the common ancestor of Pleosporales species. Per contra, each of the Pcls, P678, P807, P821 and P1307, only existed in one of the 13 Pleosporales mitogenomes. And these scarce introns had been detected in some distantly related species, such as Arthrobotrys musiformis [84] and Monilinia fructicola [85], suggesting potential horizontal gene transfer events. Although E. turcicum possessed the second largest number of introns (18) in the cox1 gene, only one less than that of Agaricus bisporus [55], there were no novel Pcls detected in it, whereas in S. lycopersici, a novel Pcl (P678) was probably found, which had been never detected in all other known Ascomycota mitogenomes. For fungal mitogenomes, group I introns were considered to be an important genetic element [82], [83], and their dynamic changes could significantly affect the size and organization of mitogenomes [29], [65].

Dothideomycetes is the largest class of Ascomycota, including a great number of economically important plant pathogens [86]. Apparent characteristics of many closely related species or species complex in the Dothideomycetes are easy to be confused, which makes it difficult to distinguish these taxa accurately only according to the morphology. Therefore, it is very urgent to find reliable molecular markers to solve the above classification problem. Although some researchers had begun to use single-gene and multi-gene phylogenetic analyses to explore the genetic relationships within the class Dothideomycetes [87], [88], it is still often a hard task to identify some complicated genera and species. In recent years, due to their advantages over nuclear genomes including the uniparental inheritance and accelerated evolutionary rate [89], mitogenomes have been widely used in phylogeny and population study of plants, animals, and some kinds of fungi [22], [24], [64], [82]. Here, we established a highly supported phylogenetic tree of 98 Ascomycota species based on BI and ML analyses of the combined mitochondrial gene set. All these selected species were well gathered or divided into different independent clades and subclades, each of which exactly corresponded to one fungal taxa (such as class and order). The two Exserohilum species were observed to be the most closely related to three Bipolaris species, which was consistent with the previous analyses of multigene phylogeny based on nuclear loci [3]. Our result indicated that mitochondrial genomes might be useful markers (such as the cox1, cox2, and cob genes) for fungal phylogenetic analysis in different taxon level. These mitochondrial barcoding genes have been also used to study biological evolution of animals [90], [91]. However, the number of known mitogenomes in Ascomycota is still limited, which restricts our more accurate reconstruction of the phylogenetic relationships between related species. A better understanding of mitogenome evolution and phylogeny within Ascomycota would require the sequencing of more members from this group.

Conclusion

In this study, the mitochondrial genomes of two Exserohilum species, E. turcicum and E. rostratum, collected from China were sequenced, assembled, annotated and compared. The mitogenome of E. turcicum was the largest (264,948 bp) among all reported Pleosporales mitogenomes, whereas the E. rostratum mitogenome was extraordinarily small (64,620 bp), actually less than a quarter of the former. Through the comparative analysis of genetic compositions, increase/decrease of the intronic regions was considered as the main contributing factor of their mitogenome expansion/contraction. Despite the huge difference in size between the two Exserohilum mitogenomes, we confirmed that they shared the same set of 13 core PCGs and two rRNAs in identical order and orientation. All the 12 core PCGs (excluding rps3) of 13 Pleosporales mitogenomes were identified to have been subject to purifying selection (Ka/Ks < 1). In addition, a total of 93 group I introns, classified into 22 Pcls, were differentially distributed in the cox1 genes of 13 Pleosporales species, of which E. turcicum possessed the second largest number (18). Loss/gain and horizontal transfer events of these introns might promote size variation of the Pleosporales mitogenomes. Mitogenome-wide phylogeny well exhibited the genetic relationships between selected ascomycetes including the two Exserohilum species, suggesting that mitochondrial genes were reliable molecular markers for phylogenetic analysis of Ascomycota. This study served as the first investigation on Exserohilum mitogenomes, which will lay the foundation for further understanding the genetic evolution, species differentiation, and ecological adaptation of Exserohilum species.

Funding

This study was supported by the (202300410200), China and Program for Innovative Research Team (in Science and Technology) in University of Henan Province (18IRTSTHN021), China.

CRediT authorship contribution statement

Qingzhou Ma: Data curation, Formal analysis, Visualization, Supervision, Writing – original draft. Yuehua Geng: Data curation, Formal analysis, Project administration, Visualization. Qiang Li: Methodology, Software. Chongyang Cheng: Formal analysis, Software. Rui Zang: Data curation, Formal analysis, Project administration. Yashuang Guo: Data curation, Formal analysis, Project administration. Haiyan Wu: Data curation, Methodology. Chao Xu: Data curation, Funding acquisition, Supervision, Writing – review & editing. Meng Zhang: Funding acquisition, Project administration, Supervision.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.