Quercus variabilis samples were collected from an ancient tree (more than 400 years old) growing in Culai Mountain National Forest Park, Shandong Province, China.

Fresh leaves were collected and immediately frozen in liquid nitrogen for transport back to the lab. The genomic DNA was then isolated using a Plant DNeasy Mini kit (Qiagen China, Shanghai, China) according to the manufacturer’s instructions. The quality and quantity of the DNA were determined using agarose gel electrophoresis and with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States). A library of short-insert-size genomic DNA fragments of length 300–400 bp was constructed according to the manufacturer’s instructions, and was sequenced on a DNBSEQ platform (Beijing Genomics Institute, Shenzhen, China) for 150 bp pair end sequencing. For long-read sequencing, a 20 kb high-fidelity (HiFi) library was constructed following the manufacturer’s protocol¹ on the PacBio Sequel II platform (Pacific Biosciences of California, Inc.). To increase continuity of the genome, a Hi-C library was constructed and sequenced on the DNBSEQ platform (BGI, Shenzhen, China).

For RNA sequencing, fresh leaves, young twigs, fruits and seeds were sampled and immediately frozen in liquid nitrogen. RNA was extracted using TRIzol reagent (Invitrogen), the genomic DNA was eliminated using DNase and the samples were then mixed for RNA sequencing. We used agarose gel electrophoresis, a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific) and an Agilent Bioanalyzer 2,100 (Agilent Technologies, Santa Clara, CA, United States) to evaluate the quality of the RNA. High-quality RNA was then used to build a cDNA library following the manufacturer’s instructions. Paired-end sequencing was performed on the DNBSEQ platform (BG I, Shenzhen, China), generating 150-bp paired-end reads.

SOAPnuke V1.6.5 (Chen et al., 2018b) was employed to filter out PCR duplications, low-quality reads (≥ 10% of nucleotides with a quality score ≤ 20 or the proportion of N is greater than 1%) and adapter sequences, with the following parameters: -n 0.01-l 20-q 0.1-i -Q 2-G -M 2-A 0.5 –d. Next, the software Jellyfish v2.1.4 (Marçais and Kingsford, 2011) was used to count k-mers of length 17–31. GenomeScope 2.0 (Ranallo-Benavidez et al., 2020) was then applied for the estimation of genome size and other features such as heterozygosity and repetition rate. Smudgeplot (Ranallo-Benavidez et al., 2020) was used for the estimation of ploidy.

The PacBio SMRT-Analysis² was used as quality control to eliminate adaptors and low-quality short reads, producing a total of 51.7 G bases (~ 72 × coverage) of PacBio HiFi data. The initial assemblies were then performed in HiFi-asm v0.15.2 (Cheng et al., 2021). To acquire high-quality, chromosome-level assemblies, Hi-C reads were compared to the contigs assembled above using Juicer (Durand et al., 2016b). Unique mapped reads with map-quality scores > 40 were subsequently used for Hi-C association chromosome assembly using the 3D-DNA pipeline (Dudchenko et al., 2017). Scaffolds were then manually checked and refined with Juicebox (Durand et al., 2016a) and visualized in Hicplotter (Akdemir and Chin, 2015). A BUSCO analysis was conducted to determine gene/genome completeness using BUSCO v4 (Simão et al., 2015) together with the embryophyta_odb10 database with 1,614 plant single-copy orthologues. BWA (Li, 2013) was used to map short reads of DNBSEQ data against the assembly. SAMtools (Danecek et al., 2021) was then employed to create a pile-up file summary of the aligned reads, and the results were imported to BCFtools (Danecek et al., 2021) for SNP and INDEL calling. The heterozygosity was then calculated as the proportion of the heterozygous sites to the total sites.

Several different methods were employed to annotate the repetitive sequences. First, Tandem Repeats Finder v4.09 (Price et al., 2005) was used for the identification of tandem repeats. Then, RepeatProteinMask v4.07 and RepeatMasker v4.07 (Chen, 2004) were used with their default parameters against RepBase v21.12 (Bao et al., 2015) to identify known repeats in a homology-based approach. Thirdly, RepeatMasker (Bedell et al., 2000) identified repeat elements with a de novo library, built in RepeatModeler (Abrusán et al., 2009) and LTR_FINDER v1.06 (Zhao and Hao, 2007).

To annotate the protein-coding genes, we combined RNA-based, homology-based and de novo methods. For the RNA-based method, we generated 63.27 million raw reads (9.49 Gb) with DNBSEQ sequencing (Supplementary Table S1). After quality control and filtering by fastp (Chen et al., 2018a), 9.48 Gb clean data were retained and aligned to scaffolds using hisat2 v2.2.1 (Kim et al., 2019). Reference genome-guided transcriptome assemblies were then constructed with StringTie v2.2.0 (Pertea et al., 2015). For homology-based predictions, the Q. variabilis genome was aligned against the Arabidopsis thaliana, Q. lobata, Q. robur and Q. suber genomes using TBLASTN v2.2.26 (Mount, 2007) with an E-value cutoff of 1e–5. Finally, GeneWise v2.4.1 (Birney et al., 2004) was employed for structural inspection of these alignments. For ab initio gene prediction, MAKER v3.01.03 (Holt and Yandell, 2011) was used to compute annotation edit distance (AED) for each protein-coding gene, based on transcript assembly from the transcriptome data, as well as from homologous annotations of the four genomes. Augustus v3.4.0 (Stanke et al., 2008; Keller et al., 2011) and SNAP (Johnson et al., 2008) were then employed for ab initio gene prediction using model training, based on coding sequences of 1,200 genes with structural integrity selected based on AED. Finally, the predictions obtained using these methods were combined using EVM v1.1.1 (Haas et al., 2008). The predicted genes were functionally annotated using seven public biological databases: NR, TrEMBL (Boeckmann et al., 2003), SwissProt (Boeckmann et al., 2003), KEGG (Kanehisa and Goto, 2000), InterPro (Zdobnov and Apweiler, 2001), KOG (Koonin et al., 2004), and GO (Consortium, 2004). Blast v2.2.26 was used for homolog searches with an E-value cutoff of 1e-5, and InterproScan v5.55 (Jones et al., 2014) was used for protein function prediction based on the conserved protein domains.

Homology-based non-coding RNA (ncRNA) was identified using Infernal v1.14 (Nawrocki and Eddy, 2013) by mapping plant small nuclear RNA (snRNA) and microRNA (miRNA) genes from the Rfam database (Kalvari et al., 2018). Transfer RNAs (tRNAs) were detected with tRNAscan-SE v1.3.1 (Lowe and Chan, 2016). BLASTN was used for the identification of ribosomal RNAs (rRNAs) by alignment with known plant rRNA sequences (Vitales et al., 2017).

OrthoFinder v2.5.41 (Emms and Kelly, 2019) was used to identify homologous gene families among the assembled genomes of Q. variabilis and 13 further representative flowering plant species (Amborella trichopoda, Arabidopsis thaliana, Castanea crenata, Castanena mollissima, Eucalyptus grandis, Juglans regia, Oryza sativa, Prunus persica, Q. lobata, Q. robur, Q. suber, Vitis vinifera and Xanthoceras sorbifolia). GO enrichment analysis was conducted using ClusterProfiler with an adjusted p value cutoff of 0.05 (Wu et al., 2021). We performed collinearity analysis of homologous gene pairs between Q. lobata, Q. robur, and Q. variabilis using MCScanX (Wang et al., 2012).

For phylogenetic analysis and estimation of species divergence time, MUSCLE (Edgar, 2004) was applied to align the amino acid sequences of single-copy orthologous genes. The concatenated amino acid sequences were further used for construction of the phylogenetic tree in IQ-TREE2 (Minh et al., 2020). MCMCTREE of PAML (Yang, 2007) was used to estimate phylogenetic dating using a BRMC method (Sanderson, 2003) with the soft fossil calibrations obtained from the TimeTree website³: split of A. trichopoda from O. sativa, 173–199 million years ago (MYA); split of X. sorbifolia from A. thaliana, 96–104 MYA; split of A. thaliana from J. regia, 107–135 MYA; split of V. vinifera from A. thaliana, 89–113 MYA. Gene families were filtered out if more than 200 genes were present in one species but only 2 or fewer in the other species. The remaining gene families were used to infer the expansions and contractions of protein family in CAFÉ v3.0 (Han et al., 2013).

Searches for putative paralogous genes were conducted for Q. variabilis and P. persica against each other using BLASTP (E-value ≤ 1e–5). Syntenic blocks were then identified using MCScanX (Wang et al., 2012) with parameters of –a –e 1e–5–s 5. Synonymous substitutions per synonymous site (Ks) values were calculated with codeml in the PAML package (Yang, 2007). For interspecific orthologues, the protein sequences of the homologous genes in Q. variabilis, Q. robur, and Q. lobata were aligned in BLASTP (E-value ≤1e-5), and the results were sorted according to their bit-scores and E-values to obtain reciprocal optimal gene pairs. Then codeml was used to calculate the Ks values of reciprocal optimal gene pairs. Finally, the Ks distributions of intraspecific paralogs and interspecific orthologues were evaluated to infer whole genome duplication (WGD) events and divergence time in the species genome.

We sequenced the Q. variabilis genome using a combination of PacBio and Hi-C technologies, and obtained a high-quality diploid reference genome (Smudgeplot, Supplementary Figure S1). A 20 kb DNA library was constructed and sequenced on a PacBio Sequel II platform, generating 51.70 Gb HiFi reads, approximately 72 × the estimated genome size (713.93 Mb; Supplementary Figure S2; Supplementary Table S1). Then, initial genome sequences spanning 796.30 Mb (327 contigs, N50 of 26.04 Mb; Supplementary Table S2) were constructed, slightly larger than the total genome size as estimated at 713.93 Mb using the 21-mer peak and distribution from DNBSEQ data (Supplementary Figure S2; Supplementary Table S1). This is perhaps due to chimerism caused by the relatively high heterozygosity (estimated to be 2.15%; Supplementary Figure S2). The contig N50 of Q. variabilis is significantly higher than that of other published congeneric species, e.g., Q. acutissima (1.44 Mb; Fu et al., 2022), Q. mongolica (2.64 Mb; Ai et al., 2020), Q. robur (0.07 Mb; Plomion et al., 2018), Q. lobata (1.9 Mb; Sork et al., 2022) and Q. suber (0.08 Mb; Ramos et al., 2018; Table 1). We next used 3D-DNA derived from the Hi-C data (Supplementary Table S1) to generate 12 pseudo-chromosomes (787.15 Mb, Supplementary Table S2), with lengths ranging from 39.05 to 97.21 Mb (Figure 1A; Supplementary Table S3). Interestingly, the number of pseudo-chromosomes of the assembled haploid is the same as that of other Quercus genomes (Q. acutissima, Q. mongolica, Q. robur, Q. lobata, and Q. suber). The chromosomal genome of Q. variabilis was characterized by 245 scaffolds, with a scaffold N50 of 64.86 Mb which is similar to that of Q. mongolica (66.7 Mb), slightly smaller than that of Q. lobata (75 Mb), but ~ 22-fold, ~ 50-fold and ~ 130-fold larger than those of Q. acutissima, Q. robur, and Q. suber, respectively (Table 1). We calculated the heterozygosity based on the 10,014,769 heterozygous sites (including SNPs and INDELs), and found that the heterozygosity of this genome was 1.26%, which was slightly lower than estimated (2.15%) due to underestimation (considering only SNPs and INDELs). We further evaluated the completeness of the genome assembly using the BUSCO.v4 plant datasets, and identified 1,587 (98.3%) of the 1,614 plant single-copy orthologues, with 1,526 (94.5%) presented as single-copy (Supplementary Table S4), a value superior to that of Q. lobata (95%), Q. robur (91%), Q. suber (95%), Q. acutissima (91%) and Q. mongolica (92.71%), indicating that our genome assembly is of high quality and nearly complete.

The total length of the repetitive sequences in the Q. variabilis genome was 538.34 Mb, covering 67.6% of the assembled genome (Supplementary Table S5). This proportion was higher than that observed in the Q. mongolica genome (435.34 Mb, ~53.75% of the genome) identified using the same process (Ai et al., 2020), and also higher than those in Q. lobata (54%; Sork et al., 2022), Q. suber (51%; Ramos et al., 2018), and E. grandis (55%; Myburg et al., 2014), which have been calculated using other processes. TEs accounted for 61.09% of the Q. variabilis genome (Supplementary Table S6). Long terminal repeat retrotransposons (LTR-RT), which often contribute to variations in genome size (Feschotte et al., 2002; Du et al., 2010), were identified as being the most abundant repeats (46.63%), followed by long interspersed nuclear elements (LINE; 8.14%) and DNA elements (17.88%; Figure 1B; Supplementary Figures S3–S13; Supplementary Table S6).

Gene models for the Q. variabilis genome were obtained using a comprehensive approach including ab initio, RNA sequence-based and homology-based predictions (Supplementary Table S7). In total, we predicted 32,466 protein-coding genes with an average gene length of 5,272.04 bp, an average coding-sequence length of 1,139.49 bp, an average exon length of 226.50 bp and an average exon number per gene of 5.03 (Supplementary Table S7). Functional annotation using the NR, Swissprot, KEGG, KOG, TrEMBL, Interpro and GO databases allowed 30,878 (95.11%) of the total 32,466 genes to be assigned putative functions (Supplementary Table S8). Of these, 52.32% (16,985) of the total genes could be functionally annotated through NR, InterPro, KEGG, SwissProt and KOG simultaneously (Figure 2A). We also predicted 12,220 rRNA, 942 tRNA, 157 miRNA, and 1,148 snRNA genes in the Q. variabilis genome (Supplementary Table S9).

Orthologous gene families were identified using the proteomes of Q. variabilis predicted in our project and those of 13 other flowering plant species, including the three congeneric species (Q. lobata, Q. robur and Q. suber; Supplementary Table S10). In total, the 29,808 Q. variabilis genes (91.81% of the total) clustered into 14,930 gene families, of which 6,245 gene families (including 11,113 Q. variabilis genes) were shared among all the 14 plant species. We also found that 964/4,854, 747/3,644, 1,821/12,372, 37/95, 891/6,186, 206/1,814, 2,486/12,784, 438/1,633, 340/891, 129/370, 4,311/11,160, 326/2,693, 814/3,015, 286/938 gene families/genes appeared to be unique to A. trichopoda, A. thaliana, C. crenata, C. mollissima, E. grandis, J. regia, O. sativa, P. persica, Q. lobata, Q. robur, Q. suber, Q. variabilis, V. vinifera, and X. sorbifolia, respectively. The gene families unique to Q. variabilis were mainly enriched in “glycine catabolic process,” “serine family amino acid catabolic process,” “organic acid catabolic process,” “oxaloacetate metabolic process,” “tricarboxylic acid cycle” and “nuclear chromosome segregation” (Figure 3; Supplementary File S1).

Phylogenetic analysis was conducted based on the 483 single-copy gene families derived from Q. variabilis and 13 further flowering plant species (Supplementary Table S10). We found that within this subclade, Q. variabilis is more closely related to Q. suber than to Q. lobata or Q. robur (Figure 2B). The divergence between Q. variabilis and Q. suber occurred at approximately 13.7 (7.1–21.3) MYA, while Q. lobata and Q. robur, which belong to a different subclade, diverged from the Q. variabilis-Q. suber subclade ~27.6 (16.4–40.0) MYA (Figure 2B).

To investigate potential WGD events in the evolutionary history of Q. variabilis, we studied the distribution of the Ks between homologous gene pairs derived from Q. variabilis, Q. lobata, Q. robur, Q. suber, and P. persica. One peak was found based on the paralogous gene pairs in Q. variabilis and P. persica (~ 1.5 Ks units), indicating a shared ancient WGD event (γ) for these two species (Murat et al., 2015; Figure 2C). The divergence between Q. variabilis and three congeneric species (Q. lobata, Q. robur and Q. suber; 0.02–0.05 Ks units) occurred later than the WGD event (Figure 2C). Further investigation of the genomic collinearity between Q. variabilis and Q. lobata, Q. robur showed a clear one-to-one syntenic relationship, and the overall gene synteny was largely conserved (Figure 4), suggesting that no large amounts of chromosome fusion or species-specific WGD events occurred after species divergence (Ai et al., 2020).

We analyzed gene family expansion and contraction based on the gene families in the 14 studied flowering plant genomes (Supplementary Table S10) using OrthoFinder. The number of expanded/contracted gene families in Q. variabilis compared with its common ancestor were 739/2,409, while in Q. suber, which is genomically the most similar to Q. variabilis, these numbers were 2,177/787 (Figure 2B). A significant number of expanded genes in Q. variabilis were enriched in “monoterpene metabolic process,” “terpene biosynthetic process,” “intrachromosomal DNA recombination,” “hydrocarbon biosynthetic process” and “oxaloacetate metabolic process” (Supplementary Figure S14; Supplementary File S2), while contracted gene families were enriched in “glutathione metabolic process,” “isoflavonoid biosynthetic process,” “toxin catabolic process,” “programmed cell death induced by symbiont” and “regulation of response to red or far red light” (Supplementary Figure S15; Supplementary File S3).