The Mitochondrial Genome of a Plant Fungal Pathogen Pseudocercospora fijiensis (Mycosphaerellaceae), Comparative Analysis and Diversification Times of the Sigatoka Disease Complex Using Fossil Calibrated Phylogenies

Mycosphaerellaceae is a highly diverse fungal family containing a variety of pathogens affecting many economically important crops. Mitochondria play a crucial role in fungal metabolism and in the study of fungal evolution. This study aims to: (i) describe the mitochondrial genome of Pseudocercospora fijiensis, and (ii) compare it with closely related species (Sphaerulina musiva, S. populicola, P. musae and P. eumusae) available online, paying particular attention to the Sigatoka disease’s complex causal agents. The mitochondrial genome of P. fijiensis is a circular molecule of 74,089 bp containing typical genes coding for the 14 proteins related to oxidative phosphorylation, 2 rRNA genes and a set of 38 tRNAs. P. fijiensis mitogenome has two truncated cox1 copies, and bicistronic transcription of nad2-nad3 and atp6-atp8 confirmed experimentally. Comparative analysis revealed high variability in size and gene order among selected Mycosphaerellaceae mitogenomes likely to be due to rearrangements caused by mobile intron invasion. Using fossil calibrated Bayesian phylogenies, we found later diversification times for Mycosphaerellaceae (66.6 MYA) and the Sigatoka disease complex causal agents, compared to previous strict molecular clock studies. An early divergent Pseudocercospora fijiensis split from the sister species P. musae + P. eumusae 13.31 MYA while their sister group, the sister species P. eumusae and P. musae, split from their shared common ancestor in the late Miocene 8.22 MYA. This newly dated phylogeny suggests that species belonging to the Sigatoka disease complex originated after wild relatives of domesticated bananas (section Eumusae; 27.9 MYA). During this time frame, mitochondrial genomes expanded significantly, possibly due to invasions of introns into different electron transport chain genes.


Introduction
Mycosphaerellaceae is a highly diverse fungal family containing endophytes, saprobes, epiphytes, fungicolous and phytopathogenic species in more than 56 genera [1,2]. Family members can cause significant economic losses to a large number of important plants including ornamentals, food crops and commercially propagated trees [3][4][5][6][7][8]. Three Mycosphaerellaceae members, Pseudocercospora eumusae, P. fijiensis, and P. musae, [1] are major pathogens of bananas and plantains. They comprise the so-called Sigatoka disease complex which is responsible for one of the most economically destructive diseases for banana shaker. Then, mycelia were harvested from the liquid using vacuum filtration. Total DNA was extracted using a previously described Cetyl Trimethylammonium Bromide (CTAB) method [33]. DNA quality and quantity were measured using a fluorometer (Qubit 3.0, Thermo Fisher Scientific, Waltham, MA, USA). Furthermore, genomic DNA was visualized on 1% agarose gel to check for any break/smear or multiple bands. Library construction was performed using Illumina platform with TruSeq DNA kit (Illumina, San Diego, CA, USA) to acquire as paired-end 2 × 150-bps, with about a 350-bp insert size. Next-generation sequencing was performed by an external service (North Carolina University, Chapel Hill, NC, USA) Hiseq 2500 system ® .
Even though the Sphaerulina musiva mitogenome was available online [8] we reassembled it using raw reads available at NCBI (SRA: SRR3927043). Our major motivation was a 9322 bp inversion detected around the 10,000 bp position of this mitogenome. This inversion was splitting the gene nad2 and we wanted to make sure this inversion was present in the S. musiva mitogenome. First, raw reads were filtered using BBtools (https://sourceforge.net/projects/bbmap/ (accessed on 10 February 2021)) in Geneious 9.1.5 [43]. Then, MITObim version 1.8 [47] used cox1 as bait to map all filtered reads that partly or fully overlap with the bait. Eventually, this leads to an extension of the reference sequence and a reduction of gaps until completion of the whole mitogenome [47]. This inversion in the Sphaerulina musiva mitogenome was found to be an artifact after reassembling raw reads.

Annotation
Mitochondrial genomes of Pseudocercospora fijiensis, P. eumusae, P. musae, Sphaerulina populicola, S. musiva and Pseudovirgaria hyperparasitica were annotated in this study using a combination of software. First, predicted ORFs were determined with a translation code for "mold mitochondrial genomes" using Geneious 9.1.5 [43]. Second, genes were identified using BLASTP version 2.4.0 [48] against the non-redundant protein database from NCBI (downloaded August and December 2016); genes were also identified using MITOS [49]. Third, protein domains and sequence patterns were searched with PFAM [50] and PANTHER 11.0 [51]. Additionally, mitogenome annotation was performed using multiple alignments among the fourteen CMPCGs using MUSCLE version 3.8.31 [52] and CLUSTAL W version 2.0 [53]. Inconsistencies regarding length and position of genes was solved paying particular attention to start and stop codons. Identified ORFs larger than 300 bp with start and stop codons that did not show results with the above-mentioned annotation strategies were considered as hypothetical proteins. Circular mitogenome maps were constructed using Geneious 9.1.5. and Geneious prime [43].

PCR Amplification of cox1 Gene Copies in P. fijiensis
A PCR assay was performed to confirm the presence of two different cox1 copies: a truncated copy (cox1_1) and a complete cox1 copy with an intron (cox1_2). Primers were designed to amplify regions between the first copy (cox1_1) and the second (cox1_2), including the exons of this last copy. First, a set of primers encompassed cox1_2 exon1 and cox1_2 exon 2. A second pair of primers encompassed cox1_1 and cox1_2 exon 2. PCR amplifications were carried out in a total volume of 10 µL, containing 20 ng genomic DNA, 0.15 µM of each primer, 1× PCR buffer (without MgCl 2 ), 0.75 mM MgCl 2 , 4 µM of each dNTP and 0.65 U recombinant Taq DNA polymerase (Thermo Fisher Scientific, Waltham, Massachusetts, USA). Cycling parameters were: 3 min at 94 • C, followed by 35 cycles of 30 s at 94 • C, 30 s at a 50 to 60 • C temperature gradient to determine annealing temperature, 1 min at 72 • C, and a final elongation step of 5 min at 72 • C. PCR products were separated by electrophoresis in a 1% (w/v) agarose gel and visualized with GelRed ® (Biotium, Fremont, CA, USA) under UV light.

RT-PCR Assays for Mitochondrial Gene Pairs of P. fijiensis
Total RNA was extracted from P. fijiensis (isolate: 081022) mycelium after fifteen days of culture using TRIzol ® (Life Technologies, Carlsbad, CA, USA) according to the manufacturer's instructions. RNA concentrations were measured at 260 nm using a NanoDrop ND-1000 UV-Vis Spectrophotometer (NanoDrop Technologies, Thermo Fisher). DNAse I (Thermo Fisher Scientific, Waltham, MA, USA) was used for cDNA synthesis from RNA as template for amplification using the Maxima First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Waltham, MA, USA) according to manufacturer's instructions. Primers were designed such that the amplified product encompassed the end of one gene and the beginning of another. We used pairs of NADH-Ubiquinone Oxidoreductase Chain 3 and 2 (nad3-nad2) and mitochondrial encoded ATP Synthase Membrane Subunits 6 and 8 (atp8-atp6). Both genes and their intergenic sequences were partially amplified. PCR products were run in a 1% agarose gel electrophoresis purified using the GFX PCR DNA and gel band purification kit, according to manufacturer's instructions (GE Healthcare, Chicago, IL, USA). Purified PCR products were sequenced using Sanger Technology at Macrogen Inc. (Seoul, Korea). All sequences are available in GeneBank: atp6-atp8 cDNA (GenBank: MN171334); atp6-atp8 DNA (GenBank: MN171335); nad2-nad3 DNA (GenBank: MN171336); nad2-nad3 cDNA (GenBank: MN171337); cob DNA (GenBank: MN171338); cob cDNA (GenBank: MN171339); nad5 cDNA (GenBank: MN171340).

Phylogenetic Analysis and Divergence Times Estimates
Until now, only nuclear markers and strict clock calibration have been used to calculate diversification times for the Sigatoka disease complex species. We aimed to compare these analyses with fossil calibrated Bayesian phylogenies and mitochondrial markers. A phylogenetic tree was reconstructed to calculate diversification times. Since most mitogenomes are uniparentally inherited we used our mitochondrial phylogeny to compare topologies with already published nuclear ones for the Sigatoka disease complex species. Didymella pinodes [56], Pseudovirgaria hyperparasítica [57], Phaeosphaeria nodorum [45] and Shiraia bambusicola [46] were used as outgroups. Ten core mitochondrial genes (cox1, cox3, atp6, cob, nad1, nad2, nad4, nad4L, nad5, nad6) were aligned one by one for all species using CLUSTAL W version 2.0 [53]. Then, all aligned genes were concatenated in a single alignment for phylogenetic reconstruction. Cox2, atp8, atp9 and nad3 were excluded from the alignment either because they could not be fully recovered in P. musae (atp9, nad3) or because they were missing in outgroups S. bambusicola (atp8 atp9), P. nodorum (atp8, atp9) and D. pinodes (atp8, atp9, cox2).
A Generalized time-reversible (GTR) model was used with an estimated gamma parameter of rate heterogeneity to build maximum likelihood (ML) trees using the Randomized Accelerated Maximum Likelihood RAxML version 8.0 [58] and PhyML version 3.0 [59] programs. One hundred bootstrapped trees were generated and used to assign bootstrap support values to the consensus trees. A Bayesian phylogeny and divergence time analysis was carried out using BEAST2 version 2.5.1 [60]. Separate partitions for each gene were created with BEAUti2 (available in BEAST2). More suitable substitution models for each gene were found using the software package jModelTest2 version 2 [61] according to the Bayesian Information Criterion (BIC) [62]. To accommodate for rate heterogeneity across the branches of the tree we used an uncorrelated relaxed clock model [63] with a lognormal distribution of rates for each gene estimated during the analyses. We also used a strict clock for further comparison of results.
The fossil Metacapnodiaceae [64] was used, assuming this to be a common ancestor of the order Capnodiales with a minimum age of 100 MYA (gamma distribution, offset 100, mean 180, maximum softbound 400). Capnodiales nodes were constrained to monophyly based on the results obtained from ML analysis. A birth/death tree prior was used to model the speciation of nodes in the topology, with gamma priors on the probability of splits and extinctions. We used vague priors on the substitution rates for each gene (gamma distribution with mean 0.2 in units of substitutions per site per time unit). All XML files used to build our Bayesian phylogenies are available at Figshare (https://doi.org/10.608 4/m9.figshare.12101055.v1 (accessed on 10 February 2021)). To ensure convergence we ran analyses five times for 50 million generations each, sampling parameters every 5000 generations, assessing convergence and sufficient chain mixing (Effective Sample Sizes > 200) using Tracer version 1.5 [65]. After removal of 20% of each run as burn-in, the

P. fijiensis Mitochondrial Genome
The whole-genome assembly of P. fijiensis had an N50 = 39,827 bp and a size of 70.53 Mb in 8040 scaffolds. Recovered scaffolds had an estimated genome coverage of 50X. To separate scaffolds belonging to P. fijiensis mitogenome, whole-genome scaffolds were mapped to a draft and unpublished P. fijiensis mitogenome (see Table S1). Scaffolds belonging to the mitogenome of P. fijiensis were recovered and assembled in a circular sequence deposited in GenBank (accession number: MK754071).
yses five times for 50 million generations each, sampling parameters every 5000 generations, assessing convergence and sufficient chain mixing (Effective Sample Sizes > 200) using Tracer version 1.5 [65]. After removal of 20% of each run as burn-in, the remaining trees were combined using LogCombiner (available in BEAST2), summarized as maximum clade credibility (MCC) trees in TreeAnnotator (available in BEAST2), and visualized using FigTree version 1.3.1 [66].

P. fijiensis Mitochondrial Genome
The whole-genome assembly of P. fijiensis had an N50 = 39,827 bp and a size of 70.53 Mb in 8040 scaffolds. Recovered scaffolds had an estimated genome coverage of 50X. To separate scaffolds belonging to P. fijiensis mitogenome, whole-genome scaffolds were mapped to a draft and unpublished P. fijiensis mitogenome (see Table S1). Scaffolds belonging to the mitogenome of P. fijiensis were recovered and assembled in a circular sequence deposited in GenBank (accession number: MK754071).
Twelve putative ORFs of unknown function were predicted to produce hypothetical proteins containing 77 to 583 amino acids, which were described as ORF1-ORF11 and ORF Cytb-like. CMPCGs covered 65.55% of the genome (including 11.77% putative ORFs) ( Figure 1A, Table 1). tRNA genes and both rnl and rns corresponded to 8.8% and 10.67% respectively ( Figure 1B, Table 1). These values were similar to those reported for other representatives of the Ascomycota phylum (Table 2). Overall, P. fijiensis mtDNA has a nucleotide composition of: 36.6% of A, 13.4% of C, 12.6% of G and 37.7% of T. GC-content was 26% with coding and non-coding parts of the genome having, on average, the same GC-percentage.  The three most frequent codons were TTA (708 counts), ATA (498 counts) and TTT (498 counts) encoded amino acids leucine, isoleucine, and phenylalanine, respectively. These amino acids have hydrophobic side chains commonly found in transmembrane helices. These three codons accounted for 6.9% of all codons in the mitogenome. Eight codons (ATA, ATT, TTA, and TTG encoding methionine; CGA, CGC, CGG encoding arginine and TAG encoding a stop codon) were under-represented, being used from one to five times each. CMPCGs started with ATA, ATG, ATT or TTA encoding methionine translation initiation codon. The preferred stop codon was TAA, present in 21 protein-coding genes; the alternative stop codon was TAG. Codon usage of the ORFs was similar to that of the protein-coding loci. Atp9 also contained the highest GC contents (40.9%) among all the CMPCGs (Table 1).
Twenty-nine repeat regions ranging from 611 bp to 128 bp were located in non-coding regions of the mtDNA, but only one was present three times (length 193 bp) ( Table 2). A search within the P. fijiensis mitogenome detected 33 SSR markers. The most common SSR were mononucleotide (30 SSR), and only three dinucleotide SSR were found (https: //doi.org/10.6084/m9.figshare.12101058 (accessed on 10 February 2021)). A difference of seven nucleotides between P. fijiensis isolate 081022 and the unpublished mitochondrial genome available in MycoCosm was found, being 99% identical (Table S1).

Presence of Truncated Conserved Mitochondrial Protein-Coding Genes (CMPCGs)
P. fijiensis mitogenome had truncated copies of four CMPCGs. They were considered truncated copies because they have PFAM and PANTHER domains and high local sequence similarity (>90%) with a complete CMPCG copy despite being an incomplete sequence without start and/or stop codons. Truncated copies included atp6 (two truncated copies and one complete gene sequence), and cox1, cob, nad2 and atp9 (one truncated copy and one complete gene sequence). PCR amplifications were performed to confirm the presence of CMPCG copies in the P. fijiensis mitochondrial genome (Figure 2). The annotation showed that cox 1 had one complete copy with an intron and a second truncated copy (see Figure 2A). Primers were designed to amplify three intra and intergenic regions in these two copies (Figure 2A-D). The amplified PCR bands were of the expected size showing that the assembly on this region was correct ( Figure 2B). The presence of an intron in the complete cytochrome b (cob) gene of P. fijiensis was also experimentally confirmed by PCR amplification ( Figure S1).

Homing Endonucleases and Introns in the P. fijiensis Mitogenome
A total of two different mobile introns were annotated in the mitogenome of P. fijiensis ( Figure 1A). All encoding introns were characterized as group I intron type, which encodes homing endonucleases (HE) [67]. They belonged to the LAGLIDADG family and were 680 bp for an intronic-ORF of 2968 bp length in nad5 and 743 bp long for an intronic-ORF of 1056 bp length in cob (Figure 1, Tables 2 and 3). 2). The annotation showed that cox 1 had one complete copy with an intron and a second truncated copy (see Figure 2A). Primers were designed to amplify three intra and intergenic regions in these two copies (Figure 2A-D). The amplified PCR bands were of the expected size showing that the assembly on this region was correct ( Figure 2B). The presence of an intron in the complete cytochrome b (cob) gene of P. fijiensis was also experimentally confirmed by PCR amplification ( Figure S1).

Figure 2.
Duplicated Protein-coding genes in the mitogenome of Pseudocercospora fijiensis. A. A truncated cox1 copy and a complete copy with an intron that are localized in tandem. B. PCR amplification including intra and intergenic regions between truncated cox1 copies. C. Pairwise alignments of truncated cox1 copies reveal a region of around 600 bp of low similarity between them (non-identical sites appear in black in the alignment). D. Primers used for the amplification of both cox1 copies.

Homing Endonucleases and Introns in the P. fijiensis Mitogenome
A total of two different mobile introns were annotated in the mitogenome of P. fijiensis ( Figure 1A). All encoding introns were characterized as group I intron type, which encodes homing endonucleases (HE) [67]. They belonged to the LAGLIDADG family and were 680 bp for an intronic-ORF of 2968 bp length in nad5 and 743 bp long for an intronic-ORF of 1056 bp length in cob (Figure 1, Tables 2 and 3).

Figure 4. Comparison of Conserved Mitochondrial Protein
Coding Genes (CMPCGs) order and orientation among selected Mycosphaerellaceae species. Gene order is not conserved across members even for closely related species except for Sphaerulina species. Colored gene pairs were always recovered as neighbors; asterisks indicate atp8 is neighbor to cox1 in P. eumusae. Ellipsis in P. musae and P. eumusae mean these genomes were both recovered in different contigs each (two and three respectively); where nad6 and atp9 could not be assembled to any contig we not complete a circular mitogenome.
In addition to these core genes, five out of seven mitogenomes also have hypothetical protein coding genes, accessory genes, additional copies of CMPCGs and genetic mobile elements (Figure 3, Table S2). Genome sizes, gene numbers and contents are heterogeneous among species (Table 2, Figure 3). All annotated genomes were found to have truncated gene duplications of some CMPCGs (Table S2). Alignments of truncated gene copies showed their sequences were not identical. However, one of the copies always had a complete coding sequence with a start and a stop codon (https://doi.org/10.6084/m9.figshare.13542191). Mitochondrial genomes from P. musae and P. eumusae also had RNA and DNA polymerases (Table S2).
The order of CMPCGs among selected Mycosphaerellaceae mitogenomes was variable, except for sister species Sphaerulina populicola and S. musiva (Figure 4). Despite this variability, gene pair order was always conserved for gene pairs nad4L-nad5, nad3-nad2 and atp8-atp6 among all mitogenomes (Figure 4). We mapped RNAseq assembled transcripts from transcriptomes available online for S. musiva (SRR1652271), P. fijiensis (SRR3593877, SRR3593879), P. eumusae (GDIK00000000.1) and P. musae (GDIN00000000.1) to P. fijiensis, S. musiva, P. musae and, P. eumusae mitochondrial genomes and found neighbor genes were always part of the same transcript for each species. This suggested that they were transcribed as a single mRNA. RT-PCR amplification and subsequent Sanger sequencing confirmed bicistronic expression for nad3-nad2, atp8-atp6 gene pairs in P. fijiensis ( Figure 5). Gene order is not conserved across members even for closely related species except for Sphaerulina species. Colored gene pairs were always recovered as neighbors; asterisks indicate atp8 is neighbor to cox1 in P. eumusae. Ellipsis in P. musae and P. eumusae mean these genomes were both recovered in different contigs each (two and three respectively); where nad6 and atp9 could not be assembled to any contig we not complete a circular mitogenome.
The order of CMPCGs among selected Mycosphaerellaceae mitogenomes was variable, except for sister species Sphaerulina populicola and S. musiva (Figure 4). Despite this variability, gene pair order was always conserved for gene pairs nad4L-nad5, nad3-nad2 and atp8-atp6 among all mitogenomes (Figure 4). We mapped RNAseq assembled transcripts from transcriptomes available online for S. musiva (SRR1652271), P. fijiensis (SRR3593877, SRR3593879), P. eumusae (GDIK00000000.1) and P. musae (GDIN00000000.1) to P. fijiensis, S. musiva, P. musae and, P. eumusae mitochondrial genomes and found neighbor genes were always part of the same transcript for each species. This suggested that they were transcribed as a single mRNA. RT-PCR amplification and subsequent Sanger sequencing confirmed bicistronic expression for nad3-nad2, atp8-atp6 gene pairs in P. fijiensis ( Figure 5). Major differences in genome size among members of Mycosphaerellaceae seem be related to the invasive presence of HEG and ORFs in some mitogenomes. Zasmi cellare and Zymoseptoria tritici have compact mitogenomes with core mitochondrial g (CMPCGs) and lack of HEG. While, in other Mycosphaerellaceae, the presence o quences containing LAGLIDADG or GIY-YIG domains related to Homing Endonuc Genes (HEGs) was observed ( Table 2). These ORFs ranged from one in Pseudocerco musae and P. eumusae to 28 in Sphaerulina populicola (Table 2). S. populicola has the la mitogenome report in this study (139606bp), suggesting HEGs might have caused mentation of CMPCGs. In almost all instances, pieces of CMPCGs were collinearly tributed and each fragment was found to be followed by an insertion of a HEG re sequence. An extreme case of HEG invasion to CMPCGs was found in cox1 of Sphaer populicola (CDS: 1542bp). This gene has twelve fragments distributed along 22070bp, of them containing a piece of cox1 followed by a HEG related domain ( Figure S3).

Inferred Mitochondrial Phylogeny and Diversification Times
A robust phylogeny with posterior probabilities greater than 0.97 was recov containing two main lineages (Pleosporales + Capnodiales). Divergence time estim using fossil calibration are shown in Figure 6 Major differences in genome size among members of Mycosphaerellaceae seemed to be related to the invasive presence of HEG and ORFs in some mitogenomes. Zasmidium cellare and Zymoseptoria tritici have compact mitogenomes with core mitochondrial genes (CMPCGs) and lack of HEG. While, in other Mycosphaerellaceae, the presence of sequences containing LAGLIDADG or GIY-YIG domains related to Homing Endonuclease Genes (HEGs) was observed ( Table 2). These ORFs ranged from one in Pseudocercospora musae and P. eumusae to 28 in Sphaerulina populicola (Table 2). S. populicola has the largest mitogenome report in this study (139,606 bp), suggesting HEGs might have caused fragmentation of CMPCGs. In almost all instances, pieces of CMPCGs were collinearly distributed and each fragment was found to be followed by an insertion of a HEG related sequence. An extreme case of HEG invasion to CMPCGs was found in cox1 of Sphaerulina populicola (CDS: 1542 bp). This gene has twelve fragments distributed along 22,070 bp, each of them containing a piece of cox1 followed by a HEG related domain ( Figure S3).

Discussion
In this study, the complete mitochondrial genome of a plant pathogenic fungus, Pseudocercospora fijiensis, was sequenced and annotated. We also used comparative analysis and fossil calibration phylogenies to further understand the evolution of Mycosphaerellaceae mitogenomes. To date, more than 700 complete fungal mitogenomes are available online, but the mitogenomes of Pseudocercospora species have not been reported in the organelle genome database of NCBI (August 2017). The mitogenome of P. fijiensis and related species provides a molecular basis for further studies on molecular systematics and evolutionary dynamics of Ascomycota fungi especially belonging to Dothideomycetes.

Discussion
In this study, the complete mitochondrial genome of a plant pathogenic fungus, Pseudocercospora fijiensis, was sequenced and annotated. We also used comparative analysis and fossil calibration phylogenies to further understand the evolution of Mycosphaerellaceae mitogenomes. To date, more than 700 complete fungal mitogenomes are available online, but the mitogenomes of Pseudocercospora species have not been reported in the organelle genome database of NCBI (August 2017). The mitogenome of P. fijiensis and related species provides a molecular basis for further studies on molecular systematics and evolutionary dynamics of Ascomycota fungi especially belonging to Dothideomycetes.
Ascomycetes mitochondrial genomes like most mitochondrial genomes along the tree of life generally consist of: two ribosomal subunits (rnl and rns), a distinct set of tRNAs and fourteen genes of the respiratory chain complexes (cox1, cox2, cox3, cob, nad1 to nad6, atp6, atp8 and atp9) [17]. The mtDNA of P. fijiensis contains the 14 mitochondrial inner membrane proteins involved in electron transport and coupled oxidative phosphorylation, as well as rnl and rns (Figure 1). These genes were also found in mitochondrial genomes of selected Mycosphaerellaceae species studied here. Additionally, DNA polymerases, RNA polymerases and Reverse transcriptases were found in Pseudocercospora musae, P. eumusae and S. populicola (Table S1). These polymerases and transcriptases might come from mitochondrial plasmids integrated into their mitochondrial genomes [21,68,69].
A variable number of open reading frames of unknown function and introns related to homing endonuclease genes (HEG), often including GIY-YIG or LAGLIDADG protein domains, were found in several Mycosphaerellaceae genomes including P. fijiensis (Table 2; Figure 1). Despite not being ubiquitous in Ascomycete mitogenomes these mobile elements are fairly common [16,17,22,67,70]. Differences in mitogenome sizes, gene order and gene duplication among Mycosphaerellaceae are attributed to heterogeneous content of accessory genes, and intron mobile sequences (Table 1, Figures 3 and 4).
The comparative study of Mycosphaerellaceae species selected mitogenomes showed that sizes and gene order are not conserved among members of the family (Figures 3  and 4). Divergence times within sister clades Pseudocercospora 13.31 MYA (9.49-17.28 MYA, 95% HPD) and Sphaerulina 13.39 MYA (9.39-17.69 MYA, 95% HPD) are roughly the same (Table 3; Figure 6). But gene order in Sphaerulina species is conserved while in Pseudocercospora species it is not (Figures 4 and 6). Nonetheless, some gene pairs were always found together, atp6-atp8, nad2-nad3, nad4L-nad5, in most studied species (Figure 4). Gene order variation among Mycosphaerellaceae could be due to mtDNA rearrangements caused by different processes, such as fusion, fission, recombination, plasmid integration and mobility [21]. Aguileta et al. (2014) compared 38 fungal mitogenomes to understand mitochondrial gene order evolution. They found evidence of gene rearrangements and a relationship with intronic ORFs and repeats. Their results support recombination in all fungal phyla. Despite rearrangements being pervasive in fungal mitogenomes, they found conserved gene pairs nad2-nad3 and nad4L-nad5 in most species [20]. Bicistronic transcription of atp6-atp8, nad2-nad3 gene pairs was confirmed experimentally in P. fijiensis using RT-PCR ( Figure 5). Maintenance of such proximity through evolution is important for mitochondria functions and modifications of such proximity can negatively affect organisms [71].
Even though sister species Sphaerulina musiva and S. populicola shared core gene content and gene order, their genomes sizes were quite different (53,234 bp and 139,606 bp). This could be due to a widespread occurrence of homing endonucleases HEG related to ORFs and accessory genes in the S. populicola mitogenome (Table S1, Figure S2). Mitogenome size variability due to the occurrence of mobile elements and accessory gene invasions was also observed in nine phylogenetically related species belonging to the genera Aspergillus and Penicilium [72]. In a phytopathogenic fungus Sclerotinia borealis mitogenome size expansion was also shown to be due to plasmid-like sequences and HEGs related to ORFs [73].
Mitochondrial gene duplication is seldom described in fungi. HEGs invasion has been previously described in fungal mtDNA, where truncated genes were ubiquitous [74]. P. fijiensis and other Mycosphaerellaceae mitogenomes have 2 to12 truncated copies of Conserved Mitocondrial Protein Coding Genes (CMPCGs) (Table S1). Truncated CMPCGs have partial sequences lacking start and/or stop codons. Truncated gene copies did not have a particular distribution pattern; they were either close to each other or dispersed in the genomes. Mardanov et al. (2014) found duplications of truncated extra copies of atp9 and atp6 in the phytopathogenic fungus Sclerotinia borealis [73]. Two incomplete copies of atp6 were also found on different strands of the mtDNA of Shiraia bambusicola (Pleosporales) [46]. However, it is not common to find highly fragmented genes such as the cox1 of Sphaerulina populicola ( Figure S3). A similar case was reported in the mitogenome of Sclerotinia borealis, where thirteen introns of cox1 and truncated copies of CMPCGs were found [73].
Fossil calibrated phylogenies for the Mycospharellaceae had later diversification times 66.66 MYA (55.47-78.27 MYA, 95% HPD) compared to previous studies 186.7-143.6 MYA [10]. The Sigatoka disease complex had an early divergent Pseudocercospora fijiensis, that splits from sister species P. musae + P. eumusae 13.31 MYA (9.49-17.28 MYA, 95% HPD); while sister species P. eumusae and P. musae split from their shared ancestor in the late Miocene 8.22 MYA (5.6-11.07 MYA, 95% HPD) ( Figure 6). Chang et al. (2016) estimated the divergence of P. fijiensis from their last common ancestor with P. musae + P. eumusae to be between 39.9-30.6 MYA and the divergence of P. musae and P. eumusae to be between 22.6-17.4 MYA. Diversification ages estimated here were based on mitochondrial markers and Bayesian analysis using both relaxed and strict clock models. These have placed all diversification times in the Mycosphaerellaceae at later times than those calculated by  (Table 3).
These differences in diversification times compared to Chang et al. (2016) could be due to different calibration points and dating methods. Chang et al. (2016) implemented r8s Likelihood methods [75] and a calibration at the origin of the Dothideomycetes crown group (394-285 MYA) using previous Bayesian estimations [76]. For this study Bayesian analysis in BEAST2 [60] using fossil calibration with a Metacapnodiaceae fossil was implemented [64]. We found that Sigatoka disease members (13.31 MYA (9.49-17.28 MYA, 95% HPD)) appeared after the genus Musa (27.9 MYA (21.5-34.4, 95% HPD)) [77]. Species member of this genus within the section Eumusa sensu latto comprise cultivated banana and are the host of the Sigatoka disease complex. This naturally prompts a further question: did the Sigatoka disease complex originate through host-tracking evolution? This hypothesis explains that a pathogen is likely to be younger than the host until changes related to the genetics of the pathogens or/and exogenous factors have observed alterations in their virulence spectra [70]. We are currently limited in terms of a good taxonomic sampling and biogeographical analysis for Pseudocercospora species in arriving at an answer to this question. A host tracking coevolution hypothesis has also been proposed for Zymoseptoria tritici (syn. Mycosphaerella graminicola) [5].

Conclusions
We successfully sequenced and analyzed mitochondrial genomes of Pseudocercospora fijiensis, P. eumusae, P. musae, Sphaerulina populicola and S. musiva. A robust mitochondrial phylogeny containing two main lineages (Pleosporales + Capnodiales) was obtained and divergence times were estimated using fossil calibration. Fossil calibrated phylogenies are reported for the first time here for fungal plant pathogens that had later diversification times for the origin of all the species involved, compared to previous studies. Genome size variation and organization among Mycosphaerellaceae could be related to the proliferation of type I and II introns, gene duplications and possible plasmid insertions, phenomena known for many fungal mitogenomes. Despite their order variability, some genes were always recovered as neighbors in all mitogenomes analyzed. Bicistronic expression for nad3-nad2, atp8-atp6 gene pairs in P. fijiensis was confirmed experimentally. Further gene editing and virulence assays will be important to shed light on fungal adaptation and more effective disease control strategies. Phylogenomic studies including a good taxonomic sampling and biogeographical analysis for Pseudocercospora species will further clarify whether the Sigatoka disease-causing species virulence flared-up after Musa domestication.
Supplementary Materials: The following are available online at https://www.mdpi.com/2075-1 729/11/3/215/s1: Table S1. Sources of genomic data used in this study. Table S2: Protein-coding genes annotated in Mycosphaerellaceae mitochondrial genomes, Figure S1: Mitochondrial gene duplication in Pseudocercospora fijiensis. A. The intron of cob. B. PCR amplification of intra and intergenic sequences among exons. C. Primers used for the amplification of cob, Figure S2 Funding: This research was funded by Instituto para el desarrollo de la Ciencia y la Tecnología "Francisco José de Caldas (Colciencias)", Colombia, grant nunmber: 221356934854; The Asociación de