Fungal isolates were isolated from TBD samples collected across Ghana as reported previously (Amoako-Attah et al., 2020). Based on that previous study (Amoako-Attah et al., 2020), one strain each from M. crinis-equi, M. palmivorus and Mar. scandens, and seven isolates from M. tenuissimus (Table 1) were chosen for whole genome Illumina sequencing and their complete mitogenomes were assembled for this study. Fungal biomass generation and DNA extraction was performed as reported previously (Amoako-Attah et al., 2020).

Genomic DNA of the above mentioned 10 isolates were sequenced using Illumina X-ten paired-end short-read technology (library preparation and sequencing done by Beijing Genome Institute, Shenzhen, China) as reported previously (Ali et al., 2020). The short reads were assembled using ‘‘SOAPdenovo’’ (¹ version: 2.01) as described previously (Ali et al., 2017). The mitogenomes were identified from the assembled contigs based on the average GC content, size and depth of coverage of the contigs. As the fungal mitogenomes are normally AT-rich (Hausner, 2003), we targeted the contigs with < 35% GC content, high coverage and >1,000 bp long, and these were BLASTn (Altschul et al., 1990) against the Mon. roreri mitogenome (Costa et al., 2012). Only one contig per isolate showed similarity on the nucleotide level to the Mon. roreri mitogenome and showed the presence of common mitochondrial genes. To confirm the circularity of the mitogenomes, contigs were searched for sequence repeats on each end and trimmed leaving only the unique sequence. The mitogenomes from all the 10 isolates satisfied that criteria.

Mitochondrial gene annotation of the 10 isolates and the publicly available Mon. roreri (GenBank no. HQ259115.1) and Mon. perniciosa (GenBank no. AY376688.1) were performed with MFannot (Valach et al., 2014) using the NCBI translation code 4 (The Mold, Protozoan, and Coelenterate Mitochondrial Code). The tRNAs and introns were re-confirmed by RNAweasel². To further confirm genes boundaries as well as intron-exon boundaries, MFannot predictions were again checked individually by aligning against their orthologous in closely related fungal species. A physical map of the mitogenome was created with OrganellarGenome-DRAW (OGDRAW) v 1.2 (Lohse et al., 2007). Functional annotation of the predicted open reading frames (ORFs) was complemented with Blast2GO Basic (Conesa et al., 2005). The mitogenomes were deposited at NCBI GenBank under the accession number MZ615343- MZ615352.

To evaluate the application of mitogenomes for fungal phylogeny, a phylogenetic tree was constructed using Amino acid sequences of the 14 conserved mitochondrial proteins of these species, including three ATP synthase subunits (atp6, atp8, and atp9), four Cytochrome c oxidase subunits (cox1, cox2, cox3, and cob), and seven NADH dehydrogenase subunits (nad1, nad2, nad3, nad4, nad4L, nad5, and nad6). Accessions of completely sequenced mitogenomes of four species in Marasmiaceae family and five species from other families of the order Agaricales were retrieved from NCBI Organelle Genome Resources website³. In addition, the mitogenomes of Rhizoctonia solani was used as outgroups. Protein sequences were combined and aligned using ClustalW2 tool (Larkin et al., 2007) under default settings and a phylogenetic tree was reconstructed using the Maximum Likelihood method based on the JTT matrix-based model using MEGA v. 6 (Tamura et al., 2013). The best fit amino acid substitution model with lowest BIC score was determined using MEGA v. 6. Bootstrap values were computed with 1000 resampling iterations using an approximate likelihood ratio test. Similarly, phylogenetic trees were also constructed using Amino acid sequences of the rps3 gene, cox1 introns and the internal and external ORF’s of Marasmiaceae species listed in Table 1. The best fit model used, as determined using MEGA v. 6, for these phylogenetic trees were as follows: rps3 of three M. tenuissimus isolates with the other five species- WAG; rps3 of seven isolates of M. tenuissimus- JTT; cox1 introns of two M. tenuissimus isolates with the other five species- GTR; cox1 introns of seven isolates of M. tenuissimus- T92; and the internal and external ORF’s of all the six Marasmiaceae species- WAG.

The complete mitogenomes of Marasmiaceae species listed in Table 1 were used for reciprocal BLASTn searches to identify regions of similarity, insertions, and rearrangements. First, to identify major rearrangements, structural comparisons between the three types of mitogenomes based on gene order (M. tenuissimus, Mar. scandens, and M. palmivorus) were generated with Circoletto (Darzentas, 2010), combining BLASTn (threshold of 1e ^––10) searches with the Circos output. Second, to identify regions of similarity and insertions, complete mitogenome alignments were carried out using the mVISTA program⁴. Default parameters were utilized to align the genomes in Shuffle-LAGAN mode and a sequence conservation profile was visualized in an mVISTA plot (Frazer et al., 2004).

To show the distribution of the GC content along the mitogenomes, GC-plots were created using DNAplotter (Release 18.1.0) with a window size of 10,000 bp, and step size of 200 bp. The pink represents below average and the green represents above average GC content for each mitogenomes.

Tandem repeats were identified using the Tandem Repeats Finder (Benson, 1999) under default parameters.

The mitogenomes of M. crinis-equi (strain: GHA76), M. palmivorus (strain: GHA12), Mar. scandens (strain: GHA19), and seven isolates of M. tenuissimus were sequenced, assembled, and annotated in this study (Table 1 and Supplementary Excel File 1). The Mon. roreri (GenBank no. HQ259115.1) and Mon. perniciosa (GenBank no. AY376688.1) mitogenomes were originally annotated manually using BLAST searches and ORF finder software (Formighieri et al., 2008; Costa et al., 2012). For comparative genomics, it is important to use the same annotation technique. Therefore, the Mon. roreri and Mon. perniciosa mitogenomes were re-annotated using the MFannot tool (Valach et al., 2014). Both the current and previous methods resulted in the same number and orientation of the primary mitochondrial genes and similar numbers of intronic and accessory ORFs (Supplementary Figure 1). Mitogenomes are known for their wide range of sizes despite their seemingly similar functional gene composition, even among closely related species (Zhang et al., 2020). There is a remarkable interspecies variation in mitogenome size within the six cacao-associated Marasmiaceae species so far sequenced, ranging from 43,121 to 51,210 bp for M. tenuissimus, to 63,005 bp (M. crinis-equi), 79,017 bp (M. palmivorus), and 86,201 bp (Mar. scandens), as compared to 93,722 bp for Mon. roreri and 109,103 bp for Mon. perniciosa (Table 1). A portion of Mon. perniciosa’s expanded genome size can be attributed to the stable insertion of a linear plasmid (Formighieri et al., 2008) between nad1 and cob in association with other rearrangements (Figure 1).

All the mitogenome sequences (M. crinis-equi, M. palmivorus, Mar. scandens, seven isolates of M. tenuissimus, Mon. roreri, and Mon. perniciosa) carry the 14 core genes (Bullerwell and Lang, 2005) involved in oxidative phosphorylation, ATP synthesis and mitochondrial protein synthesis (Figure 1 and Supplementary Excel File 1). In addition, these genomes carry genes encoding the small and large subunits of rRNA. All the rRNA and core genes were in the same orientation, except for atp8. Though it is not common in fungi for atp8 be in the opposite orientation, it seems to be a common phenomenon in the Marasmiaceae. Beside the above mentioned six Marasmiaceae species, Lentinula edodes and Omphalotus japonicus maintained the same orientation (Supplementary Figure 2). But it is not a common phenomenon for the order Agaricales or class Agaricomycetes (Supplementary Figure 2). All the genomes also carry a complete set of tRNAs and the rps3 gene (Supplementary Excel File 1). The tRNAs are generally distributed in small blocks around the genomes, the exception being Mar. scandens which has 17 of 25 tRNA copies located in a single block between the rns and cox1 genes (Figure 1). Further analysis comparing the mitogenome structure of Mar. scandens with that of M. tenuissimus and M. palmivorus indicates that the large tRNA block of Mar. scandens has resulted from recombination of at least four small tRNA blocks (Supplementary Figure 3), suggesting a recombination hotspot in the Mar. scandens mitogenome. The gene order for the core genes within the mitogenomes for the six species varied and show distinct relationships depending on the species compared. M. palmivorus has a conserved gene order aligning with Mon. roreri and Mon. perniciosa for both core genes and tRNAs (Figure 1 and Supplementary Figure 4). Similarly, M. crinis-equi and M. tenuissimus have a conserved gene order with each other but not with the other species being studied (Figure 1 and Supplementary Figure 5). When compared to the other species being studied, the Mar. scandens mitogenome sequence is unique (Figure 1). Comparing these three patterns of gene order, multiple gene block rearrangements are evident (Figure 2). All the seven M. tenuissimus isolates showed a conserved gene order despite differences in their mitogenome sizes and ORF composition (Figure 3 and Supplementary Figure 6). Though the mechanisms of gene rearrangement in fungal mitogenomes are not fully understood, these gene arrangements provide a large amount of information for understanding the evolutionary status of species (Wu et al., 2021).

A phylogenetic analysis was carried out incorporating the 14 concatenated core proteins for each of the 10 mitogenomes sequenced in this study along with additional Marasmiaceae and Agaricales mitogenomes obtained from the public domain (Figure 4). The resulting phylogenetic tree separated each species within the Marasmiaceae family with bootstrap support of more than 98%. The tree topology shows three distinct clades (each with internal nodes of 99–100% bootstrap support) with M. palmivorus, Mon. roreri, and Mon. perniciosa forming one clade, M. crinis-equi and M. tenuissimus forming another clade and Mar. scandens grouping with L. edodes and O. japonicus to form a third clade. Notably, similar grouping was also observed based on the conservation of gene order for the first two clades. But there was no conserved gene order between Mar. scandens, L. edodes and O. japonicus (Figure 1 and Supplementary Figure 2). The presence a of large tRNA block in both Mar. scandens and O. japonicus further supports the presence of a recombination hotspot possibly contributing to the non-conserved gene order in this clade. These phylogenetic findings were consistent with earlier phylogenetic observations based on ITS sequence of a large group of Marasmiaceae fungi from around the world (Amoako-Attah et al., 2020). In addition to clearly showing interspecies variation, the phylogenetic tree also showed some intraspecies variation within the seven isolates of M. tenuissimus with bootstrap support of more than 78% (Figure 4). Four of the five M. tenuissimus isolates, each carrying a single intron in the cox1 gene (Figure 3), showed 100% sequence similarity (Figure 4). M. tenuissimus isolate GHA64, also carrying a single cox1 intron and having the smallest genome, was separated from the four isolates with bootstrap support of 79% (Figure 4) and sequence similarity of 99.75% suggesting this variation was not critically linked to variation in cox1 intron composition. M. tenuissimus isolates GHA79 and GHA63, both carrying three cox1 introns, were not clonal (Figure 4). Amoako-Attah et al. (2020) also reported different morphotypes within M. tenuissimus isolates. GHA79 was reported as morphotype B and GHA63 was reported as morphotype C, further indicating the isolates are not clonal. Whether these three different phylogenetic groups within the M. tenuissimus are different cryptic species, need further investigation. The whole genome sequence analysis of these isolates which is ongoing in our lab should shed more light on this subject. The rRNA ITS sequence is the only molecular marker so far used in Marasmiaceae phylogeny (Koch et al., 2018; Amoako-Attah et al., 2020). The mitogenome sequence provide more genetic information than the ITS sequence alone, and is a better choice for analyzing the origin and phylogeny of this lesser known family of Basidiomycetes.

Despite similarities in primary gene compliments and, in specific comparisons, gene order associations, the genome lengths vary significantly both between closely related species and, for M. tenuissimus, within species. These differences in genome lengths are primarily due to the number and lengths of introns, differences in intergenic space, and differences in the size and numbers of unidentified ORFs (uORF). The numbers of introns found in the 12 genomes compared ranged from 1 to 13 with total lengths varying from 1,054 to 16,979 bp. Mon. roreri has four protein coding genes (cox1, nad5, nad1, and cob) with introns. Mon. perniciosa has 6 (cox1, nad5, nad4, nad1, cob, and cox2), Mar. scandens, M. crinis-equi and M. palmivorus have 2 (cox1, cob), and M. tenuissimus has 1 (cox1) protein coding gene(s) with introns. It is typical in fungi that cox1 carries the most intronic sequences and cob carries the second most for a specific species (Megarioti and Kouvelis, 2020). Mon. roreri has 5 cox1 introns, Mon. perniciosa 6 cox1 introns, Mar. scandens 4 cox1 introns, M. crinis-equi 5 cox1 introns and M. palmivorus 4 cox1 introns (Figure 1). M. tenuissimus isolates GHA63 and GHA79 carry 3 introns within the cox1 gene while the remaining isolates carry only one intron within the cox1 gene (Figure 3). Another discrepancy observed is the nad4 gene of M. tenuissimus isolates GHA63 which is almost double in size compared to the other isolates and carries no intron (Figure 3 and Supplementary Excel File 1). The intronic regions show a high degree of sequence variation compared to the conserved sequences of the 14 core genes (Supplementary Figures 4, 5). Variation in intron numbers and length typically contributes to the size variation of the fungal mitogenome (Hamari et al., 2002).

The total intergenic space found in the 12 genomes being compared range from 16,610 to 36,404 bp. Isolate GHA63, the largest M. tenuissimus mitogenome, has the smallest intergenic space (16,610 bp), while Mar. scandens has the largest intergenic space (36,404 bp) (Table 1). TRs are involved in expansion of intergenic space and have been associated with gene order changes (Aguileta et al., 2014). The number and types of TRs detected also vary with these species (Supplementary Excel File 1). Mon. perniciosa is the extreme with 60 TRs detected. Mon. roreri has 24 TRs while the closely related M. palmivorus has only 11 TRs, the lowest for any mitogenome studied here. Costa et al. (2012) also reported high repeat sequence in Mon. perniciosa compared to Mon. roreri. M. tenuissimus genomes that carry a single intron in cox1 have 16 to 17 TRs while the 2 genomes carrying 3 cox1 intron have 18 and 21 TRs. M. crinis-equi carries 24 TRs, including 4 within cox1 introns (Supplementary Excel File 1). Although the types of TRs vary between M. tenuissimus isolates GHA63 and 79 and M. crinis-equi, these isolates share a pattern of increased TR numbers starting after rps3 continuing to the following trnR. This pattern is less developed in the remaining M. tenuissimus isolates. Comparing the TRs between rps3 and trnR region of the seven M. tenuissimus isolates shows that the overall number of TR regions for some of the isolates are increasing and the number of repeats for one TR region is also increasing (Supplementary Excel File 2). But at the same time, GHA79 has gained unique TR regions and lost a common TR region (Supplementary Excel File 2). Comparative analysis of the mitogenomes from three isolates of L. edodes, another Agaricales also revealed variable number TR regions (Kim et al., 2019). The authors proposed that the elongation of the repeat regions occurs through reciprocal incorporation of basic repeat units. Variations in TR number are thought to be generated by slipped-strand mispairing in the mitogenome (Mundy and Helbig, 2004) and/or genetic recombination (Nishizawa et al., 2000; Wang et al., 2015). In animals and plants, the TRs are responsible for mitogenome stability by insertion and deletion of the repeat array (Lobachev et al., 1998; Albert and Sellem, 2002). In fungi, repeat rearrangement occurs mostly through genetic recombination and contributes to mitogenome evolution (Aguileta et al., 2014). Whether selection acts upon the number and size of the repeat region in the mitogenome or results from some stochastic events related to erroneous DNA replication and repair needs further investigation.

As might be expected (Zardoya, 2020), the presence of uORF outside the introns also contributes to the differences in mitogenome sizes among the species studied. The large genome of Mon. perniciosa carries 40 accessory genes, many of which are part of an inserted linear plasmid (Formighieri et al., 2008), Mon. roreri carries 24 uORFs, M. scandens 7 uORFs, and M. crinis-equi 2 uORFs, M. palmivorus 14 uORFs, and M. tenuissimus 3 to 4 uORFs depending on the isolate. Most of the uORFs show extreme sequence variation within and between species (Supplementary Figure 7 and Supplementary Excel File 3). Most exceptions where similarities are significant occur between Mon. perniciosa, Mon. roreri, and M. palmivorus.

The average GC content of the 12 mitogenomes range between 26.58 and 32.04% (Table 1) which is in the higher range of values for fungal mitogenomes (average = 24.4 ± 7.3%, estimated from all fungal mitogenomes deposited in the NCBI organelle genome database) (Franco et al., 2017). Moreover, distribution of the GC content along the mitogenomes varied greatly among the species studied here in. The three closely related species, Mon. perniciosa, Mon. roreri, and M. palmivorus, have distinct patterns of GC content but do share some similarities (Supplementary Figure 8). All three species show above average GC content in association with the rps3 and rrnL. M. palmivorus and Mon. perniciosa have significant areas of below average GC content associated with the cox1 gene, its introns, and the block of uORFs following the cox 1 gene. This block of uORFs is missing in Mon. roreri. Mon. perniciosa and Mon. roreri also have distinct areas of low GC content associated with different blocks of uORFs not found in M. palmivorus. Otherwise, the shifts in GC content are muted in M. palmivorus compared to Mon. roreri and Mon. perniciosa. As might be expected, M. tenuissimus and M. crinis-equi show similar patterns of GC content. The pattern of GC content in Mar. scandens, with is higher overall GC content (32.04%), generally lacks the extremes observed in the Mon. perniciosa, Mon. roreri, and M. palmivorus grouping. Overall, patterns of low GC content are associated with the occurrence of introns containing ORFs and groupings of uORFs. There is also a tendency for areas of tRNA blocks to show above average GC content.

Rps3 is the only ribosomal protein encoded in the fungal mitogenome and is involved in the assembly of the 37S ribosomal subunit (Seif et al., 2005). The rps3 gene is extremely diverse in location and organization within mitogenomes and has a complex evolutionary history of acquisition by group I introns, loss of the intron, and establishment of rps3 as a free-standing gene (Sethuraman et al., 2009). Korovesi et al. (2018) proposed two evolutionary routes for rps3 gene accusation into the mitogenome, as a free-standing gene in the case of Basidiomycota and most of the yeasts or an anchored gene within the omega intron in case of Pezizomycotina. Complying with the model, at least one rps3 is found in each mitogenome, existing as a free-standing gene (Figures 1, 3). The rps3 genes in Mon. perniciosa, Mon. roreri, and M. palmivorus show sequence identity (>73%) (Figure 5) and are similarly located in front of the rnL gene (Figure 1). The M. tenuissimus and M. crinis-equi rps3 genes also show sequence identity (>64%) (Figure 5A) and are similarly located prior to nad1 (Figure 3). M. tenuissimus isolate GHA64 carries 2 copies of the rps3 gene (68.9% identity), one being unique, and sharing sequence homology with the Mar. scandens rps3 (73.5% identity) but not the other isolates of M. tenuissimus (Figures 5A,B). Alignment of the rps3 gene region between M. tenuissimus isolate GHA64 and GHA63 shows an insertion of a 479 bp segment in the GHA64 mitogenome, which is the part of the second rps3 copy (Figure 5C). The phylogenetic relationships using the 14 core proteins are largely mimicked by relationships between the rps3 gene comparisons (Figures 4, 5). The rps3 gene provides an alternative to the whole mitogenome as a molecular marker for phylogenic study. The rps3 gene was used for phylogenetic analysis of Colletotrichum species (Pszczółkowska et al., 2020).

Intron dynamics play an important role in altering organization and size of fungal mitogenomes (Himmelstrand et al., 2014; Li et al., 2021). Group I and group II introns are commonly found, among other places, in organelles of higher eukaryotes (Haugen et al., 2005). Often considered selfish DNA, they can act as ribozymes catalyzing their own splicing from a precursor RNA and restoring the translational reading frames of their host genes. Group I introns, which are mobile due to the activity of HEGs, are abundant in fungal mitogenomes (Lang et al., 2007). In some instances, HEGs function as mobile elements, moving independently from their host intron to new locations (Edgell, 2009). These endonucleases, having settled within introns, provide a valuable source of genetic information. The mitogenomes of the six Marasmiaceae family members harbor varying numbers of introns, all corresponding to group I. Most of the introns found in these genomes are group 1A and 1B introns. Mon. perniciosa also carries group 1C introns and M. crinis-equi carries group 1D introns, while Mar. scandens and Mon. roreri carry both group 1C and 1D introns (Supplementary Excel File 1). Introns in mitogenomes can have direct biological consequences. Possibly, the best example of this is the involvement of group I introns blocking the main mutation involved in the resistance against Quinone outside Inhibitor (QoIs) fungicides (Grasso et al., 2006). Recently, Cinget and Bélanger (2020) identified 216 novel group 1D introns involved in QoIs resistance and hypothesized that mobility of the intron across fungal mitogenomes influences the ability to develop resistance to QoIs. Mon. roreri has 14 ORFs within 13 introns, Mon. perniciosa 10 ORFs within 13 introns, Mar. scandens 8 ORFs within 9 introns, and M. crinis-equi and M. palmivorus 6 ORFs within 7 introns. M. tenuissimus isolates GHA63 and GHA79 carry 3 ORFs within 3 introns all within the Cox1 gene while the remaining isolates carry only one ORF within one intron within the cox1 gene (Figure 3). The intronic ORFs identified regardless of species encode GIY-YIG or LAGLIDADG domain-containing HEGs (Supplementary Excel File 1).

Specific sequence comparisons can be made among the intronic ORFs found in the Marasmiaceae mitogenomes due to their limited numbers (Supplementary Figure 9). As expected, the intronic ORFs show substantial sequence differences depending on whether they are GIY-YIG or LAGLIDADG domain-containing HEGs and clusters separately (Figure 6A). For M. tenuissimus, the five isolates carrying a single cox1 intronic ORF show close similarity (>99% identity) while the 3 cox1 intronic ORFs of isolates 63 and 79 show close similarity (>99% identity) in a positionally dependent manner (Figure 6B). As seen in other comparisons, the genetic relationships among these intronic ORFs across species, in some cases, also show specificity as to position within specific core genes (Huang et al., 2021). For example, the terminal GIY-YIG containing intronic ORF in cox1 shows sequence similarity (>95%) across all species except M. tenuissimus (Figure 6A). The positional relationships among intronic ORFs are not consistent, as is expected considering the variable number of intronic ORFs found depending on the species. The terminal intronic ORF (LAGLIDADG type) of the cox1 gene in M. tenuissimus (orf322-I3) is similar to (>92% sequence similarity) the next to last ORFs (orf332-I4 and I5) in the more distantly related Mon. roreri, and Mon. perniciosa. Another example occurs between the middle LAGLIDADG type intronic ORF of cox1 in M. crinis-equi (orf388-I3) and the initial cox1 intronic ORF in Mon. roreri (orf415-I1), and Mon. perniciosa (orf422-I1). The relationships between LAGLIDADG type intronic ORFs appear more focused, involving smaller subsets within the six species. These relationships do not always involve the Moniliophthora, for example, the next to last intronic ORF (LAGLIDADG type) of the cox1 gene in M. crinis-equi (orf278-I4) and M. palmivorus (orf278-I3) have 99% sequence similarity not shared with any Moniliophthora intronic ORF. The single GIY-YIG containing intronic ORF found in cob also shows sequence similarity in M. palmivorus, Mar. scandens, M. crinis-equi, and Mon. roreri (Supplementary Figure 9). The cob genes in M. tenuissimus and Mon. perniciosa do not carry a GIY-YIG containing intron. Despite the close sequence similarities between species for many of the intronic ORFs, all species compared except M. crinis-equi have at least one intronic ORF that has minimal sequence similarity (% < 40) with any other intronic ORF. The middle cox1 intronic ORF of M. tenuissimus isolate GHA63 (orf311-I2) has minimal sequence homology with any cox1 intronic ORF identified and instead shows its closest sequence homology with a Mon. roreri cob intronic ORF (orf357). The sequence and positional similarities and dissimilarities found for many of the Marasmiaceae intronic ORFs indicates a shared ancestral origin for many and varied histories of expansion and contraction dependent on the species, which is typical for mitogenomes (Megarioti and Kouvelis, 2020; Mukhopadhyay and Hausner, 2021). The intronic ORFs also share significant homologies with mitogenome intronic ORFs of unrelated fungal species (Supplementary Excel File 2), often outside the basidiomycetes, something commonly cited as evidence of introgression between diverse species. Similar intron dynamics were observed within the Agaricales (Huang et al., 2021).