The mitochondrion is an important organelle of eukaryotic cells, involved in multiple processes that are necessary for cell survival and reproduction. It is responsible for the ATP production, the universal energy-transfer biomolecule in living organisms, through the oxidative phosphorylation pathway. In addition to respiratory metabolism and energy production function, mitochondria is also involved in other processes such as senescence during the cell cycle and the maintenance of ion homeostasis (Burger et al., 2003; Basse, 2010; Chatre and Ricchetti, 2013).

Conversely to the majority of other intracellular membrane structures, mitochondria contain their own genome that is capable of independent replication and inheritance. Nonetheless, over evolutionary time, most of the genes in the initial endosymbiont have been transferred to the nuclear genome and only a few genes are retained in mitochondrial genomes (mitogenomes). Indeed, most of the mitochondrial proteins are encoded by the nuclear genome and transported to the mitochondria (Burger et al., 2003). Within each mitochondrion, there are typically multiple copies of the genome, which replicate independently of the nuclear genome and cell cycle (Chatre and Ricchetti, 2013). Furthermore, most of the mitochondrial genomes are composed of a single circular chromosome, although some organisms contain linear mitochondrial genomes made up of multiple subunits (Burger et al., 2003).

Fungal mitogenomes are generally characterized by the presence of 14 conserved protein-coding genes involved in electron transport (cytochrome oxidase subunits 1, 2, and 3, NADH dehydrogenase subunits 1, 2, 3, 4, 4L, 5, and 6) and ATP synthesis (ATP synthase subunits 6, 8, and 9, and apocytochrome b). In addition, mitochondrial genome typically has one of the small (SSU rRNA) and large subunits ribosomal RNAs (LSU rRNA), and a set of transfer RNA (tRNA) genes (Bullerwell and Lang, 2005; Al-Reedy et al., 2012; Ferandon et al., 2013; Brankovics et al., 2018). Despite the relatively conserved gene content, genome structural changes caused by insertion of repetitive elements, gain or loss of introns, and/or gene transfers to the nucleus often lead to mitogenome size variations (Burger et al., 2003).

Introns are present in many fungal mitogenomes. These sequences can be divided into two groups: group I and group II. Introns Group I (IGI) presents six subclasses and are more abundant in mitogenomes than the group II. IGI can encode cellular proteins, such as rps3 ribosomal protein or maturases like homing endonucleases (HEGs), which are likely involved in the splicing process (Hausner, 2003). IGI are located within genes and are associated with shaping of mitogenomes, gene rearrangements and also in the host fitness (Hausner, 2003; Belfort, 2017). Introns type II are self-splicing and consist in a sequence of approximately 600 nucleotides. Usually, they are similar to IGI and can be classified into four groups (Michel et al., 1989; Hausner, 2003).

A recent study (Sandor et al., 2018) based on 20 species/groups of fungi revealed, for most of them, higher genetic diversity of nuclear genes and genomes than for the mitochondrial ones. This pattern is more similar to the observed in plants but different from most of the animals. The mechanisms responsible for the structural variations in fungal mitochondrial genomes and its relationship with fungi divergence time are still poorly understood. Nevertheless, the increasing availability of fungal mitogenome sequences in public databases can help us to understand the factors that drive fungi mitogenome variations of fungi.

The order Hypocreales (Ascomycota) is a monophyletic group with more than 2700 species of fungi, distributed in 240 genera divided in sexual (teleomorph) and asexual (anamorph) morphs (Varshney et al., 2016; Wijayawardene et al., 2018). Organisms of this order show a great variability of ecological functions and may display roles ranging from saprobes to pathogens. Special interest has been given to economically important species in the fields of pharmacology and medicine, biological control and biotechnology (Rehner and Samuels, 1995; Varshney et al., 2016), such as Trichoderma harzianum. This species is widely distributed in soil and plants due to its diversity of ecological niches (Druzhinina et al., 2011; Contreras-Cornejo et al., 2016). Furthermore, T. harzianum acts as a biofertilizer, promoting plant growth, and as a biocontrol agent during plant infections (Druzhinina et al., 2011; Harman, 2011), and has been considered an attractive alternative to the usage of chemical fungicides (Schmoll et al., 2016). Despite its economic relevance and wide distribution, the nuclear genome of T. harzianum is not yet fully assembled, and the mitochondrial genome is still not available in public databases.

Herein, we performed the deep sequencing of T. harzianum HB324 mitochondrial genome and compared it with 34 other reference species from the Hypocreales order that were available in public databases. We performed de novo structural and functional annotation of all species analyzed. The standardized annotation allowed us to perform comparative analyses revealing high variation among the mitogenomes, even within individuals from the same family. Mitochondrial genome length varied from 24,565 to 103,844 bp. Despite the size variation, the copy number, size, and structure of protein-coding genes were highly conserved, suggesting that differences in the genome length are likely related to non-coding regions, such as introns, HEGs and unidentified ORF (uORF). Introns were the most widespread non-coding element found in mitogenomes. Analyses based on sequence similarity have revealed greater sharing of intronic sequences within families. In some cases, we observed conservation even in the position of intron insertion, such as the intron within gene rrnL, present in almost all species. In contrast, species with higher divergence time showed smaller similarity among intronic sequences. We also investigated whether the genetic divergence between fungal species would be positively correlated with the size of the mitogenome and the prevalence of protein-coding and non-coding mitochondrial elements. Fungi species that evolve faster presented a higher frequency of non-coding elements on their mitogenomes. Therefore, our data suggests that the difference in the size of the mitogenomes is due to the presence of non-coding components. Besides, species that evolve faster have a greater proportion of non-coding elements, corroborating their role in the shaping of fungi mitochondrial genomes.

Comparative genomics has been used to elucidate the population structure and to understand the functions of each species in an environment (Nadimi et al., 2016). In the last years, this methodology has revealed the main genomic organization of fungal species, factors related to fungal pathogenicity, production of substances and enzymes of commercial interest, and other aspects of fungal evolution and biology (Stajich, 2017). Mitogenomes are studied for phylogenetic relationships, phylogeography and structure dynamics in the population since their genomes present conserved gene content, rapid evolutionary rate and many copies in the host‘s cell (Moore, 1997). Herein, we present a new mitochondrial genome of Trichoderma harzianum and performed comparative analysis including 34 species from Hypocreales order to investigate variations in the mitochondrial genome size and structure, intron sharing, presence of NUMTs, and the relationship among divergence time and coding and non-coding elements that compose these mitogenomes.

Trichoderma species are found in almost every environment, mainly in soil and plant root ecosystems (Harman et al., 2004) and are widely applied in agriculture as biocontrol agents, inhibiting the growth of other fungi and nematodes (Lorito et al., 2010). Furthermore, some species can stimulate plant growth and development by stablishing mutualistic beneficial interactions (Bonfante and Requena, 2011). The species T. harzianum is widely used in agricultural industries due to its mycoparasitic activity, which inhibits many plant pathogens (Lorito et al., 2010). T. harzianum HB324 was isolated by our research group from Hevea brasiliensis (rubber tree) leaves and our data have shown that this isolate can act as biocontrol agent against the phytopathogen Colletotrichum sp. (Fonseca et al., 2018). Despite the economic and ecological importance of this fungus, its mitochondrial genome has not been yet available. In the current work, the complete sequencing of mitochondrial genome of T. harzianum HB324 isolate is presented, and its mitogenome has 32,277 bp, being the second largest of the genus and contain four uORFs, five introns and six domains of HEGs.

The remarkable variation in the mitogenomes sizes of the 35 fungal species evaluated in this study (from 24 to 103 Kb) was observed at both genus and species levels: within strains of Trichoderma, Beauveria, and Fusarium, for example, it reached up to 70 Kb. These results reinforce those previously reported. There are studies describing genomes of around 18 Kb in Hanseniaspora uvarum (Pramateftaki et al., 2006) and up to 235 Kb in Rhizoctonia solani (Losada et al., 2014). Variations in mitogenome sizes at genus and species level have also been reported within species of the genera Schizosaccharomyces (Bullerwell et al., 2003) and Fusarium (Al-Reedy et al., 2012), as well as within the species of Saccharomyces cerevisiae (Wolters et al., 2015), Cordyceps militaris (Zhang et al., 2015), and Rhizophagus irregularis (Formey et al., 2012).

The wide size range of the mitogenomes evaluated in our study may be explained in part by variations in length of intergenic regions, and differences in number and length of introns. The presence or absence of introns have been widely described as the major contributing factor for such variation among fungi (Jung et al., 2010; Joardar et al., 2012; Zhang et al., 2015). Among the fungal strains evaluated in our study, F. graminearum had the highest number of introns (35), accounting for 54% of the entire mitochondrial genome. Our results are similar to those reported to Podospora anserina (Ascomycota), for which 33 mitochondrial introns were detected (Cummings et al., 1990), as well as to the Basidiomycota mushroom Agaricus bisporus, for which 45.3% of the mitochondrial genome was composed of introns (Ferandon et al., 2013). Our analysis also revealed fungal species with no intron sequences, such as Metarhizium anisopliae and as previously described for Mycosphaerella graminicola (Torriani et al., 2008).

In fungal mitogenomes two classes of introns are found, group I and/or group II. They can be distinguished from each other by their sequence, structure and splicing mechanisms. Group I has been reported as the most widespread in fungal mitogenomes and can be subdivided into IA, IB, IC1, IC2, IC3, and ID (Saldanha et al., 1993; Hausner, 2003; Lang et al., 2007). Mitochondrial introns can self-disperse in mtDNA, acting as mobile elements due to the presence of enzymes that are able to cleave double-stranded DNA to insert introns into new genomic locations (Lazowska et al., 1980; Pellenz et al., 2002). HEGs are encoded by group I and II introns or can be found freestanding. They are classified into six families, but only two are detected in fungal mitogenomes, LAGLIDADG and GIY-YIG (Stoddard, 2014).

In our study, we identified all group I subgroups of introns, except IC3, and a few sequences of the group II. The most common intron type was IB (59.6%) and IA (17.2%), while group II introns were the least frequent, and accounted for only 2.3% of the intron sequences identified. Furthermore, 450 HE sequences, 275 LAGLIDADG and 175 GIY-YIG, were detected. Our data are corroborated by many studies, such as that published by Zubaer et al. (2018), who described the presence of 81 introns, 72 of which containing HE, in the fungus Endoconidiophora resinifera. Most of the introns detected (32) were classified as group IB, while only three introns group II were encountered. Moreover, the HEG of the LAGLIDADG family was the most frequent. In Ophiocordyceps sinensis, 52 introns were found, of which only six were classified as group II, and subgroup I was recovered as the most frequent (Li et al., 2015).

Group I introns (IGI) are considered the largest class of introns. Introns from this group are mostly self-splicing or are spliced by protein factors that stabilize the intronic RNA that undergoes conformational changes, such circularization. This type of intron can also have repetitive sequences or protein coding sequences at the ends of the intron, as seen in the case of the rps3 ribosomal protein gene (Nielsen and Johansen, 2009). In addition, some introns have been described to be involved in rRNA and tRNA folding (Cao and Woodson, 2000; Rangan et al., 2004), suggesting that introns are important regulators of the expression of genes encoding proteins and ribosomes in mitogenomes (Nielsen and Johansen, 2009). IGIs may also be important for the fitness of the fungal host. Intron group I colonized by HEGs can promote RNA splicing, due to the presence of maturase which acts as a cofactor promoting the splicing between the intron and the precursor RNA (Lambowitz and Perlman, 1990; Belfort, 2003). IGIs may have their splicing rate altered due to exposure to certain environmental factors (Lambowitz et al., 1998). For example, in chloroplasts, it has been reported that type IGI splicing was activated by light (Deshpande et al., 1997). Splicing inhibitory factors have also been described, such as the molecular Flavin mononucleotide (FMN), which directly interferes with the affinity of the intron with molecules involved in the catalysis of the intron itself (Kim and Park, 2000). Since growth conditions such as temperature, luminosity, pH and salinity have already been described influencing the rate of splicing in other organisms, it was suggested by Belfort (2017) suggested that some IGIs can work as biosensors when exposed to certain environmental factors, promoting a selective advantage over the intronless fungi. Fungal species such as budding yeast and Neurospora crassa, which do not present introns have been described presenting respiratory defects (Lambowitz and Perlman, 1990; Hausner, 2003). Additionally, type I introns have been studied to be used as ribozymes that have the ability to inactivate specific regions of the host itself or viral sequences based on the messenger RNA (mRNA) cleavage (Johansen et al., 1997; Lambowitz et al., 1998).

Our data also revealed a positive correlation between the number and length of IGIs and HE. This has also been reported for the fungus Agaricus bisporus (Ferandon et al., 2013) and E. resinifera, for which most of the introns contained HE domain, and 50% of intron size was due to the presence of HEGs (Zubaer et al., 2018).

Introns were distributed all over the mitochondrial genes evaluated, and cox1 was the main introns reservoir, followed by rrnl, cox2, cob, cox3, and nad2. According to Ferandon et al. (2013) and Zubaer et al. (2018), cox1 is the main intron reservoir in fungal mitogenomes. The authors reported the presence of 19 introns in A. bisporus and 23 introns in E. resinifera, which was considered the fungal species with the largest number of introns in this gene. Introns have also been identified in other genes, such as cox2, cob, and rrnl in two species of Rhizopogon (Agaricomycotina) (Li et al., 2019). In the fungus O. sinensis (Hypocreales), the genes cox1 and rrnl were reported as the main reservoir of introns (Li et al., 2015). Moreover, in the fungus Stemphillium lycopersici the putative cox1 and cox2 gene fusion accounted for nearly half of the number of introns encountered in the mitochondrial genome (Franco et al., 2017).

It is well known that groups I and II introns can propagate as moving elements through different mechanisms: group I by HE activity and group II by reverse protein transcription (Bonen and Vogel, 2001; Hausner, 2003; Lang et al., 2007). Since introns have the ability to disperse throughout the genome, we investigated intron sharing among genes and species. According to our data, there is a high intron sharing rate among fungal species of the order Hypocreales (267 out of 349 introns were shared among species). These findings may be an indicative of horizontal intron transfer. Horizontal intron transfer mechanisms in fungal mitochondrial DNA have already been described in many species, such as Chrysoporthe (Kanzi et al., 2016), Aspergillus and Penicillium (Joardar et al., 2012) and Armillaria (Kolesnikova et al., 2019). Moreover, besides the introns sharing in species of the order Hypocreales being more frequent among those taxonomically related, it was also observed for species belonging to different families. Mardanov et al. (2014) identified nine introns in rRNA genes from Sclerotinia borealis, but only three introns sequences presented similarity with rRNA introns in the Helotiales order. Other four sequences exhibited similarity with mitochondrial genes from species that are not related to this order. In the same study, it was described that species from this order had the same intron insertion position in the cox1 gene. In our study, we also checked for the presence of intron sharing by position in mitochondrial genes. We found that the cox1 gene, followed by the cob and rrnl genes, had the highest sharing rate among species, displaying respectively 21, 8, and 5 different insertions events. Nevertheless, we also found conserved insertion positions in the species evaluated. Joardar et al. (2012) showed that the cox1 gene for Aspergillus species presented variable insertions sites. Furthermore, Guha et al. (2018) analyzed 129 fungal species from Ascomycota phylum and found 21 different sites of intron insertion in the gene cob. These authors also described that some positions are more common for insertion sites and suggested that these locations may be preferred by HEGs activity. In addition, it was shown that introns insertion site preservation could be related to the host fitness. In Saccharomyces cervisiae, for instance, the presence of introns in nuclear and mitochondrial genes helps to resist starvation conditions (Parenteau et al., 2019).

Two IGI sequences classified as type IA located in the rrnl gene were the most shared among the fungal species in our analyses and were also detected in the genes cox1 and cob. According to Zubaer et al. (2018), some introns occurring in the rrnl gene are well conserved and widespread in many fungal species of the phylum Ascomycota. In Saccharomycetales, the presence of an IGI in the gene rrnl, called omega intron, has been studied in detail (Goddard and Burt, 1999). The authors observed a cycle of invasion and degeneration for this intron, which, in order to be maintained in mitogenome continue transposing into other genes or genomes (horizontal transfer). This invasion cycle occurs through the presence of HEGs that promote intron insertion into a mitochondrial DNA region. This mechanism of invasion cycle is mainly responsible for intron diversity in mtDNAs, as well as sequence sharing among different fungal species (Goddard and Burt, 1999; Wu and Hao, 2014). IGIs may also transpose via reverse splicing of RNA. This mechanism is based on the recombination of the intron-containing rRNA molecule that is reverse transcribed into DNA and inserted by recombination processes in other sites of the mitochondrial genome (Roman and Woodson, 1995; Hausner, 2003). This mechanism can corroborate that both introns IA detected in the gene rrnl are the most shared among species and present in other genes other than rrnl. Furthermore, some introns were found exclusively in one species (such as for Metarhizium anisopliae) or in a unique gene (atp8 and atp9), not being shared with no other species or gene.

Mitochondrial genes have numerous sequences of HEGs, which are capable of self-mobilization. Since HEGs can interrupt or structurally modify mitochondrial genes (Joardar et al., 2012; Kolesnikova et al., 2019), we evaluated the duplication of mitochondrial genes to the nuclear genome. Duplication was found a common phenomenon among the mitochondrial genes of species belonging to the order Hypocreales, since among the 17 core genes, only atp8, atp9, and cox3 were not detected in the nuclear genome. Duplication of NUMTs is a common phenomenon in other kingdoms, such as in mammals (Metazoa) (Tsuji et al., 2012). In Fungi, few studies that have verified duplications of mitochondrial genes in the nuclear genome have been described (Wright and Cummings, 1983; Specht et al., 1992; Brankovics et al., 2018; Wang et al., 2018). The species Schizopphillum commune has a duplication of the atp9 gene (Specht et al., 1992), and in Podospora anserina there are mitochondrial plasmids containing cox1 and cox3 genes during its senescence (Wright and Cummings, 1983). Furthermore, in Fusarium gramineraum, a mitochondrial pangenomic analysis revealed a 3,174 bp NUMT fragment (Brankovics et al., 2018). In our gene ordering analyses, some species, such as Metacordyceps chlamydosporia, Cordyceps militaris and Fusarium circinatum did not show any duplication of mitochondrial genes, as found and described in Hirsutella thompsonii (Wang et al., 2018).

Based on the NUMT data, we also evaluated the divergence time by molecular clock analyses using mitochondrial genes that are exclusively from the mitochondrial genome. Molecular clock is used to investigate the timing of phylogenetic events, such as dating the origin of taxonomic groups or events of gene duplication, diversification or gene loss (Weir and Schluter, 2008). Herein, we estimated that the order Hypocreales probably originated around 137.39 Mya. Compared to an estimation based on 638 protein sequences in the nuclear genome, which indicates that Hypocreales originated approximately at 193 Mya (Kubicek et al., 2019). It is well known that the evolutionary rate of fungal mtDNA is close to that of plants, which have the lowest nucleotide substitution rates (Aguileta et al., 2014; Sandor et al., 2018). The evolution rate is based on the percentage of base substitutions in conserved genetic regions. Clark-Walker (1991) demonstrated that the percentage of substitution in mitochondrial genes was lower than that of nuclear genes, suggesting that fungal mitogenomes evolve slower than their respective nuclear genomes. Therefore, the time found in our analyses is expected, since mitochondrial genomes evolve slower than nuclear genomes.

Also, based in the value of d⁢S¯ of each species, it was possible to estimate the relationship among the genome size and the neutral substitution rates. The more fast evolving a mitogenome is, the larger its size and more abundant in non-coding regions, such as introns, HEGs and uORFs, while the protein coding-region are mostly unchanged. This indicates that the protein-coding genes of mitogenomes are under the influence of negative selection, since the coding regions remains stable independently of the neutral substitutions rates.

Fungal mitochondrial genomes usually contain 14 core protein-coding genes involved in the electron transport chain ATP synthase complex (Kang X. et al., 2017). In addition to these genes, they also have two ribosomal genes and rps3 protein that are also conserved. In our analysis, most of the species presented the 14 core genes widely described for fungal mitogenomes, the two ribosomal genes and rps3 protein.

The ordering of the mitochondrial genes of the evaluated species of the order Hypocreales is remarkably conserved. Nonetheless, the species Acremonium fuci, A. chrysogenum and Clonostachys rosea presented the gene cox2 displaced, while Nectria cinnabarina and Epichloe typhina displayed a rearrangement for the nad4l gene. Differences in genome ordering can be attributed due to non-homologous recombination processes, plasmid sequence integration or presence of transposable elements, such as IGI sequences and HEGs (Gordenin et al., 1993; Lavrov et al., 2002; Rocha, 2003; Repar and Warnecke, 2017).

Species of the Beauveria genera exhibited an additional copy of the gene rps3 as a freestanding gene. The rps3 gene is the only ribosomal protein encoded by fungal mtDNA and may be located as a freestanding gene or inserted within an intron of the rrnl gene (Sethuraman et al., 2009). Gene duplication may occur due to truncated sequences or loss of function. In the mitogenomes of Beauveria evaluated, rps3 copies were almost the same size. Nevertheless, only freestanding copies had coding potential, suggesting that the rps3 copy located in the rrnl gene could be a pseudogene.

In the current study we provided the complete sequencing of the mitogenome of the fungus Trichoderma harzianum, a commercially important fungal strain. Comparative analyses of 35 species of the order Hypocreales revealed a structural dynamic in the mitochondrial genome of this well-studied and diverse order of the Fungi kingdom. The genome size variability found was mainly correlated with the presence of non-coding elements, such as introns and HEGs, which could be considered the determinant elements for the shape of mitochondrial genomes, as described in other Orders. Some studies, such as the one described by Mardanov et al. (2014) verified the difference in mitogenomes size and presence of introns in the orders Peltigetales and Helotiales and showed a dynamic of acquisition and loss of introns during evolution. In our work, we estimated that mitogenomes that evolve faster, have longer length of mitogenome and non-coding region.

The in-house scripts generated for this study can be found in the Supplementary Text S1 and are available in GitHub repository (https://github.com/paulaluize/mitogenomes). The original contributions presented in the study are publicly available. This data can be found here: https://www.ncbi.nlm.nih.gov/nucleotide/; accession number MT263519.

PF, BB, VA, EA, and AG-N conceived and designed experiments. PF, RD-P, DA, DB, EA, FB, L-ED-B, and AG-N analyzed the data. PF, FB, L-ED-B, EA, and AG-N wrote the manuscript. All authors read and approved the final manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.