Strobilomyces alpinus strain H2022S10 was collected from Pudacuo Nature Reserve in Shangri-La, Yunnan, China, on 6 August 2022. The specimens were collected with the utmost care and stored at ultra-low temperature freezer at -80 °C. Identification was conducted, based on a combination of morphological characteristics and multigene phylogenetic analysis by Li-Hong Han. The specimen has been deposited in the Herbarium of the College of Biological Resources and Food Engineering at Qujing Normal University, under the accession number H2022S10.

Additionally, four unassembled sequencing datasets of Strobilomyces species were retrieved from the Sequence Read Archive database. The four species include S. echinocephalus from China (accession number SRR28540251), an unidentified species, Strobilomyces sp. 1 (SRR28540241) from the USA, S. confusus from Vietnam (SRR28540078) and another unidentified species, Strobilomyces sp. 2 (SRR28540220) from Vietnam (Tremble et al. 2024).

Whole-genomic DNA was extracted from the unpolluted context of the stipe using the CTAB method (Doyle and Dickson 1987). Libraries were constructed by MGIEasy Universal DNA Library Prep Kit v.1.0, following the standard protocol. Briefly, 1 μg of genomic DNA was randomly fragmented by Covaris and fragments averaging 200–400 bp in size were selected with MGIEasy DNA Clean beads. Subsequently, the fragments were end-repaired, 3’ adenylated and ligated with adapters. The DNA samples were amplified via PCR and the products were purified by the MGIEasy DNA Clean beads. The double stranded PCR products were heat denatured and circularised by the splint oligo sequence in MGIEasy Circularisation Module. The single strand circular DNA was formatted as the final library and assessed for quality control. Sequencing was then performed on a DNBSEQ-T7RS platform, generating 150 bp paired-end reads by GrandOmics Biosciences Co., Ltd. (Wuhan, China), yielding over 5.0 GB of raw data. The mitogenome sequences were assembled using GetOrganelle v.1.7.5 (Jin et al. 2020) with the filtered data. Annotation of the mitogenomes was primarily carried out using GeSeq (Tillich et al. 2017) and MITOS (http://mitos.bioinf.uni-leipzig.de/index.py) , following the mould mitochondrial genetic code, with manual corrections made where necessary. The graphical representation of the mitogenome was visualised using Chloroplot (Zheng et al. 2020).

The strand asymmetries in the five Strobilomyces mitogenomes were evaluated using the formulae: AT skew = [A - T] / [A + T] and GC skew = [G - C] / [G + C]. The codon usage frequency and preference in the S. alpinus mitogenome were analysed using the Sequence Manipulation Suite (https://detaibio.com/sms2/index.html), based on genetic code 4. The RNA editing sites in 15 PCGs of S. alpinus were predicted using the PmtREP programme (http://112.86.217.82:9929/#/tool/alltool/detail/336), with a threshold value of 0.2. The pairwise genetic distances amongst the 14 core PCGs, excluding rps3 due to incompleteness in the mitogenomes Strobilomyces sp. 1 and S. echinocephalus, were calculated using the Kimura-2-parameter (K2P) substitution model in MEGA X (Kumar et al. 2018). Furthermore, the non-synonymous (Ka) and synonymous (Ks) substitution rates for the core PCGs in the five Strobilomyces mitogenomes were estimated using DnaSP v.6 (Rozas et al. 2017). To identify simple sequence repeats (SSRs) within the mitogenomes, we utilised the MISA programme (Beier et al. 2017) with predefined thresholds for the minimum number of repeats: 10 for mononucleotides, five for dinucleotides, four for trinucleotides and three for tetranucleotides, pentanucleotides and hexanucleotides. For the detection of tandem repeats, we employed Tandem Repeats Finder v.4.09 software (Benson 1999) with default parameters. Additionally, we identified forward, reverse, palindromic and complementary repeats based on REPuter (Kurtz et al. 2001). We set the maximum number of computed repeats to 5,000, the minimum repeat size to 30 and the Hamming distance to 3.

To gain insights into the evolution of mitogenomes in Strobilomyces, we conducted a comparative analysis of five mitogenomes. This analysis included evaluating genome size, gene content, gene order, genetic distance and selection pressure. For a more detailed visualisation of syntenies and potential gene re-arrangement events, we performed a mitogenome collinearity analysis using Mauve (Stajich et al. 2010).

To determine the phylogenetic position of Strobilomyces within Boletales, we assembled two datasets based on nucleotide sequences. The first dataset comprised concatenated sequences of 14 core PCGs from mitogenome, excluding the rps3 gene due to incompleteness in two of the mitogenomes. The second dataset consisted of nuclear ribosomal DNA (nrDNA) sequences, specifically ETS-18S-ITS1-5.8S-ITS2-26S. We utilised 36 representative mitogenomes and 85 nrDNA sequences from Boletales taxa, along with two outgroup species from Russulales: Russula compacta and Lactarius deliciosus (Suppl. material 1: table S1). Sequence alignment was performed using MAFFT (Katoh et al. 2019) with default parameters. Phylogenetic relationships were inferred by both Maximum Likelihood (ML) and Bayesian Inference (BI) approaches, with parameters specified in our previous study (Yu et al. 2023a). The GTR+F+R3 substitution model was selected for the mtPCGs dataset, while the GTR+F+R8 model was applied to the nrDNA trees. To further assess the phylogenetic topology, we employed the BI method in MrBayes v.3.1.2 (Ronquist et al. 2012). Four Markov Chain Monte Carlo (MCMC) chains were run simultaneously for 20 million generations, sampling trees every 100 generations. In evaluating phylogenetic confidence, we considered nodes to be strongly supported by both ML bootstrap values (MLB ≥ 70%) and Bayesian posterior probabilities (BPP ≥ 0.95).

BI Bayesian Inference

K2P Kimura-2-parameter

Ka Non-synonymous substitution rate

Ks Synonymous substitution rate

Mitogenome Mitochondrial genome

ML Maximum Likelihood

nrDNA Nuclear ribosomal DNA

PCG Protein-coding gene

RSCU Relative synonymous codon usage

SSRs Simple sequence repeats

The five Strobilomyces species exhibited complete mitogenome sequences ranging from 35,618 bp in S. sp. 2 to 42,088 bp in S. alpinus (Table 1). The comparative analysis indicated that all five mitogenomes possessed circular DNA molecules. The nucleotide composition revealed a higher abundance of AT (77.3%) compared to GC (22.7%), with 37.7% A, 39.5% T, 11.7% G and 11.0% C. The negative AT skew and positive GC skew (except in S. sp. 1) indicated a higher frequency of T and G over A and C in the forward strand, respectively, which is consistent with the overall nucleotide ratio.

The mitogenomes encoded 14 core PCGs essential for energy metabolism. These include three ATP synthases, three cytochrome c oxidases, seven NADH dehydrogenases and one ubiquinol cytochrome c reductase. Additionally, the gene rps3 encodes a ribosomal protein, which is involved in translation, along with two rRNAs and 24 tRNAs, were identified. In the S. alpinus mitogenome, most genes—excluding atp6, atp8, cox1 and nine tRNA genes—are transcribed on the forward strand (Fig. 1).

Figure 1. 

Circular map of the mitogenome of S. alpinus, with genes colour-coding according to their functional groups. The inner ring illustrates the GC content.

Amongst the amino acids found in the mitogenomes of Strobilomyces species, leucine, isoleucine and serine were the most frequently utilised. Specifically, leucine (14.8%) was the most prevalent, followed by isoleucine (12.1%) and serine (8.71%). Cysteine (0.6%) and tryptophan (1.1%) were the least frequently encountered residues (Suppl. material 2: fig. S1). Notably, the codon usage patterns exhibited remarkable consistency across the Strobilomyces species.

For further analysis, we focused on the mitogenome of S. alpinus to compare the distribution and frequency of codon usage. All core PCGs and the ribosomal protein-coding gene rps3 in the S. alpinus mitogenome possessed the canonical AUG start codon and UAA stop codon. However, in the atp8 and nad4L genes from two Vietnamese samples (S. confusus and S. sp. 2), the UAG codon functioned as the stop codon.

The relative synonymous codon usage (RSCU) was calculated for the 15 PCGs in the mitogenome of S. alpinus. The RSCU values of 27 codons exceeded 1.00, with 14 codons ending in U, 12 ending in A and only one ending in G (Fig. 2). This observation indicates a strong bias towards NNA and NNT codons when available. Amongst these preferred codons, despite arginine being encoded by six different codons, only AGA (with an RSCU value 6.0) was utilised to encode this amino acid. The codon for leucine (UUA) was the second most preferred, with an RSCU value of 5.7. Seven codons (AGG, CGA, CGC, CGG, CGU, CUC and UGA) were never used, while five codons (GGC, CUG, AAG, ACG and UAG) were used in some Strobilomyces species at a very low frequency. Notably, the tryptophan codon UGA, which functions as a stop codon in the standard genetic code, was not utilised in the mitogenome of S. alpinus.

Figure 2. 

Codon usage for 15 PCGs in the mitogenome of S. alpinus. The 15 PCGs consist of 14 typical PCGs and rps3.

In our analysis of the mitogenome of S. alpinus, we predicted 41 RNA editing sites that involve the conversion of cytosine (C) to uridine (U) within 12 of the 15 PCGs (Fig. 3). Notably, no RNA editing was observed in three genes: atp6, atp8 and nad4L. The frequency of these editing events varied significantly amongst the PCGs. The nad5 gene exhibited the highest number of RNA editing sites (7), followed by the nad2 and nad4 genes, each with six sites. Amongst these editing sites, 26.8% (11) were located at the first position of the codon triplet, while 73.2% (30) occurred at the second position.

After RNA editing, the hydrophobicity of 43.9% of the amino acids remained unchanged. However, 39.0% of the amino acids were predicted to transition from hydrophilic to hydrophobic, while 17.1% were predicted to change from hydrophobic to hydrophilic. The most frequent amino acid substitution, occurring 12 times, was from alanine to valine, while the rarest substitution was from leucine to phenylalanine, occurring only once amongst all editing events.

Figure 3. 

Prediction of RNA editing site numbers and amino acid changes.

As the rps3 gene was frequently absent or incomplete in the mitogenomes, we focused on 14 core PCG genes: atp6, atp8, atp9, cob, cox1, cox2, cox3, nad1, nad2, nad3, nad4, nad5, nad6 and nad4L to calculate genetic distances and substitution rates amongst the five Strobilomyces species. Our results showed that the genetic distance ranged from 0.01 (atp9) to 0.05 (atp8) based on the K2P model (Fig. 4A), demonstrating a high degree of conservation amongst these genes. Due to multiple instances of zero Ka or Ks values in atp8 and atp9, only 12 remaining genes were included in the analysis. Amongst these 12 genes, the nad6 gene exhibited the highest Ka value, followed by nad2, while cox2 showed the lowest Ka value (Fig. 4B). The nad4 gene had the highest Ks value, whereas nad4L displayed the lowest Ks value (Fig. 4C). The overall ratio of Ka/Ks for all detected genes was below one, ranging from 0.03 (cox2) to 0.27 (nad6) (Fig. 4D). This finding indicates that these genes have undergone purifying selection throughout their evolution. It is important to note that when comparing Ka/Ks ratios, three genes—nad6, nad2 and atp6—exhibited significantly higher ratios (Ka/Ks > 0.1), while three other genes—cox2, cob and cox1—showed notably lower ratios (Ka/Ks < 0.05).

Figure 4. 

Genetic and evolutionary analyses of core genes in the five Strobilomyces mitogenomes. A Kimura-2-parameter (K2P) distance, B number of non-synonymous substitutions per non-synonymous site (Ka), C number of synonymous substitutions per synonymous site (Ks), D Ka/Ks ratio.

Figure 5. 

Number of SSRs (A) and long repeats (B) in the five Strobilomyces mitogenomes.

To gain a deeper understanding of the characteristics of repeat sequences, we identified SSRs, tandem repeats and dispersed repeats within the mitogenomes of Strobilomyces in this study. Our results showed that the number of SSRs detected ranged from 30 in S. alpinus to 49 in S. echinocephalus (Suppl. material 1: table S2 and Fig. 5A). Generally, mononucleotide and tetranucleotide repeats were the most prevalent, followed by trinucleotide repeats, while hexanucleotide repeats were the least common. Mononucleotide repeats of A/T were found to be the most abundant type of repeat. Notably, only three instances of A/T repeats, each consisting of 10 nucleotides, were found in S. alpinus. Mononucleotide repeats of C/G were exceedingly rare, with only one instance of an 11-nucleotide repeat identified in S. sp. 2. Dinucleotide repeats of AT/AT were detected between one and three times across the Strobilomyces mitogenomes, while no AG/CT repeats were observed. It is intriguing to note that, despite possessing the largest mitogenome (42,088 bp) amongst the studied species, S. alpinus had the fewest SSRs.

The analysis of interspersed repetitive sequences revealed an uneven distribution amongst the Strobilomyces mitogenomes (Fig. 5B). The number of repeats ranged from 7 in S. echinocephalus to 25 in both S. confusus and S. sp. 2. Forward repeats, which ranged from 2 to 19, were the most common type, followed by palindromic repeats, ranging from 3 to 10. Reverse and complementary repeats occurred one to two times in each mitogenome. Notably, 91% of the repetitive sequences were between 30 and 49 bp in length, while only 9% ranged from 50 to 80 bp.

Furthermore, the analysis of tandem repeats revealed their presence in each mitogenome, with quantities ranging from 8 to 18 varying between 25 and 86 bp (Suppl. material 1: table S3 and Fig. 3). Most of these tandem repeats consisted of two copies, with 71% measuring between 25 and 50 bp and 29% measuring between 51 and 86 bp (Suppl. material 1: table S3). Interestingly, the number of tandem repeats did not correlate with the size of the mitogenome.

To conduct an intensive study of the evolution of Strobilomyces mitogenomes, a comparative analysis was performed amongst five species. The Mauve alignment highlighted a conserved synteny amongst the mitogenomes, which were organised into 11 homologous regions represented by distinct coloured synteny blocks (Fig. 6A). This alignment revealed a high degree of synteny across the mitogenomes. However, two instances of gene re-arrangement were identified. Specifically, the gene linkage comprising trnG, atp9, trnD, trnW, trnK, trnQ, trnT, trnF, trnA, cox2 and rnl underwent a reverse re-arrangement in the mitogenome of Strobilomyces sp. 1 from the USA. Another re-arrangement, involving the gene linkage of trnH, trnM and atp6, was observed in two mitogenomes—S. confusus and Strobilomyces sp. 2—both originating from Vietnam (Fig. 6B).

Figure 6. 

Collinearity (A) and gene order comparison (B) amongst five Strobilomyces mitogenomes. Genes (green), rRNA (orange) and tRNA (purple) are represented in their respective colours, while gene re-arrangements across the mitogenomes are indicated by red boxes.

Two phylogenetic construction methods, BI and ML, produced consistent topological structures. However, the phylogenetic trees constructed from different datasets, specifically mtPCGs and nrDNA, exhibited distinct topologies (Figs 7, 8). The summary coalescent tree revealed major nodes with substantial statistical support.

Overall, the concatenated mtPCGs and nrDNA phylogenies exhibited congruent topologies at the family level, with a few exceptions. Species within the same family clustered together, including Russulaceae, Coniophoraceae, Gomphidiaceae, Rhizopogonaceae, Sclerodermataceae, Paxillaceae and Boletaceae. One notable exception was the family Boletinellaceae, which formed a monophyletic clade in the phylogeny of mtPCGs, but was integrated into the Boletaceae family in the phylogeny of nrDNA. The phylogenetic relationships at the family level also revealed subtle differences between the mtPCGs and nrDNA trees (Figs 7, 8). For instance, Strobilomyces species clustered together to form a sister branch alongside other Boletaceae species in the mtPCGs tree. In contrast, in the nrDNA tree, the Strobilomyces species were positioned at the base of Boletaceae. Boletaceae represents a complex clade characterised by conflicting phylogenetic markers. Despite this, fungal species within each family clustered together with high support values (Fig. 7).

Many of the recovered groups within the Boletales are consistent with the previous studies (Tremble et al. 2024). In our analysis of the mitogenome sizes, we found that the families situated at the base of the phylogenetic tree tend to possess larger mitogenomes. As species evolve, the size of their mitochondria generally decreases.

Figure 7. 

Phylogenetic analysis of Boletales species, based on concatenated nucleotide sequences of 14 typical mtPCGs. The tree presented here is the single best topology recovered from Maximum Likelihood analysis. Support values for the nodes are indicated as ML/BI with an asterisk (*) denoting a value of 1 or 100. The phylogenetic clades are colour-coded at the family level, with dark blue representing the genus Strobilomyces.

Figure 8. 

Phylogenetic relationships amongst representative specimens of Boletaceae, inferred from the nrDNA dataset using Maximum Likelihood and Bayesian Inference methods (only the ML tree is displayed). The phylogenetic clades are distinguished by different colours at the family level, with dark blue representing the genus Strobilomyces.

The authors have declared that no competing interests exist.

No ethical statement was reported.

All the fungal strains used in this study have been legally obtained, respecting the Convention on Biological Diversity (Rio Convention).

This study was funded by the National Natural Science Foundation of China (32100010 and 32060710) and Program of Doctoral Innovation Research Team from Qujing Normal University.

CL: Writing – original draft, review and editing, Conceptualisation, Methodology, Software, Formal Analysis, Funding acquisition; WYL: Investigation, Resources; LXZ: Investigation, Resources; MD: Investigation, Resources; HHC: Writing – review and editing, Methodology, Supervision; LHH: Writing – original draft, review and editing, Conceptualisation, Methodology, Investigation, Resources, Project administration, Funding acquisition. All authors have read and approved the final version of the manuscript.

Chao Liu https://orcid.org/0000-0001-6811-2218

Huan-Huan Chen https://orcid.org/0000-0001-9477-1014

Li-Hong Han https://orcid.org/0000-0002-6127-0915

All data used in this study are publicly available. The assembled mitogenome sequence of S. alpinus was submitted to the China National Center for Bioinformation with accession number C_AA098151.1 (https://ngdc.cncb.ac.cn/genbase/) and National Center for Biotechnology Information with accession number PQ568401 (https://www.ncbi.nlm.nih.gov).