Comparative Mitogenomic Analysis Reveals Intraspecific, Interspecific Variations and Genetic Diversity of Medical Fungus Ganoderma

Ganoderma species are widely distributed in the world with high diversity. Some species are considered to be pathogenic fungi while others are used as traditional medicine in Asia. In this study, we sequenced and assembled four Ganoderma complete mitogenomes, including G. subamboinense s118, G. lucidum s37, G. lingzhi s62, and G. lingzhi s74. The sizes of the four mitogenomes ranged from 50,603 to 73,416 bp. All Ganoderma specimens had a full set of core protein-coding genes (PCGs), and the rps3 gene of Ganoderma species was detected to be under positive or relaxed selection. We found that the non-conserved PCGs, which encode RNA polymerases, DNA polymerases, homing endonucleases, and unknown functional proteins, are dynamic within and between Ganoderma species. Introns were thought to be the main contributing factor in Ganoderma mitogenome size variation (p < 0.01). Frequent intron loss/gain events were detected within and between Ganoderma species. The mitogenome of G. lucidum s26 gained intron P637 in the cox3 gene compared with the other two G. lucidum mitogenomes. In addition, some rare introns in Ganoderma were detected in distinct Basidiomycetes, indicating potential gene transfer events. Comparative mitogenomic analysis revealed that gene arrangements also varied within and between Ganoderma mitogenomes. Using maximum likelihood and Bayesian inference methods with a combined mitochondrial gene dataset, phylogenetic analyses generated identical, well-supported tree topologies for 71 Agaricomycetes species. This study reveals intraspecific and interspecific variations of the Ganoderma mitogenomes, which promotes the understanding of the origin, evolution, and genetic diversity of Ganoderma species.


Sequencing, Assembly, and Annotation of Mitogenomes
Four Ganoderma specimens were collected from Sichuan Province in southwest China, which were then stored in Culture Collection Center of Sichuan Academy of Agricultural Sciences (Contact information: Wenli Huang, wenlih11@126.com), including G. subamboinense s118, G. lucidum s37, G. lingzhi s62, and G. lingzhi s74. The growth environment information and collection information of these specimens can be found in Supplementary Table S1. We extracted the genomic DNA of the four Ganoderma specimens using a fungal DNA extraction kit (Omega Bio-Tek, Norcross, GA, USA). Sequencing libraries were constructed with the genomic DNA using a NEBNext ® Ultra™ II DNA Library Prep Kit (NEB, Beijing, China) according to the instructions. We conducted whole genomic sequencing (WGS) using the Illumina HiSeq 2500 Platform (Illumina, San Diego, CA, USA). In order to process the raw data, a series of quality control steps were carried out. Among these steps, ngsShoRT was used for filtering low-quality sequences [46] and AdapterRemoval v2 was used to remove adapter reads [47]. We de novo assembled the four Ganoderma mitogenomes using NOVOPlasty v3.7.2 [48] with a k-mer size of 25. The four mitogenomes obtained were annotated using the MFannot tool [49] and MITOS [50] according to our previously described methods [51]. We drew the circular maps of the four Ganoderma mitogenomes using OGDRAW [52].

Intron Analysis
Fungal mitochondrial introns vary greatly, and usually do not contain any introns in most eukaryotic mitogenomes while basidiomycetes, including Polyporales, usually contain varying numbers of introns [58,59]. We further separated the 15 core PCG introns in Ganoderma and other Polyporales mitogenomes (the accession numbers are shown in Table S7) into position classes (Pcls) with the mitogenome of Ganoderma calidophilum [25] as a reference in accordance with reported methods [60]. First, the core PCGs excluding introns were aligned with reference mitochondrial sequences of G. calidophilum using Clustal W [61]. Each Pcl comprised introns inserted at the same position in the PCG coding regions [60]. Introns that are of the same Pcl typically have high sequence similarities and are deemed orthologous [62]. Most Pcls that differ possess low sequence similarities and non-orthologous mobile genetic elements [45]. Pcls of 15 core PCGs in Polyporales were given names based on the reference host gene's coding region insert position.

Phylogenetic Analysis
A phylogenetic tree of 71 Agaricomycetes species was constructed based on the combined mitochondrial gene set in order to explore the phylogenetic relationships of Ganoderma species and the phylogenetic status of these Ganoderma species among the Agaricomycetes class (the accession numbers are shown in Table S7). The mitogenome of Hannaella oryzae from the order Tremellales was designated as an outgroup [63]. The combined mitochondrial gene set contains 14 core protein-coding gene sequences, excluding the rps3 gene and intron regions. For alignment of single mitochondrial genes, we used MAFFT v7.037 [64]. We subsequently used SequenceMatrix v1.7.8 to concatenate the resulting alignments into a gene set [65]. We employed PartitionFinder 2.1.1 to identify the gene set's best-fit evolution and partitioning scheme models [66]. MrBayes v3.2.6 was used [67] along with the Bayesian inference (BI) method based on the combined gene set to evaluate the phylogenetic relationships between Agaricomycetes species. BI analysis involved two independent runs with three heated chains and one cold chain, each performed simultaneously for 2 × 10 6 generations. Once every 1000 generations, each run was sampled. Stationarity was assumed to have been achieved when the estimated sample size (ESS) reached more than 100, and the potential scale reduction factor (PSRF) became close to 1.0 (the closer the PSRF value was to 1, the better the convergence effect). We eliminated the first 25% of samples as burn-in. Then, we made use of the remaining trees to evaluate the Bayesian posterior probabilities (BPPs) in a 50% majority-rule consensus tree [68]. In addition, the maximum likelihood (ML) technique was employed for assessment of the phylogenetic relationships of 71 Agaricomycetes species using RAxML v8.0.0 [69] with the combined gene set. We performed an evaluation of the bootstrap values (BS) using an ultrafast bootstrap approach with 1000 replicates.

Characterization of the Four Ganoderma Mitogenomes
The four complete Ganoderma mitogenomes were all made up of circular DNA molecules that ranged in size from 50,603 to 73,416 bp (Figure 1), and the molecular structure was further verified by polymerase chain reaction amplification and pyrophosphate sequencing. The mitogenome of G. subamboinense s118 was the largest while the G. lingzhi s74 was the smallest among the four Ganoderma mitogenomes. The GC content of the four mitogenomes was within 26.33-26.56%, with a 26.41% average GC content (Table S1). The mitogenomes of G. subamboinense s118 and G. lucidum s37 exhibited positive AT skews while the two strains of G. lingzhi had negative AT skews. The GC skews of the four Ganoderma mitogenomes were all positive. In total, 23 to 28 non-intronic open-reading frames (ORFs) were detected in the 4 Ganoderma mitogenomes, which included 15 core PCGs and 8-12 free-standing ORFs. These free-standing ORFs mainly encoded DNA polymerases, RNA polymerases, and proteins with unknown functions in the four Ganoderma mitogenomes (Table S2). A total of 48 introns were detected in the 4 Ganoderma mitogenomes, with each containing 8-21 introns. Each intron contained one or two intronic ORFs, encoding LAGLIDADG homing endonucleases or GIY-YIG homing endonucleases. A total of 49 intronic ORFs were detected in the 4 Ganoderma mitogenomes. Among the 48 introns detected in the 4 Ganoderma specimens, 29 belonged to Group IB, 7 belonged to Group IA, 6 belonged to Group ID, 4 belonged to Group IC2, 1 belonged to Group II, and 1 was an unknown type. Two rRNA genes were found in all four Ganoderma mitogenomes, namely, the small subunit ribosomal RNA (rns) and the large subunit ribosomal RNA (rnl). In addition, the 4 Ganoderma mitogenomes all contained 26 tRNA genes.

Codon Usage Analysis
Up to now, nine complete mitochondrial genomes of Ganoderma species have been published in the database (the accession numbers are shown in Table S7). We further compared these 13 Ganoderma mitogenomes to reveal the conservation and variation between Ganoderma mitogenomes. Of the 15 core PCGs in the 13 Ganoderma mitogenomes we detected, 14 used ATG as start codons except for cob, which used GTG as start codons in all the 13 Ganoderma species (Table S3). The majority of the core PCGs in thee Ganoderma species used TAA as stop codons while the atp6 gene of G. calidophilum, G. applanatum, and G. subamboinense s118, and the cox2 gene of G. calidophilum used TAG as stop codons. The nad1 and nad2 genes of most Ganoderma species used TAG as stop codons.
Codon usage analysis revealed that the most frequently used codons were AAA in the G. subamboinense s118 and G. lingzhi s74 mitogenomes, TTT in G. lucidum s37, and TTA in the G. lingzhi s62 mitogenome ( Figure 2 and Table S4). In general, AAA (for lysine; Lys), TTA (for leucine; Leu), TTT (for phenylalanine; Phe), AAT (for asparagine; Asn), and ATT (for isoleucine; Ile) were the five most frequently used codons in the four newly sequenced Ganoderma mitogenomes. The high AT content in the four Ganoderma mitogenomes (average: 73.59%) was due to frequently used codons ending in A and T.

Codon Usage Analysis
Up to now, nine complete mitochondrial genomes of Ganoderma species have been published in the database (the accession numbers are shown in Table S7). We further compared these 13 Ganoderma mitogenomes to reveal the conservation and variation between Ganoderma mitogenomes. Of the 15 core PCGs in the 13 Ganoderma mitogenomes we detected, 14 used ATG as start codons except for cob, which used GTG as start codons in all the 13 Ganoderma species (Table S3). The majority of the core PCGs in thee Ganoderma species used TAA as stop codons while the atp6 gene of G. calidophilum, G. applanatum, and G. subamboinense s118, and the cox2 gene of G. calidophilum used TAG as stop codons. The nad1 and nad2 genes of most Ganoderma species used TAG as stop codons.

Repetitive Sequences Analysis
We identified 47, 9, 7, and 7 repeat sequences in the mitogenomes of G. subamboinense s118, G. lucidum s37, G. lingzhi s62, and G. lingzhi s74, respectively, by comparison of the whole mitogenomes against themselves via BLASTn analysis (Table S5). In the four Ganoderma mitogenomes, repeat sequence lengths ranged from 35 to 868 bp while the range of pairwise nucleotide similarities was 66.07-100%. The orf410, orf211 and orf168 protein-coding regions had the largest repeats, and in the intergenic region between orf211 and orf168 in the G. subamboinense mitogenome. The repetitive sequences comprised 2.33-12.12% of the whole mitogenomes of the four Ganoderma specimens. The G. subamboinense s118 mitogenome exhibited the highest proportion of repeat sequences while the lowest proportion of repeat sequences was found in the G. lingzhi s74 mitogenome.
Codon usage analysis revealed that the most frequently used codons were AAA in the G. subamboinense s118 and G. lingzhi s74 mitogenomes, TTT in G. lucidum s37, and TTA in the G. lingzhi s62 mitogenome ( Figure 2 and Table S4). In general, AAA (for lysine; Lys), TTA (for leucine; Leu), TTT (for phenylalanine; Phe), AAT (for asparagine; Asn), and ATT (for isoleucine; Ile) were the five most frequently used codons in the four newly sequenced Ganoderma mitogenomes. The high AT content in the four Ganoderma mitogenomes (average: 73.59%) was due to frequently used codons ending in A and T.

Repetitive Sequences Analysis
We identified 47, 9, 7, and 7 repeat sequences in the mitogenomes of G. subamboinense s118, G. lucidum s37, G. lingzhi s62, and G. lingzhi s74, respectively, by comparison of the whole mitogenomes against themselves via BLASTn analysis (Table S5). In the four Ganoderma mitogenomes, repeat sequence lengths ranged from 35 to 868 bp while the range of pair-wise nucleotide similarities was 66.07-100%. The orf410, orf211 and orf168 protein-coding regions had the largest repeats, and in the intergenic region between orf211 and orf168 in the G. subamboinense mitogenome. The repetitive sequences comprised 2.33-12.12% of the whole mitogenomes of the four Ganoderma specimens. The G. subamboinense s118 mitogenome exhibited the highest proportion of repeat sequences while the lowest proportion of repeat sequences was found in the G. lingzhi s74 mitogenome.
A total of ten, four, two, and four tandem repeats were found in the mitogenomes of G. subamboinense s118, G. lingzhi s74, G. lucidum s37, and G. lingzhi s62, respectively (Table  S6). The mitogenome of G. subamboinense s118 contained the longest tandem repeat sequence at 45 bp. Among the four newly sequenced mitogenomes, most of the tandem repeats were duplicated once or twice, with the highest replication number (12) in the mitogenome of G. subamboinense s118. Tandem repeat sequences made up 0.14-0.49% of the four Ganoderma mitogenomes. A total of ten, four, two, and four tandem repeats were found in the mitogenomes of G. subamboinense s118, G. lingzhi s74, G. lucidum s37, and G. lingzhi s62, respectively (Table S6). The mitogenome of G. subamboinense s118 contained the longest tandem repeat sequence at 45 bp. Among the four newly sequenced mitogenomes, most of the tandem repeats were duplicated once or twice, with the highest replication number (12) in the mitogenome of G. subamboinense s118. Tandem repeat sequences made up 0.14-0.49% of the four Ganoderma mitogenomes.

Variation, Genetic Distance, and Evolutionary Rates of Core Genes
Of the 15 core PCGs revealed herein, the nad3 gene and then the nad6 gene showed the greatest mean K2P genetic distance between the 13 Ganoderma species (Figure 3). On average, the atp8 gene exhibited the smallest K2P genetic distance between the 13 Ganoderma species, showing the atp8 gene was highly conserved. The nad3 gene exhibited the greatest mean nonsynonymous substitution rate (Ka) while the atp8 and atp9 genes had the lowest Ka values. The nad3 gene revealed the highest synonymous substitution rate (Ks) among the 15 core PCGs while that of the rps3 gene was the lowest. The Ka/Ks values of 14 out of 15 core PCGs were less than 1, demonstrating that these genes underwent purifying selection. However, the Ka/Ks values of the rps3 gene had values greater than one between some species, such as between G. calidophilum and G. lucidum s26, between G. subamboinense s118 and G. leucocontextum, and between G. meredithae and G. sinense, indicating that in some Ganoderma species, the rps3 gene was under either positive or relaxed selection. rate (Ks) among the 15 core PCGs while that of the rps3 gene was the lowest. The Ka/Ks values of 14 out of 15 core PCGs were less than 1, demonstrating that these genes underwent purifying selection. However, the Ka/Ks values of the rps3 gene had values greater than one between some species, such as between G. calidophilum and G. lucidum s26, between G. subamboinense s118 and G. leucocontextum, and between G. meredithae and G. sinense, indicating that in some Ganoderma species, the rps3 gene was under either positive or relaxed selection.

Intron Dynamic Analysis
A significant correlation was found between the number of introns and the mitogenome sizes of 13 Ganoderma species (p < 0.01). The Pearson and Spearman correlation coefficients were 0.923 and 0.906, respectively, indicating that introns are the main factor causing the size variations of Ganoderma mitogenomes ( Figure 4). A total of 235 introns were detected in the 13 Ganoderma mitogenomes we tested, with each containing 8-31 introns. The dynamic changes in the introns of Ganoderma species indicated that the loss or gain of introns occurred during the evolution of Ganoderma species. In total, 215 out of the

Intron Dynamic Analysis
A significant correlation was found between the number of introns and the mitogenome sizes of 13 Ganoderma species (p < 0.01). The Pearson and Spearman correlation coefficients were 0.923 and 0.906, respectively, indicating that introns are the main factor causing the size variations of Ganoderma mitogenomes ( Figure 4). A total of 235 introns were detected in the 13 Ganoderma mitogenomes we tested, with each containing 8-31 introns. The dynamic changes in the introns of Ganoderma species indicated that the loss or gain of introns occurred during the evolution of Ganoderma species. In total, 215 out of the 235 introns (91.49%) were distributed in core PCGs of the 13 Ganoderma species, indicating that the core PCGs were the largest host genes of introns in Ganoderma. The remaining 20 introns were distributed in rRNA genes. The introns of core PCGs were harbored in thee cox1, cox2, cox3, cob, nad1, nad2, nad3, nad4, and nad5 genes, and the cox1 gene contained the largest number of introns, with 114. Six core PCGs, including the atp6, atp8, atp9, nad4L, nad6, and rps3 genes, did not contain any intron in the thirteen Ganoderma mitogenomes. The uneven distribution of introns in different core PCGs indicated that fungal introns had host gene preference. Since most of Ganoderma introns are distributed in the core PCGs, we studied the dynamic changes in the introns of Ganoderma core PCGs. introns were distributed in rRNA genes. The introns of core PCGs were harbored in thee cox1, cox2, cox3, cob, nad1, nad2, nad3, nad4, and nad5 genes, and the cox1 gene contained the largest number of introns, with 114. Six core PCGs, including the atp6, atp8, atp9, nad4L, nad6, and rps3 genes, did not contain any intron in the thirteen Ganoderma mitogenomes. The uneven distribution of introns in different core PCGs indicated that fungal introns had host gene preference. Since most of Ganoderma introns are distributed in the core PCGs, we studied the dynamic changes in the introns of Ganoderma core PCGs. Using the insertion sites of introns in the protein-coding region of host genes as a basis, introns were categorized into different position classes (Pcls) depending on their corresponding reference genes. The 114 introns of the cox1 gene could be divided into 24 Pcls in the 13 Ganoderma mitogenomes ( Figure 5). P1305 was the Pcl that was most broadly distributed; it was distributed in 12 of the 13 Ganoderma species. P209 was also widely distributed in Ganoderma cox1 genes, which was distributed in 10 of the 13 Ganoderma species. Some rare intron was only distributed in 1 of the 13 Ganoderma cox1 genes, including P218, P309, P900, P971, and P1107. We also detected nine, three, two, three, one, one, two, and five Pcls in the cob, cox2, cox3, nad1, nad2, nad3, nad4, and nad5 genes. P426 from the nad5 gene was the most common intron in Ganoderma species, which was present in all of the 13 Ganoderma mitogenomes. P234 from the cox2 gene was also widely distributed in Ganoderma species, which was distributed in 10 of the 13 Ganoderma mitogenomes. However, P201, P597, and P600 from the cob gene; P453 and P543 from the cox2 gene; P775 from the cox3 gene; P153 and P657 from the nad1 gene; P320 from the nad3 gene; and P675 from Using the insertion sites of introns in the protein-coding region of host genes as a basis, introns were categorized into different position classes (Pcls) depending on their corresponding reference genes. The 114 introns of the cox1 gene could be divided into 24 Pcls in the 13 Ganoderma mitogenomes ( Figure 5). P1305 was the Pcl that was most broadly distributed; it was distributed in 12 of the 13 Ganoderma species. P209 was also widely distributed in Ganoderma cox1 genes, which was distributed in 10 of the 13 Ganoderma species. Some rare intron was only distributed in 1 of the 13 Ganoderma cox1 genes, including P218, P309, P900, P971, and P1107. We also detected nine, three, two, three, one, one, two, and five Pcls in the cob, cox2, cox3, nad1, nad2, nad3, nad4, and nad5 genes. P426 from the nad5 gene was the most common intron in Ganoderma species, which was present in all of the 13 Ganoderma mitogenomes. P234 from the cox2 gene was also widely distributed in Ganoderma species, which was distributed in 10 of the 13 Ganoderma mitogenomes. However, P201, P597, and P600 from the cob gene; P453 and P543 from the cox2 gene; P775 from the cox3 gene; P153 and P657 from the nad1 gene; P320 from the nad3 gene; and P675 from the nad4 gene were only detected in 1 of the 13 Ganoderma species, which were considered to be rare Pcls in Ganoderma. We found that these rare Pcls were also distributed in distant species, such as Paxillus rubicundulus [51], Agaricus bisporus [62], Rhizopogon salebrosus [31], Russula compacta [37], and Moniliophthora perniciosa [70], suggesting that potential intron transfers may take place in Ganoderma mitogenomes, or that intron insertions were convergent in distantly related species.
We observed loss or gain of introns within G. lucidum and G. lingzhi species. The G. lucidum s37 lost P926 in the cox1 gene, and P324 in the nad5 gene compared with the other two G. lucidum mitogenomes. The mitogenome of G. lucidum s26 gained P637 in the cox3 gene compared with the other two G. lucidum mitogenomes. Within the G. lingzhi species, the G. lingzhi s62 lost P209 in the cox1 gene, P539 in the nad4 gene, and P324 in the nad5 gene compared with the other two G. lingzhi mitogenomes.

Comparative Mitogenomic Analysis
Comparative genomic analysis showed that the mitogenome of G. calidophilum was the largest among the 13 Ganoderma mitogenomes while G. lingzhi and G. lucidum had the smallest mitogenome among the 13 Ganoderma species, with an average size of 56,383 and 57,304, respectively (Table S1). The mitogenome size also varied within the G. lingzhi and G. lucidum species. The GC content was also found to vary between and within Ganoderma species, indicating base composition dynamics in the evolution of Ganoderma species. Base biases were frequently detected in the mitogenomes of Ganoderma. The AT and GC skews of G. leucocontextum and AT skew of G. lucidum (KC763799) were negative (T over A; C over G) in their mitogenomes. However, there was an excess of Gs (but not Cs) and As We observed loss or gain of introns within G. lucidum and G. lingzhi species. The G. lucidum s37 lost P926 in the cox1 gene, and P324 in the nad5 gene compared with the other two G. lucidum mitogenomes. The mitogenome of G. lucidum s26 gained P637 in the cox3 gene compared with the other two G. lucidum mitogenomes. Within the G. lingzhi species, the G. lingzhi s62 lost P209 in the cox1 gene, P539 in the nad4 gene, and P324 in the nad5 gene compared with the other two G. lingzhi mitogenomes.

Comparative Mitogenomic Analysis
Comparative genomic analysis showed that the mitogenome of G. calidophilum was the largest among the 13 Ganoderma mitogenomes while G. lingzhi and G. lucidum had the smallest mitogenome among the 13 Ganoderma species, with an average size of 56,383 and 57,304, respectively (Table S1). The mitogenome size also varied within the G. lingzhi and G. lucidum species. The GC content was also found to vary between and within Ganoderma species, indicating base composition dynamics in the evolution of Ganoderma species. Base biases were frequently detected in the mitogenomes of Ganoderma. The AT and GC skews of G. leucocontextum and AT skew of G. lucidum (KC763799) were negative (T over A; C over G) in their mitogenomes. However, there was an excess of Gs (but not Cs) and As (but not Ts) in the replication leading strands of most Ganoderma mitogenomes. Different Ganoderma species contained varied numbers of non-conserved ORFs, encoding DNA polymerase, RNA polymerase, homing endonucleases, and unknown functional proteins ( Table 1). The six G. lucidum and G. lingzhi strains all contained the orf109, which encoded RNA polymerase. The orf283 encoding DNA polymerase could only be detected in G. lucidum. Some non-conserved PCGs only existed in the mitogenome of G. lingzhi, including orf115, orf250, orf744, and orf879 encoding DNA polymerases; orf202 and orf211 encoding RNA polymerases; orf118 and orf309 encoding homing endonucleases; and orf152 and orf413 encoding unknown functional proteins. In G. lucidum, strain s37 lost orf151, orf192, and orf294, which encoded unknown functional proteins, while strain s26 lost orf283, which encoded DNA polymerase. In G. lingzhi, the orf118 encoding GIY-YIG endonuclease was lost in strain s8 while the strain s26 gained the orf115 and orf250 encoding DNA polymerases, and orf202 and orf211 encoding RNA polymerases compared with other G. lingzhi species. These results indicate that loss or gain events of non-conserved PCGs have occurred in the evolution of Ganoderma. The position classes and number of introns also varied within and between Ganoderma species. In addition, 29.41% of the intronic ORFs in the G. sinense mitogenome have been lost, indicating that intron ORFs of G. sinense are undergoing contraction. The number of tRNA in Ganoderma mitogenomes varied between Ganoderma species but was conserved within Ganoderma species. All Ganoderma species contained two mitochondrial rRNA genes.
identical gene arrangements, which were arranged in the following order: cox1, nad4, atp6  rnl, cox3, nad4L, rps3, nad2, nad5, rns, nad6, atp9, nad1, cob, cox2, atp8, and nad3, while the  mitogenomes of G. subamboinense s118, G. lucidum s37, G. lingzhi s62, and G. lingzhi s74 had  another gene order: cox1, nad4, atp6, rnl, cox3, nad4L, nad5, rns, cox2, rps3, nad6, atp8, nad2  nad3, atp9, nad1, and cob, which were different from other Ganoderma species. Gene rear rangements within G. lucidum and G. lingzhi were observed, including gene transfer, in sertion, and inversion events. The findings showed that large-scale gene rearrangements took place during the evolution of Ganoderma species. Figure 6. Gene order comparation between 13 Ganoderma species. All genes are shown in order o occurrence in the mitochondrial genome, starting from cox1. Fifteen core protein-coding genes and two rRNA genes were included in the gene arrangement analysis. The phylogenetic positions of 13 Ganoderma species were established using the Bayesian inference (BI) method and maximum likeli hood (ML) method based on concatenated mitochondrial genes. Species and NCBI accession num ber used for gene arrangement analysis in the present study are listed in Supplementary Table S7. We performed phylogenetic analysis using maximum likelihood (ML) and Bayesian inference (BI) according to the combined mitochondrial gene set (14 core PCGs). These methods generated identical, well-supported tree topologies (Figure 7). All major clades within the phylogenetic tree were well supported (BPP ≥ 0.99; BS ≥ 98). According to the phylogenetic tree, the 71 Agaricomycetes species comprised 6 major clades of the orders Agaricales, Boletales, Russulales, Polyporales, Hymenochaetales, and Cantharellales (Table S7) The 24 Polyporales species could be divided into 7 groups, corresponding to the families Phanerochaetaceae, Meruliaceae, Fomitopsidaceae, Laetiporaceae, Sparassidaceae, and Polypo raceae while Taiwanofungus could not be placed in any recognized family. Ganoderma and Tremetes are closely related and form a group Polyporaceae. The phylogenetic analyses showed that G. lingzhi was a sister species to G. lucidum. Figure 6. Gene order comparation between 13 Ganoderma species. All genes are shown in order of occurrence in the mitochondrial genome, starting from cox1. Fifteen core protein-coding genes and two rRNA genes were included in the gene arrangement analysis. The phylogenetic positions of 13 Ganoderma species were established using the Bayesian inference (BI) method and maximum likelihood (ML) method based on concatenated mitochondrial genes. Species and NCBI accession number used for gene arrangement analysis in the present study are listed in Supplementary Table S7. We performed phylogenetic analysis using maximum likelihood (ML) and Bayesian inference (BI) according to the combined mitochondrial gene set (14 core PCGs). These methods generated identical, well-supported tree topologies (Figure 7). All major clades within the phylogenetic tree were well supported (BPP ≥ 0.99; BS ≥ 98). According to the phylogenetic tree, the 71 Agaricomycetes species comprised 6 major clades of the orders Agaricales, Boletales, Russulales, Polyporales, Hymenochaetales, and Cantharellales (Table S7). The 24 Polyporales species could be divided into 7 groups, corresponding to the families Phanerochaetaceae, Meruliaceae, Fomitopsidaceae, Laetiporaceae, Sparassidaceae, and Polyporaceae while Taiwanofungus could not be placed in any recognized family. Ganoderma and Tremetes are closely related and form a group Polyporaceae. The phylogenetic analyses showed that G. lingzhi was a sister species to G. lucidum.

Mitogenome Size Variation in Ganoderma
In the present study, the mitogenome size of Ganoderma species was found to vary greatly, ranging from 50,603 to 124,588 bp. The largest mitogenome G. calidophilum [25] in the genus Ganoderma is 2.46 times larger than the smallest one G. lingzhi s62. Previous studies have found that the fungal mitogenome sizes is one of the most variable groups in eukaryotes, which is mainly due to the dynamic changes in introns, the accumulation of repeat sequences, the content of non-conserved PCGs, and so on [71][72][73][74]. The number of introns in this study had a significant effect on the size variation of Ganoderma mitogenomes (p < 0.01). This suggests that introns may be the main contributing factor affecting the size variation of Ganoderma species. The mitogenomes of G. lucidum and G. lingzhi were the smallest among Ganoderma species, which indicated that the two species may have experienced intron loss events during evolution. We observed size variations of the mitogenomes within G. lucidum (7289 bp variations) and G. lingzhi (10,712 bp variations) species. Within the G. lucidum species, G. lucidum (KC763799) contained the largest mitogenome. However, the mitogenome of G. lucidum (KC763799) did not have the largest number of introns or non-conserved PCGs. Through comparative mitogenomic analysis, we found that the expansion or contraction of the intergenic region may explain the size variations within G. lucidum mitogenomes. In G. lingzhi species, the expansion of non-conserved PCGs can explain the size increase of the G. lingzhi s74 mitogenome. In general, the size variations of the Ganoderma mitogenomes are caused by several factors, including the dynamic changes in the introns, the expansion or contraction of the intergenic region, and the accumulation of non-conserved PCGs.

Mitogenome Content Evolution in Ganoderma
Since the ancestors of eukaryotes acquired mitogenomes from Alphaproteobacteria through endosymbiosis, most mitochondrial genes have been integrated into the nuclear genome during evolution [75,76]. However, most fungal mitogenomes retain a set of core PCGs and a number of non-conserved PCGs to ensure the stability of the oxidative phosphorylation process. Partial deletion of core PCGs was also observed in some fungal mitogenomes, indicating the complex evolution of fungal mitogenomes [45]. In this study, we found that all the mitogenomes of Ganoderma contained a whole set of core PCGs. The length, base composition, codon usage frequency, and start and stop codons of these core PCGs varied greatly within and between Ganoderma species, indicating that the core PCGs of Ganoderma mutated frequently. Interestingly, we found that the rps3 gene demonstrated positive selection or relaxed selection between some Ganoderma species. This phenomenon has also been observed in other fungal species [77,78], and what causes the selection pressure on the rps3 gene needs to be further verified. Different numbers of non-conserved PCGs were observed in the mitogenomes of Ganoderma, which encoded RNA polymerase, DNA polymerase, homing endonuclease, and unknown functional proteins. Two closely related Ganoderma species, G. lucidum and G. lingzhi, both contained orf109 encoding RNA polymerase, which indicated that orf109 may play an important role in the mitochondrial function of the two Ganoderma species. Interspecific differences in the non-conserved PCGs were detected between the two closely related Ganoderma species. All the three mitogenomes of G. lucidum contained orf396, orf429, and orf568, which encoded DNA polymerase, unknown functional protein, and RNA polymerase, respectively. However, G. lingzhi s26 lost the three non-conserved PCGs. All G. lingzhi contained an orf151 gene, which encoded unknown functional protein, while G. lucidum s37 lost this gene. The dynamic changes in the RNA polymerase, DNA polymerase, homing endonuclease, and unknown functional genes were frequently observed within G. lucidum and G. lingzhi species, indicating that the loss and gain of non-conserved PCGs occurred in the evolution of Ganoderma mitogenomes, even within species.

Dynamics of Introns in Ganoderma
The dynamic change in the introns is thought to be the primary factor causing the size variation of the Ganoderma mitogenomes. Compared with mitogenomes from other genera, most of the intronic ORFs in Ganoderma species were retained, and only G. sinense lost 29.41% of the intronic ORFs [25]. It was found that introns were unevenly distributed in the Ganoderma core PCGs and rRNA genes, and the cox1 gene was the largest host gene of introns in Ganoderma. Introns can be categorized into different Pcls depending on their insertion sites, and the same Pcls were considered to be homologous [79]. There was no significant correlation between the number of Pcls and the length of host genes, indicating that intron insertion sites had sequence selection preference. P209 and P1305 of the cox1 gene, P234 of the cox2 gene, and P426 of the nad5 gene were considered to be the most common Pcls in Ganoderma, which were distributed in more than 10 of the 13 Ganoderma species. However, some rare Pcls were only detected in 1 of the 13 Ganoderma species, including P201, P597, and P600 from the cob gene; P453 and P543 from the cox2 gene; P775 from the cox3 gene; P153 and P657 from the nad1 gene; P320 from the nad3 gene; and P675 from the nad4 gene. These rare Pcls were found in distant species, which indicated possible horizontal gene transfer events. In previous studies, introns were found to be horizontally transferred between different organelles [80]. We also detected intron loss and gain events in G. lucidum and G. lingzhi. The results show that the introns of Ganoderma are highly variable within and between species, and the physiological and functional effects of these intron dynamics need to be further revealed.

Gene Rearrangements and Phylogenetic Analysis
This research work found that the Ganoderma mitochondrial gene arrangement varied greatly within and between species. Some species with distant phylogenetic relationships had identical gene arrangements while some closely related species had varied gene arrangements, suggesting that the mitochondrial gene arrangement in Ganoderma species is both random and dynamic. Previous studies found that gene rearrangements occur frequently in fungal mitogenomes [81]. However, the fungal gene rearrangement mechanism remains unknown. Animal mitochondrial gene rearrangement has been investigated in depth, and researchers have designed multiple models to explore the mechanism of animal mitochondrial gene rearrangement [82,83]. The accumulation of repetitive sequences in fungal mitogenomes was believed to cause gene rearrangements in fungal mitogenomes [84]. However, the result cannot explain the intraspecific gene rearrangements of Ganoderma mitogenomes. More intraspecific mitochondrial gene rearrangements and their mechanisms need to be elucidated to understand the evolutionary patterns of fungal mitogenomes.
The Ganoderma species is widely distributed on Earth, with high species diversity [1,85]. Some Ganoderma species are pathogenic fungi of trees while some are traditional Asian medicinal materials [2,86]. The discovery of bioactive natural products from Ganoderma has attracted the enthusiasm of medical researchers all over the world and promoted the continuous expansion of the Ganoderma cultivation scale [87,88]. However, the limited and overlapped morphological features led to confusion in the accurate classification of Ganoderma species [5,19]. The phenomenon of synonyms and homonyms is very prevalent in cultivation and product development of Ganoderma species, which limits the large-scale development and utilization of this important medicinal fungus [20,21]. The mitogenome is an effective tool to understand the phylogenetic relationship of species [89][90][91]. However, the lack of a mitochondrial reference genome of Agaricomycetes, especially Ganoderma species, limits the application of the mitogenome in the classification and phylogenetic relationship analysis of Ganoderma. In this study, we newly obtained four complete mitogenomes of Ganoderma species. Based on the ML and BI phylogenetic inference methods, we obtained phylogenetic trees for 71 Agaricomycetes species with high support rates in major clades. The phylogenetic relationships among the Ganoderma species and the phylogenetic status of Ganoderma in Agaricomycetes were also revealed. The phylogenetic tree of the 24 Polyporales species based on the mitogenome is highly consistent with that based on polygenic molecular markers [92][93][94]. This study provides useful reference data for the classification and identification of Ganoderma species and promotes our understanding of intraspecific and interspecific variations of Ganoderma species.

Conclusions
In the present study, the mitogenomes of four Ganoderma specimens, including G. subamboinense s118, G. lucidum s37, G. lingzhi s62, and G. lingzhi s74, were sequenced and assembled. We performed comparative mitogenomic analysis of intraspecific and interspecific Ganoderma species to reveal the conservation and variation of Ganoderma mitogenomes. The sizes of Ganoderma mitogenomes varied within and between Ganoderma species, and the intron was the main factor contributing to the size variations of Ganoderma mitogenomes. Frequent intron loss/gain events were detected within and between Ganoderma species. The mitogenome of G. lucidum s26 gain intron P637 in cox3 gene compared with the other two G. lucidum mitogenomes. The core PCG rps3 of the Ganoderma species experienced positive or relaxed selection while non-conserved PCGs are dynamic within and between Ganoderma species. Large-scale gene arrangements were detected even within Ganoderma species. The results showed that there was high mitogenome diversity within and between Ganoderma species, regarding intron types, non-conserved genes, and gene arrangement. This study revealed intraspecific and interspecific variations of Ganoderma mitogenomes for the first time, which promotes the understanding of the evolution and genetic diversity of Ganoderma species.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/jof8080781/s1. Table S1: Comparison on mitogenomes among Ganoderma species; Table S2: Characterization of the mitogenome of four Ganoderma species; Table S3: Start and stop codon analyses of 15 core protein-coding genes in 13 Ganoderma species; Table S4: Codon usage analysis of the four Ganoderma mitogenomes; Table S5: Local BLAST analysis of the four Ganoderma mitogenomes against themselves; Table S6: Tandem repeats detected in the four Ganoderma mitogenomes using the online program Tandem Repeats Finder; Table S7: Species, and GenBank accession number used for phylogenetic analysis in this study. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: The four Ganoderma mitogenomes, including G. subamboinense s118, G. lucidum s37, G. lingzhi s62, and G. lingzhi s74, were submitted to GenBank under the accession numbers MW752412, MW752414, MW752415, and MW752413, respectively. All data generated or analyzed during this study are included in this published article and its Supplementary Materials.