Mitochondrial Genomics of Six Cacao Pathogens From the Basidiomycete Family Marasmiaceae

Thread blight disease has recently been described as an emerging disease on cacao (Theobroma cacao) in Ghana. In Ghana, thread blight disease is caused by multiple species of the Marasmiaceae family: Marasmius tenuissimus, M. crinis-equi, M. palmivorus, and Marasmiellus scandens. Interestingly, two additional members of the Marasmiaceae; Moniliophthora roreri (frosty pod rot) and Moniliophthora perniciosa (witches’ broom disease), are major pathogens of cacao in the Western hemisphere. It is important to accurately characterize the genetic relationships among these economically important species in support of their disease management. We used data from Illumina NGS-based genome sequencing efforts to study the mitochondrial genomes (mitogenomes) of the four cacao thread blight associated pathogens from Ghana and compared them with published mitogenomes of Mon. roreri and Mon. perniciosa. There is a remarkable interspecies variation in mitogenome size within the six cacao-associated Marasmiaceae species, ranging from 43,121 to 109,103 bp. The 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). Among seven M. tenuissimus mitogenomes sequenced, there is variation in size and sequence pointing to divergent evolution patterns within the species. The intronic regions show a high degree of sequence variation compared to the conserved sequences of the 14 core genes. The intronic ORFs identified, regardless of species, encode GIY-YIG or LAGLIDADG domain-containing homing endonuclease genes. Phylogenetic relationships using the 14 core proteins largely mimic the phylogenetic relationships observed in gene order patterns, grouping M. tenuissimus with M. crinis-equi, and M. palmivorus with Mon. roreri and Mon. perniciosa, leaving Mar. scandens as an outlier. The results from this study provide evidence of independent expansion/contraction events and sequence diversification in each species and establish a foundation for further exploration of the evolutionary trajectory of the fungi in Marasmiaceae family.

Thread blight disease has recently been described as an emerging disease on cacao (Theobroma cacao) in Ghana. In Ghana, thread blight disease is caused by multiple species of the Marasmiaceae family: Marasmius tenuissimus, M. crinis-equi, M. palmivorus, and Marasmiellus scandens. Interestingly, two additional members of the Marasmiaceae; Moniliophthora roreri (frosty pod rot) and Moniliophthora perniciosa (witches' broom disease), are major pathogens of cacao in the Western hemisphere. It is important to accurately characterize the genetic relationships among these economically important species in support of their disease management. We used data from Illumina NGS-based genome sequencing efforts to study the mitochondrial genomes (mitogenomes) of the four cacao thread blight associated pathogens from Ghana and compared them with published mitogenomes of Mon. roreri and Mon. perniciosa. There is a remarkable interspecies variation in mitogenome size within the six cacao-associated Marasmiaceae species, ranging from 43,121 to 109,103 bp. The 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). Among seven M. tenuissimus mitogenomes sequenced, there is variation in size and sequence pointing to divergent evolution patterns within the species. The intronic regions show a high degree of sequence variation compared to the conserved sequences of the 14 core genes. The intronic ORFs identified, regardless of species, encode GIY-YIG or LAGLIDADG domain-containing homing endonuclease genes. Phylogenetic relationships using the 14 core proteins largely mimic the phylogenetic relationships observed in gene order patterns, grouping M. tenuissimus with M. crinis-equi, and M. palmivorus with Mon. roreri and Mon. perniciosa, leaving Mar. scandens as an outlier. The results from this study provide evidence of independent expansion/contraction events and sequence diversification in each species and establish a foundation for further exploration of the evolutionary trajectory of the fungi in Marasmiaceae family.

INTRODUCTION
There are more than 1,000 species of fungi in the family Marasmiaceae, which includes the genus Marasmius Fr. (1836: 339) (Agaricales, Basidiomycota). Although most are saprobes, the family does include plant pathogens. Examples include Marasmiellus scandens on cacao and other tropical trees (Massee, 1910), Marasmius crinis-equi on cacao and the tea tree (Su et al., 2011), M. puerariae on kudzu vine (Kirschner et al., 2013), M. palmivorus on cacao and oil and coconut palms (Desjardin and Perry, 2017) and M. graminum on rice (Gaire et al., 2020). In addition, two devastating cacao diseases, frosty pod rot and witches' broom disease which occur in the Western hemisphere, are caused by the Marasmiaceae species Moniliophthora roreri and Moniliophthora perniciosa, respectively (Evans et al., 2003).
Thread blight disease (TBD) is an emerging threat to cacao production in Ghana (Amoako-Attah et al., 2016, 2020. It is caused by canopy-dwelling basidiomycetes and mostly seen as a network of mycelial rhizomorphs, often called "threads, " that grow along twigs and branches covering leaves creating distinctive leaf blight symptoms. Historically, based on morphological methods, TBD in cacao has been attributed to the causal agent Mar. scandens (syn: M. byssicolor) in Ghana and in other cacao growing areas of the world (Opoku et al., 2007;Amoako-Attah et al., 2016). However, recent morphological and molecular analyses identified four different fungal species, of the Marasmiaceae, causing TBD-like symptoms on cacao in Ghana (Amoako-Attah et al., 2020). The four species were M. crinisequi (abundant thin, black, "horse hair"-type rhizomorphs), M. tenuissimus (scattered brown or whitish to brownish-white rhizomorphs), M. palmivorus (aggregates of shiny-or silkywhite hyphae), and Mar. scandens (faint-cream or dull-white hyphae). Among the four species, M. tenuissimus was the most frequently isolated TBD-associated fungus in Ghana (Amoako-Attah et al., 2020). When considered with the Moniliophthora spp. that also cause disease on cacao, these six species within the Marasmiaceae create a fertile ground for dissection of genetic diversity in association with varied disease processes in a single crop.
Mitochondria play important roles in eukaryotic cells, mainly with respect to respiratory metabolism and energy supply (Kobayashi et al., 1996). Production of reactive oxygen species during ATP synthesis within the mitochondria makes the mitochondrial DNA more susceptible to damage and mutations, compared to nuclear DNA (Richter, 1992). At the same time, the mitochondrial genome is subject to all the mechanisms responsible for maintaining nuclear genome integrity (Kaniak-Golik and Skoneczna, 2015). Recent findings suggest that the mutation rate of mitochondrial DNA can be significantly lower in plants and fungi than in animals, with some exceptions (Torriani et al., 2014;Sandor et al., 2018). Fungal mitochondrial genomes (mitogenomes) are typically small and present in multiple copies in each cell with a highly compact gene organization (Franco et al., 2017). Mitogenomes usually harbor 14 coregenes encoding proteins involved in electron transport and oxidative phosphorylation, as well as untranslated genes of the small and large ribosomal RNA (rRNA) subunits and a set of transfer RNA (tRNA) genes (Bullerwell and Lang, 2005). Fungal mitogenomes sometimes also carry genes encoding the mitochondrial ribosomal protein S3 (rps3) and the RNA subunit of the mitochondrial RNase P (Franco et al., 2017). Beside these genes, fungal mitogenomes are also characterized by their variable number of group I and group II introns that may carry open reading frames with motifs of homing endonuclease genes (HEGs) (Bullerwell and Lang, 2005;Lang et al., 2007). HEGs are genetic mobile elements that encode site-specific DNA endonucleases and promote their own propagation (Burt and Koufopanou, 2004;Edgell, 2009) resulting in the insertion, deletion or mutation of DNA sequences (Stoddard, 2011). Such, mobile introns like HEGs represent one of the major sources of variation within fungal mitogenomes (Wu and Hao, 2014;Repar and Warnecke, 2017). Fungal mitogenomes often differ in gene order and composition, pseudogene content and length of intergenic regions due to rearrangements caused by recombination. The presence of double-stranded RNA elements and self-splicing introns can also add to fungal mitogenome variability (Kulik et al., 2020). Variations in mitogenomes can provide vital clues into the evolution, population genetics, and biology of the organisms involved and are a rich source of novel genotyping markers due to the presence of these mobile introns.
Unlike plants and animals, fungi exhibit a diversity of mitochondrial DNA inheritance patterns from strictly uniparental to biparental, or a mixture of both, as well as recombinant mitochondrial DNA genotypes (Wilson and Xu, 2012;Xu and Wang, 2015). Uniparental mitochondrial inheritance is very common in filaments basidiomycete (Xu and Li, 2015). Moreover, Mon. roreri propagates clonally (Ali et al., 2015;Díaz-Valderrama and Aime, 2016), whereas Mon. perniciosa is primarily homothallic (Kües and Navarro-González, 2010), and thus is expected to have primarily uniparental mitochondrial DNA inheritance. The mitochondrial DNA inheritance patterns of other Marasmiaceae need further confirmation.
With the emergence of easier and more affordable next generation sequencing, access to whole genome sequences has increased dramatically, along with the number of completely sequenced mitogenomes. But within the Marasmiaceae family, only two complete mitogenomes from Mon. roreri and Mon. perniciosa have been reported (Formighieri et al., 2008;Costa et al., 2012), representing a gap in terms of understanding the evolutionary and biological aspects of this large group of fungi. As emerging pathogens of cacao and other tropical crops, it is important to understand the genomics of these economically important Marasmiaceae species. Therefore, we sequenced and assembled the complete mitogenomes of M. crinisequi, M. tenuissimus, M. palmivorus, and Mar. scandens. As M. tenuissimus was the major causal agent TBD in Ghana and reported to have different morphotypes (Amoako-Attah et al., 2020), we also sequenced and assembled the complete mitogenomes of six additional M. tenuissimus isolates and performed a comparative mitogenome analysis along with Mon. roreri and Mon. perniciosa.

Fungal Isolates
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).

Sequencing and Mitogenome Assembly
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 . The short reads were assembled using "SOAPdenovo" ( 1 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.

Annotation of the Mitogenome
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 2 . To further confirm genes boundaries as well as intronexon 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.

Phylogenetic Inference
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 3 . 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.

Mitogenome Synteny and Structural Comparisons
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 4 . 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).

GC-Plots
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 (TRs)
Tandem repeats were identified using the Tandem Repeats Finder (Benson, 1999)

Overview of the Mitogenomes
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  (Formighieri et al., 2008) between nad1 and cob in association with other rearrangements (Figure 1).

Primary Mitochondrial Gene Composition and Order
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  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.

(Supplementary
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. crinisequi 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).

Phylogenetic Analysis
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 . 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.

Mitogenome Size Variation and Sources of Genome Expansion
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 FIGURE 3 | Conservation of gene order and Intraspecies variation in the mitogenome structure of seven isolates of the cacao thread blight associated pathogen Marasmius tenuissimus. Mitochondrial gene annotation was performed with MFannot using the NCBI translation code 4, and the physical map of the mitogenomes were created with OrganellarGenome-DRAW (OGDRAW) v 1.2.
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 scandens) and the Western hemisphere frosty pod rot (Mon. roreri) and witches' broom (Mon. perniciosa) pathogens of cacao along with some other members of Marasmiaceae and Agaricales. The Rhizoctonia solani was used as an outgroup. The analysis was based on the amino-acid sequences of 14 conserved mitochondrial proteins with 3,154 distinct alignment positions and 1,000 rapid bootstrap inferences. Sequences were combined and aligned using ClustalW2 tool under default setting and the phylogenetic tree was reconstructed using the Maximum Likelihood method. The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. Analyses were conducted in MEGA6 (Tamura et al., 2013). The analysis was re-run with seven M. tenuissimus and M. crinis-equi as an outgroup.
of Mon. perniciosa carries 40 accessory genes, many of which are part of an inserted linear plasmid (Formighieri et al., 2008)

GC Content
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. crinisequi 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.

The Ribosomal Protein Rps3
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 Sequences were combined and aligned using ClustalW2 tool under default setting and the phylogenetic tree was reconstructed using the Maximum Likelihood method. The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. Analyses were conducted in MEGA6 (Tamura et al., 2013). (C) Visualization of rps3 gene region sequence alignment between M. tenuissimus isolate GHA64 and GHA63. The mVISTA program (http://genome.lbl.gov/vista/mvista/submit.shtml) was used to compare the section of mitochondrial genomes with default parameters in Shuffle-LAGAN mode and a sequence conservation profile was visualized in an mVISTA plot (Frazer et al., 2004). intron in case of Pezizomycotina. Complying with the model, at least one rps3 is found in each mitogenome, existing as a freestanding 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).

Introns and Accessory Mitochondrial Genes
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 Sequences were combined and aligned using ClustalW2 tool under default setting and the phylogenetic tree was reconstructed using the Maximum Likelihood method. The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. Analyses were conducted in MEGA6 (Tamura et al., 2013). (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. crinisequi 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).

CONCLUSION
The mitogenomes of the six Marasmiaceae species compared display a wide range of genome sizes and composition. The Mon. perniciosa and Mon. roreri mitogenomes are the most complex, having the largest sizes and numbers of both intronic ORFs and uORFs. The most commonly encountered TBD-associated pathogen in Ghana, M. tenuissimus, has the smallest mitogenomes. The well-defined gene order patterns and core gene sequence similarities group M. tenuissimus with M. crinis-equi, and M. palmivorus with Mon. roreri and Mon. perniciosa, leaving Mar. scandens as an outlier among the species studied. Most of the intronic ORFs of the TBD-associated pathogens share sequence homology with intronic ORFs found in the Mon. perniciosa and Mon. roreri mitogenomes, although exceptions exist. Together, the reduced numbers of intronic ORFs found in the TBD-associated pathogens compared to the Moniliophthoras, the common homologies of most intronic ORFs among the species, and the common occurrence of uORFs lacking shared homologies provide evidence of independent expansion/contraction events and sequence diversification in each species depending on their physiological and developmental needs. The freestanding rps3 gene is a promising marker that may also be useful in inferring phylogenies within Marasmiaceae. Among the seven mitogenomes of M. tenuissimus sequenced, there is variation in both size and sequences pointing to three distinct phylactic groups within the species. This is the first mitogenome phylogeny of the genus Marasmius, and results highlight the potential for future systematics changes in Marasmius and crypticism in M. tenuissimus.

DATA AVAILABILITY STATEMENT
The mitogenome sequence data presented in the study are deposited at NCBI GenBank under the accession numbers MZ615343-MZ615352.

AUTHOR CONTRIBUTIONS
BB, SA, and LM provided intellectual and editorial comments. SA and BB conceived and designed the experiments and wrote the first draft manuscript. IA-A and EK-A isolation of fungal isolates. SA performed the experiments. SA and JS analyzed the data. All authors contributed to the manuscript revision, read, and approved the submitted version.

FUNDING
This work was supported by the U.S. Department of Agriculture.

ACKNOWLEDGMENTS
IA-A thanks International Forest Services Program for processing of a travel grant to USDA/ARS laboratory at Beltsville Agricultural Research Center-West, Beltsville, United States. References to a company or product were only for information and do not imply approval or recommendation of the company or product to the exclusion of others that may also be suitable. USDA is an equal opportunity provider and employer.