The First Five Mitochondrial Genomes for the Family Nidulariaceae Reveal Novel Gene Rearrangements, Intron Dynamics, and Phylogeny of Agaricales

The family Nidulariaceae, consisting of five genera including Cyathus, is a unique group of mushrooms commonly referred to as bird’s nest fungi due to their striking resemblance to bird’s nests. These mushrooms are considered medicinal mushrooms in Chinese medicine and have received attention in recent years for their anti-neurodegenerative properties. However, despite the interest in these mushrooms, very little is known about their mitochondrial genomes (mitogenomes). This study is the first comprehensive investigation of the mitogenomes of five Nidulariaceae species with circular genome structures ranging in size from 114,236 bp to 129,263 bp. Comparative analyses based on gene content, gene length, tRNA, and codon usage indicate convergence within the family Nidulariaceae and heterogeneity within the order Agaricales. Phylogenetic analysis based on a combined mitochondrial conserved protein dataset provides a well-supported phylogenetic tree for the Basidiomycetes, which clearly demonstrates the evolutionary relationships between Nidulariaceae and other members of Agaricales. Furthermore, phylogenetic inferences based on four different gene sets reveal the stability and proximity of evolutionary relationships within Agaricales. These results reveal the uniqueness of the family Nidulariaceae and its similarity to other members of Agaricales; provide valuable insights into the origin, evolution, and genetics of Nidulariaceae species; and enrich the fungal mitogenome resource. This study will help to expand the knowledge and understanding of the mitogenomes in mushrooms.


Introduction
Mitochondria are essential organelles in fungi, and they play crucial roles in the growth, development, and adaptation of fungal species [1,2].The fungal mitochondrial genome (mitogenome) is typically made up of a circular DNA that encodes fourteen mitochondria special protein-coding genes (PCGs), ribosomal RNAs, and transfer RNAs, which are usually found in clusters [3,4].The size of fungal genomes varies greatly [4], ranging from 19,431 bp in Schizosaccharomyces pombe [5] to 332,165 bp in Golovinomyces cichoracearum [6].Mitochondrial genome size varies considerably among and within fungal species [7,8], and the principal drivers of this variation are the presence and proliferation of introns within the protein-coding genes, as well as significant non-coding regions (NCRs) [3,9].To gain insight into the evolution of fungal mitochondrial genomes, it is -critical to investigate the dynamics of introns and determine the origin and function of these NCRs [10].The mitochondrial genomes (mitogenomes) of fungi are valuable tools for studying the evolutionary history and phylogenetic relationships due to their small sizes, conserved orthologous genes, high copy numbers, low recombination rates, and high evolutionary rates [11][12][13][14].In fact, the mitogenome is often referred to as the "second genome" of eukaryotes due to its significant contribution to understanding eukaryotic individual growth and environmental adaptation [15,16].Additionally, studying the mitogenomes of pathogenic fungi can provide insight into the mechanisms underlying their pathogenicity in plants [17,18].Despite its importance, the fungal mitochondrial genome has received relatively little attention.To date, fewer than 1300 fungal mitochondrial genomes have been deposited in the NCBI database (https://www.ncbi.nlm.nih.gov/genome/browse#!/organelles/, accessed on 1 January 2023), and this number is far less than the number of published fungal nuclear genomes and particularly the number of DNA parts.
Mushrooms, the common name for macrofungi capable of forming specific fruiting bodies, have emerged as a significant source of food and medicine for humans due to their nutritional value and rich content of bioactive compounds [19][20][21].As the importance of medicinal mushrooms for human health continues to grow and sequencing technologies become increasingly available, the genomes of valuable medicinal mushrooms are being sequenced and deciphered [22][23][24].However, the mitogenomes of these mushrooms remain relatively understudied.Among the well-known medicinal mushrooms in Chinese medicine, Cyathus species have garnered much attention in recent years for the production of cyathane diterpenoids with significant neurotrophic and neuroprotective activities [25][26][27][28][29].The family Nidulariaceae, commonly known as "bird's nest fungi", is a distinctive and fascinating group of fungi in the fungal kingdom [30], comprising five genera: Cyathus, Crucibulum, Mycocalia, Nidularia, and Nidula (Figure 1).Although the family Nidulariaceae, represented by the genus Cyathus, is classified in the order Agaricales, there are obvious morphological differences between the family Nidulariaceae and other mushrooms of the order Agaricales (Figure 1).To date, 30 mitogenomes have been reported from the order Agaricales [31], but none have been reported from the family Nidulariaceae.In addition to the obvious differences in morphology, whether there are also obvious differences in the mitogenomes of Nidulariaceae and other Agaricales mushrooms is an intriguing scientific question that deserves attention and investigation.
The present study aimed to investigate the mitogenomes of five Cyathus species, namely C. striatus 87405, C. striatus AH44044, C. jiayuguanensis 765, C. pallidus QL1, and C. stercoreus NPCB004.This is the first report of Nidulariaceae mitogenomes and allowed the characterization of these genomes, including their base composition, gene distribution, intron dynamics, codon usage bias, and PCG characteristics.Comparative analyses were performed within the order Agaricales to explore the correlation between genome size and the proportion of unannotated open reading frames, homing endonucleases, intronic regions, and intergenic regions.The results showed a positive correlation between genome size and these factors.Genetic distance and evolutionary selection pressure (Ka/Ks) analyses indicated that the nad6 and rps3 were not conserved in the order Agaricales and that rps3 showed species specificity.Furthermore, gene rearrangement analysis revealed uniformity within the genus Cyathus and heterogeneity between Cyathus and other members of the order Agaricales.Phylogenetic inference based on combined sequences of fourteen conserved proteins and four datasets (PCG, PCG12, PCGR, and PCG12R) within Agaricales further confirmed the evolutionary position of Nidulariaceae within the order and the topology of the Agaricales evolutionary tree.This study provides valuable insights into the origin, evolution, and genetics of Nidulariaceae species and closely related members within Agaricales.The reported mitogenomes provide a valuable resource for future studies of Cyathus species and for comparative mitogenomics among Basidiomycetes mushrooms.The present study aimed to investigate the mitogenomes of five Cyathus species, namely C. striatus 87405, C. striatus AH44044, C. jiayuguanensis 765, C. pallidus QL1, and C. stercoreus NPCB004.This is the first report of Nidulariaceae mitogenomes and allowed the characterization of these genomes, including their base composition, gene distribution, intron dynamics, codon usage bias, and PCG characteristics.Comparative analyses were performed within the order Agaricales to explore the correlation between genome size and the proportion of unannotated open reading frames, homing endonucleases, intronic regions, and intergenic regions.The results showed a positive correlation between genome size and these factors.Genetic distance and evolutionary selection pressure (Ka/Ks) analyses indicated that the nad6 and rps3 were not conserved in the order Agaricales and that rps3 showed species specificity.Furthermore, gene rearrangement analysis revealed uniformity within the genus Cyathus and heterogeneity between Cyathus and other members of the order Agaricales.Phylogenetic inference based on combined sequences of fourteen conserved proteins and four datasets (PCG, PCG12, PCGR, and PCG12R) within Agaricales further confirmed the evolutionary position of Nidulariaceae within the order and the topology of the Agaricales evolutionary tree.This study provides valuable insights into the origin, evolution, and genetics of Nidulariaceae species and closely related members within Agaricales.The reported mitogenomes provide a valuable resource for future studies of Cyathus species and for comparative mitogenomics among Basidiomycetes mushrooms.

Characterization of the Five Cyathus Mitogenomes
The five Cyathus species possessed circular mitogenomes of varying sizes: 115,936 bp for C. jiayuguanensis 765, 127,920 bp for C. striatus 87405, 120,589 bp for C. striatus AH440044, 114,236 bp for C. pallidus QL1, and 129,263 bp for C. stercoreus NPCB004 (Figure 2).With an average GC content of 28.39%, interspecific variation is evident-the highest being 30.12% in C. jiayuguanensis 765 and the lowest being 27.34% in C. striatus AH44044.While C. jiayuguanensis 765 exhibited negative AT skew (−0.011) and positive GC skew (0.010), the remaining four mitogenomes showed positive AT skews (0.066-0.094) and negative GC skews (−0.147 to −0.122).In total, the mitogenomes contained 43, 44, 39, 31, and 38 protein-coding genes for C. striatus 87405 and AH44044, C. jiayuguanensis 765, C. pallidus QL1, and C. stercoreus NPCB004, respectively (Table S1).Each mitogenome encoded 15 core PCGs, including five cytochrome C oxidase subunit genes, five ATP synthase subunit genes, six NADH dehydrogenase subunit genes, and one ribosomal protein S3 gene.These core PCGs, except for rps3, are conserved mitochondrial genes involved in the oxidative phosphorylation pathway [32].Furthermore, two ribosomal RNAs (rns and rnl) were present in all five Cyathus mitogenomes, with the number of tRNA genes ranging from 24 to 28 (Figure 2, Table S2).The proportions of intergenic regions, PCGs, intronic regions, and RNA regions we comparable across all five mitogenomes.Intergenic regions were the most abundant, a eraging 51.94% of the mitogenomes, while RNA regions contributed the least, averagin 6.50%.PCGs comprised an average of 30.37% of the five Cyathus mitogenomes, and i tronic regions represented 11.20%.There was significant variation in the proportions PCGs and intergenic regions between species.PCGs comprised an average of 34 S3).The majority of fungal mitogenomes contained a collection of uncharacterized open reading frames (un_ORFs), which encode either unidentified proteins or misidentified known proteins [33].Notably, five un_ORFs encoding unknown proteins were identified in C. striatus 87405 and AH44044, accompanied by several copies of genes encoding GIY and LAGLIDADG endonucleases.Similarly, C. jiayuguanensis 765 harbored un_ORFs and putative reverse transcriptase encoding genes, as well as an abundance of GIY and LAGLIDADG endonuclease gene copies.The remaining two species, C. pallidus QL1 and C. stercoreus NPCB004, exhibited a limited number of un_ORFs, alongside multiple copies of GIY/LAGLIDADG-related genes (Table S2).We also observed introns present in the mitochondrial genomes of different species within this genus, with significant differences in their distribution.Intronic analyses revealed that C. striatus 87405 and AH44044 contained more introns than other species, distributed in cob, cox1, cox3, nad1, nad5, and rnl genes, whereas C. jiayuguanensis 765 had fewer introns, distributed in atp9, cob, cox1-3, nad1, nad5, and rnl genes.Except for C. pallidus QL1, all species exhibited group II introns (Table S2).
The proportions of intergenic regions, PCGs, intronic regions, and RNA regions were comparable across all five mitogenomes.Intergenic regions were the most abundant, averaging 51.94% of the mitogenomes, while RNA regions contributed the least, averaging 6.50%.PCGs comprised an average of 30.37% of the five Cyathus mitogenomes, and intronic regions represented 11.20%.There was significant variation in the proportions of PCGs and intergenic regions between species.PCGs comprised an average of 34 2, Table S3).

Codon Usage Analysis
To analyze the codon usage patterns in different Agaricales mitogenomes, we compared 34 mitogenomes, including 29 published in the database and 5 Cyathus mitogenomes assembled by ourselves.The cox1 genes in 19 Agaricales mitogenomes utilized ATG as start codons, while others employed GTG as start codons, with Amanita muscaria, C. jiayuguanensis 765, and Stropharia rugosoannulata being exceptions, using TTG.Most of the remaining fourteen core PCGs predominantly used ATG as start codons, although the atp6 gene of Marasmius tenuissimus and the rps3 gene of S. rugosoannulata used GTG as start codons.In addition, the cox2 gene of C. striatus, the nad1 gene of Laccaria bicolor, and the nad4 and nad6 genes of S. rugosoannulata used TTG as start codons.Among the core PCGs of the 34 mitotic genomes examined, TAA was the most prevalent stop codon, followed by TAG (Table S4).
Subsequent analysis revealed considerable differences in the use of start and stop codons among Agaricales species, even among closely related species.For instance, the cox1 gene of C. jiayuguanensis 765 used TTG as the start codon, whereas C. striatus 87405 and AH44044 employed ATG as the start codon.The cox3 gene of C. jiayuguanensis 765, C. pallidus QL1, and C. stercoreus NPCB004 adopted TAG as stop codons, while that of C. striatus 87405 and C. striatus AH44044 used TAA as stop codons.Within Amanita species, the atp6 gene of A. muscaria used TAG as the stop codon, while A. phalloides and A. thiersii utilized TAA as stop codons (Table S4).
Codon usage has been shown to directly affect the rate and energy required for translation.Genes in mitochondria are believed to favor specific codons to save time and conserve energy for cell growth [34].Our in-depth analysis revealed that the most commonly used codons of five Cyathus species mitogenomes were UUA (for leucine; Leu), followed by AAA (for arginine; Arg).The codon usage bias for different amino acids is essentially the same among the five Cyathus species; for example, UGU was the most frequently used codon for cysteine, and GAA was the preferred codon for glutamic acid (Figure 3A, Table S5).Moreover, the codon usage bias of the Cyathus mitogenomes is mainly consistent with that of Agaricales species mitogenomes (Figure 3B, Table S5), which suggests that mitogenomes within the same order tend to share similar codon usage biases.Some exceptions may be attributed to the diversity of un_ORFs in mitotic genomes.

Molecular-Signature-Based Assessment of the Coding Potential of PCGs
AT-rich mitochondrial genomes have been demonstrated to preferentially use AT rich codons to encode their genes [35].To evaluate this observation in 34 Agaricales specie we analyzed the correlation between AT content, GC content, and the usage frequency o AT-rich and GC-rich codons.Our results revealed a significant positive correlation b tween GC/AT content and the usage frequency of GC-rich/AT-rich codons (correlation co efficient = 0.85, p-value = 1.60 × 10 −10 ) (Figure 4A, Table S6).
We then calculated the GC content of the CDS and intronic regions, and the resul demonstrated that the GC content of the intronic regions was significantly different from that of the CDS by the Wilcoxon test (p-value = 0.022).In mammals, GC content can b used as a splicing marker because of exons' low GC content and intronic sequences' hig

Molecular-Signature-Based Assessment of the Coding Potential of PCGs
AT-rich mitochondrial genomes have been demonstrated to preferentially use AT-rich codons to encode their genes [35].To evaluate this observation in 34 Agaricales species, we analyzed the correlation between AT content, GC content, and the usage frequency of AT-rich and GC-rich codons.Our results revealed a significant positive correlation between GC/AT content and the usage frequency of GC-rich/AT-rich codons (correlation coefficient = 0.85, p-value = 1.60 × 10 −10 ) (Figure 4A, Table S6).
over, distinct replication processes on the forward and reverse strands can result in different mutation rates, affecting codon usage [38].Our t-SNE findings also imply that un_ORFs and conserved genes might experience different selective pressures.Considering the limited number of tRNA gene copies for each amino acid in mitogenomes, codon optimization is expected to enhance mRNA stability and translation rate.Hence, we recommend that molecular biology experiments be carried out to explore the structure and function of un_ORFs further.

Repetitive Sequence Analysis
The mitogenome of organisms is known to contain various types of repetitive sequences, including tandem and interspersed repeats [39,40].We analyzed the mitogenomes of five Cyathus species to identify and characterize their repetitive sequences.A total of 30 repetitive sequences were present in these five Cyathus mitogenomes, with the number of repeats ranging from 11 to 13.The size of these sequences varied from 39 to 624 bp, with the largest repeat (624 bp) being identified in C. jiayuguanensis 765.This repeat sequence was located in the intergenic region between the genes cox3 and orf275, as well as in the partial region of the rns gene.Another large repeat sequence of 502 bp was identified in the partial region of the fourth exon and third intron of the nad3 gene of C. jiayuguanensis 765.Pairwise nucleotide identities between the five Cyathus mitogenomes ranged from 80.54% to 100%, indicating a high degree of conservation among these species.The repetitive sequences accounted for 0.78% to 3.10% of the five Cyathus mitogenomes.The highest proportion of repeat sequences was found in C. jiayuguanensis 765, while C. stercoreus NPCB004 had the lowest content of repetitive sequences (Table S8).
A total of 1228 tandem repeats were found in the five Cyathus mitogenomes, with the number of repeats ranging from 177 to 313.The longest tandem repeat sequence of 952 bp was identified in C. jiayuguanensis.Most of the tandem repeats were duplicated once or We then calculated the GC content of the CDS and intronic regions, and the results demonstrated that the GC content of the intronic regions was significantly different from that of the CDS by the Wilcoxon test (p-value = 0.022).In mammals, GC content can be used as a splicing marker because of exons' low GC content and intronic sequences' high GC content [36].Since the GC content of intronic regions (28.35 ± 4.06%) is higher than that of CDS (26.46 ± 2.31%), our results suggest that the GC content of intronic regions may act as a potential splicing signal in Agaricales mitogenomes (Figure 4B, Table S6).
The t-SNE results visually showed that the relative synonymous codon usage (RUSU) value of core PCGs differed from that of un_ORFs.Approximately 19.25% of the core protein-coding genes and 17.45% of the unidentified genes were present on the reverse strand (Table S7), which might contribute to the differential RSCU values between the two groups tested.Interestingly, the un_ORFs formed only one cluster, whereas the core PCGs formed two clusters, one of which partially overlapped with the un_ORFs cluster (Figure 4C).This suggested the presence of variation in codon usage bias within the same group.A correlation between gene composition in mitogenomes and location on the strand has been reported in several mammals [37], hinting at a similar pattern in Agaricales species.Moreover, distinct replication processes on the forward and reverse strands can result in different mutation rates, affecting codon usage [38].Our t-SNE findings also imply that un_ORFs and conserved genes might experience different selective pressures.Considering the limited number of tRNA gene copies for each amino acid in mitogenomes, codon optimization is expected to enhance mRNA stability and translation rate.Hence, we recommend that molecular biology experiments be carried out to explore the structure and function of un_ORFs further.

Repetitive Sequence Analysis
The mitogenome of organisms is known to contain various types of repetitive sequences, including tandem and interspersed repeats [39,40].We analyzed the mitogenomes of five Cyathus species to identify and characterize their repetitive sequences.A total of 30 repetitive sequences were present in these five Cyathus mitogenomes, with the number of repeats ranging from 11 to 13.The size of these sequences varied from 39 to 624 bp, with the largest repeat (624 bp) being identified in C. jiayuguanensis 765.This repeat sequence was located in the intergenic region between the genes cox3 and orf275, as well as in the partial region of the rns gene.Another large repeat sequence of 502 bp was identified in the partial region of the fourth exon and third intron of the nad3 gene of C. jiayuguanensis 765.Pairwise nucleotide identities between the five Cyathus mitogenomes ranged from 80.54% to 100%, indicating a high degree of conservation among these species.The repetitive sequences accounted for 0.78% to 3.10% of the five Cyathus mitogenomes.The highest proportion of repeat sequences was found in C. jiayuguanensis 765, while C. stercoreus NPCB004 had the lowest content of repetitive sequences (Table S8).
A total of 1228 tandem repeats were found in the five Cyathus mitogenomes, with the number of repeats ranging from 177 to 313.The longest tandem repeat sequence of 952 bp was identified in C. jiayuguanensis.Most of the tandem repeats were duplicated once or twice in the five Cyathus mitogenomes, with the highest number of duplications (22) observed in the C. pallidus QL1 mitogenome.The proportions of tandem repeat sequences gradually decreased in the mitogenomes of C. jiayuguanensis, C. stercoreus NPCB004, C. pallidus QL1, C. striatus 87405, and AH44044 and were 19.95%, 15.69%, 15.37%, 11.46%, and 9.76%, respectively (Table S9).Our results suggest that repetitive sequences, especially tandem repeats, are common features of Cyathus mitogenomes.

Genetic Distance and Evolutionary Rates of PCGs
Within the 15 detected core PCGs, the rps3 gene exhibited the highest median K2P genetic distance in the 34 Agaricales mushrooms, preceded by the nad6 and nad3 genes.This was an indication that the three genes may have diverged significantly in the course of evolution.Conversely, the nad4L gene had the lowest median K2P distance (Figure 5, Table S10), suggesting high conservation, probably due to its role as a minimal multi-pass membrane protein required for catalysis in the mitochondrial membrane respiratory chain NADH dehydrogenase (Complex I).
twice in the five Cyathus mitogenomes, with the highest number of duplications (22) S9).Our results suggest that repetitive sequences, especially tandem repeats, are common features of Cyathus mitogenomes.

Genetic Distance and Evolutionary Rates of PCGs
Within the 15 detected core PCGs, the rps3 gene exhibited the highest median K2P genetic distance in the 34 Agaricales mushrooms, preceded by the nad6 and nad3 genes.This was an indication that the three genes may have diverged significantly in the course of evolution.Conversely, the nad4L gene had the lowest median K2P distance (Figure 5, Table S10), suggesting high conservation, probably due to its role as a minimal multi-pass membrane protein required for catalysis in the mitochondrial membrane respiratory chain NADH dehydrogenase (Complex I).
For the 15 scrutinized core PCGs, the rps3 gene had the highest observed Ka value, while the nad4L gene demonstrated the lowest Ka value, parallel to the K2P patterns.The synonymous substitution rate (Ks) of the cox1 gene was the highest, while that of the rps3 gene was the lowest among the 34 Agaricales species.The Ka/Ks values of 15 core PCGs were <1, indicating that these genes were under purifying selection pressure.In general, a higher rps3 ratio is more likely with greater species diversity in the Ka/Ks calculation.Intriguingly, the average Ka/Ks value for rps3 from the 34 Agaricales species was 0.82 (Figure 5, Table S10), which showed that the rps3 gene may not have been under positive selection and that it is still in the process of phylogenetic differentiation.In conclusion, our findings illuminate the evolutionary dynamics of core PCGs in Agaricales species and emphasize the pivotal role of purifying selection in maintaining structural integrity and biological function.For the 15 scrutinized core PCGs, the rps3 gene had the highest observed Ka value, while the nad4L gene demonstrated the lowest Ka value, parallel to the K2P patterns.The synonymous substitution rate (Ks) of the cox1 gene was the highest, while that of the rps3 gene was the lowest among the 34 Agaricales species.The Ka/Ks values of 15 core PCGs were <1, indicating that these genes were under purifying selection pressure.In general, a higher rps3 ratio is more likely with greater species diversity in the Ka/Ks calculation.Intriguingly, the average Ka/Ks value for rps3 from the 34 Agaricales species was 0.82 (Figure 5, Table S10), which showed that the rps3 gene may not have been under positive selection and that it is still in the process of phylogenetic differentiation.In conclusion, our findings illuminate the evolutionary dynamics of core PCGs in Agaricales species and emphasize the pivotal role of purifying selection in maintaining structural integrity and biological function.
Position class (Pcl) is a term used to define the exact location of the coding region in the cox1 gene [41], typically used to characterize the positional information of introns contained in the cox1 gene.Two hundred thirty-six introns detected in the cox1 genes of 34 Agaricales mushrooms were categorized into 43 Pcls based on the sequence comparison of the cox1 gene with a reference species, Ganoderma calidophilum [42].Discrepancies in the class and quantity of introns in 34 species indicated the occurrence of gain or loss of introns.Pcls occurring in more than or equal to 20% of the 34 species were defined as common insertion sites, while the rest were classified as rare insertion sites.This study identified 15 common Pcls and 28 rare Pcls, of which P1305 was the most universally distributed intron detected in 20 of the 34 species.P612 ranked second and was found in 19 of the 34 species.Rare Pcls (P262, P547, P1193, P1458, P864, P193, P196, P278, P941, P385, P387, P1122, P371, P537, P668, and P821) were present in only 1 of the 34 mitogenomes (Figure 6).Interestingly, several rare Pcls, especially P941, P278, and P369, were identified in species with distant evolutionary relationships, such as Rhodotorula mucilaginosa, G.meredithae, Fomitiporia mediterranea, and G. calidophilum from Basidiomycota, suggesting possible gene transfer events.Position class (Pcl) is a term used to define the exact location of the coding region in the cox1 gene [41], typically used to characterize the positional information of introns contained in the cox1 gene.Two hundred thirty-six introns detected in the cox1 genes of 34 Agaricales mushrooms were categorized into 43 Pcls based on the sequence comparison of the cox1 gene with a reference species, Ganoderma calidophilum [42].Discrepancies in the class and quantity of introns in 34 species indicated the occurrence of gain or loss of introns.Pcls occurring in more than or equal to 20% of the 34 species were defined as common insertion sites, while the rest were classified as rare insertion sites.This study identified 15 common Pcls and 28 rare Pcls, of which P1305 was the most universally distributed intron detected in 20 of the 34 species.P612 ranked second and was found in 19 of the 34 species.Rare Pcls (P262, P547, P1193, P1458, P864, P193, P196, P278, P941, P385, P387, P1122, P371, P537, P668, and P821) were present in only 1 of the 34 mitogenomes (Figure 6).Interestingly, several rare Pcls, especially P941, P278, and P369, were identified in species with distant evolutionary relationships, such as Rhodotorula mucilaginosa, G.meredithae, Fomitiporia mediterranea, and G. calidophilum from Basidiomycota, suggesting possible gene transfer events.

Comparative Analysis of Cyathus Mitogenomes
A synteny analysis based on five Cyathus mitogenomes was conducted to investigate the sequence identity of the complete mitogenome.The results revealed that there were substantial regions of identical sequences in the mitogenomes of C. striatus 87405 and AH44044, as well as many regions of high sequence identity between C. pallidus QL1 and C. stercoreus NPCB004 (Figure 7A).Conversely, fewer identical sequences were identified in C. jiayuguanensis compared to the other four Cyathus mitogenomes (Figure 7A).The analysis indicated that C. pallidus QL1 and C. stercoreus NPCB004 share close intragroup evolutionary relationships, whereas C. jiayuguanensis is more distantly related.Furthermore, the protein sequence identity of the fourteen concatenated conserved genes between C. striatus 87405 and AH44044 exceeded 99%, while it was 97.47% between C. pallidus QL1 and C. stercoreus NPCB004.In contrast, C. jiayuguanensis displayed significantly lower sequence identity with other species, ranging from 81.24% to 84.92% (Figure 7B, Table S11).
in C. jiayuguanensis compared to the other four Cyathus mitogenomes (Figure 7A).The analysis indicated that C. pallidus QL1 and C. stercoreus NPCB004 share close intragroup evolutionary relationships, whereas C. jiayuguanensis is more distantly related.Furthermore, the protein sequence identity of the fourteen concatenated conserved genes between C. striatus 87405 and AH44044 exceeded 99%, while it was 97.47% between C. pallidus QL1 and C. stercoreus NPCB004.In contrast, C. jiayuguanensis displayed significantly lower sequence identity with other species, ranging from 81.24% to 84.92% (Figure 7B, Table S11).
Cyathus jiayuguanensis 765 harbored some un_ORFs homologous to those in C. striatus, including orf334, orf353, and orf358, which encode LAGLIDADG homing endonucleases.Two pairs of homologous genes (orf261 and orf261, orf389 and orf366) were distributed in C. pallidus QL1 and C. stercoreus NPCB004, all encoding LAGLIDADG homing endonucleases.Interestingly, homologous LAGLIDADG homing endonucleases, orf 273 and orf261, all existed in the cob gene of C. striatus 87405 and AH44044, C. pallidus QL1, and C. stercoreus NPCB004 (Figure 7C, Table S12).Therefore, it is speculated that these two un_ORFs are commonly present in the genus Cyathus.In addition, we observed varying counts of un_ORFs in different species:  Cyathus jiayuguanensis 765 harbored some un_ORFs homologous to those in C. striatus, including orf334, orf353, and orf358, which encode LAGLIDADG homing endonucleases.Two pairs of homologous genes (orf261 and orf261, orf389 and orf366) were distributed in C. pallidus QL1 and C. stercoreus NPCB004, all encoding LAGLIDADG homing endonucleases.Interestingly, homologous LAGLIDADG homing endonucleases, orf 273 and orf261, all existed in the cob gene of C. striatus 87405 and AH44044, C. pallidus QL1, and C. stercoreus NPCB004 (Figure 7C, Table S12).Therefore, it is speculated that these two un_ORFs are commonly present in the genus Cyathus.In addition, we observed varying counts of un_ORFs in different species: 4 in C. striatus 87405, 2 in C. striatus AH44044, 17 in C. jiayuguanensis 765, 9 in C. pallidus QL1, and 16 in C. stercoreus NPCB004.These results suggest that C. jiayuguanensis 765 may have undergone accelerated evolutionary changes compared to other analyzed Cyathus species, offering valuable insights into the evolutionary relationships and divergence among Cyathus species at the mitogenome level.

Comparative Analysis of Agaricales Mitogenomes
A comparative analysis of the mitochondrial genomes of 34 Agaricales species, including five Cyathus species, was conducted to understand their mitochondrial features better.The mitogenome sizes showed considerable variation, ranging from 43,328 bp (Asterophora parasitica) to 256,807 bp (Clavaria fumosa) (Figure 8, Table S1).The analysis encompassed five genomic components, namely un_ORFs, intronic regions, intergenic regions, RNA regions, and core PCGs, to determine their contribution to mitogenomic differences.The relative content of each component displayed substantial disparity among the mushroom species studied.Core PCGs and RNA regions occupied nearly 50% of mitogenomes in mushrooms such as Asterophora parasitica, Schizophyllum commune, and M. tenuissimus, while accounting for only one-eighth or less in others, including C. striatus 87405, Omphalotus japonicus, and Clavaria fumosa.Intergenic regions also contributed to mitogenome size variation, constituting 46% of C. jiayuguanensis 765 and 28% of Leucoagricus naucinus mitogenomes.Meanwhile, un_ORFs percentages ranged from 5% to 30% across species (Figure 8, Table S1).
Previous studies have demonstrated that repetitive sequences, particularly in noncoding regions, facilitate the recombination of mitogenomes, thereby increasing the likelihood of gene shuffling events [4].Therefore, the most plausible explanation for this phenomenon is that C. jiayuguanensis possessed the greatest proportion of repetitive (3.10%) and tandem repeat (19.95%) sequences relative to other Cyathus species (Tables S8 and S9).
Previous studies have demonstrated that repetitive sequences, particularly in noncoding regions, facilitate the recombination of mitogenomes, thereby increasing the likelihood of gene shuffling events [4].Therefore, the most plausible explanation for this phenomenon is that C. jiayuguanensis possessed the greatest proportion of repetitive (3.10%) and tandem repeat (19.95%) sequences relative to other Cyathus species (Tables S8 and S9).
Fungal mitochondrial genomes typically exhibit relatively high variation rates, characterized by disparities in ribosomal RNA sequences and specific protein-coding genes across different families.In light of this, we performed a phylogenetic analysis using the PCG, PCG12, PCGR, and PCG12R datasets derived from 34 Agaricales species, as opposed to the amino acids of conserved genes in the 111 Basidiomycete mitogenomes.After subjecting the datasets to heterogeneity analysis and base substitution saturation testing, we constructed eight phylogenetic trees.The Cyathus clade displayed identical phylogenetic topologies (C.jiayuguanensis + ((C.pallidus + C. stercoreus) + (C.striatus 87405 + C. striatus AH44044))) across all four datasets, strongly supported by high BPP (1) and BS (100) values.However, the phylogenetic relationships among Nidulariaceae, Strophariaceae, and Hydnangiaceae remained unresolved, as three ML trees and one BI tree supported (Nidulariaceae + other species) + (Strophariaceae + Hydnangiaceae), while other trees concurred with the tree based on 111 species and supported other species + (Nidulariaceae + (Strophariaceae + Hydnangiaceae)).Furthermore, discrepancies in topologies involving V. volvacea and L. sordlda were observed across eight trees, with (((((L.decastes + L. shimeji) + A. parasitica) + M. boudieri) + T. bakamatsutake + T. matsutake) + L. sordlda) + V. volvacea being the most dominant topology and possessing the highest node support (Figure 12, Tables S14 and S15).Based on these results, we confirm the evolutionary position of five Cyathus species and believe that the integrated mitochondrial gene data can serve as a robust and reliable genetic marker to characterize the phylogeny of Basidiomycetes.Fungal mitochondrial genomes typically exhibit relatively high variation rates, characterized by disparities in ribosomal RNA sequences and specific protein-coding genes across different families.In light of this, we performed a phylogenetic analysis using the PCG, PCG12, PCGR, and PCG12R datasets derived from 34 Agaricales species, as opposed to the amino acids of conserved genes in the 111 Basidiomycete mitogenomes.After sub- phariaceae + Hydnangiaceae)).Furthermore, discrepancies in topologies involving V. volvacea and L. sordlda were observed across eight trees, with (((((L.decastes + L. shimeji) + A. parasitica) + M. boudieri) + T. bakamatsutake + T. matsutake)) + L. sordlda) + V. volvacea being the most dominant topology and possessing the highest node support (Figure 12, Tables S14 and S15).Based on these results, we confirm the evolutionary position of five Cyathus species and believe that the integrated mitochondrial gene data can serve as a robust and reliable genetic marker to characterize the phylogeny of Basidiomycetes.

Discussion
Compared to what is observed in other eukaryotes, the size of fungal mitogenomes is usually not conserved, and the length of mitogenomes varies considerably even among species of the same genus [3], which may be caused by intron insertions, repetitive sequence distribution, and non-conserved genes caused by the horizontal transfer [47].The comparative analysis of the mitogenomes of 34 species of Agaricales showed that they exhibited some variation in size, ranging from 43,328 bp to 256,807 bp.The mitogenomes of five Nidulariaceae species were in the intermediate size range, smaller than those of A. bisporus and Amanita thiersii and larger than those of Armillaria sinapina and Tricholoma

Discussion
Compared to what is observed in other eukaryotes, the size of fungal mitogenomes is usually not conserved, and the length of mitogenomes varies considerably even among species of the same genus [3], which may be caused by intron insertions, repetitive sequence distribution, and non-conserved genes caused by the horizontal transfer [47].The comparative analysis of the mitogenomes of 34 species of Agaricales showed that they exhibited some variation in size, ranging from 43,328 bp to 256,807 bp.The mitogenomes of five Nidulariaceae species were in the intermediate size range, smaller than those of A. bisporus and Amanita thiersii and larger than those of Armillaria sinapina and Tricholoma bakamatsutake.Correlation analysis revealed that the sizes of the un_ORFs, HEGs, intronic region, and intergenic region (Figure 9A-D) all contributed significantly to the mitogenome length (correlation coefficients > 0.9), which suggested the presence of a large number of intron drift events and the acquisition and deletion of un_ORFs in these Agaricales species.Furthermore, comparative analyses revealed that the GC skew of the 34 species varied considerably.With the exception of Nidulariaceae, most of the species had positive GC skew, which reflected obvious interspecific differences.In animal, plant, and some fungal genomes, GC skew usually carries information on the location of outgoing transcriptional starts, whereas in the bacterial genome, GC skew is closely related to the positioning of the leading and lagging strands [48].However, there are few studies on mitochondrial GC skew in fungi.Considering the differences in the affinity of DNA double strands for RNA polymerase and related transcription factors, it is speculated that there is a certain relationship between the difference in the GC skew of fungal mitochondrial DNA and the transcription efficiency of mRNA.Therefore, the intrinsic link between GC-skew differences and replication and transcription for species of the order Agaricales, especially species of the family Nidulariaceae, needs to be further investigated and validated.
Correlation analysis based on the GC content and GC-rich codons showed that the coding genes showed a preference for GC-rich codons among the 34 Agaricales species (Figure 4A).Although the GC content of the 34 species varied widely, from 19.74% to 31.73% (Table S1), the codon bias of these Agaricales species showed some conservation (Figure 3B).The reason for this bias may be that changes in GC content occur more in non-coding regions than in coding regions.Furthermore, among the five Nidulariaceae species, not only are codon preferences highly conserved, but also encoded tRNA types are completely consistent.All Nidulariaceae species contain tRNAs for 20 natural amino acids, with trnM, trnL, and trns all having two copies, and all of their trnM anticodons are CAU (Figures S1-S7).Given the correlation between codon bias and tRNA abundance [49,50], it is speculated that Nidulariaceae species share similar gene expression patterns.
Mitochondrial gene rearrangements are commonly studied in animals, and these studies suggest that gene rearrangements may be the result of tandem repeat sequence duplication and random loss [51,52].However, the mechanism of fungal mitochondrial gene rearrangement is less clear.Some studies have shown that the accumulation of tandem repeats and the degree of dispersion of homing endonucleases and tRNAs will lead to changes in gene arrangement [4].Gene rearrangement analyses have identified large-scale gene rearrangement events in 17 genes in the order Agaricales species, and it is thought that such gene rearrangements are caused by these two types of factors.However, in the genus Cyathus, the gene order of C. stercoreus, C. pallidus, C. striatus 87405, and AH44044 is highly conserved.In contrast, the gene order of C. jiayuguanensis 765 differs from that of other Cyathus species.The reason may be that C. jiayuguanensis 765 contains the highest proportion of repetitive sequences (Tables S8 and S9).This finding suggests the complexity of gene rearrangement in the fungal mitogenome.
A phylogenetic tree constructed from the amino acid sequences of fourteen PCGs in the mitogenomes of 111 Basidiomycete species revealed well-supported evolutionary divergence (Figure 11).In this tree, C. striatus 87405 and AH44044 form one of the smallest sister groups, while C. stercoreus and C. pallidus form the other.This result is consistent with the findings of the synteny analysis (Figure 7A) for the five Nidulariaceae species.In the synteny analysis, C. striatus 87405 and AH44044 share a large number of sequence identity regions, and C. stercoreus and C. pallidus also share a large number of highly similar regions (Figure 7A).The consistent results obtained with various analytical methods indicate the dependability of the analytical results.Due to the degeneration of codons, amino acid-based phylogenetic trees cannot accurately capture differences at the codon level.Therefore, to gain a more comprehensive understanding of the phylogenetic topology of the Nidulariaceae fungi, we rationally adopted the Agaricales dataset to construct eight additional phylogenetic trees.Among them, the Agaricales species share a fixed topology with high support.However, it is worth noting that among the eight trees, half of the trees show that the Nidulariaceae species form a monophyletic group with the species of Strophariaceae and Hydnangiaceae, while the other half show that the species of Nidulariaceae form a parallel group with the species of the Strophariaceae and Hydnangiaceae.Furthermore, the topology of other members of the order Agaricales shows differences in different evolutionary trees, such as the evolutionary relationship between V. volvacea and L. sordlda.Five trees together support the topology of (((((L.decastes + L. shimeji) + A. parasitica) + M. boudieri) + T. bakamatsutake + T. matsutake) + L. sordlda) + V. volvacea (Figure 12), and these results were consistent with those of the phylogenetic tree based on fourteen conserved proteins (Figure 11).However, the evolutionary relationship between V. volvacea and L. sordlda was not the case in the other three phylogenetic trees.These findings provide new insights into the phylogenetic relationships of Agaricales fungi and highlight the importance of considering codon-level changes when constructing phylogenetic trees.

Sequence Analysis of Mitogenomes
In this extensive study, we analyzed the mitogenomes of 34 distinct Agaricales species, focusing on features such as nucleotide composition, codon usage bias, and evolutionary pressures.We assessed the GC content, GC skew, and AT skew utilizing a Python script, where GC content is defined as (G + C)/(G + C + A + T), GC skew as (G − C)/(G + C) and AT skew as (A − T)/(A + T).We determined the RSCU values of mitochondrial PCGs, core PCGs, and un_ORFs using an online tool (http://cloud.genepioneer.com:9929/#/login, accessed on 25 March 2023).We performed t-distributed stochastic neighbor embedding (t-SNE) analyses based on RSCU values with the Rtsne R package (v0.15)(https://github.com/jkrijthe/Rtsne,accessed on 25 March 2023).
To investigate the presence of intra-genomic duplications of large fragments or interspersed repeats within the five Cyathus mtDNAs, we performed BLASTn analysis (e-value 10 −10 ) by comparing the mitogenomes against themselves.We utilized Tandem Repeat Finder (v4.07b) (https://tandem.bu.edu/trf/submit_options, accessed on 26 March 2023) with default parameters to identify tandem repeats longer than 10 bp.
For 34 Agaricales species examined, the lengths of core PCGs, RNA regions, intergenic regions, intronic regions, HEGs, and un_ORFs were compared.A stacked bar graph illustrating the length and compositional distribution of the mitogenomes was generated using the ggplot2 package (v.3.4.2) (https://github.com/tidyverse/ggplot2,accessed on 25 March 2023).If a gene was situated within an intron, the intronic sequence length was calculated by subtraction of the size of the gene.Mitogenome size and the associations between its six compositions (core PCGs, RNA regions, intergenic regions, intronic regions, homing endonuclease genes HEGs, and un_ORFs) were evaluated using the cor.test() function in R. Previous reports were consulted for intron analysis of the cox1 genes of 34 species [41,42].

Conclusions
In this study, we de novo assembled and annotated five mitochondrial genomes from the genus Cyathus (C.striatus 87405, C. striatus AH44044, C. jiayuguanensis 765, C. pallidus QL1, and C. stercoreus NPCB004), generating high-quality data for future research.Comparative analysis of the mitogenomes of 34 Agaricales species revealed conservation in codon usage preferences and Ka/Ks ratios for most conserved genes.However, significant differences in mitochondrial genome length and intron insertion sites were observed.These genome length divergences were attributed to un_ORFs, HEGs, and intronic and intergenic regions.The Agaricales species exhibited a plethora of homing endonucleases and diverse intron sites, indicating that large-scale intron gain or loss events have occurred during evolution.In addition, gene rearrangements were evident among the 34 species, indicating substantial gene order variation, even within the same genera.Phylogenetic trees constructed using amino acids from fourteen conserved genes and four datasets (PCG, PCG12, PCGR, and PCG12R) revealed the evolutionary position of Cyathus fungi within Agaricales and confirmed their monophyly.In particular, subspecies 87405 and AH44044 of C. striatus showed high intersubspecific homology and sequence identity, whereas C. pallidus QL1 and C. stercoreus NPCB004 showed significant interspecific similarity.These results are consistent with the affinities shown in the five phylogenetic trees.
In conclusion, the mitogenome analyses not only reflect the presence of gene rearrangements and abundant intron dynamics in the mitogenome of the genus Cyathus, but also clarify the phylogenetic position of Cyathus within the family Nidulariaceae.These results provide valuable information for species identification in the genus Cyathus and offer new perspectives for understanding the evolutionary status of Cyathus species.

Figure 1 .
Figure 1.Comparison of the fruiting body morphology of Nidulariaceae and other Agaricales outside Nidulariaceae.

Figure 1 .
Figure 1.Comparison of the fruiting body morphology of Nidulariaceae and other Agaricales outside Nidulariaceae.

Figure 2 .
Figure 2. Circular maps of the mitogenomes of five Cyathus species.Genes are shown by blocks different colors.The doughnut chart represents different compositions of the total mitochondri genomes of the five Cyathus species.From the inner doughnut to the outer ring, C. stercore NPCB004, C. pallidus QL1, C. jiayuguanensis 765, C. striatus 87405, and C. striatus AH44044 are s quenced.The symbol * indicates that the corresponding gene contains an intron (introns).

Figure 2 .
Figure 2. Circular maps of the mitogenomes of five Cyathus species.Genes are shown by blocks of different colors.The doughnut chart represents different compositions of the total mitochondrial genomes of the five Cyathus species.From the inner doughnut to the outer ring, C. stercoreus NPCB004, C. pallidus QL1, C. jiayuguanensis 765, C. striatus 87405, and C. striatus AH44044 are sequenced.The symbol * indicates that the corresponding gene contains an intron (introns).

Figure 3 .
Figure 3. (A) Stacked column plots of the RSCU of the mitochondrial genomes of five Cyathus sp cies.(B) Heatmap of the RSCU of the mitogenomes of 34 Agaricales species.

Figure 3 .
Figure 3. (A) Stacked column plots of the RSCU of the mitochondrial genomes of five Cyathus species.(B) Heatmap of the RSCU of the mitogenomes of 34 Agaricales species.

Figure 4 .
Figure 4. (A) Pearson correlation matrix between the usage of AT-rich codons, GC-rich codons, AT content, and GC content in 34 Agaricales species (B) GC content of mitochondrial CDS and intronic regions from 34 Agaricales species.(C) T-SNE analysis of the RSCU in mitochondrial core PCGs and un_ORFs.

Figure 4 .
Figure 4. (A) Pearson correlation matrix between the usage of AT-rich codons, GC-rich codons, AT content, and GC content in 34 Agaricales species (B) GC content of mitochondrial CDS and intronic regions from 34 Agaricales species.(C) T-SNE analysis of the RSCU in mitochondrial core PCGs and un_ORFs.

Figure 5 .
Figure 5. Genetic analysis of 15 protein-coding genes conserved in the 34 Agaricales mitogenomes.K2P, Kimura 2-parameter distance; Ka, the mean number of non-synonymous substitutions per non-synonymous site; Ks, the mean number of synonymous substitutions per synonymous site.

Figure 5 .
Figure 5. Genetic analysis of 15 protein-coding genes conserved in the 34 Agaricales mitogenomes.K2P, Kimura 2-parameter distance; Ka, the mean number of non-synonymous substitutions per non-synonymous site; Ks, the mean number of synonymous substitutions per synonymous site.

Figure 6 .
Figure 6.Intron insertion sites in the cox1 gene of 34 Agaricales species.The phylogenetic tree was constructed using amino acids from fourteen conserved PCGs.

Figure 6 .
Figure 6.Intron insertion sites in the cox1 gene of 34 Agaricales species.The phylogenetic tree was constructed using amino acids from fourteen conserved PCGs.
4 in C. striatus 87405, 2 in C. striatus AH44044, 17 in C. jiayuguanensis 765, 9 in C. pallidus QL1, and 16 in C. stercoreus NPCB004.These results suggest that C. jiayuguanensis 765 may have undergone accelerated evolutionary changes compared to other analyzed Cyathus species, offering valuable insights into the evolutionary relationships and divergence among Cyathus species at the mitogenome level.

Figure 10 .
Figure 10.Gene order analyses for 34 Agaricales mitogenomes.The same gene is represented by the same background color.All genes (fifteen core PCGs and two ribosomal RNAs) are shown in the order of their appearance in the mitogenome, starting with cox1.

Figure 10 .
Figure 10.Gene order analyses for 34 Agaricales mitogenomes.The same gene is represented by the same background color.All genes (fifteen core PCGs and two ribosomal RNAs) are shown in the order of their appearance in the mitogenome, starting with cox1.

23 Figure 11 .
Figure 11.Phylogeny of 111 Basidiomycetes species based on amino acids of fourteen conserved PCGs.Values before the slash at the clades indicate Bayesian posterior probabilities (BPPs), while those after the slash at the clades indicate bootstrap (BS).The asterisk at the clade indicates a BPP value of 1 and a BS value of 100.The area where the Cyathus species are located is marked with an orange background.

Figure 11 .
Figure 11.Phylogeny of 111 Basidiomycetes species based on amino acids of fourteen conserved PCGs.Values before the slash at the clades indicate Bayesian posterior probabilities (BPPs), while those after the slash at the clades indicate bootstrap (BS).The asterisk at the clade indicates a BPP value of 1 and a BS value of 100.The area where the Cyathus species are located is marked with an orange background.