High-quality assembly of the Chinese white truffle genome and recalibrated divergence time estimate provide insight into the evolutionary dynamics of Tuberaceae

A high-quality assembly and gleba-associated microbiome of T. panzhihuanense

Using k-mer-based approaches, we estimated a genome size of 109 Mb, compatible with estimates of other Tuber species (Murat et al. 2018). As expected for the haploid state of the gleba fruiting body, the genome was completely homozygous (Supplemental Fig. S1A). The assembled genome spans 119 Mb with an N50 value of 7.3 Mb, a consensus quality value (QV) of 57.72, and complete single-copy BUSCO scores of 96.4%. (Supplemental Table S1; Supplemental Fig. S1B).

Among the 29 assembled contigs, contigs 1–17 encompass almost the entire nuclear genome (Fig. 1A). Contig 18 represents the 302 kb mitochondrial genome, for which we obtained a complete and circularized sequence, whereas contigs 19–29 contain only 45S rDNA arrays. We found an enrichment of TTAGGG telomere repeats, typical of many eukaryotes, including true truffles (Martin et al. 2010), toward both ends of all 17 majors nuclear contigs, except for contigs 16 and 17, which show enrichment only at one end (Supplemental Fig. S2). These two contigs carry 45S rDNA loci at the ends lacking telomeric repeats (see Results section “45S rDNA genes are organized as a single tandem array with only one diverging locus dispersed throughout the genome”), suggesting that they may represent a single chromosome. These findings indicate that our assembly is close to chromosome level, with a haploid chromosome number of 1n = 16.

Figure 1.

A highly compartmentalized genome colonized by Chromoviridae-related Gypsy transposons. (A) From outer to inner circles: distribution of transposons, genes, TE hotspots, and TE coldspots across the 17 main contigs of the T. panzhihuanense nuclear genome. (B) Correlation between TE density against gene density across 50 kb nonoverlapping genomic windows; rho = Spearman's rank correlation coefficient. (C) Phylogenetic tree of representative Gypsy RT protein segments longer than 100 amino acids. Green circles highlight crown nodes of identified families with bootstrap support values of 75 or more. (D) Genomic occurrence of the four identified Gypsy clades in the T. panzhihuanense genome.

Open in new tabDownload PowerPointLink to figure

On the 17 nuclear contigs, we predicted 8701 protein-coding genes (PCGs) with a complete and single-copy BUSCO score slightly higher than the assembly score (Supplemental Table S1). The T. panzhihuanense genome encodes for 363 secreted proteins (SPs), of which 173 are small SPs (SSPs), and a reduced set of 198 carbohydrate-active enzymes (CAZys) (Supplemental Table S2).

Microbial profiling (Supplemental Fig. S3A) and Blobtools analyses (Supplemental Fig. S3B) for T. panzhihuanense, along with other publicly available truffle fruiting bodies, reveal the exclusive presence of ascocarp-associated bacteria of the gleba. Bacteria of the genus Bradyrhizobium were ubiquitously detected across all species and samples (Supplemental Fig. S3A). Their association with Tuber spp. mycelia and fruiting bodies was already described and suggested to be potentially involved in nitrogen fixation (Monaco et al. 2022; Graziosi et al. 2024).

A highly compartmentalized genome dominated by retrotransposons

Manual curation of repeat libraries is considered an important step for accurate TE annotation in nonmodel species (Goubert et al. 2022; Peona et al. 2024). To improve TE annotation in T. panzhihuanense within a reasonable time frame, we manually curated the automatically generated libraries by selecting the most represented consensus sequences and/or those containing detectable protein fragments. Although the curated consensus sequences represent only 18% of the partially curated library, we found that they account for 54% of the total TE content, demonstrating that our selection strategy successfully captured most of the TE complement. Although the total TE content did not change much when using the raw versus the partially curated TE library, resulting in both instances >60%, we observed a drastic reduction in the number of unclassified elements using the former (Supplemental Fig. S4A,B). Additionally, the curation process yielded significantly longer consensus sequences compared with the uncurated ones (Supplemental Fig. S4C), suggesting a more accurate reconstruction of their complete structure (Peona et al. 2024). Finally, we managed to recover multiple consensus sequences belonging to the Plavaka group of the CMC/CACTA superfamily that were misclassified as LTR/Pao by the automatic annotation of RepeatModeler. Closer inspection of these consensus sequences revealed the presence of 2 bp target site duplications and terminal inverted repeats characterized by the motif 5′-CATAC-3′, which are features typical of Plavaka transposons previously identified in other fungal species (Supplemental Fig. S4D; Iyer et al. 2014; Hiltunen et al. 2021).

Regarding the general TE composition, Gypsy long terminal repeats (LTRs) and long interspersed nuclear element (LINE) Tad1 greatly dominate the genomic TE landscape of the Chinese white truffle (Supplemental Fig. S4B). TEs and genes have a strong and significant negative correlation between their densities (Spearman's rho = −0.92; P < 0.001) (Fig. 1A,B). Moreover, the majority of the assembly is composed of either significantly TE-enriched (TE hotspots; 51% of the genome; mean TE content = 96%) or TE-depleted genomic regions (TE coldspots; 32% of the genome; mean TE content = 15%) (Supplemental Table S3), with TE hotspots significantly longer than TE coldspots (Welch two-sample t-test, P < 0.001) (Supplemental Fig. S5). TE hotspots are richer in long retrotransposon insertions, whereas DNA elements are predominant transposons in TE coldspots (Supplemental Fig. S6). Moreover, the latter genomic compartment contains almost all PCGs (7047), with a gene density of 0.32, 2.6 times higher than the genome-wide estimation. These results are mirrored by the distance of TEs from the closest gene (Supplemental Fig. S7A), with DNA transposons more closely associated with genes compared with LINEs and LTRs (Kruskal–Wallis rank-sum test, P < 0.001; pairwise Wilcoxon rank-sum test with Bonferroni correction, P < 0.001), likely owing to the existence of numerous short nonautonomous and/or highly degenerated versions (Supplemental Fig. S7B,C). The classification of almost the entire genome into these two compartments, along with their striking differences in TE composition, strongly suggests a highly compartmentalized genomic organization. Effector-like SSPs showed no significant enrichment within TE hotspots (sample size = 8701; X² = 1.6749, P = 0.25).

Chromoviridae-related Gypsy elements dominate the repeat landscape in the Chinese white truffle genome

Because of the significant contribution of Gypsy elements to the high TE content in the T. panzhihuanense genome, as well as in other true truffles (Payen et al. 2016; Murat et al. 2018b), we further investigated their phylogenetic placement and domain structure. Based on network clustering (Supplemental Fig. S8A) and phylogenetic analyses, we identified 10 evolutionary distinct families (Fig. 1C; Supplemental Fig. S8B). These families belong to four distinct clades: Skipper-Like, TCN-Like, Tmt1, and a rich and diversified clade that we named Prometheus (Fig. 1C; Supplemental Fig. S9). The two former clades correspond to two different elements of the Tmt6 clade founded by Payen et al. (2016), whereas the Prometheus elements correspond to the Tmt3 clade. The Tmt2, Tmt4, and Tmt5 clades identified in the work of Payen et al. (2016) seem to be absent from T. panzhihuanense. Similarly to Tmt1 and TCN-like elements, Prometheus LTRs possess a chromodomain (CHD) and belong to the Chromoviridae Gypsy branch (Supplemental Fig. S10). On the other hand, the identified Skipper-like element, despite its placement with high support in a sister relationship with the reference Skipper transposon (bootstrap = 88) (Supplemental Fig. S9), is lacking the characteristic CHD domain (Supplemental Fig. S10; Marín and Lloréns 2000). We found that Prometheus is the most abundant Gypsy clade (Fig. 1D) followed by Tmt1, TCN-like, and Skipper-like (Fig. 1D). Although most of the Gypsy insertions are composed of degenerated and/or fragmented elements (Supplemental Fig. S11), we also observed 3066 insertions that almost perfectly match LTR regions of their parent consensus sequence, resembling solo-LTR elements that arise from nonallelic homologous recombination (NAHR) between the two flanking LTRs (Kent et al. 2017).

45S rDNA genes are organized as a single tandem array with only one diverging locus dispersed throughout the genome

The ribosomal-only contigs 19 to 29 (Fig. 2A) display a conserved structural pattern characterized by 45S rDNA genes (18S–5.8S–28S) flanked by a large array of tandemly repeated elements (Fig. 2B,C). We could identify other characteristic motifs: a 160 bp sequence at the 3′ end of each rDNA unit with strong similarity to LTR regions of Prometheus Gypsy transposon, and, at the 5′ end before the main tandemly repeated region, a 100 bp palindromic structure followed by a smaller GC-rich tandem repeat (Fig. 2C). Beside these 45S rDNA-only contigs, we identified four complete and two partial 45S loci with the same structure at one end of contig 16 and contig 17, as well as an additional one isolated along contig 3, 2.8 Mb away from the closest contig end (Fig. 2A). On the other hand, 5S rDNA genes were only identified outside 45S rDNA loci, with a total of 28 copies scattered along multiple contigs (Fig. 2A).

Figure 2.

Genomic features of rDNA loci. (A) Genomic locations of nuclear rDNA genes. The size of the contigs is proportional to their length, but different scales are used for contigs ctg_1 to ctg_17 and for contigs ctg_19 to ctg_29. The length of each contig is reported below its name. (B) Self-alignment of ctg_19 as an example of the 45S rDNA locus structure. (C) Simplified structure of 45S rDNA loci; the length of the blocks is not proportional to their actual length. (IR) Inverted repeat, (DR) direct repeat, (IGS) intergenic spacer, and (ITS1 and ITS2) internal transcribed spacers. (D) Alignment of 45S rDNA loci from 18S to 28S. Red lines highlight nucleotides that differ from the consensus sequence. The names of each locus reflect the naming system used in panel A. (E) Same as D, but for 5S rDNA genes.

Open in new tabDownload PowerPointLink to figure

We aligned all 45S rDNA loci and all 5S rDNA genes to evaluate a possible intragenomic variation. 45S rDNA loci identified from contig 16 to contig 29 are highly homogenized, without any insertions or deletions and only seven positions harboring variants unique or shared across paralog copies (Fig. 2D). Instead, the intergenic spacer (IGS) shows variable lengths ranging from 6312 bp to 8510 bp, mainly owing to different copy numbers in the main tandem array (Fig. 2B). Contrary to other 45S loci, rDNA paralog genes found isolated on contig 3 are more diverging (Fig. 2D). Compared with the 18S_10 and 28S_10 rDNA loci, the most diverging regions are the two ITSs, both showing 96% of identity, whereas 18S, 28S, and 5.8 show an identity of 99%, 99%, and 98%, respectively. This locus is flanked by two IGS regions significantly shorter than those previously described (lengths of ∼3600 and ∼1600 bp) (Supplemental Fig. S12A,B), with one characterized by the insertion of a nonrepetitive genomic region of 480 bp of unknown origin, followed by a solo-LTR (Supplemental Fig. S12C,D). On the other hand, 5S rDNA genes present lower sequence similarity levels that can drop down to 80.8% of identity between 5S_14 and 5S_22 paralogs, with a mean identity of 93.6% (Fig. 2E).

Fossil-calibrated and genome-scale divergence time reveals a late Cretaceous origin and Paleogene radiation of Tuberaceae

To provide a robust timeline for Tuberaceae emergence and diversification as well as a backbone for all comparative genomic analyses, we first produced a recalibrated divergence time among a set of 32 Ascomycetes, including eight Tuberaceae (Fig. 3). The species were selected to maximize the number of possible calibration points, allowing us to use the highest number of fossils with clear phylogenetic placement (10 fossil calibrations) to date for an Ascomycota divergence time estimation, assembling a phylogenomic supermatrix of 1057 complete and single-copy BUSCO genes.

Figure 3.

Late Cretaceous emergence of Tuberaceae and a Paleogene diversification of true truffles. Phylogenetic relationships and divergence time estimations obtained with 1057 complete and single-copy BUSCO genes. The species tree was constructed using an amino acid supermatrix, and both protein (black) and nucleotide (red) data were used for divergence time estimation. All nodes received maximum support values from the species tree inference, except for the split between Arthrobotrys and Pezizomycetes. The outgroup Schizosaccharomyces osmophilus was removed for visualization purposes. Numbers within circles represent the calibration points described in Supplemental Table S10. (Ne) Neoproterozoic, (Cam) Cambrian, (Or) Ordovician, (Si) Silurian, (De) Devonian, (Car) Carboniferous, (Pe) Permian, (Tr) Triassic, (Ju) Jurassic, (Cr) Cretaceous, and (Pa) Paleogene.

Open in new tabDownload PowerPointLink to figure

Generally, we did not observe over- or underestimation of divergence times of Ascomycota, and for most nodes, the use of nucleotide or protein data leads to confidence intervals largely overlapping (Fig. 3), providing evidence on the robustness of our estimates. Such estimates place the origin of the Tuberaceae crown group near the end of the Cretaceous (mean: ∼76 million years ago [MYA]), between previous estimations of ∼140 MYA (Bonito et al. 2013; Murat et al. 2018b) and ∼45 MYA (Miyauchi et al. 2020). The dawn of diversification within the genus Tuber is instead the Paleogene, ∼56 MYA, which sees the emergence of major subgeneric lineages: the Aestivum clade at ∼34 MYA and the Melanosporum clade at ∼26 MYA. The estimated split between T. panzhihuanense and T. borchii (Puberulum clade), its closest relative in the data set, is estimated to have occurred between ∼26 and ∼30 MYA.

High and variable TE content in Tuberaceae

We estimated the TE content of the 13 included Pezizales species using automatically de novo generated libraries. Although we did not manually curate these elements, as previously shown, the overall TE content and its composition did not change significantly when using the raw and partially curated repeat library in the Chinese white truffle genome (Supplemental Fig. S4B), making our results suitable for a genome-wide comparison of the TE content between different species (Carrasco-Valenzuela et al. 2025). First, we explored the presence of possible biases in assembly sizes (as a proxy for genome size) and TE estimations owing to different levels of assembly contiguities measured in terms of N50 values, and we did not identify any significant relationship (all P > 0.05) (Supplemental Table S4). Instead, we detected a significant positive correlation between assembly size (proxy for genome size) and both total TE content and Gypsy content, also when correcting for shared evolutionary histories using phylogenetic-independent contrasts (PICs) (Fig. 4A,B; Supplemental Tables S5, S6). Moreover, for total TE, DNA, and non-Gypsy LTR content, we found strong and significant phylogenetic signals, with Pagel's λ values approaching one, whereas no such relationship was observed for Gypsy elements (Supplemental Table S6). Indeed, Gypsy LTRs contribute very differently to the overall repeat content, even between closely related Tuber species such as T. melanosporum and T. indicum (Fig. 4C). The four identified clades account for most of the Gypsy content in Tuberaceae, with CHD-related elements being the main contributors (Fig. 4D; Supplemental Table S5), whereas the relative proportions of the different clades vary greatly among species (Fig. 4D). Lineage-specific activity of Gypsy clades in T. panzhihuanense was confirmed by repeat landscape profiles using a neutral substitution rate of 3.36 × 10⁻³ per million years (Fig. 4E). We did not produce repeat landscape profiles for other species owing to the significantly lower quality of their assemblies; however, an evolutionary scenario mainly involving lineage-specific amplifications was supported by Gypsy phylogenetic analyses with the presence of highly dense species-specific clades that are evolutionarily distant between each other (Fig. 4F).

Figure 4.

Lineage-specific proliferation of different Gypsy clades drives genome expansion in Tuberaceae. Correlations between total TE content (A) and Gypsy content with assembly size (B), used as a proxy for genome size after correcting for shared evolutionary history using phylogenetic-independent contrasts (PICs). (C) Assembly size and TE content for 13 Pezizales species used in comparative genomic analyses. (D) Relative contribution of identified clades to the total Gypsy content of Tuberaceae. Others refer to Gypsy clades that were not identified and characterized in T. panzhihuanense. (E) Repeat landscape profile describing the activity over time of Gypsy clades in the T. panzhihuanense genome after transforming the CpG-corrected Kimura distance into millions of years using a substitution rate of 3.36 × 10⁻³ per million year. The green box highlights the mean divergence time between T. panzhihuanense and T. borchii based on nucleotides (lower bound) and amino acids (upper bound). (F) Clade-specific phylogenetic trees of all Gypsy RT segments with a length of at least 100 amino acids extracted from all Tuberaceae genomes. Different colors highlight different species. All species abbreviations are reported in Figure 3.

Open in new tabDownload PowerPointLink to figure

Members of all four Gypsy clades were identified also within Morchellaceae, Gyromitra esculenta Pers. ex Fr., and Pyronema confluens Tul. & C. Tul. (Supplemental Fig. S13; Supplemental Table S5), suggesting their ancestral presence in Pezizomycetes.

Synteny conservation between true truffles and Morchellaceae despite high gene-family turnover rates and TE accumulation

Murat et al. (2018b) found a high rate of gene-family turnover within Tuberaceae together with an unexpected high level of microsynteny conservation within them. Here, we took advantage of the increased number of truffle genomes together with our high-quality T. panzhihuanense assembly to further test these findings and their relationship with TE accumulation. Because no gene annotation was publicly available for Verpa conica (O.F. Müll.) Sw. and G. esculenta, we performed a de novo annotation following the same pipeline used for T. panzhihuanense (for summary statistics, see Supplemental Table S7). We estimated, compared with other analyzed species, a fourfold higher gene-family (CAFE λ) turnover rate in Tuberaceae (0.163 vs. 0.036) (Fig. 5A) with their early evolution mainly impacted by gene losses. However, we also identified 27 significantly expanded gene families in their stem branch (CAFE Viterbi P < 0.05) (Fig. 5A).

Figure 5.

High gene-family turnover rates within Tuberaceae and synteny conservation with Morchellaceae. (A) Gene-family expansions and contractions along the Pezizales phylogeny. At each node, the top number indicates the number of gene-family expansions, and the bottom number indicates the number of contractions. The numbers in parentheses represent the significantly expanded and contracted gene families. (λ) Birth–death model rates fitted with CAFE for Tuberaceae and all other Pezizales species. Species abbreviations are shown in Figure 3. (B) Syntenic relationships between the T. panzhihuanense and M. sextelata genomes. The length of syntenic blocks is proportional to their actual length in base pairs. (C,D) Examples of syntenic genomic regions between T. panzhihuanense and M. sextelata involved in gene loss and small-scale rearrangements. (Private genes) Genes without any known orthologous relationship. Blocks have different scales and are shown at the bottom. (E) Cumulative distribution of the length of syntenic blocks in T. panzhihuanense and M. sextelata. (F,G) Comparison of TE and gene coverage (% of base pairs) across syntenic and nonsyntenic genomic regions in T. panzhihuanense, calculated over nonoverlapping genomic windows of 50 kb. Nonsyntenic genomic regions were found to be significantly more TE-rich and less gene dense (Wilcoxon rank-sum test, P < 0.01). (H) TE coverage (% of base pairs) across exons, introns, and 5 kb of gene-flanking regions for genes belonging to significantly expanded gene families identified in the branch leading to Tuberaceae diversification, as well as in the T. panzhihuanense terminal branch, compared with all other genes. Genes belonging to significantly expanded gene families were found to be significantly more TE-rich across all genomic compartments (Wilcoxon rank-sum test, P < 0.01).

Open in new tabDownload PowerPointLink to figure

Gene-based synteny analysis between T. panzhihuanense and the 16 longest scaffolds of T. magnatum confirms a notable conservation of mid-scale synteny within true truffles (Supplemental Fig. S14A). Moreover, two T. magnatum scaffolds align with two different contigs of T. panzhihuanense. Considering that all involved contigs possess telomeric repeats at both ends, these rearrangements might reflect chromosomal fusions or fissions. Also, the comparison of T. panzhihuanense with the long-read-based high-quality Morchella sextelata M. Kuo assembly revealed relatively high levels of mesosynteny conservation between the two species (Fig. 5B; Supplemental Fig. S14B). Indeed, 61% of the T. panzhihuanense genome and 70% of single-copy orthologs results were syntenic with synteny blocks that can reach up to 2 Mb of length (mean = 269.5 kb) (Supplemental Fig. S15A,B). Synteny blocks are markedly elongated in T. panzhihuanense relative to M. sextelata (Wilcoxon rank-sum test, P < 0.001) (Fig. 5E) exhibiting up to 20-fold length increases. TE accumulation in intergenic regions emerged as the main driver of synteny block and genome size expansion. Although M. sextelata shows slightly longer and more numerous introns than T. panzhihuanense (Wilcoxon rank-sum test, P < 0.001) (Supplemental Fig. S16A,B), intergenic regions are on average about 4.8 times longer in T. panzhihuanense (Wilcoxon rank-sum test, P < 0.001) (Supplemental Fig. S16C), with TE accumulation having a much greater impact on their length than on introns (Supplemental Fig. S16D,E). Syntenic genomic regions are slightly but significantly depleted in TEs (TE density: 0.59 vs. 0.65; Wilcoxon rank-sum test, P < 0.01) (Fig. 5F) and are almost two times more gene-dense compared with nonsyntenic ones (gene density: 0.159 vs. 0.08; Wilcoxon rank-sum test, P < 0.01) (Fig. 5G). However, 59% of the TE hotspots are completely contained within synteny blocks and flanked by syntenic genes (Fig. 5C,D). Visual inspection of some of these regions revealed some notable cases of loss of CAZy-encoding genes in the T. panzhihuanense lineage related to small scale rearrangements (e.g., GH37 and GH39) (Fig. 5C) and TE accumulation (e.g., PL1_2) (Fig. 5D).

Finally, we found that genes belonging to significantly expanded families in the branch leading to Tuberaceae diversification as well as in T. panzhihuanense terminal branch are significantly more TE-rich (Wilcoxon rank-sum test, all P < 0.01), in their exons, in their introns, and within 5 kb of flanking sequences (Fig. 5H). Additionally, genes included in TE hotspots are enriched in paralogous copies of these fast-evolving families (X² test, sample size = 8701, degree of freedom = 1, X² = 149.1907, P < 0.01), with 4.7 times more of these genes than expected by chance. These findings led us to hypothesize that genes located within TE hotspots are likely to be members of multicopy gene families. Indeed, 54% of the gene content of TE hotspots belongs to families with more than one gene compared with 18% of the genome.

Gene-family expansions at the root of Tuberaceae are related to the establishment of an ECM lifestyle

Gene duplication can lead to neofunctionalization and protein diversification (Copley 2020). We therefore tested whether gene duplication contributed to the transition from a saprotrophic to an ECM lifestyle in the Tuberaceae stem branch. Murat et al. (2018b) proposed that conserved ancient gene networks underlie the development and functioning of Tuberaceae ectomycorrhizae. Accordingly, we expect orthologous ECM-induced genes to have retained their function throughout the diversification of the clade. We therefore identified gene families associated with ECM (ECM-induced families) by intersecting our gene-family set with genes found to be upregulated in T. magnatum ECM compared with free-living mycelium by Murat et al. (2018b). In total, we identified 328 gene families containing T. magnatum ECM-induced orthologs, of which 18 are expanded, and five of these are significantly expanded in the branch of interest (Fig. 6A; Supplemental Table S8; Supplemental Fig. S15). We then applied a Fisher's exact test to assess the hypothesis that ECM-induced gene families are significantly overrepresented among expanded and significantly expanded gene families in the stem branch of Tuberaceae. Contingency tables were constructed by counting the number of ECM-induced and non-ECM-induced gene families in the expanded/significantly expanded categories versus all other gene families having at least one gene in one Tuberaceae species (i.e., present at the root of Tuberaceae). In both instances, we found a significant overrepresentation of ECM-induced gene families of 2.5-fold and 3.7-fold, respectively, among expanded (P < 0.001) and significantly expanded gene families (P = 0.013).

Figure 6.

Ancestral gene duplications of ECM-induced gene families. (A) CAFE results for three ECM-induced gene families that are significantly expanded on the stem branch of Tuberaceae. For each node, inferred and observed gene-family counts are shown on internal and terminal branches, respectively. Green indicates expansions, red indicates contractions, and significant changes are marked with an asterisk. Numbers in brackets indicate node IDs and refer to Supplemental Table S8. Species abbreviations are given in Figure 3. (B) InterProScan domain annotations of representative T. panzhihuanense proteins from each of the three gene families. For OG0000075, we used a reference protein from T. magnatum because no family members were detected in T. panzhihuanense. Symbols correspond to those in A. (TPR) Tetratricopeptide repeat. Supplemental Figure S17 shows the results for the two others significantly expanded ECM-induced gene families, OG0000009 and OG0000138, encoding ankyrin repeat–containing proteins and AAA ATPases, respectively.

Open in new tabDownload PowerPointLink to figure

Significantly expanded gene families of particular interest are the OG0000007, OG0000028, and OG0000075 (Fig. 6A,B). The former encodes for NOD-like receptors (NLRs), and a total of 187 genes are clustered in this family, with only one ortholog in P. confluens and one in V. conica. Morchellaceae have apparently lost this gene family, whereas Tuberaceae possess a variable but high number of paralogous copies. The other two gene families encode for putative Ras-like and But2 C-terminal-like proteins, respectively. For the former family, Tuberaceae species encodes between six and 30 genes, whereas other Pezizomycetes encode between three and zero. For the latter family, Tuberaceae species encode at least five gene copies, except for T. panzhihuanense, which appears to have lost this gene family. In other analyzed species, there are between one and three copies. All Tuberaceae But2 C-terminal-like genes encode for SSP. The two remaining significantly expanded ECM-induced gene families OG0000009 and OG0000138 encode, respectively, ankyrin repeat–containing proteins and AAA ATPases (Supplemental Fig. S17).

Sampling and DNA/RNA sequencing

The sample used for genome sequencing was collected in 2020 by the Jinsha river basin (river section of Zhongba village, Renhe district, city of Panzhihua, China). The fruiting body was isolated and taxonomically identified based on morphology and molecular characteristics and sent to Personalbio for whole-genome sequencing (WGS). WGS was performed both on a PacBio Sequel II obtaining 38+ Gb of HiFi reads and on an Illumina NovaSeq 6000 platform. RNA-seq data used for genome annotation were obtained from three replicates of a mix of 100 ECM individuals isolated from one P. massoniana–T. panzhihuanense mycorrhizal seedlings. Samples were sent to Novogene Technology for library construction and sequencing on an Illumina NovaSeq 6000 platform in 150 PE mode. For detailed sample identification, mycorrhizal synthesis preparation, and library construction, see the Supplemental Methods.

Genome survey and assembly

Illumina reads were quality checked with FastQC v0.11.7 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) and trimmed with Trimmomatic v0.39 (-phred33 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36) (Bolger et al. 2014). Trimmed reads were used with KAT v2.4.2 (Mapleson et al. 2017) hist to produce k-mer frequency histogram. Genome size and heterozygosity were estimated with GenomeScope 2.0 (Ranallo-Benavidez et al. 2020). PacBio subreads were processed to obtain final circular consensus sequencing (CCS) reads using PacBio CCS tool v6.4.0 (‐‐min-rq 0.99). The genome was assembled with Hifiasm v0.16.1 (Cheng et al. 2021) under default parameters, and Blobtools v1.1 (Laetsch and Blaxter 2017) was used to identify contigs deriving from fruiting body–associated microorganisms. To obtain coverage information, we mapped HiFi reads back to the assembly with minimap2 v2.24-r1122 (-ax map-hifi) (Li 2018), whereas for taxonomic annotation, we blasted all contigs against the NCBI nt database (BLASTN -max_target_seqs 10 -max_hsps 1 -evalue 1 × 10⁻²⁵). Reads mapped to identified microorganism-associated contigs were excluded, and the genome was reassembled with Hifiasm. TTAGGG telomeric repeats were identified with tidk (Brown et al. 2025) under default settings. The quality and completeness of the final assembly were assessed with BUSCO v5.2.2 (fungi_odb10 reference database) (Mosè et al. 2021) and Merqury (Rhie et al. 2020) with a k-mer size of 21. Metagenomic taxonomic profiling was performed using raw short reads from our study, along with publicly available sequencing data from the gleba of Tuber brumale, T. indicum, and T. melanosporum. A T. borchii isolate was used as a control. Paired-end sequencing reads were classified using Kraken2 v2.1.2 (Wood et al. 2019) with a custom Kraken database built from merged bacterial and fungal genomes downloaded from NCBI RefSeq. For details about Kraken2 classification, see the Supplemental Methods.

Repeat annotation

For repeat annotation, we produced a starting raw repeat library with RepeatModeler2 v2.0.4 (Flynn et al. 2020) with the LTR structure extension. Raw consensus sequences were used for a first genome annotation with RepeatMasker v4.1.2-p1 (Smit et al. 2013–2015) in sensitive mode (-s). Based on RepeatMasker results, we isolated all consensus sequences with at least 50 instances in the genome and/or characterized by protein fragments based on Blastx (E-value 1 × 10⁻⁵) results against the RepeatPeps database from the RepeatMasker package. Selected sequences were subjected to manual curation and classification based on structural features (presence and length of target site duplications, presence of terminal inverted repeats or long terminal repeats, and TE-related domains identified through the NCBI Conserved Domain Database) following a “blast–extend–extract” process as described by Goubert et al. (2022) and Peona et al. (2024). Curated and uncurated consensus were merged and redundancy removed following the 80–80 rule (i.e., requiring a minimum 80% identity along 80% of the shortest sequence) (Wicker et al. 2007; Goubert et al. 2022) with cd-hit-est (-c 0.8 -n 5 -aS 0.8 -g 1 -G 0 -t 1). This partially curated library was used for repeat annotation and soft-masking the genome with RepeatMasker in sensitive mode. We postprocessed RepeatMasker results with the parseRM.pl script (https://github.com/4ureliek/Parsing-RepeatMasker-Outputs/blob/master/parseRM.pl) to obtain genome-wide repeat estimations and RepeatCraft (Wong and Simakov 2019) in loose mode to defragment closely spaced repeat loci and obtain a final TE annotation.

Gene annotation

Quality of RNA-seq reads was assessed with FastQC v0.11.7 and trimmed with Trimmomatic v0.39 (-phred33 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36). We assembled the transcriptome with Trinity v2.1.1 (Grabherr et al. 2011) under default parameters.

For gene annotation, we used Funannotate v1.8.16 (https://doi.org/10.5281/zenodo.1134477) on the soft-masked version of the genome performing an ab initio gene prediction with AUGUSTUS (Stanke et al. 2008), SNAP (Korf 2004), GlimmerHMM (Majoros et al. 2004), CodingQuarry (Testa et al. 2015), and GeneMark-ES (Lomsadze et al. 2005). As external evidence, together with the assembled transcriptome, we selected a high-quality and consistent set of proteins composed by four Ascomycota RefSeq proteomes available on the NCBI GenBank database (GCF_000151645.1: T. melanosporum; GCF_003444635.1: Morchella importuna; GCF_020137385.1: M. sextelata; and GCF_024521635.1: Tricharina praecox) and the Swiss-Prot database (The UniProt Consortium 2023). We excluded other available Tuber spp. gene annotation available on NCBI as the gene annotation is not part of the RefSeq collection (https://www.ncbi.nlm.nih.gov/refseq/annotation_euk/process/).

Functional annotation of predicted genes was performed with InterProScan V5.70 (Jones et al. 2014) to obtain Pfam domains (Mistry et al. 2021) and InterPro entries. CAZymes were separately annotated with dbCAN3 (Zheng et al. 2023). rDNA genes were annotated using RNAmmer (Lagesen et al. 2007). SPs were annotated following the method of Pellegrin et al. (2015). We considered SP proteins as SSPs when their length was equal to or lower than 300 aa, as in the work by Pellegrin et al. (2015).

TE genomic distribution analyses

To test the interplay between gene density and TE densities, we subdivided the T. panzhihuanense nuclear genome into nonoverlapping windows of 50 kb using BEDTools makewindows (Quinlan and Hall 2010) and calculated the proportion of base pairs occupied by genes and transposons for each interval.

To detect genomic regions significantly enriched (TE hotspots) and depleted (TE coldspots) in TEs, we applied a binomial test on the number of base pairs occupied by transposons in genomic sliding windows of 20 kb with a window step size of 5 kb. For each window, we computed the TE coverage using BEDTools coverage and compared the observed number with the genome-wide average estimation computed across all windows. For each genomic interval, we separately tested for both significant greater and lower TE content compared with genome-wide estimation as alternative hypotheses, generating two sets of intervals representing genomic regions significantly enriched and depleted in TEs, respectively. P-values were corrected for multiple testing with the Benjamini–Hochberg procedure and a false-discovery rate (FDR) cut-off of 0.01. Because applying a sliding window approach may result in overlapping genomic intervals being annotated as both TE-rich (TE hotspots) and TE-depleted (TE coldspots), we used BEDTools subtract reciprocally on both interval sets to remove overlapping regions. Genes completely falling within TE-enriched or TE-depleted windows were extracted using BEDTools intersect (-f 1).

Identification, classification, and annotation of LTR Gypsy families

We reconstructed T. panzhihuanense Gypsy families from the evolutionary relationships of RT domains of single insertions. RT nucleotide segments >400 nt, and the corresponding protein sequences were retrieved via BLASTX (E-value ≤ 1 × 10⁻⁵) of all LTR insertions against RT domains from GypsyDB (Llorens et al. 2011) and confirmed via BLASTP (E-value ≤ 1 × 10⁻⁵) against the RepeatPep library, retaining only best hits to Gypsy elements.

For phylogenetics, confirmed RT protein sequences were clustered with CD-HIT (-c 0.8), aligned with MAFFT v7.520 (G-INS-i) (Katoh and Standley 2013), and trimmed with trimAl (-gt 0.8) (Capella-Gutiérrez et al. 2009). A maximum likelihood (ML) tree was inferred with IQ-TREE v2.2.2.6 (Minh et al. 2020), using ModelFinder (Kalyaanamoorthy et al. 2017) for model selection and 1000 ultrafast bootstrap replicates (Hoang et al. 2018) for assessing nodal support.

To complement the phylogenetic analyses, we built a network from all-to-all BLASTN (E-value ≤ 1 × 10⁻⁶) of RT nucleotide sequences; alignments <300 bp were excluded. Bit scores were used as edge weights to detect communities with greedy_modularity_communities function of the NetworkX Python package (Hagberg et al. 2008). Gypsy communities were mapped onto the RT phylogeny, and families were defined as well-supported (bootstrap ≥ 75) divergent clades, splitting communities when necessary.

To build consensus sequences representative of the identified families, we recovered all insertions from which the RT segments used in the phylogenetic analyses were derived and aligned them separately with MAFFT v7.520. From these alignments, we generated initial consensus sequences with CIAlign (Tumescheit et al. 2022), refined them using a “blast–extend–extract” process as previously described, and finally assessed their completeness with TE-aid (Goubert et al. 2022).

We classified curated Gypsy families relying on phylogenetic relationships with previously described elements. RT protein domains extracted from each consensus were added to a set of reference Gypsy elements, including known Gypsy clades described in GypsyDB and elements from Novikova et al. (2010), Riccioni et al. (2008), and Payen et al. (2016). To identify representative elements of the six Gypsy lineages described by Payen et al. (2016), we isolated the RT segments from all Gypsy sequences mined by the authors based on BLASTX alignment coordinates, as previously described. All mined Gypsy sequences were then clustered at 80% identity at the protein level using CD-HIT. All sequences were aligned with MAFFT G-INS-I, and gappy positions were removed with trimAl (-gt 0.8). We inferred a ML tree with IQ-TREE together with ModelFinder and 1000 UltraFastBootstrap replicates. Each Gypsy family was assigned to the closest known clade if the bootstrap value was 75 or more.

Classified Gypsy families were used in an additional RepeatMasker analysis in sensitive mode, and results were postprocessed to obtain their genome-wide estimations and repeat landscapes, as previously described. We used the parseRM.pl script to produce a repeat landscape describing the activity of transposons through absolute time using a Tuber-specific neutral substitution rate (see the subsection “Phylogenomics and divergence time estimation”).

To identify solo-LTRs arising from NAHR events, we blasted back (BLASTN, E-value 1 × 10⁻⁵) all insertions against their source consensus sequence and identified alignments that reciprocally cover both the insertion and one of the two LTR regions of the source consensus for at least 90% of their length. If the alignment reciprocally covers the insertion and the whole consensus sequence for 90% of their length, we considered the elements as full-length, and all other instances were considered as degenerated and/or fragmented elements.

Comparative genomics analyses