life
Article
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
Juliana E. Arcila-Galvis
1
, Rafael E. Arango
2,3
, Javier M. Torres-Bonilla
2,3,4
and Tatiana Arias
1,
*
,†
Citation: Arcila-Galvis, J.E.; Arango,
R.E.; Torres-Bonilla, J.M.; Arias, T.
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. Life 2021, 11,
215. https://doi.org/10.3390/
life11030215
Academic Editor: Andrea Luchetti
Received: 20 January 2021
Accepted: 8 February 2021
Published: 9 March 2021
Publisher’s Note: MDPI stays neutral
with regard to jurisdictional claims in
published maps and institutional affil-
iations.
Copyright: © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).
1
Corporación para Investigaciones Biológicas, Comparative Biology Laboratory, Cra 72A Medellín, Antioquia,
Colombia; juearcilaga@unal.edu.co
2
Escuela de Biociencias, Universidad Nacional de Colombia-Sede Medellín, Cl 59A Medellín, Antioquia,
Colombia; rarango@cib.org.co (R.E.A.); javier.torres@colmayor.edu.co (J.M.T.-B.)
3
Corporación para Investigaciones Biológicas, Plant Biotechnology Unit,
Cra 72A Medellín, Antioquia, Colombia
4
Colegio Mayor de Antioquia, Grupo Biociencias, Cra 78 Medellín, Antioquia, Colombia
* Correspondence: tatiana.arias48@tdea.edu.co; Tel.: +57-300-845-6250
† Current address: Tecnológico de Antioquia, Cl 78B Medellín, Antioquia, Colombia.
Abstract:
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.
Keywords:
banana; diversification times; mitochondrial genome; Mycosphaerellaceae; plant pathogens;
Pseudocercospora; sigatoka disease
1. 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
–
8
]. Three My-
cosphaerellaceae 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
Life 2021, 11, 215. https://doi.org/10.3390/life11030215 https://www.mdpi.com/journal/life
Life 2021, 11, 215 2 of 19
growers [
9
,
10
]. Diseases caused by these three pathogens induce plant physiological
alterations including a reduction in photosynthetic capacity, crop yield, and fruit qual-
ity [
9
]. The Sigatoka disease complex causal agents form a robust clade, with P. fijiensis
diverging earlier (39.9–30.6 MYA) than sister species P. eumusae and P. musae (22.6–17.4
MYA) [
10
–
12
]. Among them, Pseudocercospora fijiensis (teleomorph Mycosphaerella fijiensis)
is the causal agent of black leaf streak disease (BLSD; aka Black Leaf Spot Disease), one the
most damaging and costly diseases for banana and plantain worldwide [13].
Fungal mitochondrial genomes (mitogenomes) are circular or linear, usually AT en-
riched and range in size from 1.1 kb (Spizellomyces punctatus) [
14
] to 272 kb (Morchella
importuna) [
15
]. Their size variation is mostly due to the presence or absence of accessory
genes including RNA and DNA polymerases, reverse transcriptases and transposases,
mobile introns, and size variation in intergenic regions [
16
,
17
]. In spite of the variation in
size, their core gene content is largely conserved, even though their relative gene order is
highly variable, both between and within the major fungal phyla [
18
–
20
]. Mitogenomes
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 [
16
,
21
–
25
].
Typically, group-II introns contain ORFs that code for reverse-transcriptase-like proteins. In
contrast, group-I introns encode proteins with maturase and/or endonuclease activity [
16
].
Because of the limited comparative analysis of complete fungal mitogenome sequences,
it has been difficult to estimate the timeframes and molecular evolution associated with
mitochondrial genes or genomes [26].
Mitochondria have proven to be useful in evolutionary biology and systematics be-
cause they contain their own genome capable of independent replication, uniparental
inheritance [
27
], near absence of genetic recombination, and uniform genetic backgrounds
for some species [
28
]. Attempts to determine a time frame for fungal evolution are ham-
pered by the lack of reliable fossil records. Hence, so far studies have focused on relating
rates of DNA base substitutions and molecular clocks [
29
], based on the assumption that
mutation rates of nuclear genes are similar to their counterparts in organisms with datable
fossils [
19
]. Mitochondria plays a major role in fungal metabolism and fungicide resistance
but until now only two annotated mitogenomes have been published in Mycosphaerel-
laceae (Zasmidium cellare and Zymoseptoria tritici) [
30
,
31
]. Sigatoka disease comparative
mitogenome studies will provide answers on the evolution and adaptation of these plant
pathogenic fungi.
This study aimed to: (i) sequence and characterize the complete mitogenome of
Pseudocercospora fijiensis; (ii) compare mitogenomes of P. fijiensis with closely related species
P. eumusae, and P. musae (causal agents of Sigatoka), and species with publicly available
high throughput data such as Sphaerulina musiva and S. populicola (causal agents of leaf spot
and canker diseases in poplar); and (iii) estimate timeframes and mitochondrial molecular
evolution using fossil records to calibrate the Sigatoka disease complex phylogeny. We
found that in mitogenomes analyzed herein, there were differences in content of free-
standing and intronic Homing Endonuclease Genes (HEGs), genes coding for hypothetical
proteins, and accessory genes such as DNA/RNA polymerases, reverse transcriptases and
transposases. This work contributes to the understanding of mitogenome organization in
Mycospharellaceae. In addition, new fossil calibrations for the Sigatoka’s complex species
and mitochondrial comparative analysis aid in our understanding of the tempo and mode
of evolution of these plant fungal pathogens.
2. Materials and Methods
2.1. Fungal Strain, DNA Extraction, and Library Construction and Sequencing
P. fijiensis isolate 081022 was obtained from naturally infected banana leaves coming
from a commercial plantation located in Carepa, Antioquia, Colombia. Taxonomic affil-
iation has been confirmed based on both morphological criteria and Polymerase Chain
Reaction (PCR) [
32
]. For DNA extraction mycelia from 7-day old culture were transferred
to potato dextrose broth and incubated for 5–7 days at room temperature in a rotary
Life 2021, 11, 215 3 of 19
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
®
.
2.2. Sequence Sources, Data Filtering and Assemblies
Eleven mitochondrial genome (mitogenomes) sequences were used for this study.
Seven belonging to Mycosphaerellaceae: Pseudocercospora fijiensis (syn. Mycosphaerella
fijiensis), Pseudocercospora eumusae (syn. Mycosphaerella eumusae), Pseudocercospora musae (syn.
Mycosphaerella musicola), Sphaerulina musiva (syn. Septoria musiva), Sphaerulina populicola
(syn. Septoria populicola), Zasmidium cellare, and Zymoseptoria tritici (syn. Mycosphaerella
graminicola). One species is from Capnodiales: Pseudovirgaria hyperparasitica and three are
from Pleosporales, the sister group of Capnodiales: Didymella pinodes (syn. Mycosphaerella
pinodes), Parastagonospora nodorum (syn. Phaeosphaeria nodorum), and Shiraia bambusicola
(Table S1). Mitogenomes were obtained either from our own sequencing data, or sequence
data available at GenBank [
34
], RefSeq [
35
] or MycoCosm [
36
]. Authors, seq ID and
databases are listed in Table S1.
Read quality was assessed using FastQC v. 0.11.5 [
37
] for the P. fijiensis isolate 081022
raw reads recovered here. Low-quality reads and/or bases were trimmed using Trim-
momatic version 0.36 [
38
]. First, we de novo assembled whole DNA using Spades 3.9.0
(parameter “-careful”) [
39
] at different k-mer sizes (k = 61, 71, 81, and 91). The assem-
bly with the highest N50 and assembly size was scaffolded by SSPACE version 3.0 [
40
].
Remaining gaps between scaffolds were closed using GapFiller version 1.10 [
41
] and a
final genome assembly was evaluated by REAPR version 1.0.18 [
42
]. Scaffolds from the
whole genome sequencing assembly were mapped to a draft and an unpublished P. fijiensis
mitogenome available in MycoCosm [
36
] using Geneious 9.1.5 [
43
]. Mitogenomes were also
filtered from de novo whole-genome assemblies for Pseudocercospora musae and P. eumusae
available online [
8
,
10
]. To separate mitochondrial contigs or scaffolds from the nuclear
contigs or scaffolds, we used BLASTn [
44
] and the Electron Transport Chain Conserved
Mitochondrial Protein Coding Genes (CMPCGs) compiled from published mitogenomes of
Zasmidium cellare, Zymoseptoria tritici, Didymela pinodes, Phaeosphaeria nodorum and Sharaia
bambusicola as queries [30,31,45,46].
Even though the Sphaerulina musiva mitogenome was available online [
8
] we reassem-
bled 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.
Annotated mitochondrial genomes filtered from whole-genome assembly projects for
Pseudocercospora musae and P. eumusae or reassembled for Sphaerulina musiva are available
in Figshare (dataset: https://doi.org/10.6084/m9.figshare.12101058.v1 (accessed on 10
February 2021)).
Life 2021, 11, 215 4 of 19
2.3. 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].
2.4. 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.
2.5. Transcriptome de novo Assembly
RNA-seq raw reads of S. musiva (SRR1652271) and P. fijiensis (SRR3593877, SRR3593879)
were downloaded from the European Bioinformatics Institute (EMBL EBI) database. Reads
were quality filtered and trimmed using BBDuck from BBtools (https://sourceforge.net/
projects/bbmap/ (accessed on 10 February 2021)) before carrying out transcriptome
de novo assemblies with Trinity version 2.3.1 [
54
]. The P. eumusae (GDIK00000000.1)
and P. musae assembled transcriptomes (GDIN00000000.1) were also downloaded from
GeneBank. RNA-seq Geneious 9.1.5. plugins were used to map the assembled transcripts
to mitogenomes of either S. musiva, P. fijiensis, P. eumusae or P. musae paying particular
attention to gene pairs atp6-atp8, nad2-nad3.
2.6. 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 manu-
facturer’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
Life 2021, 11, 215 5 of 19
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).
2.7. Identification of Repetitive Elements
Repetitive sequences in mitogenomes of Mycosphaerellaceae were identified and
annotated using Geneious Primer Tandem Repeats Finder and using a minimum repeat
length of 100, excluding repeats up to 10 bp longer [
43
]. Simple sequence repeat (SSR)
markers and loci were identified using the MicroSAtellite Identification tool (MISA) [
55
]
(https://doi.org/10.6084/m9.figshare.12101013 (accessed on 10 February 2021)).
2.8. Phylogenetic Analysis and Divergence Times Estimates
Until now, only nuclear markers and strict clock calibration have been used to calcu-
late 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 mi-
togenomes 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 us-
ing 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 Ran-
domized 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
![Figure 1. Circular map of the mitogenome of Pseudocercospora fijiensis. Genes are visualized by arrows, which are presented in a clockwise direction (forward). (A) Annotated introns, hypothetical Open Reading Frames (ORFs), group I mobile introns (Homing Endonuclease Genes (HEGs)) and genic regions; clusters of tRNA and rRNA are also indicated here. (B) Annotated tRNA and rRNA regions. Green color arrows: protein-coding genes; yellow color arrows: Coding regions CDS and introns; orange arrows: hypothetical ORFs; blue arrows: genes of large and small ribosomal subunits; bright yellow arrows: tRNAs genes; grey arrow: HEGs. Circular mitogenomes were generated by using the Geneious Primer [43].](/figures/figure-1-circular-map-of-the-mitogenome-of-pseudocercospora-hfku58rg.webp)
