Many studies have shown that fungal mitogenomes have great structural plasticity (Fonseca et al., 2020; Araújo et al., 2021; Kulik et al., 2021). In this study, the length difference among mitogenomes ranges from 12 (R. allomycis Chytridiomycota – NC_003061) to 272 kbp (Ophiocordyceps camponoti Ascomycota – CM022976). Since the two mitogenomes are from different and distant phyla, we assessed whether mitogenome length would be related to phyla classification (Figure 2A). In general, Ascomycota (63,888 ± 35,298 bp) and Basidiomycota (77,394 ± 49,928 bp) mitogenomes have quite variable lengths. However, Cryptomycota (25,401 bp) and Mucoromycota (61,936 ± 14,895 bp) have smaller and more homogeneous sizes. This result may be related to the availability of data, since more than 90% are composed of mitogenomes from Ascomycota and Basidiomycota, while Mucoromycota and Cryptomycota represent 3.38% and 0.13% of mitogenomes, respectively.

As the length of mitogenomes is variable and is related to variations in both protein-coding and accessory element regions (Paquin et al., 1997; Bullerwell, 2003; Al-Reedy et al., 2012; Araújo et al., 2021), we re-annotated all mitogenomes to facilitate the comparison between the features present in each genome. After the annotation, we estimated the presence/absence of mitochondrial genes that are considered essential in fungal mitogenomes. Several previous studies have already demonstrated the presence of atp, cox, cob, nad, rns, and rnl in fungal mitogenomes (Song et al., 2020; Kulik et al., 2021). According to our analyses, although most of the mitogenomes have all the previously expected core genes, there was not a single gene shared in all evaluated fungal species. Figure 2B displays the frequency of each gene per phylum. In this figure, the phylum Cryptomycota is not shown, as it has only one representative (Paramicrosporidium saccamoebae – CM008827) containing all genes except for nad6. Genes of atp, cox, and cob families are the most frequent among fungal mitogenomes. Homologs of nad, rns and rnl are the least found within Ascomycota. It is worth noting that 149 mitogenomes are from the genus Saccharomyces, which does not have nad homologs (Ruan et al., 2017). In addition, some S. cerevisiae mitogenomes also have no homologs of atp and cob (Figure 3). Many studies have already reported that the coding content of the genus Saccharomyces is well conserved (Freel et al., 2015; Ruan et al., 2017). No member of Saccharomycetaceae were shown to have nad homologs, and studies suggested that nad homologs were lost after the divergence of this family (Dujon, 2010).

As a controversy of the relationship between size and presence of core genes, the mitogenome of the yeast Malassezia furfur (CP046241) is 43 kb long but does not show any of the core genes. This mitogenome has 47 ORFs, of which 22 have unknown protein domains, and the others have domains from the superfamilies haloacid dehydrogenase (HAD), Glutenin high molecular weight, among others (Supplementary File 2). The species M. furfur is a basidiomycotan yeast of medical importance as it can cause skin and blood infections (Theelen et al., 2018), and the mitogenome used in this study was previously used for taxonomic classification of the species (Sankaranarayanan et al., 2020). The genes found in this mitogenome, despite having mitochondrial-related functions, are not commonly found in mitogenomes, but in nuclear genomes (Theelen et al., 2018). Another species lacking many mitochondrial genes is Rozella allomyces which has the smallest mitogenome size in our dataset (12 kbp), and lacks nine core genes (atp6, nad1-6, and rnl).

During evolution, many mitochondrial genes migrated to the nuclear genome. In fungi, the majority of genes related to mitochondrial function are found in the nuclear genome (Bolender et al., 2008). The main theory for the mechanism involved in the escape of DNA to the cytosol and, consequently, to the nucleus states that these genetic fragments are transferred by mobile elements that use the non-homologous (NHEJ) machinery for integration (Berg and Kurland, 2000; Tsuji et al., 2012). In S. cerevisiae, transfer of mtDNA to the nuclear genome has already been observed when the cell had mutations in specific nuclear genes, depending on the structure of the mitogenome and availability of sugars to allow fermentation in the medium, which could lead to the high variation of phenotypes to anaerobic environments (Shafer et al., 1999; Peter et al., 2018). The absence of core genes in the species evaluated may be an indication that these genes migrated to the nuclear genome. Studies exploring the nuclear genome of these species must be carried out to confirm the presence of genes of mitochondrial origin in the nuclear genome.

After assessing the number and length of the following accessory elements: introns, HEGs, and uORFs, in all the fungal mitogenomes filtered by CD-HIT, we characterized the number of HEGs found per genome and grouped by phylum. Many studies have been discussing the presence and influence of non-coding elements in mitogenomes of different fungal groups (Fonseca et al., 2020; Araújo et al., 2021; Kulik et al., 2021). HEGs are classified as selfish elements that can transpose by breaking the DNA at specific sites, which consequently generates gene conversion events. These accessory elements can be found in introns or as free-standing, which ensures their dissemination and fixation in the mitogenome and may even alter the reading frame of protein-coding genes (Hausner, 2003; Barzel et al., 2011; Stoddard, 2014).

In our analysis, we identified Ascomycota species that can have up to 100 HEG elements, such as Nemania diffusa (NC_049077). The other phyla exhibited a more homogeneous number of HEGs, with a mean of six HEGs per mitogenome (Figure 4A). HEGs have not been identified in many mitogenomes, such as Acremonium fuci (NC_029851) and Hyaloraphidium curvatum (NC_003048). Currently, six families of HEGs are known, with only two families found in fungal mitogenomes, LAGLIDADG and GIY-YIG (Hausner, 2003; Megarioti and Kouvelis, 2020). We identified a total of 4,432 LAGLIDADG and 1,247 GIY-YIG across all mitogenomes. Interestingly, we found 15 H-N-H elements, showing that these HEGs occur in the mitogenome of fungi, although they are rare (Hafez and Hausner, 2012; Megarioti and Kouvelis, 2020). H-N-H are usually found as free-standing or in group I introns (IGI) in Bacteria and their phages (Hafez and Hausner, 2012). All three aforementioned types of HEGs are able to cleave double-stranded DNA at specific sites, facilitating the mobility of the elements within the genome (Arbuthnot, 2015). Since numerous HEG elements were found in most mitogenomes of fungi, we asked whether there is a correlation between the mitogenome total length, and the total length of the fraction composed by HEGs. We found a significant positive Pearson correlation (R = 0.37; p < 2.2e-16) suggesting that the abundance of HEGs length is one of the factors responsible for the observed variation in fungal mitogenome size (Figure 4B). However, our data indicate that they may not be the only factor, leading us to investigate whether other non-coding elements would also be responsible for the variation in size and diversity in the composition of mitogenomes.

In total, 5,516 intronic sequences were identified. Ascomycota and Basidiomycota had 3,976 (72.1%) and 1,191 (21.59%) of all introns identified, respectively. Moreover, Ascomycota mitogenomes have an average of nine introns per mitogenome while the other phyla averaged less than five (Figure 5A). Mitogenomes with a high number of introns were identified, such as Pyrrhoderma noxium (Basidiomycota – CM008263) with 39 introns, Phycomyces blakesleeanus (Ascomycota- NC_027730) with 38 introns, Parasitella parasítica (Mucoromycota – NC_024944) with 22 introns, and Rhizophydium sp. (Chytridiomycota – NC_003053) with 17 introns. We found that the presence and global length of introns is correlated with the size of each mitogenome (R = 0.48, p < 2.2e-16) (Figure 5B). Many introns contain domains of HEGs, and since these elements can transpose to other parts of the genome, we investigated whether the length of introns was related to the HEGs length (Figure 5C), and found a correlation R = 0.86, which suggests that HEGs may be partly responsible for intron size, such as proposed previously (Edgell, 2009).

We assessed the percentage of HEGs within introns as well as the opposite (Supplementary Table 4). We found that 43.46% of all HEGs are within intronic regions (31.54% of GIY-YIG; 33.33% of HNH; and 46.99% of LAGLIDADG), and 48.23% of the analyzed introns overlap with HEGs coordinates, and, most frequently, LAGLIDADG (9.00% overlap GIY-YIG; 0.22% overlap H-N-H; and 39.01% overlap LAGLIDADG). In fact, these results agree with the literature, which shows that 30% of IGI contain internal ORFs that might represent HEGs (Chevalier and Stoddard, 2001; Hafez and Hausner, 2012). Most probably, introns originated before HEGs as first suggested by Loizos et al. (1994), providing a “safe haven” for these endonucleases to get established, although this theory does not explain why there is still many free-standing HEGs (Loizos et al., 1994). Edgell (2009) states that both HEGs and introns are selfish, independent elements, but they often benefit from being in close genomic proximity (e.g., HEGs are established by introns, and introns can use the endonuclease properties of HEGs to spread more easily) (Edgell, 2009). The association of HEGs with introns allows their persistence in the genome, since if they inserted themselves into genes, they would be subject to negative selective pressure to eliminate the selfish elements from the genome (Stoddard, 2014).

Generally, IGI is more common in fungi while IGII is more common in plants and most of them are self-splicing (Hausner et al., 2014; Fonseca et al., 2020). In this study, we performed the classification of identified introns and found six types (IA, IB, IC1, IC2, ID, and II). Introns from group IB were the most common among mitogenomes (562 or 76% of mitogenomes), while group IC1 was the least frequent (31 or 4.19% of mitogenomes). Moreover, when IGI is colonized by HEG domains, it can split a gene region due to the HEG cutting function (Lambowitz and Perlman, 1990; Belfort, 2003). Although IGII are usually infrequent, they were identified in 238 (32.21%) mitogenomes (Figure 5D). IGII introns have catalytic (ribozyme) and intron-encoded-protein sites, which allow their own splicing and proliferation in the genome (Lambowitz and Zimmerly, 2011). A previous study carried out by our group showed that intronic sequences can be shared among distinct species and at different genomic positions, indicating that the presence of self-splicing introns associated with maturases, such as HEGs, may be responsible for the transfer and sharing between genomes (Fonseca et al., 2020).

When we evaluated the presence of uORFs in mitogenomes, species of Ascomycota had the highest number of uORFs. For example, Morchella importuna (NC_012621) had 102 uORF elements and a mitogenome length of 107 kbp (Figure 6A). Mitogenome and uORF lengths present a weak positive linear relationship (R = 0.28, p = 1.2 e-14) (Figure 6B), and a moderate positive linear correlation with introns length (Figure 6C). These uORF elements have been identified in different fungal mitogenomes, and some may be shared among species, as recently demonstrated (Fonseca et al., 2020; Araújo et al., 2021; Kulik et al., 2021), which suggests that uORFs within mitogenomes of fungi may play a significant role in shaping mitochondrial genome length. This is in agreement with a recent investigation of mitochondrial landscape in eukaryotes, which has found a positive correlation between mitochondrial genome length and the prevalence of unassigned regions (Formaggioni et al., 2021). Nonetheless, differences in the literature suggest that there are probably inter-phyla variations. For instance, the linear correlation coefficient is greater than 0.50 in both the Agaricomycetes class and the Hypocreales order (Fonseca et al., 2020; Araújo et al., 2021), a higher value than the one obtained in this study. As aforementioned, uORFs are not homogeneously present in all phyla analyzed, and therefore they may not influence mitogenome length at the same magnitude for different groups.

In our attempt to identify uORFs, we found that 64.49% of them have known protein domains (Supplementary Table 5) and most of them were sequences with HEGs domains (69.05%). We also identified some uORFs with protein domains related to plasmid genes: dpo (5.60%) and rpo (1.05%). A total of 35.51% of the uORFs lacked any similarity with known protein domains and may represent unidentified protein-coding genes. Nevertheless, by implementing a computational technique that assess coding potential based on sequence intrinsic features, 43.86% of uORFs were predicted to be coding (Supplementary Table 6). Since we were working with a kingdom-wide dataset, it was out of the scope of our study to conduct genus- or species-level transcriptomic analyses to verify if those uORFs are transcribed. When possible, this method can help shed light on whether uORFs are in fact transcribed or not, as demonstrated by Torriani et al. (2014) for some Rhynchosporium species (Torriani et al., 2014).

Mitochondrial plasmids are widely found in Fungi and some mitochondrial genes, such as dpo and rpo, are commonly associated with plasmid integration events (Formighieri et al., 2008; Férandon et al., 2013; Formaggioni et al., 2021). In our study, we identified dpo homologs in at least one species of six different phyla, with Basidiomycota and Ascomycota showing the highest numbers of sequences – 227 and 116, respectively (Supplementary Table 5). On the other hand, rpo genes were only annotated in Ascomycota (46 sequences), Basidiomycota (23), and Chytridiomycota (6). In total, amongst all mitogenomes analyzed, it was possible to observe that 10.74% of Ascomycota fungi presented a dpo and/or rpo homolog, in contrast with 57.14% of Basidiomycota species. For the other phyla, the lack of representatives makes difficult to draw meaningful conclusions. For instance, most Mucoromycota species in our dataset, have a dpo homolog; however, considering that we are working with only 24 Mucoromycota species, it is not a fair comparison with the most represented phyla Ascomycota and Basidiomycota.

When analyzing similarity among annotated dpo sequences (Figure 7), some form phyla-specific clades, whereas others seem to be shared among species of different phyla. A similar pattern was observed for rpo homologs (Figure 8). One possible explanation for this finding is horizontal gene transfer (HGT). As discussed by Beaudet et al. (2013), mitochondrial plasmids of fungi and plants exhibit many similarities, indicating that HGT events between organisms of those kingdoms may occur (Beaudet et al., 2013). According to Handa (2008) a plasmid from Brassica has an ORF with similarity to rpo from several organisms, including fungi (Handa, 2008). Moreover, since different species of fungi often share the same habitat, the hypothesis that HGT events between fungi occur is often discussed by many studies that identified highly similar dpo and/or rpo homologs among phylogenetically distant fungi (Wang et al., 2008; Xavier et al., 2012; Andrade and Goes-Neto, 2015).

Although most mitogenomes are circular (Kulik et al., 2021), some genomes were annotated as linear. We wondered whether the conformation of the genome would affect its length, so we compared these values using Wilcoxon’s statistics and found no statistical difference (Figure 9A). Even though, the mitogenomes conformation had no effect in mitogenome length, fungal linear mitogenomes have distinct features when compared with circular mitogenomes, such as the presence of an invertron (terminal inverted repeats) and of one to six uORFs that may produce DNA or RNA polymerases (Handa, 2008). In our analysis, no invertron were detected in mitogenomes. The presence of linear mitogenomes has already been described in different eukaryotes, including fungi. One of the hypotheses is that linearization occurs due to the integration of plasmids containing genes that code for dpo (Williamson, 2002; van de Vossenberg et al., 2018). In our analysis, 37 (5.22%) linear genomes were identified and seven of them have dpo homologs. Examples are Candida subhashii and Saccharomycopsi malanga, with two dpo each. Despite the presence of these genes in some linear mitogenomes, the vast majority did not present dpo, suggesting that other mechanisms may be involved in mitogenome linearization.

Another little explored feature is the usage of different genetic codes in fungal mitogenomes. Among the 685 mitogenomes evaluated, the use of four different genetic codes was detected: code 1 (The Standard Code), code 3 (The Yeast Mitochondrial Code), code 4 (The Mold, Protozoan, and Coelenterate Mitochondrial Code and the Mycoplasma/Spiroplasma Code) and code 16 (Chlorophycean Mitochondrial Code)⁶. Overall, 62% of all mitogenomes use genetic code 4 (the UGA codon codes for Tryptophan), belonging to Ascomycota, Basidiomycota, Monoblepharidomycota, Cryptomycota, and Blastocladiomycota (Figure 9B). Most species of Mucoromycota are coded by code 1, while many species of Chytridiomycota are coded by code 16 (the TAG codon codes for Leucine) (Figure 9B and Supplementary Table 7). These results defy the classical ‘Frozen Accident Theory’ that proposes the immutability of the standard genetic code to avoid lethal genetic alterations (Crick, 1968). According to our genetic table analysis, most fungal species likely use the UGA codon (stop codon in standard genetic code 1) for translating Trp, similar to what was previously reported for UAA (Knight et al., 2001). Nonetheless, there is evidence that, in fact, a minority of fungal mitochondria uses UGA to code Tryptophan, showing that the annotation of fungal mitogenomes needs to be widely updated (Nibert, 2017). There are multiple theories to explain the differences of genetic code usage, but the one that better fits the fungal scenario may be related to the simplification of the mitochondrial genome and translation machinery (‘genome streamlining’ hypothesis), that would make the replication of mitogenomes easier, less energy costly and would also slowly simplify the mitogenome length (Andersson and Kurland, 1995; Knight et al., 2001; Swire et al., 2005).

We then compared the length of mitogenomes with different genetic codes (Figure 9C). Our results showed that fungi with genetic code 3 have larger mitogenomes than fungi with other genetic codes. According to NCBI Taxonomy⁷, the genetic code 3 classically represents yeasts. These unicellular fungi do not necessarily need mitogenomes for energy production, and they might have more accessory elements or non-coding regions, which contributes to the genomic size increase (Ruan et al., 2017). Many studies have demonstrated that the length of mitogenomes is generally affected by the presence of HEGs (Fonseca et al., 2020; Araújo et al., 2021), and we detected the same trend in our data when comparing genomes with genetic code 3 and 4 (Figure 9D).

We also observed a difference between the GC content of mitogenomes at phylum level (Figure 9E). Usually, fungal mitogenomes have a low GC content, and in our dataset the mean GC content was ∼25%. Ascomycota, Mucoromycota and Basidiomycota showed the highest variation. Some outliers can also be observed regarding GC content, such as M. furfur (Basidiomycota – CP046241) with 68% GC, some species of the genus Candida ranging from 13–53% (Ascomycota – NC_014337 and NC_022174), and the species Glomus cerebriforme with 47% GC (Mucoromycota – NC_022144). Some studies suggest that low GC content restricts mutations in mitochondrial genes (Grivell, 1989; Lang et al., 1999). Thus, mitogenomes that have a high GC content may be subject to higher mutational variations.