The complete mitochondrial genomes of two model ectomycorrhizal fungi (Laccaria): features, intron dynamics and phylogenetic implications

Journal Pre-proofs The complete mitochondrial genomes of two model ectomycorrhizal fungi (Laccaria): features, intron dynamics and phylogenetic implications Qiang Li, Luxi Yang, Dabing Xiang, Yan Wan, Qi Wu, Wenli Huang, Gang Zhao PII: DOI: Reference:

S0141-8130(19)33128-9 https://doi.org/10.1016/j.ijbiomac.2019.09.188 BIOMAC 13469

To appear in:

International Journal of Biological Macromolecules

Received Date: Revised Date: Accepted Date:

28 April 2019 10 July 2019 20 September 2019

Please cite this article as: Q. Li, L. Yang, D. Xiang, Y. Wan, Q. Wu, W. Huang, G. Zhao, The complete mitochondrial genomes of two model ectomycorrhizal fungi (Laccaria): features, intron dynamics and phylogenetic implications, International Journal of Biological Macromolecules (2019), doi: https://doi.org/ 10.1016/j.ijbiomac.2019.09.188

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier B.V.

The complete mitochondrial genomes of two model ectomycorrhizal fungi (Laccaria): features, intron dynamics and phylogenetic implications Running title: Complete mitogenome of two Laccaria species

Qiang Lia,b,c, Luxi Yanga,b, Dabing Xianga,b, Yan Wana,b, Qi Wua,b, Wenli Huang*c, Gang Zhao*a,b

a College of Pharmacy and Biological Engineering, Chengdu University, Chengdu, Sichuan, China;

b Key Laboratory of Coarse Cereal Processing, Ministry of Agriculture and Rural Affairs, Chengdu, Sichuan, China;

c Biotechnology and Nuclear Technology Research Institute, Sichuan Academy of Agricultural Sciences, Chengdu, Sichuan, China

*Corresponding author:

Wenli Huang and Gang Zhao

E-mail: [email protected]; [email protected];

Phone: +86-028-84616653;

1

*Present address: Key Laboratory of Coarse Cereal Processing, Ministry of Agriculture and Rural Affairs, 2025 # Chengluo Avenue, Chengdu 610106, Sichuan, China.

Abstract Laccaria amethystine and L. bicolor have served as model species for studying the life history and genetics of ectomycorrhizal fungi. However, the characterizations and variations of their mitogenomes are still unknown. In the present study, the mitogenomes of the two Laccaria species were assembled, annotated, and compared. The two mitogenomes of L. amethystine and L. bicolor comprised circular DNA molecules, with the sizes of 65,156 bp and 95,304 bp, respectively. Genome collinearity analysis revealed large-scale gene rearrangements between the two Laccaria species. Comparative mitogenome analysis indicated the introns of cox1 genes in Agaricales experienced frequent lost/gain eveants, which promoted the organization and size variations in Agaricales mitogenomes. Evolutionary analysis indicated the core protein-coding genes in the two mitogenomes were subject to strong pressure of purifying selection. Phylogenetic analysis using the Bayesian inference (BI) and Maximum likelihood (ML) methods based on a combined mitochondrial gene set resulted in identical and well-supported tree topologies, wherein the two Laccaria species were most closely related to Coprinopsis cinerea. This study severed as the first study on the mitogenomes of Laccaria species, which 2

promoted a comprehensive understanding of the genetics and evolution of the model ectomycorrhizal fungi.

Key words: Laccaria; mitochondrial genome; intron; gene rearrangement; phylogenetic analysis

1. Introduction The genus Laccaria, belonging to Hydnangiaceae, Basidiomycota, is widely distributed in the world and contains around 75 recognized species [1]. Laccaria species is an important group of ectomycorrhizal fungi (ECMF), which plays an important role in maintaining forest ecosystem and promoting natural carbon and nitrogen cycle [2]. Laccaria spp. have been found from both temperate and tropical regions of the world, and form mutualistic symbiotic relationships with a variety of vascular plant families, such as Betulaceae, Salicaceae, Fagaceae, Pinaceae, and Myrtaceae [1, 3]. The adaptation of Laccaria genus to this mutualistic lifestyle was thought to be accomplished by a series of genetic evolution and adaptation, including carbohydrate-active enzymes [4], transcription factors [5, 6], aquaporins [7, 8] and so on [9, 10]. Two species of Laccaria, L. amethystine and L. bicolor, have served as model species for studying the life history and genetics of ECMF [11]. Their genomes have been published for revealing ectomycorrhizal lifestyles [12, 13]. However, up to now, the mitochondrial genomic characterizations and variations of Laccaria species

3

are still unknown [14], which limits our comprehensive understanding of the genetic characteristics of the model ECMF.

Traditional classification of Laccaria species was mainly based on morphological characteristics, including the color, size, shape of the basidiocarp, size and shape of basidiospore, and the number of sterigmata per basidium. However, less differentiable morphological features and the overlap of similar character states made it difficult to accurately distinguish species from the Laccaria genus [1]. The introduction of molecular markers has greatly promoted the development of taxonomy and biogeography of Laccaria genus, leading to the revision of some species and the identification of new species [15]. Phylogenetic trees based on four-gene nucleotide sequence datasets revealed an increase in Laccaria’s known diversity and the origin of the most recent common ancestor to Laccaria [11]. Cho et al. [1] systematically revised the Laccaria genus in Korea using polygenic phylogenetic tree. However, no phylogenetic study of Laccaria species based on mitochondrial gene set has been conducted.

As the “second genome” of eukaryotes, mitochondrial genome plays an important role in regulating the growth and development, aging, and death of eukaryotes, as well as coping with oxidative stress [16, 17]. In addition, the independent origin, rapid mutation rate, and several available molecular markers make mitochondrial genome a powerful tool for taxonomy, population genetics and evolution studies [18, 19]. However, so far, less than 100 mitochondrial genomes of Basidiomycota have been available in the public database 4

(https://www.ncbi.nlm.nih.gov/genome/browse#!/organelles/), which was less than one tenth of the available genomes. Deficiencies in the availability of fungal mitogenomes greatly limited our understanding of the characteristics, variations and evolution of mitogenomes from the largest mushroom-forming group (Basidiomycota). Previous studies have shown that the mitogenome of fungi varied greatly in gene content, gene order, genome structure, intron, repeat sequence and so on [20, 21], even among fungi with close phylogenetic relationships [22-24]. As mobile genetic elements in fungal mitogenomes, introns could affect the structure, size and GC content of fungal mitogenomes [25, 26]. Two types of introns were detected in fungi, namely the Group I and the Group II, which harboured homing nuclease genes [27]. The homing endonuclease can identify specific sites of the mitogenome and insert into it [28]. Mitgenome characteristics, variation and intron dynamics were helpful to reveal the genetic and evolutionary processes of fungal species [29]. It was reported that the mitogenomes of Basidiomycota contained 14 core protein coding genes (atp6, atp8, atp9, cob, cox1, cox2, cox3, nad1, nad2, nad3, nad4, nad4L, nad5, and nad6) for energy metabolism and one rps3 gene for transcriptional regulation [30, 31]. So far, no mitogenome has been reported in the Laccaria genus, even in the Hydnangiaceae family.

In the present study, the complete mitogenomes of two model ectomycorrhizal Basidiomycetes, L. bicolor and L. amethystine, were assembled and annotated. The two mitogenomes were compared to identify variations and similarities in genome organization, genome size, gene content, and gene order. The dynamic changes of 5

introns in the cox1 gene of Agaricales were also revealed. In addition, phylogenetic relationships among various Basidiomycota species were analyzed based on the combined mitochondrial gene sets. The mitogenomes of the two Laccaria species further our understanding of the taxonomy, evolutionary biology, and genetics of this important ectomycorrhizal genus.

2. Materials and Methods

2.1 Mitochondrial genome assembly and annotations Fruiting bodies of L. amethystina were obtained from our 90 mushroom sequencing project [32]. After sequencing, about 15 G of raw data were obtained for genome and mitogenome assembly. The raw sequencing data of L. bicolor was downloaded from the Joint Genome Institute (JGI) database [12]. Clean reads were obtained from the raw sequencing data through a series of quality control steps. First, adapter reads were removed using the AdapterRemoval v 2 [33], and then low-quality sequences were filtered using our own compling pipeline. The two Laccaria mitogenomes were assembled with the obtained clean reads using the SPAdes 3.9.0 software [34]. MITObim V1.9 [35] was used to fill gaps between contigs. Complete mitogenomes of Laccaria were first annotated according to our previously described methods [30, 31]. Briefly, the protein-coding genes (PCGs), rRNA genes and tRNA genes of the two mitogenomes were initially annotated using the MFannot [36] and MITOS [37], based on the genetic code 4. PCGs were then modified or predicted with the NCBI Open Reading Frame Finder [38], and further annotated by BLASTP 6

searches against the NCBI non-redundant protein sequence database [39]. tRNA genes were also predicted with tRNAscan-SE v1.3.1 [40]. Graphical maps of the two complete mitogenomes were drawn with OGDraw v1.2 [41].

2.2 Sequence analysis of the Laccaria mitogenomes

DNASTAR Lasergene v7.1 (http://www.dnastar.com/) was used to analyze the base compositions of the two Laccaria mitogenomes. Strand asymmetries of the two mitogenomes were assessed using the following formulas: AT skew = [A - T] / [A + T], and GC skew = [G - C] / [G + C] [42]. The Sequence Manipulation Suite [43] was used to analyze codon usages within the two mitogenomes, based on the genetic code 4. DnaSP v6.10.01 [44] was used to calculate the synonymous (Ks) and nonsynonymous (Ka) substitution rates for core PCGs in the two mitogenomes. MEGA v6.06 [45] was used to calculate the overall mean genetic distances between each pair of the 14 core PCGs (atp6, atp8, atp9, cox1, cox2, cox3, nad1, nad2, nad3, nad4, nad4L, nad5, nad6, and cob), and rps3 gene, using the Kimura-2-parameter (K2P) substitution model.

2.3 Repetitive element analysis To identify intra-genomic duplications of large fragments or interspersed repeats throughout the two mitogenomes, we conducted BLASTN searches of each mitogenome against itself [46] using an E-value of <10−10. The Tandem Repeats Finder [47] was used to detect tandem repeats (>10 bp in length) in the two mitogenomes with default parameters. Repeated sequences were also searched by 7

REPuter [48] to identify forward (direct), reverse, complemented, and palindromic (revere complemented) repeats.

2.4 Comparative mitogenomic analysis and intron analysis The genome sizes, GC content, base composition, start and stop codon, gene numbers, and intron numbers were compared among different Agaricales species to assess variations and conservation among Agaricales mitogenomes. Genome colinearity of the two Laccaria species and 5 closely related species was analyzed with Mauve v2.4.0 [49]. Group I introns of cox1 genes in 17 Agaricales mitogenomes reported were classified into different position classes (Pcls) according to the method described by Férandon et al. [50]. 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 in letter according to the similarity with the described Pcls [50].

2.5 Phylogenetic analysis In order to investigate the phylogenetic status of the two Laccaria species among Basidiomycota phylum, we constructed a phylogenetic tree of 59 species based on the combined mitochondrial gene set (14 core PCGs + rps3 +2 rRNA genes) [30]. Annulohypoxylon stygium from the Ascomycota phylum was used as the outgroup [51]. We first used the MAFFT v7.037 software [52] to align individual mitochondrial gene, and then concatenated these alignments in SequenceMatrix v1.7.8 [53] to obtain 8

a combined mitochondrial gene set. DAMBE 7 [54] was used to test the saturation of mitochondrial DNA sequences, and it was determined that the mitochondrial DNA sequence were suitable for phylogenetic analysis. A preliminary partition homogeneity test was carried out to detect potential phylogenetic conflicts between different genes. Best-fit models of evolution and partitioning schemes for the gene set were determined according to PartitionFinder 2.1.1 [55]. Phylogenetic trees were constructed using both bayesian inference (BI) and maximum likelihood (ML) methods.㻌MrBayes v3.2.6 [56] was used for the BI analysis, and RAxML v 8.0.0 [57] was used to perform the ML analysis .

2.6 Data availability

The complete mitogenomes of L. amethystine and L. bicolor were deposited in the GenBank database under the accession numbers MK697669 and MK697670, respectively.

3. Results

3.1 Features and protein-coding genes of Laccaria mitogenomes The mitogenomes of L. amethystine and L. bicolor were all composed of circular DNA molecules, with the sizes of 65,156 bp and 95,304 bp, respectively (Fig. 1). The GC content of L. bicolor mitogenome was 4.93% higher than that of L. amethystine, reaching 28.32% (Table S1). Both the GC skew and AT skew of the L. amethystine mitogenome were positive, while both were negative in the L. bicolor mitogenome.

9

A total of 23 protein-coding genes (PCGs) were found in the L. amethystine mitogenome, including 6 open reading frames (ORFs) in introns and 17 non-intronic ORFs (Table S2). Twenty-nine PCGs were detected in the L. bicolor mitogenome, of which 9 were located in introns. Non-intronic PCGs in the two mitogenomes were mainly core genes for energy metabolism, including atp6, atp8, atp9, cob, cox1, cox2, cox3, nad1, nad2, nad3, nad4, nad4L, nad5, and nad6. In addition, both Laccaria species contained an rps3 gene for transcriptional regulation. The L. amethystine and L. bicolor mitogenomes contained 2 and 3 PCGs with unknown functions, respectively. The L. bicolor mitogenome also contained 1 DNA polymerase 2 gene and 1 DNA-directed DNA polymerase gene. Intronic ORFs in the two Laccaria mitogenomes encoded homing endonucleases from different families, such as GIYYIG endonuclease and LAGLIDADG endonuclease.

3.2 rRNA genes and tRNA genes The mitogenomes of L. amethystine and L. bicolor both contained two rRNA genes, the large subunit ribosomal RNA (rnl) gene and the small subunit ribosomal RNA gene (rns) (Table S2). The length and organization of rRNA genes in the two mitogenomes varied. The length of rns gene in the L. amethystine mitogenome was 1,936 bp, while it was 2,137 bp in the L. bicolor mitogenome. One intron was detected in the rnl gene of L. amethystine mitogenome, and two introns were found in L. bicolor rnl gene. The length of rnl gene in the L. amethystine mitogenome was 3,415 bp, while it was 3,540 bp in the L. bicolor mitogenome.

10

Both the L. amethystine and L. bicolor mitogenomes contained 25 tRNA genes, which encoded 20 standard amino acids (Table S2). The two mitogenomes both contained two tRNAs that code for arginine, serine, and leucine with different anticodons and three tRNAs coding for methionine with the same anticodons. All tRNAs were folded into classical cloverleaf structures (Fig. S1), with the length of individual tRNAs ranging from 71 bp to 87 bp. The length variations of tRNA genes in the two mitogenomes were mainly due to the expansion or contraction of extra arms. The length expansions of extra arms in trnl, trns, and trnY genes contributed to their lengths over 80 bp. Of the 25 tRNAs detected in the two mitogenomes, 15 contained sites that varied between the two mitogenomes. The most common variable site was located on the D arm (seven sites varied between the two mitogenomes), indicated that RNA genes in the two Laccaria mitogenomes varied greatly.

3.3 Codon usage analysis ATG was the common start codon in most of the core PCGs from the 17 Agaricales species we tested, except in the cox1 gene, which used GTG as the start codon in most of the Agaricales species (Table S3). TAA was the most commonly used stop codon of core PCGs in Agaricales species, followed by TAG. The nad1 gene of L. bicolor used TTG as the start codon, while ATG was used as start codons in L. amethystine and other Agaricales mitogenomes. GTG was used as start codons for cox1 gene in the two Laccaria species. The stop codons of core PCGs in the two Laccaria species were all TAA.

11

Codon usage analysis indicated that the most frequently used codons in the two mitogenomes were AAA (for lysine; Lys), TTT (for phenylalanine; Phe), TTA (for leucine; Leu), AAT (for asparagine; Asn), and ATT (for isoleucine; Ile) (Fig. S2). The high frequency of A and T used in codons contributed to the high AT content of the two Laccaria mitogenomes (average: 74.15%).

3.4 Intergenic sequence and mitoenome composition The mitogenomes of the two Laccaria species both contained two overlapping regions, one of which was located between nad5 and nad4L (1 bp) genes, and the other was located across the neighboring nad3 and nad2 genes (1 bp) (Table S2). A total of 32,083 bp and 53,482 bp intergenic sequences were detected in the mitogenomes of L. amethystine and L. bicolor, respectively. The length of individual intergenic sequence in the two mitogenomes ranged from 73 bp to 2,835 bp, with the longest intergenic sequence located between cox1 and trnM genes in the L. bicolor mitogenome.

Intergenic regions accounted for the largest proportion of the two Laccaria mitogenomes, comprising 49.24% - 56.12% of the total size, indicating that the mitogenomes of the two Laccaria species had relatively relaxed structures (Fig. 2). Protein coding regions were the second largest part in the two mitogenomes, accounting for 26.08% - 21.92% of the two Laccaria mitogenomes. Intronic regions accounted for 13.60% - 14.04% of the two Laccaria mitogenomes. RNA coding regions were the smallest part of the two mitogenomes, accounting for only 7.91% 12

11.98% of the two mitogenomes. L. bicolor was 30,148 bp larger than the mitogenome of L. amethystine, of which the intergenic region contributed the most to the expansion of the mitogenome in L. bicolor, accounting for 70.98% of the mitogenome expansion, followed by the intronic region, which contributed 15.00% of the mitogenome expansion in L. bicolor. The protein coding region contributed 12.94% to the expansion of L. bicolor mitogenome. The contribution rate of RNA region to the expansion of the L. bicolor mitogenome was 1.08%.

3.5 Repeat sequences of the two Laccaria mitogenomes By comparing the mitogenomes of two Laccaria species with themselves through BLASTN search, we identified 11 repetitive sequences in the L. amethystine mitogenome and 740 in the L. bicolor mitogenome (Table S4). The length of these repeat sequences ranged from 28 bp to 453 bp, with pair-wise nucleotide similarities ranging from 80.80% to 100%. The longest repeat sequence was located in the intergenic region between orf281 and trns, and also in the intergenic region between trnE and orf255 in the L. bicolor mitogenome. The longest repeat sequence in the L. amethystine mitogenome was located between orf224 and trnC, and also between cox3 and trnE, with a length of 335 bp. Repeat sequences accounted for 2.22% of the L. amethystine mitogenome. While in the mitogenome of L. bicolor, repetitive sequences increase to 41.10% of the total mitogenome.

A total of 53 and 160 tandem repeats were detected in the L. amethystine and L. bicolor mitogenomes, respectively (Table S5). The longest tandem sequence was 13

located between orf255 and cox3 in the L. bicolor mitogenome, with the length of 164 bp. Heptanucleotide repeats were the most frequent repeat loci in the two Laccaria mitogenomes. Most tandem repeat sequences were repeated 2 - 4 times in the two Lacccaria mitogenomes, with the highest replication number (37) observed in the L. bicolor mitogenome. Tandem repeat sequences accounted for 3.79% and 9.61% of the L. amethystine and L. bicolor mitogenomes, respectively. Using the REPuter, we identified 6 complemented, 18 forward, 12 palindromic, and 14 reverse repeats in the mitogenome of L. amethystine, accounting for 3.02% of the entire mitogenome (Table S6). Repeats identified by REPuter accounted for 4.38% of the L. bicolor mitogenome.

3.6 Gene rearrangement in Laccaria mitogenomes Of the 15 core PCGs and 2 rRNA genes detected, most genes varied in relative positions between the two Laccaria mitogenomes and 5 closely related mitogenomes in Agaricales, indicating that the gene order of Agaricales were highly variable (Fig. 3). Agaricus bisporus and Coprinopsis cinerea have close phylogenetic relationships with the two Laccaria species. However, frequent gene rearrangements lead to significant differences in the mitochondrial gene order between them. Eleven of the 17 genes tested showed variations in location between the two Laccaria species, which indicated that gene rearrangements frequently occurred during the evolution of Laccaria species.

14

A total of 20 homologous regions were detected between mitogenomes of the 7 closely related species (Fig. 4). Genome collinearity analysis showed that the relative positions of these homologous regions varied greatly between different species. A total of 18 homologous regions were detected in mitogenomes of the two Laccaria species, of which 13 varied in the relative positions between the two mitogenomes. The results further indicated that large-scale gene rearrangements occurred in Laccaria species, as well as in the Agaricales order during the evolutionary process.

3.7 Variation, genetic distance, and evolutionary rates of core genes Among the 15 core PCGs detected, 4 genes varied in length between the two Laccaria mitogenomes, including cox1, nad1, nad2 and rps3 genes (Fig. 5). Among them, the rps3 gene had the greatest variation between the two mitogenomes, with 54 amino acid variations. GC contents of most core PCGs varied between the two mitogenomes, except that of the atp9 gene. The GC content of atp9 gene was the highest among the core PCGs detected, with an average of 37.84%, while that of atp8 gene was the lowest, with an average of 16.04%. AT skews of most core PCGs were negative, except the rps3 gene. GC skews of the core PCGS in the two mitogenomes were highly variable, of which atp6, atp8, nad2, nad3, nad4, nad6, and rps3 genes were negative and others were positive. Among the 25 tRNAs detected in the two Laccaria mitogenomes, 1 varied in length between them, namely trnM. The GC content of tRNA genes in the mitogenomes of the two Laccaria species also varied greatly. The highest GC content was found in trnP gene, averaging 49.33%, while the lowest GC content found in trnI, with 29.58%. 15

Among the 15 core PCGs detected, rps3 had the highest K2P genetic distance between the two mitogenomes, followed by nad2 and atp9 (Fig. 6). The genetic distance of nad4L was the lowest among the 15 core PCGs between the two mitogenomes, indicating that this gene was highly conservative between Laccaria species. The rps3 gene had the highest synonymous substitution (Ks) rate, while cox3 had the lowest substitution rate among the 15 core PCGs detected. The highest nonsynonymous substitution rate (Ka) was in rps3, while atp6, atp9 and nad4L exhibited the lowest Ka values among the 15 PCGs. The Ka/Ks values for all 15 PCGs were less than 1, indicating that these genes have been subjected to purifying selection.

3.8 Intron dynamics of cox1 gene in Agaricales

A total of 172 introns were detected in the mitogenomes of 17 Agaricales species tested, which were distributed in cox1, cox2, cob, nad1, nad5, and rnl genes (Fig. 7). The number of introns in each species ranged from 0 to 46. The dynamic changes of the number and location of introns promoted the organization and size variations in Agaricales mitogenomes. A total of 94 introns were detected distributing in the cox1 gene of Agaricales, accounting for 54.65% of the total introns, of which 2 introns belonged to the Group II. All Group I introns of the cox1 gene could be classified into different position classes (Pcls) according to their insertion positions in the coding region of cox1 gene. Introns belonging to the same Pcl had high sequence similarities and often contained orthomologous intronic ORFs. Different Pcls showed low sequence similarities with nonhomologous intronic ORFs. Ninety-two Group I introns of cox1 genes in the 17 Agaricales species were assigned into 24 Pcls and named 16

alphabetically according to methods described by Férandon et al [50]. Among the 17 Agaricales species, A. bisporus contained the most Group I introns (45), followed by Pleurotus platypus (9) and Pleurotus ostreatus (9). While the cox1 genes of Coprinopsis cinerea, Lyophyllum decastes, and Schizophyllum commune did not contain any intron. Giant variations of intron number in the cox1 genes of Agaricales species indicated that the introns of cox1 genes in Agaricales had been subjected large scale lost/gain events during evolution. Of the 24 Pcls detected, Pcl K was the most common intron classes in the cox1 genes of Agaricales, which could be detected in 12 of the 17 Agaricales species. Pcl AI was also a widely distributed intron class in Agaricales, found in 10 of the 17 Agaricales species. Pcl B, H, T, U, V, Z, AG could only be found in one of the 17 Agaricales species. However, these Pcls could be found in distant species from order Rhizophydiales, Polyporales, Hypocreales, Harpellales, even from Charales [50], suggesting that potential horizontal gene transfer might occur between different species. Three and four introns were detected in cox1 genes of L. amethystine and L. bicolor, respectively, which belonged to 7 intron classes. All the Pcls of cox1 genes in L. amethystine and L. bicolor were different, while same Pcls could be found in distant species from other genus, suggesting that intron loss/acquisition events occurred during evolution of Laccaria.

3.9 Comparative mitogenome analysis and phylogenetic analysis

The sizes of the two Laccaria mitogenomes were in the middle among the 17 Agaricales mitogenomes detected (Table S1), which were larger than C. cinerea (42,448 bp) [58], L. decastes (50,643 bp) [59], S. commune (49,704 bp) [60], P. 17

citrinopileatus (60,694 bp) [22], and H. russula (55,769 bp) [61], and smaller than A. bisporus (135,005 bp) [50], L. edodes (121,394 bp) [62], and M. perniciosa (109,103 bp) [63]. The GC content of the two Laccaria mitogenomes were lower than A. bisporus (29.07%) [50], L. edodes (30.70%) [62], and M. perniciosa (31.89%) [63], and higher than L. decastes (22.27%) [59], Tricholoma matsutake (20.69%) [64], and Flammulina velutipes (16.50%) [65]. Of the 17 Agaricales tested, only L. bicolor, C. cinerea, and L. decastes had negative GC skews, while the others were positive. AT skews of Agaricales mitogenomes varied. Eight and ten introns were detected in the L. amethystine and L. bicolor mitogenomes, respectively, which were smaller than the 46 introns of the A. bisporus mitogenome that made A. bisporus the largest mitogenome in Agaricales. All the 17 Agaricales mitogenomes contained 2 rRNA genes. The number of tRNA genes in the 17 Agaricales mitogenomes ranged from 24 to 35.

Phylogenetic analysis using Bayesian inference (BI) and Maximum likelihood (ML) 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 (BPP ≥0.92; BS ≥ 82). Based on the phylogenetic analysis, the 58 Basidiomycota species could be divided into 13 major clades, corresponding to the orders Agaricales, Auriculariales, Russulales, Polyporales, Cantharellales, Tremellales, Trichosporonales, Microstromatales, Ustilaginales, Tilletiales, Microbotryales, Sporidiobolales, and Pucciniales. The 17 Agaricales species could be divided into five groups, wherein the first comprised only 18

one species form the Hygrophorus genus (Hygrophorus russula), the second group comprised four species within the Pleurotus genus, the third group was recovered as (Schizophyllum commune + (Flammulina velutipes + (Lentinula edodes + (Moniliophthora perniciosa + Moniliophthora roreri)), the forth group was recovered as (Tricholoma matsutake + (L. decastes + L. shimeji)), and the fifth group was recovered as (Agaricus bisporus + (Coprinopsis cinerea+ (L. amethystine + L. bicolor)). The phylogenetic analyses indicated that L. amethystine was a sister species to L. bicolor. The analyses also indicated that the Laccaria genus showed close relationships with A. bisporus [50], C. cinerea [58], T. matsutake [64], and Lyophyllum genus [59].

4. Discussion 4.1 Size variations of mitogenomes in Agaricales The fungal mitogenomes are among the most variable mitogenome of eukaryotic groups [66]. The smallest known mitogenome of fungi was the Spizellomyces punctatus mitogenome (1.14 kb) [67], while the largest mitogenome in fungi was from Rhizoctonia solani (235.84 kb) [68]. The size variations of fungal mitogenomes were mainly due to the accumulation of repetitive sequences, the change of introns, the difference of intergenic sequences and the horizontal transfer of genes [21, 68]. In the present study, great variations in the sizes of the Agaricales mitogenomes were observed. A. bisporus was found containing the largest mitogenome in Agaricales[50], while the mitogenomes of C. cinerea [58] and S. commune [60] were the smallest. The variation of intron was considered to be one of the main factors 19

resulting in the variation of mitogenome size in Agaricales. The mitogenome of A. bisporus contained the most introns (46), while only 1 or 0 introns were detected in C. cinerea and S. commune mitogenomes, respectively. We also found that the mitogenome sizes differed greatly between the two Laccaria species. The expansion of mitogenome in L. bicolor was mainly due to the increase of intergenic region, followed by the intronic region. In conclusion, the both the intergenic region and intronic region contributed to the expansion and contraction of mitogenomes in Laccaria or Agaricales [59].

4.2 Variations of mitogenome content in Laccaria species

With the rapid development of next-generation sequencing technology, more and more fungal mitogenomes have been obtained, which promotes our understanding of the content variations in fungal mitogenomes [69, 70]. During the evolution of fungi, most mitochondrial genes have been transferred to nuclear genes, which was considered had many advantages [71]. However, most fungal species retained 14 core PCGs for energy metabolism and one rps3 gene for transcriptional regulation [72-74]. These genes could be found in all published Basidiomycota mitogenomes and were regarded as the core genes of Basidiomycota species. In addition, Basidiomycota mitogenomes also contained 20-35 tRNA genes for amino acid transport. In the present study, we found that the core PCGs of Laccaria species differed in base composition and gene length, suggesting that the core PCGs of Laccaria species undergo genetic differentiation during evolution. In addition, the genetic distances of the core PCGs between the two Laccaria species were also varied, indicating that 20

different core PCGs of Laccaria species exhibited different evolutionary rates. The Ka/Ks values for all core PCGs in the two Laccaria species were less than 1, indicating that they were subjected to strong pressure of purifying selection. In addition, the tRNA genes had also undergone tremendous base variations between the two Laccaria species. Base variations of tRNA genes were reported could affect the efficiency of amino acid transport in some eukaryotic organisms [75, 76]. The effect of tRNA variations in Laccaria species on amino acid transport and protein synthesis needed further study. In addition to core PCGs, several non-intronic ORFs were found in the two Laccaria mitogenomes, which encoded DNA polymerase and some proteins with unknown function. It was suggested that there were still many unknown proteins in fungal mitogenome that need to be further elucidated, which could promote further understanding of the function and evolution of fungal mitogenomes.

4.3 Mitogenome rearrangement in Laccaria and its closely related species The gene arrangement of mitogenome could be used as a reference for analyzing the phylogenetic status and evolution of species [77]. With more and more mitogenomes being obtained, gene rearrangements have been frequently found in mitogenomes of animals, plants and fungi [78]. The arrangement of animal mitogenome has been extensively studied, and a series of models have been proposed to explain gene rearrangements in animal mitogenome [79, 80]. However, the mitochondrial gene rearrangement in fungi has been less studied compared with that of animals, which was mainly due to the insufficient number of available fungal mitogenomes. It has been reported that the mitochondrial gene order of fungi varied 21

greatly among different species, even among species from the same genus [22, 78]. In the present study, we found that large-scale gene rearrangements occurred between the two Laccaria species and their closely related species, involving core PCGs, nonintronic ORFs, tRNAs and rRNA genes. Eleven of the 17 detected genes changed their locations between the two Laccaria mitogenomes. Aguileta et al [78] reported that the accumulation of repetitive sequences in fungal mitogenomes could contribute to recombination of mitogenomes, leading to rearrangements of mitochondrial genes in fungi. In the present study, large-scale repetitive sequences were observed accumulated in the L. bicolor mitogenome, accounting for 41.10% of the whole mitogenome, which may be the main factor leading to gene rearrangement in the Laccaria mitogenome. Further studies are needed to reveal the mechanism of organization variation in mitogenomes of Laccaria species and other macrofungi.

4.4 Dynamic changes of cox1 gene introns in Agaricales Introns were considered as mobile genetic elements in mitogenome [25, 81]. The dynamic changes of introns were one of the main factors affecting the size and organization of fungal mitochondrial genomes [82, 83]. Up to now, introns have been detected in several genes of fungi, including cox1, cox 2, cob, nad1, nad5, and rnl genes, among which cox1 gene was the main host gene of fungal introns [84]. In the present study, 54.65% of the introns were detected in the cox1 gene of Agaricales, so the dynamic changes of introns in cox1 gene would exert significant influence on the size and organization of fungal mitogenomes. Ninety-two introns in cox1 genes of 17 Agaricales species were classified into 24 Pcls in this study. There were significant 22

differences in the types and quantities of Pcl between different species of Agaricales, suggesting that frequent intron loss/gain events occurred during the evolution of Agaricales. Homologous introns could not be found between the two Laccaria species, but found from distant species, suggesting that possible horizontal gene transfer events occurred. In conclusion, introns were highly variable in Laccaria genus, as well as in Agaricales, which promoted the organization and size variations of Agaricales mitogenomes and also increased the diversities of Agaricales mitogenomes.

4.5 Phylogenetic analysis of Basidiomycota based on mitochondrial genes

Limited identifiable morphological features and the overlap of similar character states made it difficult to distinguish species and species complexes in Laccaria genus accurately [1]. The introduction of molecular markers has promoted the development of taxonomy, population genetics and biogeography of the Laccaria genus. Up to now, rDNA ITS, RPB2, and EF1α genes have been used to analyze the phylogenetic relationships among Laccaria species [15, 85]. Mitogenome has been widely used in phylogenetic studies of many eukaryotic groups, such as animals and plants, due to its advantages [86]. The mitogenome has become a powerful tool for studying population genetics, taxonomy and genetics [87, 88]. However, as the model group of ectomycorrhizal fungi, no complete mitogenome has been reported in the Laccaria genus or even in the Hydnangiaceae family, which limits the comprehensive understanding of the heredity and evolution of the model ectomycorrhizal fungi. In this study, we constructed a phylogenetic tree of 58 Basidiomycota species based on a 23

combined mitochondrial gene set, revealing the phylogenetic status of the Laccaria genus in Basidiomycota. Using the BI and ML methods, we obtained identical and well-supported tree topologies, indicating that mitochondrial genes can be used as reliable molecular markers for studying phylogenetic relationships of Basidiomycota.

5. Conclusion In the present study, the mitogenomes of two Laccaria species were assembled and compared. The sizes of mitogenomes in the two Laccaria species and other Agaricales species varied greatly. Intergenic region and intronic region were considered to be the main factors leading to mitogenome size variations in Agaricales. Comparative mitogenome analysis indicated the gene contents, base compositions, gene lengths, tRNAs, and rRNAs varied greatly between the two mitogenomes. In addition, all core PCGs of the two Laccaria species were subjected to strong pressure of purifing selection. Gene collinearity analysis revealed that large-scale gene rearrangements occurred between Laccaria species and five closely related species, involving PCGs, tRNAs and rRNAs. The introns of cox1 gene in Agaricales experienced lost/gain events, which promoted the organization and size variations of Agaricales mitogenomes. Phylogenetic analysis showed that mitochondrial genes could be used as a reliable tool to analyze phylogenetic relationships of Basidiomycota. This study is the first to report complete mitogenomes of model ectomycorrhizal fungi (Laccaria spp.), which will further promote the understanding of the genetics, evolution and taxonomy of the ectomycorrhizal fungi.

Conflict of Interest: The authors declare that they have no conflict of interest. 24

References

[1] H.J. Cho, M.S. Park, H. Lee, S.Y. Oh, A.W. Wilson, G.M. Mueller, Y.W. Lim, A systematic revision of the ectomycorrhizal genus Laccaria from Korea, Mycologia 110(5) (2018) 948-961. [2] M.A. Wadud, K. Nara, C. Lian, T.A. Ishida, T. Hogetsu, Genet dynamics and ecological functions of the pioneer ectomycorrhizal fungi Laccaria amethystina and Laccaria laccata in a volcanic desert on Mount Fuji, Mycorrhiza 24(7) (2014) 551-63. [3] A. Ramos, V.M. Bandala, L. Montoya, A new species and a new record of Laccaria (Fungi, Basidiomycota) found in a relict forest of the endangered Fagus grandifolia var. mexicana, MycoKeys (27) (2017) 77-94. [4] C. Veneault-Fourrey, C. Commun, A. Kohler, E. Morin, R. Balestrini, J. Plett, E. Danchin, P. Coutinho, A. Wiebenga, R.P. de Vries, B. Henrissat, F. Martin, Genomic and transcriptomic analysis of Laccaria bicolor CAZome reveals insights into polysaccharides remodelling during symbiosis establishment, Fungal Genet Biol 72 (2014) 168-181. [5] Y. Daguerre, E. Levati, J. Ruytinx, E. Tisserant, E. Morin, A. Kohler, B. Montanini, S. Ottonello, A. Brun, C. Veneault-Fourrey, F. Martin, Regulatory networks underlying mycorrhizal development delineated by genome-wide expression profiling and functional analysis of the transcription factor repertoire of the plant symbiotic fungus Laccaria bicolor, BMC Genomics 18(1) (2017) 737. [6] J. Labbe, C. Murat, E. Morin, G.A. Tuskan, F. Le Tacon, F. Martin, Characterization of transposable elements in the ectomycorrhizal fungus Laccaria bicolor, PLoS One 7(8) (2012) e40197. [7] A. Navarro-RoDenas, H. Xu, M. Kemppainen, A.G. Pardo, J.J. Zwiazek, Laccaria bicolor aquaporin LbAQP1 is required for Hartig net development in trembling aspen (Populus tremuloides), Plant Cell Environ 38(11) (2015) 2475-86. [8] H. Xu, A. Navarro-Rodenas, J.E. Cooke, J.J. Zwiazek, Transcript profiling of aquaporins during basidiocarp development in Laccaria bicolor ectomycorrhizal with Picea glauca, Mycorrhiza 26(1) (2016) 19-31. [9] A. Flores-Monterroso, J. Canales, F. de la Torre, C. Avila, F.M. Canovas, Identification of genes differentially expressed in ectomycorrhizal roots during the Pinus pinaster-Laccaria bicolor interaction, Planta 237(6) (2013) 1637-50. [10] L. Shi, J. Wang, B. Liu, K. Nara, C. Lian, Z. Shen, Y. Xia, Y. Chen, Ectomycorrhizal fungi reduce the light compensation point and promote carbon fixation of Pinus thunbergii seedlings to adapt to shade environments, Mycorrhiza 27(8) (2017) 823-830. [11] A.W. Wilson, K. Hosaka, G.M. Mueller, Evolution of ectomycorrhizas as a driver of diversification and biogeographic patterns in the model mycorrhizal mushroom genus Laccaria, New Phytol 213(4) (2017) 1862-1873. [12] F. Martin, A. Aerts, D. Ahren, A. Brun, E.G. Danchin, F. Duchaussoy, J. Gibon, A. Kohler, E. Lindquist, V. Pereda, A. Salamov, H.J. Shapiro, J. Wuyts, D. Blaudez, M. Buee, P. Brokstein, B. Canback, D. Cohen, P.E. Courty, P.M. Coutinho, C. Delaruelle, J.C. Detter, A. Deveau, S. DiFazio, S. Duplessis, L. Fraissinet-Tachet, E. Lucic, P. Frey-Klett, C. Fourrey, I. Feussner, G. Gay, J. Grimwood, P.J. Hoegger, P. Jain, S. Kilaru, J. Labbe, Y.C. Lin, V. Legue, F. Le Tacon, R. Marmeisse, D. Melayah, B. Montanini, M. Muratet, U. Nehls, H. Niculita-Hirzel, M.P. Oudot-Le Secq, M. Peter, H. Quesneville, B. Rajashekar, M. Reich, N. Rouhier, J. Schmutz, T. Yin, M. Chalot, B. Henrissat, U. Kues, S. Lucas, Y. Van de Peer, G.K. Podila, A. Polle, P.J. Pukkila, P.M. Richardson, P. Rouze, I.R. Sanders, J.E. Stajich, A. Tunlid, G. Tuskan, I.V. 25

Grigoriev, The genome of Laccaria bicolor provides insights into mycorrhizal symbiosis, Nature 452(7183) (2008) 88-92. [13] A. Kohler, A. Kuo, L.G. Nagy, E. Morin, K.W. Barry, F. Buscot, B. Canback, C. Choi, N. Cichocki, A. Clum, J. Colpaert, A. Copeland, M.D. Costa, J. Dore, D. Floudas, G. Gay, M. Girlanda, B. Henrissat, S. Herrmann, J. Hess, N. Hogberg, T. Johansson, H.R. Khouja, K. LaButti, U. Lahrmann, A. Levasseur, E.A. Lindquist, A. Lipzen, R. Marmeisse, E. Martino, C. Murat, C.Y. Ngan, U. Nehls, J.M. Plett, A. Pringle, R.A. Ohm, S. Perotto, M. Peter, R. Riley, F. Rineau, J. Ruytinx, A. Salamov, F. Shah, H. Sun, M. Tarkka, A. Tritt, C. Veneault-Fourrey, A. Zuccaro, C. Mycorrhizal Genomics Initiative, A. Tunlid, I.V. Grigoriev, D.S. Hibbett, F. Martin, Convergent losses of decay mechanisms and rapid turnover of symbiosis genes in mycorrhizal mutualists, Nat Genet 47(4) (2015) 410-5. [14] E.M. Sheedy, A.P. Van de Wouw, B.J. Howlett, T.W. May, Mitochondrial microsatellite markers for the Australian ectomycorrhizal fungus Laccaria sp. A (Hydnangiaceae), Appl Plant Sci 2(3) (2014). [15] E.M. Sheedy, A.P. Van de Wouw, B.J. Howlett, T.W. May, Multigene sequence data reveal morphologically cryptic phylogenetic species within the genus Laccaria in southern Australia, Mycologia 105(3) (2013) 547-63. [16] G.S. Gorman, A.M. Schaefer, Y. Ng, N. Gomez, E.L. Blakely, C.L. Alston, C. Feeney, R. Horvath, P. YuWai-Man, P.F. Chinnery, R.W. Taylor, D.M. Turnbull, R. McFarland, Prevalence of nuclear and mitochondrial DNA mutations related to adult mitochondrial disease, Ann Neurol 77(5) (2015) 753-9. [17] A. Latorre-Pellicer, R. Moreno-Loshuertos, A.V. Lechuga-Vieco, F. Sanchez-Cabo, C. Torroja, R. AcinPerez, E. Calvo, E. Aix, A. Gonzalez-Guerra, A. Logan, M.L. Bernad-Miana, E. Romanos, R. Cruz, S. Cogliati, B. Sobrino, A. Carracedo, A. Perez-Martos, P. Fernandez-Silva, J. Ruiz-Cabello, M.P. Murphy, I. Flores, J. Vazquez, J.A. Enriquez, Mitochondrial and nuclear DNA matching shapes metabolism and healthy ageing, Nature 535(7613) (2016) 561-5. [18] F. Delsuc, M.J. Stanhope, E.J. Douzery, Molecular systematics of armadillos (Xenarthra, Dasypodidae): contribution of maximum likelihood and Bayesian analyses of mitochondrial and nuclear genes, Mol Phylogenet Evol 28(2) (2003) 261-75. [19] Y. Nie, L. Wang, Y. Cai, W. Tao, Y.J. Zhang, B. Huang, Mitochondrial genome of the entomophthoroid fungus Conidiobolus heterosporus provides insights into evolution of basal fungi, Appl Microbiol Biotechnol 103(3) (2019) 1379-1391. [20] Y. Deng, Q. Zhang, R. Ming, L. Lin, X. Lin, Y. Lin, X. Li, B. Xie, Z. Wen, Analysis of the Mitochondrial Genome in Hypomyces aurantius Reveals a Novel Twintron Complex in Fungi, Int J Mol Sci 17(7) (2016). [21] H. Salavirta, I. Oksanen, J. Kuuskeri, M. Makela, P. Laine, L. Paulin, T. Lundell, Mitochondrial genome of Phlebia radiata is the second largest (156 kbp) among fungi and features signs of genome flexibility and recent recombination events, PLoS One 9(5) (2014) e97141. [22] Q. Li, C. Chen, C. Xiong, X. Jin, Z. Chen, W. Huang, Comparative mitogenomics reveals large-scale gene rearrangements in the mitochondrial genome of two Pleurotus species, Appl Microbiol Biotechnol 102(14) (2018) 6143-6153. [23] Q. Li, Q. Wang, X. Jin, Z. Chen, C. Xiong, P. Li, Q. Liu, W. Huang, Characterization and comparative analysis of six complete mitochondrial genomes from ectomycorrhizal fungi of the Lactarius genus and phylogenetic analysis of the Agaricomycetes, Int J Biol Macromol 121 (2019) 249-260. [24] Q. Li, Q. Wang, C. Chen, X. Jin, Z. Chen, C. Xiong, P. Li, J. Zhao, W. Huang, Characterization and comparative mitogenomic analysis of six newly sequenced mitochondrial genomes from ectomycorrhizal fungi (Russula) and phylogenetic analysis of the Agaricomycetes, Int J Biol Macromol 26

119 (2018) 792-802. [25] Z. Hamari, A. Juhasz, F. Kevei, Role of mobile introns in mitochondrial genome diversity of fungi (a mini review), Acta Microbiol Immunol Hung 49(2-3) (2002) 331-5. [26] J.M. Burke, Molecular genetics of group I introns: RNA structures and protein factors required for splicing--a review, Gene 73(2) (1988) 273-94. [27] A. Zubaer, A. Wai, G. Hausner, The mitochondrial genome of Endoconidiophora resinifera is intron rich, Sci Rep 8(1) (2018) 17591. [28] P.Q. Anziano, D.K. Hanson, H.R. Mahler, P.S. Perlman, Functional domains in introns: trans-acting and cis-acting regions of intron 4 of the cob gene, Cell 30(3) (1982) 925-32. [29] C.W. Basse, Mitochondrial inheritance in fungi, Curr Opin Microbiol 13(6) (2010) 712-9. [30] Q. Li, M. Yang, C. Chen, C. Xiong, X. Jin, Z. Pu, W. Huang, Characterization and phylogenetic analysis of the complete mitochondrial genome of the medicinal fungus Laetiporus sulphureus, Sci Rep 8(1) (2018) 9104. [31] Q. Li, C. Chen, C. Xiong, X. Jin, Z. Chen, W. Huang, Comparative mitogenomics reveals large-scale gene rearrangements in the mitochondrial genome of two Pleurotus species, Appl Microbiol Biotechnol (2018). [32] H. Li, S. Wu, X. Ma, W. Chen, J. Zhang, S. Duan, Y. Gao, L. Kui, W. Huang, P. Wu, R. Shi, Y. Li, Y. Wang, J. Li, X. Guo, X. Luo, Q. Li, C. Xiong, H. Liu, M. Gui, J. Sheng, Y. Dong, The Genome Sequences of 90 Mushrooms, Sci Rep 8(1) (2018) 9982. [33] M. Schubert, S. Lindgreen, L. Orlando, AdapterRemoval v2: rapid adapter trimming, identification, and read merging, BMC Res Notes 9 (2016) 88. [34] A. Bankevich, S. Nurk, D. Antipov, A.A. Gurevich, M. Dvorkin, A.S. Kulikov, V.M. Lesin, S.I. Nikolenko, S. Pham, A.D. Prjibelski, A.V. Pyshkin, A.V. Sirotkin, N. Vyahhi, G. Tesler, M.A. Alekseyev, P.A. Pevzner, SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing, J Comput Biol 19(5) (2012) 455-77. [35] C. Hahn, L. Bachmann, B. Chevreux, Reconstructing mitochondrial genomes directly from genomic next-generation sequencing reads--a baiting and iterative mapping approach, Nucleic Acids Res 41(13) (2013) e129. [36] M. Valach, G. Burger, M.W. Gray, B.F. Lang, Widespread occurrence of organelle genome-encoded 5S rRNAs including permuted molecules, Nucleic Acids Res 42(22) (2014) 13764-77. [37] M. Bernt, A. Donath, F. Juhling, F. Externbrink, C. Florentz, G. Fritzsch, J. Putz, M. Middendorf, P.F. Stadler, MITOS: improved de novo metazoan mitochondrial genome annotation, Mol Phylogenet Evol 69(2) (2013) 313-9. [38] N.R. Coordinators, Database resources of the National Center for Biotechnology Information, Nucleic Acids Res (2017). [39] A.J. Bleasby, J.C. Wootton, Construction of validated, non-redundant composite protein sequence databases, Protein Eng 3(3) (1990) 153-9. [40] T.M. Lowe, P.P. Chan, tRNAscan-SE On-line: integrating search and context for analysis of transfer RNA genes, Nucleic Acids Res 44(W1) (2016) W54-7. [41] M. Lohse, O. Drechsel, S. Kahlau, R. Bock, OrganellarGenomeDRAW--a suite of tools for generating physical maps of plastid and mitochondrial genomes and visualizing expression data sets, Nucleic Acids Res 41(Web Server issue) (2013) W575-81. [42] J. Wang, L. Zhang, Q.L. Zhang, M.Q. Zhou, X.T. Wang, X.Z. Yang, M.L. Yuan, Comparative mitogenomic analysis of mirid bugs (Hemiptera: Miridae) and evaluation of potential DNA barcoding 27

markers, PeerJ 5 (2017) e3661. [43] P. Stothard, The sequence manipulation suite: JavaScript programs for analyzing and formatting protein and DNA sequences, Biotechniques 28(6) (2000) 1102, 1104. [44] J. Rozas, A. Ferrer-Mata, J.C. Sanchez-DelBarrio, S. Guirao-Rico, P. Librado, S.E. Ramos-Onsins, A. Sanchez-Gracia, DnaSP 6: DNA Sequence Polymorphism Analysis of Large Data Sets, Mol Biol Evol 34(12) (2017) 3299-3302. [45] J. Caspermeyer, MEGA Evolutionary Software Re-Engineered to Handle Today's Big Data Demands, Mol Biol Evol 33(7) (2016) 1887. [46] Y. Chen, W. Ye, Y. Zhang, Y. Xu, High speed BLASTN: an accelerated MegaBLAST search tool, Nucleic Acids Res 43(16) (2015) 7762-8. [47] G. Benson, Tandem repeats finder: a program to analyze DNA sequences, Nucleic Acids Res 27(2) (1999) 573-80. [48] S. Kurtz, J.V. Choudhuri, E. Ohlebusch, C. Schleiermacher, J. Stoye, R. Giegerich, REPuter: the manifold applications of repeat analysis on a genomic scale, Nucleic Acids Res 29(22) (2001) 4633-42. [49] A.C. Darling, B. Mau, F.R. Blattner, N.T. Perna, Mauve: multiple alignment of conserved genomic sequence with rearrangements, Genome Res 14(7) (2004) 1394-403. [50] C. Ferandon, S. Moukha, P. Callac, J.P. Benedetto, M. Castroviejo, G. Barroso, The Agaricus bisporus cox1 gene: the longest mitochondrial gene and the largest reservoir of mitochondrial group i introns, PLoS One 5(11) (2010) e14048. [51] Y. Deng, T. Hsiang, S. Li, L. Lin, Q. Wang, Q. Chen, B. Xie, R. Ming, Comparison of the Mitochondrial Genome Sequences of Six Annulohypoxylon stygium Isolates Suggests Short Fragment Insertions as a Potential Factor Leading to Larger Genomic Size, Front Microbiol 9 (2018) 2079. [52] K. Katoh, J. Rozewicki, K.D. Yamada, MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization, Brief Bioinform (2017). [53] G. Vaidya, D.L. Lohman, R. Meier, SequenceMatrix: concatenation software for the fast assembly of multi̺gene datasets with character set and codon information, Cladistics 27(2) (2011) 171-180. [54] X. Xia, DAMBE7: New and Improved Tools for Data Analysis in Molecular Biology and Evolution, Mol Biol Evol 35(6) (2018) 1550-1552. [55] R. Lanfear, P.B. Frandsen, A.M. Wright, T. Senfeld, B. Calcott, PartitionFinder 2: New Methods for Selecting Partitioned Models of Evolution for Molecular and Morphological Phylogenetic Analyses, Mol Biol Evol 34(3) (2017) 772-773. [56] F. Ronquist, M. Teslenko, P. van der Mark, D.L. Ayres, A. Darling, S. Hohna, B. Larget, L. Liu, M.A. Suchard, J.P. Huelsenbeck, MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space, Syst Biol 61(3) (2012) 539-42. [57] A. Stamatakis, RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies, Bioinformatics 30(9) (2014) 1312-3. [58] J.E. Stajich, S.K. Wilke, D. Ahren, C.H. Au, B.W. Birren, M. Borodovsky, C. Burns, B. Canback, L.A. Casselton, C.K. Cheng, J. Deng, F.S. Dietrich, D.C. Fargo, M.L. Farman, A.C. Gathman, J. Goldberg, R. Guigo, P.J. Hoegger, J.B. Hooker, A. Huggins, T.Y. James, T. Kamada, S. Kilaru, C. Kodira, U. Kues, D. Kupfer, H.S. Kwan, A. Lomsadze, W. Li, W.W. Lilly, L.J. Ma, A.J. Mackey, G. Manning, F. Martin, H. Muraguchi, D.O. Natvig, H. Palmerini, M.A. Ramesh, C.J. Rehmeyer, B.A. Roe, N. Shenoy, M. Stanke, V. Ter-Hovhannisyan, A. Tunlid, R. Velagapudi, T.J. Vision, Q. Zeng, M.E. Zolan, P.J. Pukkila, Insights into evolution of multicellular fungi from the assembled chromosomes of the mushroom Coprinopsis cinerea (Coprinus cinereus), Proc Natl Acad Sci U S A 107(26) (2010) 11889-94. 28

[59] Q. Li, Q. Wang, X. Jin, Z. Chen, C. Xiong, P. Li, J. Zhao, W. Huang, Characterization and comparison of the mitochondrial genomes from two Lyophyllum fungal species and insights into phylogeny of Agaricomycetes, Int J Biol Macromol 121 (2019) 364-372. [60] C.E. Bullerwell, G. Burger, B.F. Lang, A novel motif for identifying rps3 homologs in fungal mitochondrial genomes, Trends Biochem Sci 25(8) (2000) 363-5. [61] Q. Li, Q. Wang, X. Jin, Z. Chen, C. Xiong, P. Li, J. Zhao, W. Huang, The first complete mitochondrial genome from the family Hygrophoraceae (Hygrophorus russula) by next-generation sequencing and phylogenetic implications, Int J Biol Macromol 122 (2019) 1313-1320. [62] R. Yang, Y. Li, X. Song, L. Tang, C. Li, Q. Tan, D. Bao, The complete mitochondrial genome of the widely cultivated edible fungus Lentinula edodes, Mitochondrial DNA Part B 2(1) (2017) 13-14. [63] E.F. Formighieri, R.A. Tiburcio, E.D. Armas, F.J. Medrano, H. Shimo, N. Carels, A. Goes-Neto, C. Cotomacci, M.F. Carazzolle, N. Sardinha-Pinto, D.P. Thomazella, J. Rincones, L. Digiampietri, D.M. Carraro, A.M. Azeredo-Espin, S.F. Reis, A.C. Deckmann, K. Gramacho, M.S. Goncalves, J.P. Moura Neto, L.V. Barbosa, L.W. Meinhardt, J.C. Cascardo, G.A. Pereira, The mitochondrial genome of the phytopathogenic basidiomycete Moniliophthora perniciosa is 109 kb in size and contains a stable integrated plasmid, Mycol Res 112(Pt 10) (2008) 1136-52. [64] H. Yoon, W.S. Kong, Y.J. Kim, J.G. Kim, Complete mitochondrial genome of the ectomycorrhizal fungus Tricholoma matsutake, Mitochondrial DNA A DNA Mapp Seq Anal 27(6) (2016) 3855-3857. [65] Y.J. Park, J.H. Baek, S. Lee, C. Kim, H. Rhee, H. Kim, J.S. Seo, H.R. Park, D.E. Yoon, J.Y. Nam, H.I. Kim, J.G. Kim, H. Yoon, H.W. Kang, J.Y. Cho, E.S. Song, G.H. Sung, Y.B. Yoo, C.S. Lee, B.M. Lee, W.S. Kong, Whole genome and global gene expression analyses of the model mushroom Flammulina velutipes reveal a high capacity for lignocellulose degradation, PLoS One 9(4) (2014) e93560. [66] D. Formey, M. Moles, A. Haouy, B. Savelli, O. Bouchez, G. Becard, C. Roux, Comparative analysis of mitochondrial genomes of Rhizophagus irregularis - syn. Glomus irregulare - reveals a polymorphism induced by variability generating elements, New Phytol 196(4) (2012) 1217-27. [67] L. Forget, J. Ustinova, Z. Wang, V.A. Huss, B.F. Lang, Hyaloraphidium curvatum: a linear mitochondrial genome, tRNA editing, and an evolutionary link to lower fungi, Mol Biol Evol 19(3) (2002) 310-9. [68] L. Losada, S.B. Pakala, N.D. Fedorova, V. Joardar, S.A. Shabalina, J. Hostetler, S.M. Pakala, N. Zafar, E. Thomas, M. Rodriguez-Carres, R. Dean, R. Vilgalys, W.C. Nierman, M.A. Cubeta, Mobile elements and mitochondrial genome expansion in the soil fungus and potato pathogen Rhizoctonia solani AG-3, FEMS Microbiol Lett 352(2) (2014) 165-73. [69] Z. du Toit, M. du Plessis, D.L. Dalton, R. Jansen, J. Paul Grobler, A. Kotze, Mitochondrial genomes of African pangolins and insights into evolutionary patterns and phylogeny of the family Manidae, BMC Genomics 18(1) (2017) 746. [70] X. Liang, X. Tian, W. Liu, T. Wei, W. Wang, Q. Dong, B. Wang, Y. Meng, R. Zhang, M.L. Gleason, G. Sun, Comparative analysis of the mitochondrial genomes of Colletotrichum gloeosporioides sensu lato: insights into the evolution of a fungal species complex interacting with diverse plants, BMC Genomics 18(1) (2017) 171. [71] K.L. Adams, J.D. Palmer, Evolution of mitochondrial gene content: gene loss and transfer to the nucleus, Mol Phylogenet Evol 29(3) (2003) 380-95. [72] I.G. Johnston, B.P. Williams, Evolutionary Inference across Eukaryotes Identifies Specific Pressures Favoring Mitochondrial Gene Retention, Cell Syst 2(2) (2016) 101-11. [73] P. Bjorkholm, A. Harish, E. Hagstrom, A.M. Ernst, S.G. Andersson, Mitochondrial genomes are 29

retained by selective constraints on protein targeting, Proc Natl Acad Sci U S A 112(33) (2015) 1015461. [74] J.F. Allen, Why chloroplasts and mitochondria retain their own genomes and genetic systems: Colocation for redox regulation of gene expression, Proc Natl Acad Sci U S A 112(33) (2015) 10231-8. [75] L. Lin, P. Cui, Z. Qiu, M. Wang, Y. Yu, J. Wang, Q. Sun, H. Zhao, The mitochondrial tRNA(Ala) 5587T>C and tRNA(Leu(CUN)) 12280A>G mutations may be associated with hypertension in a Chinese family, Exp Ther Med 17(3) (2019) 1855-1862. [76] Y. Ding, Y.S. Teng, G.C. Zhuo, B.H. Xia, J.H. Leng, The Mitochondrial tRNAHis G12192A Mutation May Modulate the Clinical Expression of Deafness-Associated tRNAThr G15927A Mutation in a Chinese Pedigree, Curr Mol Med (2019). [77] B.P. Tang, Z.Z. Xin, Y. Liu, D.Z. Zhang, Z.F. Wang, H.B. Zhang, X.Y. Chai, C.L. Zhou, Q.N. Liu, The complete mitochondrial genome of Sesarmops sinensis reveals gene rearrangements and phylogenetic relationships in Brachyura, PLoS One 12(6) (2017) e0179800. [78] G. Aguileta, D.M. de Vienne, O.N. Ross, M.E. Hood, T. Giraud, E. Petit, T. Gabaldon, High variability of mitochondrial gene order among fungi, Genome Biol Evol 6(2) (2014) 451-65. [79] Y. Xia, Y. Zheng, R.W. Murphy, X. Zeng, Intraspecific rearrangement of mitochondrial genome suggests the prevalence of the tandem duplication-random loss (TDLR) mechanism in Quasipaa boulengeri, BMC Genomics 17(1) (2016) 965. [80] D.V. Lavrov, J.L. Boore, W.M. Brown, Complete mtDNA sequences of two millipedes suggest a new model for mitochondrial gene rearrangements: duplication and nonrandom loss, Mol Biol Evol 19(2) (2002) 163-9. [81] W.W. Fan, S. Zhang, Y.J. Zhang, The complete mitochondrial genome of the Chan-hua fungus Isaria cicadae: a tale of intron evolution in Cordycipitaceae, Environ Microbiol 21(2) (2019) 864-879. [82] J. Repar, T. Warnecke, Mobile Introns Shape the Genetic Diversity of Their Host Genes, Genetics 205(4) (2017) 1641-1648. [83] A.V. Mardanov, A.V. Beletsky, V.V. Kadnikov, A.N. Ignatov, N.V. Ravin, The 203 kbp mitochondrial genome of the phytopathogenic fungus Sclerotinia borealis reveals multiple invasions of introns and genomic duplications, PLoS One 9(9) (2014) e107536. [84] C. Ferandon, J. Xu, G. Barroso, The 135 kbp mitochondrial genome of Agaricus bisporus is the largest known eukaryotic reservoir of group I introns and plasmid-related sequences, Fungal Genet Biol 55 (2013) 85-91. [85] L. Vincenot, F. Popa, F. Laso, K. Donges, K.H. Rexer, G. Kost, Z.L. Yang, K. Nara, M.A. Selosse, Out of Asia: Biogeography of fungal populations reveals Asian origin of diversification of the Laccaria amethystina complex, and two new species of violet Laccaria, Fungal Biol 121(11) (2017) 939-955. [86] N. Corradi, L. Bonen, Mitochondrial genome invaders: an unselfish role as molecular markers, New Phytol 196(4) (2012) 963-5. [87] P. Johri, G.K. Marinov, T.G. Doak, M. Lynch, Population genetics of Paramecium mitochondrial genomes: recombination, mutation spectrum, and efficacy of selection, Genome Biol Evol (2019). [88] L. Wang, S. Zhang, J.H. Li, Y.J. Zhang, Mitochondrial genome, comparative analysis and evolutionary insights into the entomopathogenic fungus Hirsutella thompsonii, Environ Microbiol 20(9) (2018) 3393-3405.

30

31

Figure Legends

Figure 1. Circular maps of the mitochondrial genomes of two Laccaria species. Genes are represented by different colored blocks. Colored blocks outside each ring indicate that the genes are on the direct strand, while colored blocks within the ring indicates that the genes are located on the reverse strand.

Figure 2. The protein-coding, intronic, intergenic, and RNA gene region proportions of the entire mitochondrial genomes of the two Laccaria species. The bottom panel shows the contribution of different gene regions to the expansion of the L. bicolor mitogenome.

Figure 3. Gene order comparison between 7 Agaricales mitogenomes. Genes are represented with different color blocks. All genes are shown in order of occurrence in the mitochondrial genome, starting from cox1. Fourteen core protein coding genes, one rps3 gene, and two rRNA genes were included in the gene arrangement analysis. Species and NCBI accession number used for gene arrangement analysis in the present study are listed in Supplementary Table S7.

Figure 4. Co-linearity analysis of 7 Agaricales mitogenomes. Twenty homologous regions were detected amomg the 7 mitogenomes. The sizes and relative positions of the homologous regions varied between the mitogenomes. 32

Figure 5. Variation in the length and base composition of each of 15 protein-coding genes (PCGs) and 25 tRNA genes between two Laccaria mitogenomes. a, PCG length variation; b, GC content of the PCGs; c, AT skew; d, GC skew; e, lengths of shared tRNA genes ; d, GC content of shared tRNA genes.

Figure 6. Genetic analysis of 15 protein coding genes conserved in two Laccaria mitogenomes. K2P, the Kimura-2-parameter distance; Ka, the mean number of nonsynonymous substitutions per nonsynonymous site; Ks, the mean number of synonymous substitutions per synonymous site.

Figure 7. Pcl information of cox1 gene of the 17 Agaricales species. The same Pcl (orthologous intron) is represented by the same letter. The phylogenetic positions of 17 Agaricales species were established using the Bayesian inference (BI) method and Maximum Likelihood (ML) method based on 15 concatenated mitochondrial core proteins and 2 rRNA genes. The II in the figure above shows that the intron belongs to the group II intron. Species ID are shown in Supplementary Table S7.

Figure 8. Molecular phylogeny of 58 Basidiomycota species based on Bayesian inference (BI) and Maximum likelihood (ML) analysis of 15 protein coding genes and two rRNA genes. Support values are Bayesian posterior probabilities (before slash) and bootstrap (BS) values. Species and NCBI accession numbers for genomes used in the phylogenetic analysis are provided in Supplementary Table S7.

Supplementary figures

33

Figure S1. Putative secondary structures of the 25 tRNA genes identified in the mitochondrial genomes of two Laccaria species. Residues conserved across the two mitochondrial genomes are shown in green, while variable sites are shown in red. All genes are shown in order of occurrence in the mitochondrial genome of L. amethystine, starting from trnN.

Figure S2. Codon usage in the mitochondrial genomes of two Laccaria species. Frequency of codon usage is plotted on the y-axis. a, L. amethystine; b, L. bicolor.

34

35

36

37

38

39

40