Comparative mitogenome analysis of two ectomycorrhizal fungi (Paxillus) reveals gene rearrangement, intron dynamics, and phylogeny of basidiomycetes

In this study, the mitogenomes of two Paxillus species were assembled, annotated and compared. The two mitogenomes of Paxillus involutus and P. rubicundulus comprised circular DNA molecules, with the size of 39,109 bp and 41,061 bp, respectively. Evolutionary analysis revealed that the nad4L gene had undergone strong positive selection in the two Paxillus species. In addition, 10.64 and 36.50% of the repetitive sequences were detected in the mitogenomes of P. involutus and P. rubicundulus, respectively, which might transfer between mitochondrial and nuclear genomes. Large-scale gene rearrangements and frequent intron gain/loss events were detected in 61 basidiomycete species, which revealed large variations in mitochondrial organization and size in Basidiomycota. In addition, the insertion sites of the basidiomycete introns were found to have a base preference. Phylogenetic analysis of the combined mitochondrial gene set gave identical and well-supported tree topologies, indicating that mitochondrial genes were reliable molecular markers for analyzing the phylogenetic relationships of Basidiomycota. This study is the first report on the mitogenomes of Paxillus, which will promote a better understanding of their contrasted ecological strategies, molecular evolution and phylogeny of these important ectomycorrhizal fungi and related basidiomycete species.


INTRODUCTION
Most boreal and north temperate forest tree species form mutualistic symbiotic relationships with ectomycorrhizal fungi (ECMF), which are of great significance to the balance and stability of forest ecosystems, carbon and nitrogen cycles, and the growth and reproduction of the fungi and host tree.
The genus Paxillus is an important ectomycorrhizal fungal group, widely distributed in the northern hemisphere. Paxillus species form ectomycorrhizal associations with various host plant species, including coniferous and hardwood trees (Hedh et al. 2009). It was reported that Paxillus can enhance the abiotic stress tolerance and promote the growth of host plants (Franco et al. 2014;Li et al. 2012;Ma et al. 2014). Paxillus involutus, the most studied species in the genus, has become one of the model organisms for revealing the evolution, life-style and physiology of ECMF (Rineau et al. 2012;Shah et al. 2013).
Studies of Paxillus nuclear genomes have found that species of this genus have a reduced complement of genes encoding plant cell wall-degrading enzymes (PCWDEs) and retain a unique array of PCWDEs to adapt to mutualistic symbiotic life-styles (Kohler et al. 2015). Some nuclear genes have also evolved to better adapt to such life-styles (Kohler et al. 2015;Rajashekar et al. 2007). However, up to now, their mitochondrial genetic information has not yet been elucidated, limiting our overall understanding of the genetics and evolution of Paxillus species.
The genus Paxillus was divided into two groups in Europe, the P. involutus complex and P. rubicundulus. Paxillus rubicundulus occurs mainly in wetland habitats and only forms a relationship with Alnus species (Jargeat et al. 2016), while the P. involutus complex occurs in more diverse habitats and forms mutualistic relationships with various tree species. Therefore, host specificity could be used as a useful feature to distinguish P. rubicundulus (Jargeat et al. 2016), but accurately identifying species in the P. involutus complex remains difficult (Jargeat et al. 2014). Phylogenetic studies have revealed four genetic lineages within the P. involutus complex (Jargeat et al. 2014), and in addition there are some species whose taxonomic status remains confusing. This presents difficulties for efficient utilization and research into those Paxillus species. Mitochondrial genes have been reported as a powerful tool for the study of phylogeny due to several advantages (Delsuc et al. 2003;Nie et al. 2019), but up to now no mitochondrial genome (mitogenome) from the genus Paxillus, or even the family Paxillaceae, has been available.
The mitogenome of eukaryotic organisms was reported to have been obtained from the common ancestral Alphaproteobacteria though endosymbiosis (Lang et al. 1999). In basidiomycete species, they generally contain 15 core protein coding genes (PCGs), including 14 genes for energy metabolism (atp6, atp8, atp9, cob, cox1, cox2, cox3, nad1, nad2, nad3, nad4, nad4L, nad5, and nad6) and 1 rps3 gene for transcriptional regulation. In addition, 20-36 tRNA genes and 2 rRNA genes were also detected in basidiomycete mitogenomes (Li et al. 2018a;Li et al. 2018d). The mitogenome has become a powerful tool for the study of phylogeny due to several available molecular markers and uniparental inheritance (Andersen and Balding 2018;Wang et al. 2018). With the rapid development of next-generation sequencing technology, more and more nuclear and mitochondrial genomes of eukaryotic organisms have been obtained (Tajima et al. 2016). However, compared with their animal counterparts, the mitogenomes of fungi have been less studied (Sandor et al. 2018), especially that of basidiomycetes in the basidiomycete class Agaricomycetes (the largest mushroom-forming fungal group). The limited studies that have been carried out on fungal mitogenomes have, however, shown that these vary greatly in genome size, gene content, gene arrangement, intron number, and content of repeat sequences, even between closely related species (Li et al. 2018c;Li et al. 2019b).
In this study, the mitogenomes of two Paxillus species, P. involutus and P. rubicundulus, were assembled and annotated in order to: (1) reveal features of Paxillus mitogenomes and the variations or similarities between those of the two species in genome size, gene content, gene arrangement, and repeat sequence; (2) demonstrate any dynamic changes of introns and gene rearrangement in the two Paxillus mitogenomes compared with other basidiomycete mitogenomes to contribute to a more comprehensive understanding of their size and organization variation in basidiomycetes; and (3) analyze the phylogenetic relationships of Basidiomycota based on the combined mitochondrial gene set. The mitogenomes of the two Paxillus species could further our understanding of the taxonomy, evolution, and potential molecular mechanisms driving host specificity in an important ectomycorrhizal.

Mitochondrial genome assembly and annotations
The raw sequencing data of Paxillus involutus and P. rubicundulus were downloaded from the Joint Genome Institute (JGI) database (Kohler et al. 2015). A series of quality control steps were conducted to obtain clean reads from the raw sequencing data. First, we removed adapter reads using AdapterRemoval v2 (Schubert et al. 2016), and then ngsShoRT (Chen et al. 2014) to filter low-quality sequences. The Paxillus mitogenomes were assembled with the resulting clean reads using SPAdes 3.9.0 (Bankevich et al. 2012), and gaps between contigs were filled using MITObim V1.9 (Hahn et al. 2013). The complete mitogenomes obtained were first annotated according to our previously described methods (Li et al. 2018a;Li et al. 2018d). That is, the protein-coding genes (PCGs), rRNA genes and tRNA genes of the mitogenomes were initially annotated using MFannot (Valach et al. 2014) and MITOS (Bernt et al. 2013), based on the genetic code 4. PCGs were then predicted or modified with the NCBI Open Reading Frame Finder (Coordinators 2017), and further annotated by BLASTP searches against the NCBI non-redundant protein sequence database (Bleasby and Wootton 1990). tRNA genes were also predicted and identified with tRNAscan-SE v1.3.1 (Lowe and Chan 2016). Graphical maps of the two complete mitogenomes were drawn with OGDraw v1.2 (Lohse et al. 2013).

Repetitive element analysis
To identify if there were interspersed repeats or intragenomic duplications of large fragments throughout the two mitogenomes, we conducted BLASTN searches (Chen et al. 2015) of each mitogenome against itself using an E-value of < 10 − 10 . Tandem repeats (> 10 bp in length) in the two mitogenomes were detected using the Tandem Repeats Finder (Benson 1999) with default parameters. Repeated sequences in the two mitogenomes were also searched by REPuter (Kurtz et al. 2001) to identify forward (direct), reverse, complemented, and palindromic (reverse complemented) repeats. To identify any gene fragments transferring between nuclear genomes and mitochondrial genomes of the two Paxillus species, we performed BLASTn searches of the two mitogenomes against their previously published nuclear genomes (JOMD00000000.1 and JMDR00000000.1).

Comparative mitogenomic analysis and intron analysis
Up to now, only two mitogenomes from Boletales have been published (Li et al. 2019a). We compared the genome sizes, base compositions, gene numbers, intron numbers, and gene arrangements among different Boletales species to assess variations and conservativeness among the mitogenomes. Group I introns in cox1 genes of the 61 published Basidiomycota mitogenomes were classified into different position classes (Pcls) according to the method of Férandon et al. (Ferandon et al. 2010). Each Pcl was constituted by introns inserted at the same position in the coding region of the cox1 gene. The same Pcl from different species usually has a high sequence similarity, and contains orthologous intronic ORF. The Pcls of cox1 gene were named by letters according to the similarity with the described Pcls (Ferandon et al. 2010). The conservativeness and variations of sequences around intron insertion sites (− 15 bp -15 bp) were evaluated using WebLogo (Crooks et al. 2004).

Phylogenetic analysis
In order to investigate the phylogenetic status of the two Paxillus species within the Basidiomycota, we constructed a phylogenetic tree of 62 species based on the combined mitochondrial gene set (15 core PCGs + 2 rRNA genes) (Li et al. 2018d). Annulohypoxylon stygium (Ascomycota) was used as outgroup (Deng et al. 2018). Individual mitochondrial genes were first aligned using MAFFT v7.037 (Katoh et al. 2017), and these alignments were then concatenated in SequenceMatrix v1.7.8 (Vaidya et al. 2011) to obtain a combined mitochondrial gene set. Potential phylogenetic conflicts between different genes were detected by a partition homogeneity test; PartitionFinder 2.1.1 (Lanfear et al. 2017) was used to determine best-fit models of evolution and partitioning schemes. Both Bayesian Inference (BI) and Maximum Likelihood (ML) methods were used to construct phylogenetic trees. MrBayes v3.2.6 (Ronquist et al. 2012) was used to perform the BI analysis, and RAxML v 8.0.0 (Stamatakis 2014) the ML analysis.

Data availability
The complete mitogenomes of Paxillus involutus and P. rubicundulus were deposited in GenBank under accession numbers MK993563 and MK993564, respectively.

Genome features and PCGs of Paxillus mitogenomes
The mitogenomes of both species were both composed of circular nucleotide molecules, of 39,109 bp and 41, 061 bp, respectively (Fig. 1). The average GC content of the two mitogenomes was 21.39%. Both species contained negative AT skews and positive GC skews (Table  S1).
There were 19 and 17 PCGs detected in the P. involutus and P. rubicundulus mitogenomes, respectively (Table S2), of which 78.95-82.35% were located in the direct strand. Both the mitogenomes contained 15 core PCGs, including 14 core PCGs for energy metabolism and one rps3 gene for transcriptional regulation. In addition, the P. involutus mitogenome contained a PCG located in the second intron of the cox1 gene, which encoded LAGLIDADG endonuclease. Paxillus rubicundulus contained a PCG encoding DNA polymerase. Three and one PCGs with unknown functions were detected in the P. involutus and P. rubicundulus mitogenomes, respectively. P. involutus and P. rubicundulus contained a homologous PCG showing no significant similarity with known protein sequences in the public databases checked; this encoded 732 amino acids in the P. involutus mitogenome and 792 amino acids in that of P. rubicundulus.

rRNA genes and tRNA genes
Both the Paxillus mitogenomes contained two rRNA genes, the large subunit ribosomal RNA (rnl) and small subunit ribosomal RNA (rns) genes (Table S2). The length of the rRNA genes varied between the two: there were five nucleotide length variations between the two rnl genes, and nine between the two rns genes.
The two mitogenomes both contained 25 tRNA genes, which encode 20 standard amino acids (Table S2). Both contained two tRNAs that code for arginine and serine with different anticodons, two coding for leucine, and three coding for methionine with the same anticodons. All tRNAs in the two mitogenomes were folded into a classical cloverleaf structures (Fig. S1), with each 71-86 bp long. Three tRNA genes, including trnY, trnS and trnL contained large extra arms, which resulted in their tRNA gene length exceeding 85 bp. Of the 25 tRNA genes shared by the two Paxillus mitogenomes, 10 contained sites that varied between the two. Thirteen variable sites were detected in the 25 tRNAs between the two species, of which 53.85% occurred on the acceptor arm, suggesting that the acceptor arm was highly variable in the two Paxillus mitogenomes.

Intergenic sequence and mitogenome composition
Both the mitogenomes contained two overlapping regions (Table S2), one of which was located between the nad2 and nad3 genes (13 bp). Another overlapping region in P. involutus was located between cox3 and orf732 (46 bp), and a further one in P. rubicundulus was across the neighboring cox3 and orf792 genes (88 bp). The two Paxillus mitogenomes lacked the common overlapping region found in some other Basidiomycota mitogenomes located between the nad4L and nad5 genes. A total of 11,295 bp and 13,790 bp intergenic sequences were detected in the P. involutus and P. rubicundulus mitogenomes, respectively, with each intergenic sequence ranging from 9 bp to 3537 bp. The longest intergenic sequence detected was between cob and cox3 in the P. rubicundulus mitogenome.
Protein coding regions accounted for the largest proportion (40.92-43.60%) of the two Paxillus mitogenomes, followed by the intergenic region, which accounted for 20.88-33.58% of the whole mitogenomes (Fig. 2). RNA regions, including rRNA and tRNA, accounted for 19.41-20.38% of the two mitogenomes. The intronic region accounted for the smallest proportion of the two mitogenomes, with only 6.09-7.13%. Compared with the P. involutus mitogenome, the size of that of P. rubicundulus was expanded by 1952 bp, to which the intergenic region contributed 127.82% of the size expansion, while the other regions contracted in the P. rubicundulus mitogenome relative to rgose in that of P. involutus, indicating that the intergenic region was the most important factor contributing to the dynamic changes of mitogenomes in the two Paxillus species.

Repetitive elements in the two Paxillus mitogenomes
A total of 53 and 47 tandem repeats were detected in the P. involutus and P. rubicundulus mitogenomes, respectively (Table S3). The longest tandem repeat sequence was located in the P. rubicundulus mitogenome, with a length of 49 bp. Most of the tandem repeats contained 2-4 copies in the two mitogenomes. Using REPuter, we detected 1 complemented, 24 forward, 4 palindromic, and 21 reverse repeats in the P. involutus mitogenome, accounting for 3.13% of the whole mitogenome (Table S4). Repeated sequences identified by REPuter accounted for 3.25% of the P. rubicundulus mitogenome.
In order to detect whether there were gene fragments that naturally transferred between the mitochondrial and nuclear genomes of the two Paxillus species, we compared the two newly sequenced mitogenomes with their published nuclear genomes (Kohler et al. 2015) by BLASTN. Eleven repetitive fragments were identified between the mitochondrial and nuclear genomes of P. involutus, ranging from 33 bp to 3202 bp in length (Table S5). The pair-wise nucleotide similarities of these repetitive fragments ranged from 85.31-100%. The repetitive sequences detected in the P. involutus mitogenome accounted for 10.64% of the entire mitogenome. A total of 62 repetitive gene fragments were detected between the P. rubicundulus mitochondrial and nuclear genomes, with pair-wise nucleotide similarities ranging from 87.93-100%. The largest repeat fragment in the P. rubicundulus mitogenome reached 2367 bp, which contained a partial sequence of the rnl gene, the whole sequence of the cox2 gene, and the intergenic sequences around them. The repetitive gene fragments between P. rubicundulus mitochondrial and nuclear genomes accounted for 36.50% of the entire mitogenome.

Mitochondrial gene arrangement in Basidiomycota species
The mitochondrial gene arrangement in Basidiomycota varied greatly at the order level (Fig. 3). Of the 13 orders in Basidiomycota with data, no identical arrangement of mitochondrial genes was found between any two orders, suggesting that large-scale rearrangements of Basidiomycota mitochondrial genes had occurred in the course of evolution. At the level of family, only some species from Paxillaceae and Rhizopogonaceae were found to have the same gene order. Gene orders in other families were highly variable. At the generic level, we found large-scale gene rearrangement events in aevearl genera, including Cantharellus, Ganoderma, Lyophyllum, Pleurotus, and Russula. However, the gene arrangement of the two species of Paxillus studied was found to be conservative.
Variation, genetic distance, and evolutionary rates of core PCGs Among the 15 core PCGs detected, three genes showed length variations between the two Paxillus species, including cox3, nad2 and rps3 (Fig. 4). Among them, cox3 and nad2 genes showed two amino acid length variations between the two, while rps3 had one amino acid length variation between them. The GC content of 13 of the 15 core PCGs differed between the two Paxillus species, except for that of the atp6 and nad1 genes.
Among the 15 core PCGs tested, atp9 had the highest GC content and atp8 the lowest. Most AT skews of core PCGs were negative, except for the rps3 gene. Eleven of the 15 core PCGs showed positive GC skews, while atp6, atp8, nad2 and nad6 had negative GC skews. Among the 15 core PCGs detected, rps3 contained the largest K2P genetic distance between the two species (Fig. 5), followed by cox3 and nad2, which suggested that these genes were largely differentiated in the Paxillus species. The K2P genetic distance of atp8 was the lowest among the 15 core PCGs between the two mitogenomes, indicating that this gene was highly conserved in the two species. The cox3 gene exhibited the highest synonymous substitution (Ks) rate, while nad4L had the lowest substitution rate among the 15 core PCGs detected. The highest non-synonymous substitution rate (Ka) was in cox3, while atp9, cob and nad3 exhibited the lowest Ka values among the 15 core PCGs. The Ka/Ks values for most core PCGs were less than 1, indicating that these genes were subjected to purifying selection. However, the Ka/Ks value of nad4L gene was far greater than 1, indicating that this gene has been subjected to strong positive selection in the two Paxillus species. A nonsynonymous mutation was detected in the nad4L gene, of which the 46th nucleotide 'G' in P. involutus mutated into nucleotide 'A' in P. rubicundulus.

Intron dynamics of cox1 gene in Basidiomycota
According to previous studies (Ferandon et al. 2013), introns can be divided into different Pcls according to their precise insertion sites in the protein coding region of PCGs. Introns belonging to the same Pcls were considered as homologous and often contained homologous intronic open reading frames (ORFs). Introns belonging to different Pcls were considered as non-homologous introns, and showed low sequence similarities. In the present study, 367 introns were detected in cox1 genes A total of 32 position classes (Pcls) of group I introns were detected in the 61 species (Fig. 6), including five novel Pcls that had not been identified in previous studies, and 27 known Pcls. Pcls K was the most common intron class, which occurred in 36 of the species. Pcls P, AD, AI, N, and S were present in more than 30% of the species tested, and considered as the common Pcls in Basidiomycota. Pcls E was detected only in the cox1 gene of Ganoderma leucocontextum (MH252534) (Li et al. 2019e). However, it was present in the cox1 gene of the distantly related aquatic species Blastocladiella emersonii (Tambor et al. 2008) which belongs to the phylum Blastocladiomycota, a situation not readily explained as horizontal gene transfer between fungi which grow in such different environments seems improbable. In the two Paxillus species, we detected three group I introns, which belonged to different Pcls, indicating great variation of introns between them. The conservativeness and variation of nucleic acid sequences (30 bp) around insertion sites of different Pcls was also assessed (Fig. 7). We found that the insertion site preferences of different Pcls were quite different, while the nucleic acid sequences around the same Pcls were relatively conservative. Over 13 of the 30 nucleic acid loci detected were highly consistent around the same Pcls. Of the 32 group I intron Pcls we detected, 28 were inserted downstream of base T, while 11 Pcls were inserted behind GGT. In addition, 14 Pcls of 32 group I introns were located before base A, indicating a preference of intron insertion sites. These results could be useful for accurately identifying intron classes and annotating mitochondrial genes with multiple introns.

Comparative mitogenome analysis and phylogenetic analysis
Compared with previously published Rhizopogon mitogenomes, the two Paxillus species contained smaller mitogenomes (Table S1). That of P. involutus was 97.16% smaller than the largest mitogenome, R. vinicolor. The two Paxillus species contained a higher proportion of protein coding regions, RNA regions, and intergenic regions, but a lower proportion of intronic regions relative to the Rhizopogon species. Amongst the four Boletales species we tested, they all contained negative AT skews and positive GC skews. The two Paxillus species contained less PCGs, intron and intronic ORFs than Rhizopogon species. However, the four species from Boletales had the same number of rRNA genes and tRNA genes. Phylogenetic analysis using Maximum likelihood (ML) and Bayesian Inference (BI) methods based on the combined mitochondrial gene set (14 core PCGs + rps3 + 2 rRNA genes) yielded identical and well-supported tree topologies (Fig. 8). All major clades within the trees were well supported (BS ≥ 97; BPP ≥0.98). According to the phylogenetic tree, the 61 Basidiomycota species included could be divided into 13 major clades, corresponding to the orders Agaricales, Boletales, Cantharellales, Microbotryales, Microstromatales, Polyporales, Pucciniales, Russulales, Sporidiobolales, Tilletiales, Tremellales, Trichosporonales, and Ustilaginales (Table S6). Boletales had a close phylogenetic relationship with Agaricales and Russulales, and the four Boletales species could be divided into two groups, with P. involutus as a sister species to P. rubicundulus.

Genome size variation in Boletales
The mitochondrial genome size of fungi is highly variable, the largest being in Rhizoctonia solani (235.85 kb; (Losada et al. 2014) and the smallest in Spizellomyces punctatus (1.14 kb; (Forget et al. 2002). Intronic and intergenic regions, the content of repeat sequence and horizontally transferred genes have been found to be the main factors contributing to size dynamic changes in fungal mitogenomes (Losada et al. 2014;Salavirta et al. 2014). So far, only four complete mitogenomes from Boletales have been published, and these also varied in size. The two Rhizopogon species contained larger mitogenomes than those from Paxillus, was mainly due to their larger intronic regions. In Paxillus, the intergenic region was the main factor contributing to the size variation of the mitogenomes. Intronic and intergenic regions together contribute most to variations in mitogenome size in Boletales, as well as in other Basidiomycota so far examined (Ferandon et al. 2010).

Gene content variation in Paxillus species
Since eukaryotic organisms first acquired a mitogenome from endosymbiotic bacteria, most mitochondrial genes have been transferred to the nuclear genome during evolution, which is considered to have many advantages (Adams and Palmer 2003). However, most fungal mitochondria have so far been found to retain 14 core PCGs for energy metabolism and one rps3 gene for transcriptional regulation (Allen 2015;Bjorkholm et al. 2015). We found that the core PCGs in the two Paxillus species varied in base composition and length. In addition, the two core PCGs had different genetic distances between the species, indicating that different genes evolve at different rates. The nad4L gene is an important functional protein for ATP synthesis in mitochondria, which plays an important role in maintaining cell homeostasis and regulating celluar responses to the environment (Kazama et al. 2001). Interestingly, we found that the nad4L gene had experienced a strong positive selection between the two Paxillus species, which might be related to their host specificity. In addition, several PCGs coding homing endonucleases, DNA polymerase, and several unknown ORFs, were found in the species, indicating that some PCGs remain to be revealed. The tRNA and rRNA genes shared by the two Paxillus species also underwent base and length variations. Base variation in tRNA has been thought to affect the efficiency of amino acid transport and protein synthesis in animals (Ding et al. 2019;Lin et al. 2019). Further studies are needed to reveal the effect of tRNA variation on the growth and development of basidiomycete species.

Potential gene transfer between nuclear genome and mitogenome
The synergistic effects of nuclear and mitochondrial genes maintain the growth and development of eukaryotes (Baris et al. 2017;Latorre-Pellicer et al. 2016). During the evolution of eukaryotes, most mitochondrial genes have been transferred to the nucleus (Adams and Palmer 2003). However, some nuclear genes have also been found transferred to the mitogenome, and this has affected the evolution of eukaryotes (Li et al. 2019f;Srirattana and St John 2018). In our study, 10.64 and 36.50% of gene fragments were found potentially transferred between nuclear and mitochondrial genomes of the two Paxillus, respectively, which is a higher proportion than reported from Ganoderma lucidum (Li et al. 2013) and Pleurotus eryngii (Li et al. 2018a). Our result indicates that frequent natural gene transfer has occurred in the  two Paxillus species. The impact of this transfer on the evolution and function of the Paxillus mitogenomes, however, remains unknown.

Gene rearrangement of basidiomycete species
The evolutionary rate of mitochondrial gene in fungi has been reported as intermediate between that of animals and plants (Aguileta et al. 2014;Barr et al. 2005). As more mitochondrial genomes have been sequenced, gene rearrangements found to be widespread in animals, plants, and fungi (Park et al. 2016;Sankoff et al. 1992). Various gene rearrangement models have been proposed to reveal gene rearrangement events in animals (Lavrov et al. 2002;Xia et al. 2016), but mitochondrial arrangements in fungi have been much less studied. We found that the mitochondrial gene arrangement in Basidiomycota was highly variable. At the order level, no two orders had the same gene order. In addition, we also observed large-scale gene rearrangements at the family or even generic level (Li et al. 2018a;Li et al. 2018b;Li et al. 2019b;Li et al. 2019c).
The Basidiomycota mitogenome contains more repetitive sequences than animals, which we consider to be one of the main factors resulting in the high variation seen in the mitochondrial genome; this result is consistent with the findings of previous studies (Aguileta et al. 2014).

Dynamic changes of introns in cox1 gene of Basidiomycota
As a mobile genetic element (Burke 1988;Repar and Warnecke 2017), introns have been considered one of the main factors contributing to the variations of mitogenome organization and size (Hamari et al. 2002;Kanzi et al. 2016;Repar and Warnecke 2017). In Basidiomycota, the cox1 gene was found to harbour the most introns, with the largest proportion belonging to group I . We detected 32 Pcls of the group I introns and 3 Pcls of the group II introns amongst the 61 species for which data were available. Among these intron Pcls, five were newly identified here, indicating that the phylum contains diverse intron types. Pcl K was the most widely distributed Pcl in Basidiomycota, and so was presumed to be derived from the phylum's ancestors. Pcl E was only detected in Ganoderma leucocontextum (Li et al. 2019e), but this has also been reported in the distantly related Blastocladiella emersonii (Tambor et al. 2008), something suggestive of horizontal gene transfer, but difficult to imagine in view of the different ecologies of the two phyla. Different Basidiomycota species, even from within the same genera, could contain different numbers and Pcls of introns, indicating frequent intron loss or acquisition events during the evolution of Fig. 7 The insertion sites of different Pcls in the coding region of 61 Basidiomycete cox1 genes and the conservative analysis of the sequence around the insertion sites (− 15 bp -15 bp). The symbols "+ 1" and "+ 2" indicate that the insertion of the intron occurs inside the indicated codon: between the nt 1 and nt 2 of this codon for "+ 1" and between the nt 2 and nt 3 for "+ 2" Li et al. IMA Fungus (2020) 11:12 Page 10 of 15 the phylum. Introns contained homing endonucleases that can insert into specific sites in the mitogenome could potentially be used as gene editing tools (Chi et al. 2018;Zubaer et al. 2018). We assessed the conservativeness and variability of nucleic acid sequences around insertion sites of introns from basidiomycetes where information was available which showed that the insertion sites had certain base preferences. Our results provide a reference source for the accurate annotation of basidiomycete mitogenomes and also indicate a potential application of homing endonucleases in gene editing.

Phylogenetic relationships of Basidiomycota based on mitochondrial genes
Because of their independence from the nuclear genome, their rapid evolutionary rate, and the presence of many molecular markers, the mitogenome has become a powerful tool for studying phylogeny, population genetics, and taxonomy (Cameron 2014;Li et al. 2015;Li et al. 2019d). With the rapid development of next generation sequencing technology, more and more mitogenomes have been acquired, which has promoted the development of eukaryotic systems biology (Johri et al. 2019;Wang et al. 2018). However, the mitogenomes of fungi have been less studied than those of animals. Up to now, few mitochondrial molecular markers have been employed in studying their population genetics or phylogeny . We obtained identical and wellsupported tree topologies by BI and ML revealing phylogenetic relationships amongst 61species of Basidiomycota, based on the combined mitochondrial gene set. This demonstrates that the mitochondrial gene is a reliable molecular marker for the analysis of phylogenetic relationships of basidiomycetes and perhaps other fungal groups More mitogenomes are, however, needed to comprehensively assess the origin and evolution of the Basidiomycota.

CONCLUSION
In this study, the mitogenomes of two Paxillus species were assembled, annotated, and analyzed. Comparative mitogenome analysis indicated that the intronic and intergenic regions were the main factors contributing to mitogenome size variations in Boletales. The genome content, gene length, base composition, tRNAs, and rRNAs also varied greatly between the two Paxillus species studied. In addition, the nad4L gene was found to have been subjected to strong positive selection in the two species, which we suggest may be related to their contrasting ecological strategies. Potential gene fragment transferance between mitochondrial and nuclear genomes was also detected in the two Paxillus species. Large-scale gene rearrangements were also found to have occurred in basidiomycete species for which data were avaliable, with introns experiencing frequent loss or gain events, resulting in large variations of mitogenome organization and size in Basidiomycota. In addition, the insertion sites of introns in basidiomycete mitogenomes were found to have a certain base preference. Phylogenetic analysis revealed that mitochondrial genes were reliable tools for analyzing phylogenetic relationships within Basidiomycota. As the first report on the mitogenome from the genus Paxillus or even the suborder Paxilineae of Boletales, this study contributes to our understanding of the evolution, genetics and phylogeny of both Paxillus species and the phylum Basiiomycota as a whole.
Additional file 1 Fig. S1. Putative secondary structures of the 25 tRNA genes identified in the mitogenomes of two Paxillus species. Residues conserved across the two mitogenomes are shown in green, while variable sites are shown in red. All genes are shown in order of occurrence in the mitogenome of P. involutus, starting from trnM.
Additional file 2 Table S1. Comparison on mitogenomes between four Boletales species. Table S2. Gene features and organization of the Paxillus mitogenomes. Table S3. Tandem repeats detected in the mitogenomes of two Paxillus species using the Tandem Repeats Finder. Table S4. Distribution of repeat loci in the two Paxillus mitogenomes searched by REPuter. Table S5. Local BLAST searches between the mitochondrial and the nuclear genomes of two Paxillus species. Table  S6. Species and GenBank accession number used for phylogenetic analysis in this study.