Comparative genomics in plant fungal pathogens (Mycosphaerellaceae): variation in mitochondrial composition due to at least five independent intron invasions

Fungi provide new opportunities to study highly differentiated mitochondrial DNA. Mycosphaerellaceae is a highly diverse fungal family containing a variety of pathogens affecting many economically important crops. Mitochondria plays a major role in fungal metabolism and fungicide resistance but up until now only two annotated mitochondrial genomes have been published in this family. We sequenced and annotated mitochondrial genomes of selected Mycosphaerellaceae species that diverged ∼66 MYA. During this time frame, mitochondrial genomes expanded significantly due to at least five independent invasions of introns into different electron transport chain genes. Comparative analysis revealed high variability in size and gene order among mitochondrial genomes even of closely related organisms, truncated extra gene copies and, accessory genes in some species. Gene order variability was common probably due to rearrangements caused by mobile intron invasion. Three three cox1 copies and bicistronic transcription of nad2-nad3 and atp6-atp8 in Pseudocercospora fijiensis were confirmed experimentally. Even though we found variation in mitochondrial genome composition, there was no evidence of hybridization when comparing nuclear and mitochondrial dataset sets for fungal plant pathogens analyzed here. Disentangling the causes of variation in mitochondrial genome composition in plant pathogenic fungal move us closer to understanding the molecular mechanisms responsible for vital functions in fungi ultimately aiding in controlling these diseases.


Introduction
From their origin as an early alpha proteobacterial endosymbiont to their current state as cellular organelles, large scale genomic reorganization has taken place in the mitochondria of all main eukaryotic lineages (Aguileta et al., 2014). Although fungi are a lineage more closely related to animals (Clark-Walker, 1992;Hausner, 2003), mt genomes in these organisms show signals of recombination, a characteristic that is more similar to plant mtDNAs (Clark-Walker, 1992;Bernt, Braband, et al., 2013;Gualberto et al., 2014;Lavrov and Pett, 2016;Ladoukakis and Zouros, 2017). Fungi have large intergenic regions and intron characteristics that are more similar to plant mitochondrial genomes (Zhang, 1995;Gualberto et al., 2014). Fungal mt genomes range in size from 1.14 kb (Spizellomyces punctatus) to 236 kb (Rhiczotonia solani) with an average of 62.68 kb (Organelle Genome Resource Database of NCBI). They can be circular or linear, are usually AT enriched, and their size variation is mostly due to presence or absence of accessory genes, mobile introns and size variation in intergenic regions. Their core gene content is largely conserved, but their relative gene order is highly variable, both between and within the major fungal phyla (Wolf and Giudice, 1988;Paquin et al., 1997;Aguileta et al., 2014). Fungal mitochondria have introns and intronic open reading frames (ORFs) classified as group I and group II introns, which differ in their sequence, structure and splicing mechanisms (Kawano, Takano and Kuroiwa, 1995;Burger, Gray and Lang, 2003;Hausner, 2003;Franz Lang, Laforest and Burger, 2007;Alverson et al., 2010;Bernt, Braband, et al., 2013).
Typically, group-II introns contain ORFs that code for reverse-transcriptase-like proteins. In contrast, group-I introns can encode proteins with maturase and/or endonuclease activity (Hausner, 2003). Because of the limited comparative analysis of complete fungal mt genome sequences, it has been difficult to estimate the timeframes and molecular evolution associated with mt genes or genomes (Torriani et al., 2014).
The high copy number of mitochondrial genomes per cell facilitates the recovery of complete sequences by highthroughput sequencing technologies. Also, fungal inheritance of mitochondria is mostly uniparental and characterized by non-random segregation (Griffiths, 1995;Basse, 2010). Because of the limited availability of complete mt genome sequences from fungal pathogen species, it has been difficult to estimate the timeframes associated with mt genome evolution and intron fixation or invasion of introns into mt genes or genomes.
Comparative analysis of nuclear genomes from Mycospharellaceae and particularly the Sigatoka disease complex already provided new insights into the evolutionary mechanisms associated with nuclear genome evolution (Chang et al., 2016). Comparisons between nuclear and mitochondrial genomes are useful to identify hybridization, introgression or incomplete lineage sorting among other evolutionary events (Depotter et al., 2018;Gonthier et al., 2018). A series of reports of emerging fungal hybrids, all involving plant pathogens was reviewed by Brasier (2000). New fungal hybrids are reported on a regular basis (Depotter et al., 2018;Giordano et al., 2018). Thus, investigating the evolutionary history of both nuclear and mitochondrial genomes is a matter of interest for fungal pathogens. Up until now no studies compared the mitochondrial genomes of Mycosphaerellaceae species using a timeframe.
Diseases caused by these three pathogens induce plant physiological alterations including a reduction in photosynthetic capacity, crop yield, and fruit quality . Symptoms are very similar for all species: P. musae causes yellow Sigatoka, producing mild symptoms in the form of necrotic yellow spots in leaves; P. fijiensis causes black Sigatoka which is more agressive and produces black necrotic leaf lesions, and P. eumusae disease is even more severe than black Sigatoka but induces very similar symptoms. Despite such differences in virulence, the three species share a common hemibiotrophic lifestyle and disease cycle on their host (Chang et al., 2016). Sigatoka disease complex causal agents formed a robust clade, with P. fijiensis diverging earlier (39.9 -30.6 MYA) than the sister species P. eumusae and P. musae (22.6 -17.4 MYA) (Ohm et al., 2012;Arango Isaza et al., 2016;Chang et al., 2016). Evolutionary mechanims of mt genome evolution in the Sigatoka disease complex could be better understood using fossil calibrated phylogenies and relaxed clocks for comparative analysis of these pathogens.
Fungal plant pathogens can rapidly develop molecular mechanisms of resistance to commercially available fungicides such as the case of resistance to QoI fungicides by punctual mutations in the cytochrome b mitochondrial gene (Sierotzki et al., 2000;Ma and Michailides, 2005;Cañas-Gutiérrez et al., 2006.Therefore, increasing the number of fungicide applications necessary to control the disease (Cools and Fraaije, 2013;Kitchen et al., 2016). The extensive use of chemical fungicides increase the production cost and is harmful to both environment and crop workers (De Lapeyre et al., 2010;Churchill, 2011). Research related to the pathogen vital functions, such as energy production, must be intensified in order to develop new and more effective disease control strategies for food safety, human health and environment.
Here, we sequenced and annotated mitochondrial genomes of selected Mycosphaerellaceae species that diverged ~66 MYA. We performed a comparative analysis of mitochondrial genomes of Pseudocercospora eumusae, Pseudocercospora fijiensis, and Pseudocercospora musae (causal agents of Sigatoka); Sphaerulina musiva and Sphaerulina populicola (causal agents of leaf spot and canker diseases in poplar); Zasmidium cellare (ethanol metabolizing, wine cellar mold) and Zymoseptoria tritici (wheat plant pathogen causing Septoria leaf blotch), based on its gene organization and genetic content. Comparative analyses showed that the major sources of mitochondrial genome size variation were differences in the number and size of introns in the core mt protein-encoding genes and differences in content of free standing and intronic Homing Endonuclease Genes (HEGs), hypothetical proteins and accessory genes. Although the functions and evolutionary origins of HEGs and accessory genes are unclear, their presence illustrates the plasticity of fungal mt genomes (Torriani et al., 2014). No evidence of incongruences was found between previously published nuclear phylogenies and our mt topology. New fossil calibrations for the Sigatoka's complex species and mt comparative analysis aid in our understanding of the tempo and mode of evolution of these plant fungal pathogens.
We sequenced a total DNA sample of Pseudocercospora fijiensis isolate 081012 from Carepa, Antioquia, Colombia at North Carolina University using 100bp paired-end reads and Hiseq 2500 system® (Illumina, Albany, New York). Read quality was assessed by FastQC (Andrews S., 2010) and low quality reads and/or bases were trimmed using Trimmomatic (Bolger, Lohse and Usadel, 2014). High quality reads were assembled using Spades (Bankevich et al., 2012) at different k-mer sizes (k = 61, 71, 81, and 91) and an assembly with the highest N50 value and assembly size was selected. The chosen assembly was scaffolded by SSPACE (Boetzer et al., 2011), remaining gaps in the scaffolds were closed by GapFiller (Nadalin, Vezzi and Policriti, 2012) and a final genome assembly was evaluated by REAPR (Hunt et al., 2013). Recovered contigs were then mapped to P. fijiensis mitochondrial genome available in MycoCosm (Grigoriev et al., 2014) using Geneious 9.1.5 (Kearse et al., 2012).
Draft mitochondrial genome sequence was deposited in GenBank (accession number: MK754071).
Mitochondrial contigs for Pseudocercospora musae, P. eumusae and Sphaerulina populicola were filtered from whole genome de novo assembled contigs (Dhillon et al., 2015;Chang et al., 2016). A BLASTn (Camacho et al., 2009) was ran using a BLAST database generated with the list of assembled contigs, and using as query a multi-FASTA with Electron Transport Chain genes (ETC) from published mitochondrial genomes. A reference based assembly was also performed using previously assembled contigs, and whole genomes as queries (Table 1) and ETC genes of previously published mitogenomes as reference. Contigs recovered with the above-mentioned approaches where reassembled. Reference based assembly, and de novo assembly were performed using the Geneious software (Kearse et al., 2012). An iterative mapping approach of clean reads of Sphaerulina musiva mitochondrial genome (SRA: SRR3927043) was performed against a reference cox1 to generate an assembly of the mitochondrial genome. Not annotated drafts of mitochondrial genomes of Pseudocercospora fijiensis, P. eumusae, P. musae, S. Populicola and Pseudovirgaria hyperparasitica were directly downloaded from databases (Table 1).

Annotation
Draft mitochondrial genomes of Pseudocercospora fijiensis, P. eumusae, P. musae, Sphaerulina populicola, S. musiva and Pseudovirgaria hyperparasitica were annotated using a combination of programs as follows: Predicted ORFs were determined with a translation code for "mold mitochondrial genomes" using Geneious 9.1.5 (Kearse et al., 2012). Genes were identified with BLASTP (Johnson et al., 2008) against the non-redundant protein database at NCBI and with the MITOS web server (Bernt, Donath, et al., 2013). Protein domains and sequence patterns were searched using PFAM (Finn et al., 2014) and PANTHER (Mi et al., 2017). Annotations obtained for the aforementioned approaches were contrasted and inconsistencies regarding length and position of genes were solved through multiple alignment of related gene sequences using MUSCLE (Edgar, 2004) and CLUSTALW (Larkin et al., 2007). We considered as hypothetical proteins those ORFs larger than 300bp that did not show results with the above-mentioned annotation strategies.

PCR amplifications of cox1 gene copies in P. fijiensis
A PCR assay to confirm the presence of three different copies for cox1 in P. fijiensis mitochondrial genome was performed. Primers were designed to amplify regions between cox1_1 and cox1_2, between cox1_2 and cox1_3 and between cox1_1 and cox1_3 ( Figure 1). PCR amplifications were carried out in a total volume of 10μl, containing 20ng genomic DNA, 0.15μM of each primer, 1x PCR buffer (without MgCl 2 ), 0.75mM MgCl 2 , 4μM of each dNTP and 0.65 U recombinant Taq DNA polymerase (ThermoFisher, MA, USA). Cycling parameters were: 3 min at 94°C, followed by 35 cycles of 30s at 94°C, 30s 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® under UV light.

RT -PCR assays for mitochondrial gene pairs of P. fijiensis
Total RNA was extracted from P. fijiensis mycelium after fifteen days of culture using TRIzol® reagent (Life Technologies) according to the manufacturer's instructions. RNA concentration was measured at 260 nm, using a NanoDrop ND-1000 UV-Vis Spectrophotometer (NanoDrop Technologies, Thermo Fisher). DNAse I (Thermo Scientific) was used to treat RNA and cDNA molecules were synthetized using Maxima First Strand cDNA Synthesis Kit (Thermo Scientific) according to manufacturer's instructions and used as template for amplification.
Primers were designed such that the amplified product encompass the end of one gene and the beginning the other for gene pairs nad3-nad2 and atp8-atp6 and amplifying partially both genes and their intergenic sequences. PCR amplification was performed and PCR products were analyzed by 1% agarose gel electrophoresis and gel purified using the GFX PCR DNA and gel band purification kit, according to the manufacturer's instructions (GE Healthcare). All purified PCR products were sequenced by Sanger Technology at Macrogen Inc (http://www.macrogen.com).

Phylogenetic analysis and divergence times estimates
To analyze the evolution of plant fungal pathogens through time and across sampled species, we first reconstructed a phylogenetic tree, then calculated reliable diversification times and last mapped mt molecular features such as size, gene order and others, in an evolutionary context. Since most mt genomes are uniparentally inherited we used our mt phylogeny to compare topologies with already published nuclear ones for the Sigatoka disease complex species (Arango Isaza et al. 2016;Chang et al. 2016). Up until now only nuclear markers and strict clock calibration have been used to calculate diversification times for the Sigatoka disease complex species. We aim to compare these analyses with fossil calibrated Bayesian phylogenies and mitochondrial markers.
A bayesian phylogeny and divergence time analysis were carried out using BEAST2 version 2.5.1 (Bouckaert et al., 2014). Separate partitions for each gene were created with BEAUti2 (BEAST 2 package). More suitable substitution models for each gene were found using the software package jModelTest2 (Guindon, Gascuel and Rannala, 2003;Darriba et al., 2012) according to the Bayesian Information criterion (BIC) (Schwarz, 1978). To accommodate for rate heterogeneity across the branches of the tree we used an uncorrelated relaxed clock model (Drummond et al., 2006) with a lognormal distribution of rates for each gene estimated during the analyses. We also used a strict clock to further compare results. The fossil Metacapnodiaceae (Schmidt et al., 2014) was used assuming to be a common ancestor of the order Capnodiales with minimum age of 100 MYA (gamma distribution, offset 100, mean 180, maximum soft bound 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 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). To ensure congruence we ran the 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 1.5 (Rambaut and Drummond, 2009). After removal of 20 % of each run as burn in the remaining trees were combined using LogCombiner (part of the BEAST package), summarized as maximum clade credibility (MCC) trees in TreeAnnotator (part of the BEAST package), and visualized using FigTree. (Rambaut, 2006) 3. Results

Inferred mt phylogeny and diversification times
A comparative framework was used to understand mitochondrial evolution in fungal plant pathogens. Since most mt genomes are uniparentally inherited we built a mt phylogeny in order to compare topologies with already published nuclear ones for the Sigatoka disease complex species. Our mitochondrial topology supports previously published nuclear phylogenies (Arango Isaza et al. 2016;Chang et al. 2016) and all nodes were well supported ( Figure 6). Ideally, phylogenies should be built with single-copy genes to avoid paralogy (Guindon, Gascuel and Rannala, 2003). However, most core mitochondrial genes for species included in this phylogeny have more than one copy. ML phylogenies were built for each gene including all copies per specie. All gene copies were monophyletic in all cases (Appendix A3-14), indicating that duplicated genes are still poorly differentiated and homology among copies from same species can still be detected. One complete copy with start and stop codons was selected to build a species tree. A phylogeny using only a single copy rns gene was also used. The same topology was always recovered (Appendix A15). Finally, we used ten of fourteen core genes for divergence time estimation because when using only the rns gene, the 95% highest posterior density (HPD) intervals were significantly broader.
A robust phylogeny with posterior probabilities greater than 0.97 was recovered, containing two main linages (Pleosporales + Capnodiales). Divergence time estimates using fossil calibration are shown in Figure 6 (Figure 6).

Mitochondrial genome content and genome organization
We compared mitochondrial genomes of Pseudocercospora fijiensis, P. eumusae, P. musae, Sphaerulina populicola, and S. musiva, with previously published mitochondrial genomes of Zymoseptoria tritici and Zasmidium cellare. They showed a homogeneous GC content between 27% and 32%, high variability in genome size, ranging from 23743 bp in "cellar mold" Zasmidium cellare to 136906 bp in poplar pathogen Sphaerulina populicola (Appendix 1). Annotated genes of seven mitochondrial genomes analyzed included a set of conserved mitochondrial protein coding genes, namely the subunits of the electron transport chain (ETC) complex I (nad1, nad2, nad3, nad4, nad4L, nad5 and nad6), complex III (cox1, cox2 and cox3), complex IV (cob) and ATP synthase subunits (atp6, atp8 and atp9). A small ribosomal subunit RNA (rns), a large ribosomal subunit rRNA (rnl), and a set of transfer RNAs (tRNAs). Five of seven genomes contained hypothetical proteins, accessory proteins, additional copies of genes and genetic mobile elements (Table 3). ETC gene order among Mycosphaerellaceae mitochondrial genomes vary greatly, except for sister species Sphaerulina populicola and S. musiva which showed the same gene order (Figure 3). Despite this variability, gene pair order was conserved for nad4L-nad5, nad3-nad2 and atp8-atp6 among all mitochondrial genomes (Figure 3). We mapped RNA seq assembled transcripts to P. fijiensis, S. musiva, P. musae and, P. eumusae mitochondrial genomes and found that these neighbor genes were contained in the same transcript suggesting that they were transcribed as a single mRNA.
All annotated genomes were found to have duplications of some ETC genes. They were considered as gene duplications because they presented the same PFAM and PANTHER domains and high local sequence similarity (>90%). Duplicated gene alignments showed that gene copies were not identical in length ( Figure 1C), and only one copy had a complete coding sequence. The mitochondrial genome of P. fijiensis has the greater number of duplicated genes, with five ETC genes having two to three copies (Tables 2 -3). Duplicated genes were found either close to each other, or disperse in the genomes.

Experimental confirmation of gene copies, bicistronic transcription and intron content in P. fijiensis mt genome
We used P. fijiensis to prove experimentally the three most characteristic features of Mycosphaerellacea genomes found through bioinformatics analysis: Presence of ETC copies, bicistronic transcription of some genes and intron presence in ETC genes. First, PCR amplifications were performed to confirm the presence of three cox 1 copies in P. fijiensis mitochondrial genome. Primers were designed to amplify regions among the three different copies of cox1 ( Figure 1D). We chose to amplify this region because it was the only case of tandem duplication among all genomes. Other gene copies in P. fijiensis as well as in other species were found sparse. The three regions amplified with bands of the expected size, confirming that the different gene copies exist ( Figure 1A and B).
Second, in P. fijiensis mitochondrial genome there are two genes with introns (cob and nad5), so we confirmed the presence of cob intron experimentally but amplification of nad5 intron was unsuccessful because the PCR product exceeds 3Kb and there were few optimal primer design sites surrounding the targeted region (Appendix 2). Finally, RT-PCR amplification and subsequent Sanger sequencing confirmed polycistronic expression for nad3-nad2, atp8-atp6 gene pairs in P. fijiensis (Figures 3-4).

Homing endonucleases invasion and putative horizontal gene transfer
Differences in genome sizes among members of the Mycosphaerellaceae were attributed to the presence of different amounts of accessory genes, gene copies and intron mobile sequences (Figure 2, Tables 2-3). Accessory genes, gene copies and intron mobile sequences are specie specific. A commonly observed feature in Mycosphaerellaceae mitochondrial genomes was an invasive presence of sequences containing LAGLIDADG or GIY-YIG domains related to homing endonuclease genes (HEG) (Chevalier and Stoddard, 2001). These sequences caused fragmentation in ETC genes. In almost all instances, pieces of an ETC genes were collinearly distributed and each fragment was found to be followed by an insertion of an HEG related sequence. An extreme case was found in cox1 in Sphaerulina populicola (CDS: 1542bp), which has twelve fragments distributed along

Discussion
Our principal objective was to investigate the evolutionary history of mitochondrial genome evolution in fungal plant pathogens such Mycosphaerellaceae. We found unequivocal evidence of at least five independent horizontal gene transfer events related to variation in genome composition in this family. Even though we found variation in mitochondrial genome composition there was no evidence of hybridization when comparing nuclear and mitochondrial dataset sets for fungal plant pathogens analyzed here.
Genome variation in size and organization is not correlated with phylogenetic relationships, nor does it carry phylogenetic signal based on our current sampling. Ascomycete's mitochondrial genomes like most mitochondrial genomes along the tree of life generally consist of a highly conserved core of protein encoding genes: 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) (Clark-Walker, 1992). These were also the genes found in mitochondrial genomes of Mycosphaerellaceae. Besides this common set of genes, a variable number of free standing open reading frames of unknown function and a variable number of introns related with homing endonucleases (HEF), often including GIY-YIG or LAGLIDADG protein domains were found. These ORFs ranged from one in Pseudocercospora musae and P. eumusae to 36 in Sphaerulina populicola (Table 2).
Despite not being ubiquitous in Ascomycete mitochondrial genomes, these ORFs are common (Clark-Walker, 1992;Adams and Palmer, 2003;Burger, Gray and Lang, 2003;Hausner, 2003;Lodish, 2003;Bernt, Braband, et al., 2013;Jelen et al., 2016). HEG invasion has been previously described in fungi, where fragmented genes are the rule (Dassa et al., 2009). However, it is not common to find highly fragmented genes as the cox1 of Sphaerulina populicola ( Figure 5). A similar case was reported in the mitochondrial genome of Sclerotinia borealis, where researchers found thirteen introns of cox1 and truncated copies of ETC genes (Mardanov et al., 2014).
Differences in mitochondrial genomes sizes among members of Mycosphaerellaceae are attributed to heterogeneous content of accessory genes, gene copies and intron mobile sequences in some of the genomes ( Figure 2, Tables 2-3). These elements are species specific, indicating that they have been acquired through horizontal transfer through at least five independent events in their evolutionary history ( Figure 5). In Pseudocercospora and Sphaerulina species we found duplicated genes, ORFs related to HEGs and hypothetical proteins, all contributing to the variation in size. Mitochondrial genomes from P. musae, P. eumusae , and S.
Populicola had also accessory genes including RNA and DNA polymerases, reverse transcriptases and transposases. Even for sister species such as Sphaerulina musiva and S. populicola who share core gene content (ETC genes) and gene order, genomes sizes differ greatly due to a higher content of homing endonucleases HEG related ORFs and the presence of accessory genes in S. populicola (Appendix 1). Invasion of mobile elements and accessory genes as main source of size variability among the genomes of phylogenetically related organisms was observed before in a study where nine mitochondrial genomes of Aspergillus and Penicilium were compared (Joardar et al., 2012). In the phytopathogenic fungus Sclerotinia borealis, they also found expansion of the genome due to plasmid like sequences integrated in the mitochondrial DNA and HEG related ORFs (Mardanov et al., 2014).
Mitochondrial genomes of S. musiva, S. populicola, P. fijiensis, P. eumusae and P. musae have partial extra copies of ETC protein coding genes (Table 3) and tRNAs, while the fully annotated published Mycosphaerellaceae genomes of Zasmidium cellare and Zymoseptoria tritici do not possess any duplication (Torriani et al., 2008;Goodwin et al., 2016). Mitochondrial gene duplication is seldom described in fungi. Maradanov et al. (2014) found duplications of some mitochondrial regions, resulting in the appearance of truncated extra copies of atp9 and atp6 in phytopathogenic fungus Sclerotinia borealis (Mardanov et al., 2014).
Two incomplete copies of atp6 were found on different strands of the mitogenomic DNA of Shiraia bambusicola, a Pleosporales member used as outgroup in this study (Shen et al., 2015). Three hypotheses can be drawn regarding mitochondrial gene duplications: First, copies could be undergoing subfunctionalization (Conant and Wolfe, 2008) where recombination events of fission and fusion are the drivers of such a process (Burger, Gray and Lang, 2003), however we did not see two functional copies of any ETC genes in any genome. . Second, they can be truncated copies accounting for pseudogenes (Conant and Wolfe, 2008). We therefore suggest that additional Mycospharellaceae fungi need to be available before stronger conclusions can be drawn about the proximal causes of gene duplication. Further gene editing and virulence assays to test these hypotheses can shed light on fungal adaptation. Mycosphaerellaceae gene order is not conserved among mitochondrial genomes (Figure 3). Even though divergence times within sister clades Pseudocercospora 13.31  MYA, 95% HPD) are roughly the same (Figure 6), gene order is highly conserved in the Sphaerulina species but not in the Pseudocercospora ones ( Figure 2). Nonetheless, some pairs of genes are always found together atp6-atp8, nad2-nad3, nad4L-nad5. Here, bicistronic transcription of atp6-atp8, nad2-nad3 through RT-PCR was confirmed (Figures 3-4). The presence of mitochondrial transcripts containing more than one ORF make sense due to mt endosymbiotic origin from an α -proteobacteria (Schäfer, Hansen and Lang, 2005;Barbrook et al., 2010;Chrzanowska-Lightowlers, 2015). We attribute differences in gene order to genomic rearrangements, which could be caused by different processes, such as fusion, fission, recombination, plasmid integration and mobility among mitochondrial genomes (Kawano, Takano and Kuroiwa, 1995). Indeed, plasmid integration appears to have a role in recombination thanks to the presence of long inverted repeats that allows homologous recombination (Salavirta et al., 2014). Since most mt genomes are uniparentally inherited we built a mt phylogeny in order to compare topologies with already published nuclear ones for the Sigatoka disease complex species (Arango Isaza et al. 2016;Chang et al. 2016). When phylogenetic comparisons between mitochondrial and nuclear data sets in fungi differ, hybridization, introgression or incomplete lineage sorting among other events could explain incongruences (Depotter et al. 2018;Giordano et al. 2018 (Chang et al., 2016). This suggests that a maximum of 78 MY elapsed between the split of Mycospharellaceae species and the at least five independent intron invasion events into their mt genomes. The most parsimonious explanation for these independent intron invasions is horizontal transfer from donor species. Fungal species can directly exchange genetic material through mycelial or conidial anastomoses if they share the same ecological niche (Torriani et al., 2014). Another possibility is indirect transmission of genetic material through virus that can move DNA or RNA from one fungal species to the other (Rosewich and Kistler, 2000) As far as the Sigatoka disease complex species diversification we found an early divergent Pseudocercospora fijiensis that splits from the sister species P. musae + P. eumusae 13.31 MYA (9.49-17.28, 95% HPD); while the sister species P. eumusae and P. musae split from their shared ancestor in the late Miocene 8.22 MYA (5.6 -11.07, 95% HPD) ( Figure 6). Chang et al. (2016) using nuclear phylogenies estimated the splitting of P. fijiensis from their last common ancestor with P. musae + P. eumusae to be between 39.9-30.6 MYA and the splitting of P. musae and P. eumusae to be between 22.6 -17.4 MYA. Our results naturally prompt a further question: Was the origin of the Sigatoka disease complex through host-tracking evolution? This hypothesis explain that a pathogen is likely to be younger that the host until changes related to the genetics of the pathogens or/and exogenous factors have prompted the observed alterations in their virulence spectra and the sudden flare up and over dominance of one species over the others (Stukenbrock and Mcdonald, 2008). We are currently limited in terms of a good taxonomic sampling and biogeographical analysis for Pseudocercospora species to answer this question.
Strong differences in diversification times compared to Chang et al. (2016) were found due to differences in calibrations points and methods. Chang et al. (2016) used r8s Likelihood methods (Sanderson, 2003) and a calibration at the origin of the Dothideomycetes crown group (394-285 MYA) using previous bayesian estimations (Gueidan et al., 2011). Here bayesian analysis in BEAST2 (Bouckaert et al., 2014) and a fossil calibration using a Metacapnodiaceae (Schmidt et al., 2014) were used, assuming the common ancestor of the order Capnodiales to be constrained to a minimum age of 100 MYA (gamma distribution, offset 100, mean 180, maximum soft bound 400).
Diversification ages estimated using mitochondrial markers and bayesian analysis using both relaxed and strict clock models were younger than those calculated by Chang et al. (2016) using nuclear markers, penalized maximum likelihood analysis and strict clock, except for Sphaerullina. Divergence times were also estimated using strict clock to validate differences between mean node ages using relaxed lognormal clock versus strict clock.
For a strict clock, 95% highest posterior density (HPD) intervals were significantly broader (27-131 MY) in comparison to (6-23 MY) and mean node ages were among 9 -43 MY older (Appendix A16, Table 4). Giving the low taxonomic sampling of this study is possible that diversification times can change when adding more taxa to the phylogeny. However, we estimate this change will not be dramatic.
Mitochondrial genomes in Mycosphaerellaceae are highly variable in size and gene order, due to at least five independent intron invasion and horizontal gene transfer during their evolutionary history. Presence of duplicated genes, homing endonuclease genes (HEG) invasion and accessory genes was observed. Despite their order variability some genes were always recovered as neighbors in all mt genomes analyzed showing polycistronic transcription due to their endosymbiotic origin from an α -proteobacteria. Further gene editing and virulence assays will be important to shed light on fungal adaptation and more effective disease control strategies. Our mt phylogeny did not differ to previously published nuclear topologies. Fossil calibrated phylogenies reported for first time here for fungal plant pathogens had earlier diversification times for the origin of all the species involved compared to previous studies. Phylogenomic studies including a good taxonomic sampling and biogeographical analysis for Pseudocercospora species will further clarify whether the Sigatoka disease causing species virulence flare up after Musa domestication.

Data Availability Statement:
The annotated mitochondrial genome of Pseudocercospora fijiensis , atp6-atp8 and nad2-nad3 sequences will be accessible at GeneBank upon paper acceptance . Phylogeny xml archive and other mitochondrial annotated genomes are available at https://sites.google.com/site/tatianatatianaarias/publications/arcila-et-al-supplementary-materials