The complete reference genome for grapevine ( Vitis vinifera L.) genetics and breeding

1Since the first human genome was published in 2000, hundreds of reference genomes have successively been assembled in a variety of species. A reference genome is essential for biological and genetic studies. Thus, acquiring a high-quality genome has persistently been pursued. Despite this, there are many missing segments due to highly repetitive sequences clustered across the genome, especially three representative regions: telomere, centromere, and ribosome DNA (rDNA).

2The centromere, which hosts CENPA/CENH3-variant nucleosomes and where the kinetochore forms and attaches to spindle microtubules, plays an essential role during cell division. It consists of alpha satellites, highly repetitive DNA sequences. The alpha satellite is composed of monomeric DNA repeats known as higher-order repeats (HORs), which contain arranged monomers that range from 100 to 200 bp. Despite their conserved function across species, their structure and sequence can change rapidly within and between species, and diverse organizations can be observed from one species to another. Nevertheless, centromeres show concerted evolution within genomes. Currently, the centromere remains mostly unknown to researchers.

3Telomeres are mostly unknown as well. They are composed of tandem repeats of relatively conserved microsatellite sequences located at the ends of chromosomes in eukaryotes. Telomeres are important for protecting chromosome terminal sequences during cell division. Ribosomal DNA (rDNA) is one of the most abundant repetitive elements in a genome, and plays an essential role in ribosome formation while driving cell growth and cell proliferation.

4Because of the missing information on previously assembled genomes, the investigation of centromeres, telomeres, and rDNA has been extremely limited in the past two decades. Fortunately, benefiting from the improvement of sequencing technologies and computational algorithms, genome assembly has ushered in a new era: that of telomere-to-telomere (T2T) sequencing. Compared with fragmented genomes, a T2T genome has fewer or no gaps at all. It is based on third-generation sequencing platforms, including PacBio high-fidelity long reads (HiFi), ultra-long Oxford Nanopore Technologies (ONT), and Hi-C data. Moreover, the T2T genome includes nearly complete information on the telomere, centromere, and rDNA regions. Promisingly, the T2T genome allows us to access these regions, opening a window into understanding the structure of these regions and the function of genes in these regions. Since the first complete human X chromosome was published in 2020, T2T assembly has quickly become a research hotspot. In plants, the first T2T genome was reported in Arabidopsis thaliana in 2021. At present, T2T genome assemblies have been obtained in several species, such as rice, banana, and watermelon, fascinating researchers into genomic structure and function and their relation to crop breeding traits.

5The grapevine (Vitis vinifera ssp. vinifera), a fruit tree that originated in the Near East, is one of the most widely cultivated and economically valuable crops worldwide. Domesticated grapes often have highly heterozygous genomes, which greatly impedes the acquisition of high-quality genomes. For instance, ~15% of genes are hemizygous in the ‘Chardonnay’ genome.31 Fortunately, the PN40024 genotype, a highly homozygous cultivar derived from selfing of cv. ‘Helfensteiner’, became the reference genome of grapevine, first obtained in 2007 (8X), and was the first fruit crop to be sequenced. Subsequently, several updated versions have been released: the 12X.v2 version and its upgraded annotation VCost.v3 in 2017, and the PN40024.v4.1 version in 2021. The grape gene reference catalogue now includes a full correspondence between all of their annotation versions. In addition, fragmented genome assemblies of various grape cultivars have been produced in recent years, such as those for ‘Black Corinth’, ‘Cabernet Franc’, ‘Cabernet Sauvignon’, ‘Carménère’, ‘Chardonnay’, ‘Merlot’, and ‘Nebbiolo’. As the grapevine is a representative dicotyledonous plant among fruit trees, its high-quality genome will greatly facilitate research on gene function, genetic structure, and evolution of Vitis and eudicot species.

6Despite the great number of grape genome sequences available, these genome assemblies are incomplete in repetitive regions, centromeres, and telomeres. Here we generated a T2T-level gap-free grape genome of the PN40024 reference and aimed to address four main analyses. The application of third-generation sequencing and assembly technologies to high-fidelity long reads has contributed to gap-free genome assemblies. Thus, our first question was to see whether we could complete the grape reference genome using these new sequencing and assembly approaches. Second, as studies on the centromere, telomere, and rDNA have long been neglected, we analyzed the features, structure, and distribution of these regions based on the assembled gapless grape genome. Third, the annotation of transposable elements (TEs) and genes in highly repetitive regions was improved based on the T2T genome, which could further improve our understanding of their biological functions, especially those of gene clusters. Finally, the PN40024 genome is almost fully homozygous, but some sites remain heterozygous after nine generations of selfing. It is worthwhile to investigate the genomic distribution and genetic effects of such heterozygous sites.

7PN40024, a highly homozygous inbred line originating from ‘Helfensteiner’, was used for T2T genome assembly. In total, 21 Gb (21 024 461 524 bp, ~42× coverage) HiFi reads were generated by the PacBio platform. For the preliminary assembly, hifiasm was used to assemble the HiFi reads. We then used MUMmer and the 12X.v0 genome version (V. vinifera genome assembly 12X.v0; https://www.ncbi.nlm.nih.gov/data-hub/genome/GCF_000003745.3) to order the 38 contigs into 19 chromosomes (Fig. 1). Only one gap was left after initial assembly into contigs (Supplementary Data Fig. S1). After filling the gap with continuous long reads of PN40024.v4, a gap-free PN_T2T genome was finally generated (494.87 Mb), being 69 Mb longer than 12X.v0 (426.18 Mb, Table 1) using the same statistical method. The k-mer metric was used to evaluate genomic homozygosity, estimated at 99.8% (Supplementary Data Fig. S2A–D). BUSCO (Benchmarking Universal Single-Copy Orthologs) was used to evaluate genomic completeness; 98.5% of the core conserved plant genes were found complete in the genome assembly (Supplementary Data Fig. S2E), which is 4.8% more than in 12X.v0 (93.7%, Table 1).

8Compared with the 12X.v0 genome, a substantial improvement of several metrics was observed in our PN_T2T assembly. The contig N50 length of PN_T2T was ~250 times higher than that of 12X.v0 (25.93 Mb versus 102 kb), and all 9429 gaps in 12X.v0 and 3391 gaps present in PN40024.v4 were filled in the PN_T2T genome (Table 1, Supplementary Data Table S1, Fig. 1A). As shown in Fig. 1C, 28 gaps in 12X.v0 were filled in PN_T2T, the largest gap being 16 951 bp in the 1-Mb syntenic region on chromosome 18 (Fig. 1C). Orientation errors in 12X.v0 were also corrected, such as inversions and translocations compared with PN_T2T (Fig. 1A, Supplementary Data Fig. S3). For example, two large inversions, which were located surrounding the centromere of chromosome 3 and at the ends of chromosome 5, with the length of 4.9 and 1.9 Mb, were observed between two versions of the assembly, respectively (Fig. 1A and B, Fig. S8). Moreover, 19 centromeres and 36 out of the 38 telomeres were detected on the PN_T2T genome assembly, except one telomere on chromosome 15 and one telomere on chromosome 17. A total of 37 534 genes and 41 064 transcripts were annotated, among which 24 526 (86.01%), 27 696 (78.83%), and 27 717 (78.75%) were shared with older versions PN40024.v2.1 (https://phytozome-next.jgi.doe.gov/info/Vvinifera_v2_1) and PN40024.v4.1 (https://grapedia.org/genomes/), respectively (Supplementary Data Table S2). A total of 5472 (14.58%) genes were not found to correspond in any of the three versions. A total of 97.9% of completely assembled genes was assessed by the BUSCO analysis, and structural domains were detected in 35 508 sequences out of 40 307 unique sequences (88.1%), while PN40024.v4.1 has 38 364 unique sequences, and 29 688 sequences were detected with structural domains (77.4%, Supplementary Data Table S2).

9Based on the species-specific pan-TE database constructed by RepeatModeler2, the repeats were detected with the pipeline shown in Fig. 2A. Finally, 66.47% of our gap-free grape genome was marked as repetitive sequences (Fig. 1D). As a comparison, 62.47% of the repetitive sequences were identified in the 12X.v0 genome using the same pipeline (Supplementary Data Table S3). Among the repeats predicted in the PN_T2T genome, the largest portion comprises TEs (63.90%), with a total length of 316 Mb (59.96% and 292 Mb in 12X.v0). The TEs mainly consisted of the long terminal repeat (LTR) type (47.54%), predominantly Gypsy (20.22%) and Copia (19.67%) elements. In total, we detected 276 rDNA sequences, representing 0.019% of the genome.

10To access the telomeric and centromeric regions in PN_T2T, we identified the telomeres and centromeres using the pipeline described in Fig. 2A. For telomeres, we checked the 150-kb sequences at both ends of each chromosome, and the length of the telomere repeat unit was set to range from 5 to 12 bp. Finally, the telomere repeat unit (TTTAGGG/CCCTAAA) was detected, which was the most abundant in the genome and carried by all chromosomes. The same telomere repeat unit was reported in grapes by Melters et al. [11] and Castro et al.. We further predicted the telomeres in 36 out of 38 telomeres in the PN_T2T genome, except the short arms of chromosome 15 and chromosome 17 (Figs 1A and 2B, Supplementary Data Table S4). Among them, the longest telomere (31 kb) was in the short arm of chromosome 8, with 4479 repeats, while the shortest telomere (1260 bp) was in the long arm of chromosome 7, with only 180 repeats.

11To detect centromeric regions, we scanned candidate repeats from 30 to 500 bp along the genome. Tandem Repeats Finder (TRF) found 470 different repeat units in the PN_T2T genome. The 107-bp repeats were the most abundant unit in the whole genome, which had 182 620.5 (copies ≥2) repetitions accounted for ~3.95% of the total genome sequence, followed by 321 bp (2.45%), 214 bp (1.94%), and 135 bp (1.05%) (Fig. 3A). Interestingly, we found the sequences of 214- and 321-bp repeat units consisted of two and three copies of the 107-bp repeat unit, respectively. The TE analyses also supported the centromeric feature of this 107-bp repetitive unit (Fig. 2). Thus, the centromeres were recognized mainly based on 107-bp repeats, and localized on all 19 chromosomes (Figs 1A and 2B, Supplementary Data Table S5). As shown in Fig. 3B, the total length of 107-bp repeats varied from 1.4 kb to 3.8 Mb, but the sequences of the 107-bp repeats were highly conserved among chromosomes (Fig. 3C). The 107-bp repeats were the most abundant in all chromosomes, except chromosomes 3, 14, and 18 (Fig. 3D–H, Supplementary Data Table S6). We found that the 187-bp was the main repeat unit in chromosome 14 and was scattered throughout the whole chromosome, and that 51-, 56-, 105-, and 107-bp repeat units were highly overlapped and enriched in the centromere, which showed a core region in the chromosome through IGV visualization (Supplementary Data Fig. S4). The centromeric repeat unit in chromosome 3 was the 135-bp repeat and its integer multiples (270 and 405 bp). For chromosome 18, 66 bp and its integer multiple 132 bp were the main repeat units (Supplementary Data Fig. S4).

12To locate the centromeric repeats, we further examined the relationship between TEs and centromeres. LTR retrotransposons or centromeric retrotransposons (CRs) were usually mixed with tandem repeats and enriched in plant centromeric regions. We found (Fig. 4A) that the genes and TE repeats, such as LTR (Gypsy and Copia), DNA TE (MULE-MuDR), and RC (Helitron), had a low density in the special region when the enormous centromeric tandem repeats enriched in the chromosome were viewed in Integrative Genomics Viewer (IGV) (Fig. 2, Supplementary Data Fig. S4). We then inferred the region with centromeric repeats and low TE density as the centromeres after zooming one by one (Supplementary Data Fig. S4, Supplementary Data Table S5). The pattern of 107 bp was the target, which was highly linked with the centromeric region in grapes. However, there were likely different repeat units and patterns that appeared on chromosomes 3, 14, and 18 (Fig. 3F–H). The scattering of transposons and the distribution of the centromere showed that specific sequence-defined repeat superfamilies were correlated or anticorrelated, to various levels, with centromeric proximity (Figs 2B and 4A), forming density gradients that are the main chromosome-scale repeat-associated features, presumably reflecting overall chromatin structure (Supplementary Data Fig. S4).

13To detect the captured genes, we then screened all genes in these regions in the highly linked centromeric region. Interestingly, we found 343 genes (Supplementary Data Tables S7 and S8) captured in the centromeres, which included 179 genes with Uni-Prot ID through BLASTP. Through GO (Gene Ontology) functional annotation, 12 genes were enriched in protein binding (molecular function, MF), such as VviAMP1 (Uni-Prot ID Q9M1S8), involved in ethylene, gibberellin, and abscisic acid signaling pathways. In addition, we found 10 genes enriched in the cellular component (CC) of the cytosol, mitochondrion and cytoplasm, including auxin transport protein VviBIG (Uni-Prot ID Q9SRU2), which influences general growth and development in plants; fumarate hydratase 1 VviFUM1 (Uni-Prot ID P93033), which catalyzes the active of mitochondrial Krebs cycle-associated enzyme; and 6-phosphogluconate dehydrogenase, decarboxylating 2 VviPGD2 (Uni-Prot ID Q9FWA3), which plays a key role in the development of the male gametophytes and the interaction between the pollen tube and the ovule. Moreover, RNA modification, protein autophosphorylation, DNA integration, DNA recombination, and photomorphogenesis appeared enriched while exploring biological process (BP) related terms (Fig. 4C).

14To infer the gene clusters in the grapevine genome, protein-to-protein alignments among the PN40024 protein-coding genes exposed a rich panoply of duplication structures in terms of genomic positions and functions. Prominent and complex tandem-like blocks of high-similarity genes could be seen via visualizations of all–versus–all alignments (Supplementary Data Fig. S5). We found a total of 377 gene clusters in the grapevine reference genome (Supplementary Data Table S9). These duplications often involved local rearrangements and could extend to megabases with dozens to hundreds of genes involved (Fig. 5). On chromosome 16 (23–27 Mb), there were 599 enriched-domain genes mainly including WAKs (Wall associated receptor kinase galacturonan binding), PPR repeat, Leucine-rich, ABC transporter, Integrase domain, Peptidase family, Protein kinase and Reverse transcriptase (Fig. 5A). On chromosome 18 (25–36 Mb), there were 1 237 genes enriched in domains mainly including Integrase domain, C-JID domain, NB-ARC domain, Leucine rich repeat, Multicopper oxidase, Reverse transcriptase, Terpene synthase, and TIR. Our results show that many of the strongly enriched structural domains are part of the structural domains of plant disease resistance genes (R genes), including NB-ARC, TIR, and structures identified by the Colis database. We analyzed the domain architecture of our 41 064 PN_T2T PCGs and identified 3 381 possible R genes. Collectively, these R genes and gene clusters in grapes highlight a tremendous opportunity for exploring plant defense mechanisms.

15Based on the PN_T2T genome assembly, the resequencing data of four PN40024 clones were downloaded from NCBI and analyzed. A total of 244 215 SNPs were detected, among which 208 330 SNPs (85.3%) were shared in all four samples while the other 35 886 SNPs were only present in one to three samples (Fig. 6A). Interestingly, we found nine hotspots of heterozygous SNPs on chromosomes 1, 2, 3, 4, 7, 10, 11, and 16 (Fig. 5A, Supplementary Data Fig. S6). To further investigate the highly heterozygous regions, we examined the top 5% heterozygosity windows and identified a total of nine large continuous fragments (chromosome 1, 1.1–1.3 Mb; chromosome 2, 4.2–7.2 Mb; chromosome 3, 9.4–9.9 Mb; chromosome 4, 21.8–22.9 Mb; chromosome 7, 15.3–26.2 Mb; chromosome 10, 0.7–6.5 Mb, 17.6–18.3 Mb; chromosome 11, 7.1–7.8 Mb; chromosome 16, 13.0–13.5 Mb). The GO enrichment analysis of the genes in these regions showed that the most significantly enriched terms were response to water deprivation, protein phosphorylation, cell division, response to oxidative stress, and response to salt stress, which were closely associated with key physiological activities in plants (Supplementary Data Tables S10 and S11, Fig. 6C, Supplementary Data Fig. S7). We further phased these nine hotspots of heterozygous regions on the PN_T2T reference genome (Supplementary Data File 2).

16A complete reference genome is essential for crop genetic studies and breeding purposes. The latest version of the PN40024.v4 assembly improved the reference resource by including long-read sequences and by gathering a gold-standard annotation. Nevertheless, these previous versions still possessed thousands of gaps and lacked repetitive regions, centromeres, and telomeres, all of which limited access to variants within these regions. On occasion such unreachable regions underlie quantitative trait loci (QTLs) for important agronomic traits, such as berry color and sex determination on chromosome 2 and disease resistance on chromosome 14. A full reference genome has therefore great potential to reveal the missing heritabilities of important polygenic agronomic traits, increasing genetic gain in grapevine breeding.

17More and more investigations suggest the important functions of gene clusters, with a total of 377 gene clusters being detected in PN_T2T. The grapevine genome is also widely used in studies of plant evolution and comparative genomics because of its important phylogenetic position in the evolution of eudicots. The T2T version could be widely used in plant evolutionary genomics, especially the repetitive sequences, centromeres, and telomeres. The T2T gap-free reference genome has incorporated gene annotations of previous versions with more accurate TE annotation (up to ~67% of the genome), which will be an important resource for grapevine functional genomics and breeding.

18The centromeric region ranges from kilobases to gigabases in length, including >90% tandem repeats. The centromere is among the last great unknowns in genomics, since it was inaccessible by previous sequencing technologies. Assemblies often collapse due to the highly repetitive nature of the centromeric region. We assembled and annotated centromeres for all 19 chromosomes of the grapevine genome (Fig. 1). Most of the chromosomes have a single centromere while others could have multiple centromeric regions—the so-called holocentromere. On chromosomes 16 and 18 we found tandem repeats in many regions, while on other chromosomes only a single peak was detected (Fig. 2B), suggesting that the structure of the centromeric region might be more complicated and requires further investigation.

19In the PN40024 grapevine reference genome there are three major repetitive patterns across the 19 chromosomes, suggesting different chromosomal evolutionary histories (Fig. 3D–H). On chromosomes 3, 14, and 18, we found 135-, 56-, and 66-bp tandem repeats, respectively (Supplementary Data Fig. S4), while on other chromosomes the major unit of tandem repeats was 107 bp (Figs 3D–H and 4B and D). The evolutionary histories of the centromeres of each grapevine chromosome are still an open question to be addressed with all Vitis genomes. Previous comparative genomic analyses suggested that the centromere is conservative among closely related species with a constant number of chromosomes. Transformation of centromeric structures occurs during chromosome division and fusion when the number of chromosomes changes throughout evolution. The muscadine grape (Vitis rotundifolia) has 20 chromosomes, with chromosomes 7 and 13 collinear with subgenus Vitis chromosome 7, which is associated with a chromosome fusion event. Only one centromeric region is left on chromosome 7 in our grapevine reference genome (Fig. 2B, Supplementary Data Fig. S4), suggesting one centromere was lost during the evolution of the genus Vitis.

20Centromeric architecture shaped the content within the genome, population genetic diversity within species, and genetic differentiation among species. Population genetic analyses have previously revealed that the genetic variants in the centromeric region are highly linked, with much lower genetic diversity compared with chromosome arms. The centromeres capture tens to thousands of genes that are highly linked to the centromeric tandem repeats. These genes, along with the centromeric region, are functional as supergenes. In total, we found 343 captured genes (Supplementary Data Table S7) in the centromeric region in the grapevine reference genome. Interestingly, the genes are mainly involved in the ethylene, gibberellin, and abscisic acid signaling pathways.

21The current plant used to build the grapevine reference genome originated from the ‘Helfensteiner’ cultivar selfed for nine generations, which resulted in a 99.8% homozygous genome (Supplementary Data Fig. S2A–D). The remaining heterozygous sites are still of interest as they could represent hotspots of required heterozygosity, with lethal consequences if found in the homozygous state. Thus, we collected Illumina resequencing reads for four clones of PN40024 maintained in different international laboratories. Interestingly, the heterozygous SNPs and structure variants (SVs) were enriched in specific regions when mapped to PN_T2T. In total, we found 208 330 heterozygous SNPs shared by the four samples, and 35 886 SNPs specific to one to three samples. The former is more likely the original variant of PN40024 after nine generations of selfing while the latter could be somatic variants generated during distribution and tissue culture in the different laboratories. Interestingly, we found that hotspots of common variants were enriched in central biological processes, including the oxidation–reduction process and protein phosphorylation. The hotspots on chromosome 2 also covered the sex-determination QTL region (Fig. 6), which complicated the mining of the sex-determination genes, because the candidate genes were not present in the old version of the reference genome. It has been reported that, during the clonal reproduction of fruit trees, such heterozygous deleterious variants accumulate in the genome. The clonal processes hide recessive deleterious variants, including small SNPs and indels and large structural variants, in a heterozygous state. Strong inbreeding depression has been commonly observed in clonal crops, including potato, cassava, citrus, and grapevine, since the strongly deleterious variants in these genomic regions have been exposed to lethal or strong recessive selection during selfing cycles. In grapevine breeding, inbreeding and outcrossing depression were commonly detected because the hidden heterozygous recessive deleterious variants that increased during clonal propagation were exposed during sexual reproduction.

22Altogether, and still acknowledging all previous sequencing efforts, our work represents the completion of a full T2T sequence of the grape reference genome. This assembly, together with the previous manually curated annotation, currently being transferred into PN_T2T, should represent the gold standard for the grapevine community. In line with this forecast, the T2T assembly and its updated annotation are available for download at the Grape Genomics Encyclopedia (GRAPEDIA; https://grapedia.org/), where it will be used along with different application program interfaces, including gene cards, transcriptomic data visualizations, and software for variation–gene expression–phenotype associations.

23PN40024 is a line that belongs to one of the near-homozygous lines originally derived from the ‘Helfensteiner’ cultivar by successive selfing steps, estimated to be close to 97% homozygosity as tested by SSR markers. We got this inbred material from INRAE under a Material Transfer Agreement (MTA) and transplanted it in the greenhouse belonging to AGIS (Agricultural Genomics Institute at Shenzhen, Chinese Academy of Agricultural Sciences, Shenzhen, China) for subsequent experiments. Young leaves and ovules from PN40024 were flash-frozen in liquid nitrogen. Genomic DNA and RNA were isolated using the DNeasy Plant Mini Kit (Qiagen) following the manufacturer’s instructions. For PacBio HiFi sequencing, two single-molecule real-time cells were sequenced on a PacBio Sequel II platform, and a total of 21 Gb of HiFi reads was generated using CCS (https://github.com/PacificBiosciences/ccs) with the default parameter for the sequenced accessions. For RNA-seq, 10 μg of poly(A) mRNA that isolated from total RNA was used for preparing Illumina RNA-seq libraries for each sample. These libraries were then sequenced using the Illumina HiSeq™ 2000 system in accordance with the manufacturer’s instructions.

24Initially, the PN40024 genome was assembled by incorporating PacBio single-molecule real-time long-read sequences. Reads generated by the PacBio Sequel II platform were self-corrected, trimmed, and assembled by hifiasm, using default parameters (https://github.com/chhylp123/hifiasm). The initial output of hifiasm (v.0.13) yielded the p_ctg draft assembly. Genome heterozygosity was estimated using a k-mer-based approach by GenomeScope 2.0; it was estimated to be close to 99.8% homozygosity (Supplementary Data Fig. S2A–D). Then, homology-based scaffolds were generated with MUMmer (v.4.0.0) ‘scaffold’, using the 12X.v0 reference genome (Supplementary Data Fig. S3). By applying MUMmer tools, we ordered and oriented the contig-level assemblies into 19 chromosomes, and joined the adjacent contigs to generate a scaffold with 100 N. Finally, we adjusted the assembly manually through aligning the genome sequencing data from the previous version of PN40024, which was mapped to the genome assembly by minimap2 (v.2.21) and visualized in IGV (v.2.12.3) software to observe whether the gap regions were supported by reads (Supplementary Data Fig. S1). Filling and closing of the gaps with the selected and assigned contigs were performed by mapping the 50-bp sequences around the gap to continuous long reads of PN40024.v4 and obtaining the gapless T2T PN40024 assembly for all 19 grape chromosomes. The assembly was inspected based on BUSCO completeness and the duplication score. For the phasing of highly heterozygous regions, minimap2 was used to align all reads to the PN_T2T assembly. The primary contigs assembled by hifiasm and ragtag were used to phase these contigs into two haplotypes.

25We have used a self-developed method for genome annotation. The putative genes were first searched for by using transcripts and Uni-Prot as evidence. A preliminary gene model was then built for the putative genes and further search was performed using AUGUSTUS (v.3.4.0). All the found putative genes fragments were then filtered, including genes involving duplicated regions, genes with coding sequence lengths shorter than 90 and genes not supported by any evidence. We attempted to complement missing genes and the complete genes were subjected to alternative splicing analyses. Finally, all the results were examined by hidden Markov models downloaded from the Pfam database to obtain the final gene models. Interproscan (v.5.56–89.0) was used for function annotation for our assembly, and Pfam (v.34.0) and Coils (v.2.2.1) were used for the identification of structural domains (https://github.com/unavailable-2374/Genome-Wide-Annotation-Pipeline).

26The primary repeat analysis is outlined in Fig. 2A and began with the construction of a pan-Vitis database of repeat families by RepeatModeler (open-2.0.3) and a series of scripts, which was then applied with RepeatMasker (open-4.1.2). For building this pan-Vitis repeat database we downloaded 17 Vitis genomes from NCBI, then used RepeatModeler2 to identify TE families. After that, we got 17 consensus fasta files of TE families and by removing the single-copy and failed annotations we aggregated these files. We used NCBI-BLAST+2.9.0 to remove some redundant sequences (−i 80%, −l 80%). Next, we got the final file of repeat identity, then used deepTE with the Plant model to classify the unclassified repeat elements. Finally, the repetitive sequence of the complete reference genome was annotated by RepeatMasker.

27To compare previous versions of the grapevine genome with PN_T2T, we aligned the genomes using minimap2 and indexed the alignment BAM file using SAMtools (minimap2 -ax asm5 -t 4 --eqx A.fa B.fa | samtools sort -O BAM - > A_B.bam, samtools index A_B.bam). Next, to detect structural variations between genomes, we needed to find synteny and structural rearrangements between the genomes. For this, we used SyRI (syri -c A_B.bam -r A.fa -q B.fa -F B --prefix A_B). Finally, Plotsr was used to generate the graph (plotsr --sr A_Bsyri.out --sr B_Csyri.out --sr C_Dsyri.out --genomes genomes.txt -o output_plot.pdf, https://github.com/schneebergerlab/plotsr). MUMmer (v.4.0.0) was used to compare the 12X.v0 genome with the reference genome PN_T2T using whole-genome alignments [70]. First, we aligned the two genome sequences using nucmer (nucmer --mum) and then filtered one-to-one alignments with a minimum alignment length of 10 000 bp (delta-filter -i 95 -l 10 000).

28SAMtools (v.1.7) was used to extract the sequence of chromosome 18 (25.0–26.0 Mb) in 12X.v0 and align the sequence in PN_T2T. The gap information was detected with a python script (getgaps.py) and finally we used LINKVIEW2 (https://github.com/YangJianshun/LINKVIEW2) to visualize the alignment results.

29The telomere repeat units were explored by using the TIDK (v.0.2.0) (https://github.com/tolkit/telomeric-identifier) with options tidk explore -f genome.fa --minimum 5 --maximum 12 -o tidk_explore -t 2 --log --dir telomere_find --extension TSV. Then the whole genome was searched using the following parameters: tidk search -f genome.fa -s TTTAGGG -o tidk_search --dir telomere_find. Finally, we completed the rapid statistics of telomeres based on the TIDK plot and used the R script to visualize the telomere peak.

30For centromere annotation, TRF (v.4.09) was used to finish tandem repeat annotation with the parameters trf genome.fa 2 7 7 80 10 50 500 -f -d -m, and then we merged the results of annotation by using TRF2GFF (https://github.com/Adamtaranto/TRF2GFF). To complete the data statistics and visualization, we used information extracted by using the awk command in the Linux system and analyzed the results in IGV (v.2.12.3). We used four softwares to show more details about the centromeric region: Iqtree (v. 2.1.4-beta) was used to achieve the phylogenetic tree (options: -m GTR + I + G -bb 1000 -bnni -alrt 1000); itol (v.6) was used to visualize the phylogenetic tree; GeneDoc (v.2.7.0) (https://github.com/karlnicholas/genedoc) was used to achieve multiple sequence alignment; and R script was used to plot the data statistics and typeset details.

31To detect the functions of the genes captured in the centromeric regions, we downloaded the protein sequence library of Swiss-Prot (2022/08/30, https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/) for a local blast. After this, we extracted all the protein sequences of PN_T2T blasted by diamond (v.2.0.15) (parameter: -k 1 -e 0.00001, https://github.com/python-diamond/Diamond). We further uploaded the Swiss-Prot ID to DAVID (https://david.ncifcrf.gov/tools.jsp) and completed GO enrichment and annotation. Finally, data visualization was completed by our R scripts.

32To define the clustered genes in the reference genome, protein sequences were extracted using gffread and then filtered by e-value <1e-5 and similarity >30% using BLASTP for all-versus-all alignments. The filtered alignment results were combined with functional annotations to filter out alignment results that did not share the same structural domains. Finally, we determined the presence of gene clusters by identifying three consecutive identical Pfam accessions (https://www.ebi.ac.uk/interpro/entry/pfam/#table), using such Pfam accessions as seeds, and going up and down 30 genes to find genes with the same Pfam accessions. In total, 377 gene clusters were found (Supplementary Data Table S8).

33Four resequencing samples were downloaded from the NCBI database (SRR6156373, SRR8835144, SRR8835157, SRR8835168) and mapped to the newly assembled PN_T2T genome for SNP calling. Quality-controlled reads were mapped to the genome using bwa (v.0.7.15) with the default parameters. SAMtools (v.1.4) and GATK (v.4.1.8) were used for sorting and indexing the bam file with no duplicates. The gvcf files were combined in GATK and were used to join calling SNPs across all samples. To obtain high-quality SNPs, we performed strict filtering of the SNP calls based on the following criteria: (i) SNPs with more than two alleles were removed in all samples in vcftools with parameters --min-alleles 2 --max-alleles 2; (ii) we removed the SNPs with quality scores (GQ) <30 (--minGQ 30) and missing rate 0 (--max-missing 1); (iii) SNPs with minor allele frequencies (MAFs) ≥.01 to remove the invariable sites.

34This work was supported by the National Natural Science Fund for Excellent Young Scientists Fund Program (Overseas) to Y.Z., the National Key Research and Development Program of China (grant 2019YFA0906200), the Agricultural Science and Technology Innovation Program (CAAS-ZDRW202101), the Shenzhen Science and Technology Program (grant KQTD2016113010482651), and the BMBF-funded de.NBI Cloud within the German Network for Bioinformatics Infrastructure (de.NBI). We thank Bianca Frommer, Marie Lahaye, David Navarro-Payá, Marcela K. Tello-Ruiz, and Kapeel Chougule for their help in analyzing the RNA-seq data and in running the gene annotation pipeline. This study is also based upon work from COST Action CA17111 INTEGRAPE and from the COST Innovators Grant GRAPEDIA (IG17111), supported by COST (European Cooperation in Science and Technology). JTM is supported by PID2021-128865NB-I00 and RYC-2017-23 645 grants from AEI (Spain).

35Y.Z. conceived and designed the project with H.X., Z.C., and C.R. The PN40024 sample was provided by Z.C. under an MTA signed with INRAE. X.S., W.L., X.X., and Z.M. performed the tissue culture of the sample in the greenhouse. X.S., X.W., H.X., N.W., F.Z., H.X., H.Z. and Y.W. performed the bioinformatic analyses. A.V., K.A., D.H., J.G., J.T.M., D.W., Z.L., X.L., and W.L. performed the gene annotation. Y.P., S.H., Z.L., W.L., X.W., Y.F., Y.W., H.W. and C.L. assisted in bioinformatics analyses. X.S., S.C., X.W., H.X., and Y.Z. wrote the manuscript with comments and input from all authors.

36All PacBio sequence data have been deposited in the NCBI Sequence Read Archive under project number PRJNA882193 and the National Genomics Data Center (NGDC) Genome Sequence Archive (GSA) (https://ngdc.cncb.ac.cn/gsa/), with BioProject number PRJCA012093. The assembly and annotation as well as the sequences of centromeres and heterozygous regions have been deposited in zenodo: https://zenodo.org/record/7751391#.ZBgVmcJBy3A. The assembly and its annotation will be also hosted in the GRAPEDIA portal (https://grapedia.org/).

37All the scripts and pipelines used in this study have been archived in GitHub: https://github.com/zhouyflab.

38The authors declare no conflict of interest.

39Supplementary data is available at Horticulture Research online.