A total of 77 grape rootstocks were used as experimental martial, collected from the Zhengzhou Fruit Research Institute, CAAS, Zhejiang Academy of Agricultural Sciences, Shandong Agricultural University, Zhangjiagang Shenyuan Grape Technology Co. Ltd., and the Changli Fruit Institute, Hebei Academy of Agricultural and Forestry Sciences ( Tables S1-4 ), including 3 samples of Muscadinia and 74 samples belonging to 12 species and their hybrids of Euvitis, from North American species (V. riparia Michx., V. rupestris Scheele., V. berlandieri Planchon., V. labrusca Linn., V. champini Planchon), and Asian species (V. vinifera Linn., V. amurensis Rupr., V. dacidii Foëx., V. heyneana Roem. and Schult., V. pseudoreticulata W.T. Wang., V. adstricta Hance., V. heyneana subsp. Ficifolia). These resources have also been used to a certain extent in the breeding of grape rootstocks.

The genomic DNA was extracted with the traditional CTAB method from the grape leaves, and the DNA concentration and integrity were detected by agarose gel electrophoresis. The DNA (5 μg) was separated into ~500 bp fragments, and these fragments were used for library construction with the NEBNex DNA Library Prep Reagent Set for Illumina (BioLabs). The sequences were sequenced on an Illumina HiSeq PE150 platform at Beijing Novogene Bioinformation Technology Co. Ltd. The whole genome of each sample was re-sequenced at an average sequencing depth of ~10×.

Burrows-Wheeler Aligner (BWA) software (Version: 0.7.10-r789) was used for comparison with the reference genome, with the command ‘mem -t 4 -k -M’ (Li and Durbin, 2009). The single nucleotide polymorphism (SNP) calling was performed on a population scale with a Bayesian approach as implemented in the package SAMtools (Version: 1.3.1) and GATK (version 3.7) (Li et al., 2009). We retained high-quality SNPs, with the minor allele frequency ≥ 0.05, mapping quality ≥ 20, and missing rates < 0.20 for subsequent analysis after filtration. SNP and Indel annotations were performed using the package ANNOVAR (Version: 2015-12-14) (Wang et al., 2010), which includes filter-based annotation, gene-based annotation, region-based annotations, and other functionalities. According to the annotation, SNPs were categorized in exonic regions, intronic regions, upstream and downstream regions, splicing sites, and intergenic regions. SNPs in coding exons were further grouped into synonymous SNPs, non-synonymous SNPs, stop gain, and stop loss.

An individual-based neighbor-joining (NJ) tree was constructed to clarify the phylogenetic relationship based on the p-distance using TreeBeST (version 1.9.2) software (Vilella et al., 2009) and visualized with MEGA (version 6.0) (Tamura et al., 2013). Principal component analysis (PCA) was conducted using Genome-wide Complex Trait Analysis (GCTA, version: 1.25.3) software (Yang et al., 2011) to evaluate the genetic structure, with the significance level of the eigenvector determined using the Tracey-Widom test. The genetics structure was examined with an expectation-maximization algorithm as implemented using ADMIXTURE (Version: 1.3) (Alexander et al., 2009). We predefined the number of genetic clusters K from 2 to 8 with 10,000 iterations for each run to explore the convergence of individuals.

The pattern of linkage disequilibrium (LD) was compared using SNPs. The degree of linkage disequilibrium coefficient (r² ) between pairwise SNPs was calculated to estimate LD decay with PopldDecay software (Version: v3.31) (Zhang et al., 2019). ‘-n -dprime -minMAF 0.05’ was set in the program. The average r² was calculated in a 500-kb window and averaged in the whole genome, and LD decay figures were drawn with an R script. The IBD was calculated for all of the pairwise comparisons among the 77 rootstocks with PLINK (Purcell et al., 2007). Pairs of accessions were considered to be genetically identical if IBD > 95%, and the cut-off value for first-degree relatives (IBD value ≥ 0.45) was determined according to the distribution pattern of all of the pairwise IBD values. The network images based on IBD were constructed using Cytoscape (Version 3.6.0). Pairwise F_ST values (Weir and Cockerham, 1984) were computed in the same windows to measure the population differentiation between groups.

The phenotypic traits of 77 grape rootstock samples were investigated, and their germplasm resources of resistance to phylloxera, root knot nematode, salt tolerance, cold tolerance, drought tolerance, and waterlogging tolerance were measured. Combined with these phenotypic data and based on SNP markers, a mixed linear model (MLM) in GEMMA software was used to analyze the resistance traits, such as phylloxera, root knot nematodes, salt, cold, drought, and waterlogging, in the GWAS analysis. We screened the potential candidate SNPs based on the associated significance (P-value).

About 645 billion genome sequencing data were generated from the 77 grape rootstocks at an average depth of ~15.5×, from which a total of 4,269,582,820 clean reads and 3,943,015,251 mapped reads were acquired. The mapping rate of these raw reads to the V. vinifera reference genome (BioProject: PRJNA673645) was 92.19% ( Tables S1, 1 ), and the estimated error rate was 0.03% ( Tables S1, 2 ). The quality of the sequencing was high (Q20 ≥ 94.52, Q30 ≥ 88.27), and the GC distribution was normal. A total of 2,805,889 SNP sites were identified in this research ( Tables S1–4 ).

According to our findings, the 77 rootstock samples were spread over the phylogenetic tree’s ( Figure 1 ) ten main branches and belonged to five distinct populations: the Chinese and American group, V. riparia group, V. berlandieri group, V. rupestris group, and highly hybrid group ( Figure S1 ; Table S1-5 ). Most of the individuals in Population 1 (P1) are variants of V. riparia and their hybrids with other populations, such as R19, R2, R3 (V. berlandieri× V. riparia), R16, R17 (V. rupestris× V. riparia), R22, R8 (V. amurensis× V. riparia), R10, and R11 (V. riparia). This group also contains one variety of V. champini (R4) and one variety of V. rupestris (R15). P1 contains three groups (group1, group2, and group3) that are closely related to each other, indicating that V. champini, V. rupestris, and V. riparia have a close genetic relationship; group4 and R36 from group8 are the offspring of V. berlandieri× V. riparia and are classified as Population 2 (P2); Population 3 (P3) includes group5 and group6, which represent V. rupestris and its hybrid offspring with V. berlandieri; group5 belongs to V. rupestris; and group6 belongs to V. berlandieri× V. rupestris, and P4 represents a relatively complex group with relatively distant genetic relationships between V. labrusca (R56) and its hybrid offspring with other populations of R57, R45, R46, R47, and R48, V. champini (R52 and R58), the offspring of V. riparia (R50, R51, R53, and R54), and individuals of unknown parental origin (R55); and P5 has obvious regional characteristics and is divided into V. rotundifolia of American origin and the wild resources of China.

The evaluation of their above genetic similarity and clustering relationships through PCA indicated that the proportions of the total variance of PC1, PC2, and PC3 were 14.45%, 8.72%, and 4.87%, respectively, furthermore a total of 28.04% of all the genetic differences were observed. The first principal component (PC1) separated groups 2, 3, 4, 5, 6, 9, and 10, which had differences at the parental and geographic levels ( Figure 2 ). Group 5 was clustered in group 2 and clearly separated from other groups. Group 4 and group 6 clustered together due to the common ancestor i.e. V. berlandieri. The group1, group2, and group3 descendants of V. riparia were scattered and be separated by PC2. The dispersion between the individuals in group7 and group8 was more apparent, it shows that genetic differentiation existed in each group, with large genetic differences and rich genetic diversity. Group9 and group10 were situated far away from the other populations. The results of PCA were also consistent with the phylogenetic clusters.

Further pedigree analysis within grape rootstock groups based on identity-by-descent was also conducted. the majority of the first-degree relationships were among accessions in the same major Vitis categories ( Figure 3 ). Group2, group4, group5, and group6 are closely related to each other, forming a compact independent cluster individually, which is the result of the hybridization of V. rupestris, V. riparia, and V. berlandieri. Group3 and group9 have the common parent V. amurensis, and are weakly related, however, P5, with obvious regional characteristics, formed scattered clusters.

In plant-forward genetics studies, determining the linkage disequilibrium (LD, represented as r2) pattern is critical. The population size of the species, fertility (selfing or hybridization), selection pressure, and the rate of recombination all affect the level of LD, for which r² between the paired SNP sites were calculated using the available genome-wide SNPs. It was found that the attenuation of LD was 4.6 kb in P1, 6.1 kb in P2, 23.2 kb in P3, and 4 kb in P4, reaching half of the maximum r². These parameters are substantially greater than Liang et al. (2019) and Kui et al. (2020) identification of 2.9 kb of European wild grapes and 350 bp of cultivated grapes. The LD attenuation distance of V. amurensis from China and V. rotundifolia from the United States was 1.7 kb, which is smaller than that of wild European grapes. The r² was calculated uniformly for all of the rootstock species resources, and the results showed that an increase of 2 kb physical distance could reach half of the maximum r², and an increase of 5.1 kb attenuated to 0.1 ( Figure 4 ). The rate of LD decay differed in the different populations. The attenuation speeds of P4 and P5 were the same, whereas the remainder were in the order of P2 > P1 > P3. The LD of P4 and P5 had the lowest attenuation distance, which might be due to the significant genetic variation among these two groups. The V. rupestris population P3 had the slowest attenuation rate, the most domestication, and the most intense selection. The average coefficient size of LD at a distance of 50 kb in the genome was.41, but was reduced to 0.37 when reaching 100 kb.

Combined analysis with ADMIXTURE ( Figure S3 , S4 ) was carried out, and when K=2 was selected, group9 could be clearly distinguished from the other populations. Group4 and group6 are both hybrid descendants of V. berlandieri, and their composition was the same and could be distinguished when K=4. The composition of group2 and group5 was the same, which suggested that V. rupestris and V. riparia had a common ancestor, and their differentiation was observed until K=5. When K=3, group10 and group3 could be differentiated; when K=4, group8 was highly polymorphic and had a distant genetic distance from the other populations; when K=5, group7 was differentiated; and when K=6-8, more groups appeared in group7, 8, and 9.

The fixation index (F_ST ) can measure the degree of difference across populations by reflecting the allelic heterozygosity level of a population. The Fast of group4 and group6 was about 0.25, indicating that the genetic differentiation of these two groups was relatively large, whereas the other groups were greater than 0.25, indicating that there was strong genetic differentiation among these populations. As indicated in Figure 5 , the highest Fst value appeared between group2 and group5, indicating that there was relatively high genetic differentiation between V. rupestris and V. riparia, whereas the lowest level of genetic differentiation was observed in group4 and group6. Group4 was closer to group6 and relatively far from the group3, which was consistent with the results of the PCA and phylogenetic analysis. Overall, the Fst between the populations was greater than exceeded0.25, indicating that the high level of genetic diversity among the grape rootstocks. Counting of the overlap number indicated that the largest number of gene fragments was enriched in group3 vs. group4, distributed on 19 chromosomes. The Gene Ontology (GO) analysis terms were mainly enriched in biological process and molecular function, including nucleobase-containing compound, DNA replication, phosphatase activity, phosphoric ester hydrolase activity, spliceosome and galactose metabolism, and purine metabolism. The genes of V. riparia × V. berlandieri and V. rupestris × V. berlandieri were distributed on 19 chromosomes after selective elimination, especially chr18, and they were mainly involved in biological process and cellular components, and the pathways mainly included biosynthesis of amino acids and carbon metabolism. Group2 vs. group5 represents V. rupestris and its hybridization with V. riparia, and the hybridization resulted in differences in genes on 11 chromosomes, which mainly affected molecular function, as well as inositol phosphate metabolism, carbon metabolism, and biosynthesis of amino acids pathway. Group5 vs. group6 represents V. rupestris and its hybridization with V. berlandieri, which caused genetic differences on six chromosomes, mainly affecting the biological process, RNA transport, mRNA surveillance pathway, and nucleotide excision repair pathways.

To analyze the degree of differentiation between Chinese and other populations, we performed selection and removal analysis of group 9 and other groups, and the degree of differentiation was in the following order: group 8, group 3, group 7, group 4, group 6, group 1, group 5, and group 2 ( Figure 6 , S5 ). Most individuals in group9 belonged to V. amurensis, which had a relatively low degree of differentiation between group 3 and group 8 (which also contain components of V. amurensis), while a higher degree of differentiation was detected in the seedlings of V. riparia (group 2) and V. rupestris (group 5). The biological process category had the most significant difference in GO functions between group 9 and the other groups. Group 9 had more overlaps with group 1 and group 2, among which group 2 was relatively evenly distributed on 19 chromosomes ( Figure 6A ). It is worth mentioning that V. rotundifolia from the United States belonged to a different subgenus from the other populations and exhibited no genetic differentiation.

GWAS analysis were carried out through filtering and screening the genomic data to assess the resistance of grape rootstocks to phylloxera, root-knot nematodes, salt, drought, cold and waterlogging traits, and the locations of 631, 13, 9, 2, 810, and 44 SNP loci, related to 3312, 134, 75, 13, 5133, and 337 genes were identified ( Figure 7 ; Table S2 ). Among them, the genes related to resistance to phylloxera, waterlogging and cold were scattered on each chr. Drought-related genes were distributed on chr5 and chr10, with 9 and 4 genes mapped, respectively. Furthermore, the genes related to salt stress were distributed on chr14, chr03, chr10, chr13, chr5, chr9, and chr12. And most genes related to root-knot nematodes were distributed on chr13 and chr14, with 32 and 26 genes mapped respectively ( Figure 8 ; Table S2 ). We discovered that the association analysis result for salt was the closest to the predicted P-value of the non-associated null hypothesis using a quantile-quantile (Q-Q) plot, showing that the features were not produced by group stratification. The P-value of phylloxera and cold differed the most from the expected P-value, indicating that population stratification and individual kinship had a great influence on the association analysis and could easily result in false-positive associations, so the most SNP sites are located.

A total of 13 SNP loci were related to grape rootstock resistance to root-knot nematode in a GWAS analysis, with chr5 accounting for roughly 20% of the total. Furthermore, three were located on serine/threonine-protein kinases containing LRR structure, two PP2 proteins, one cysteine protease XCP1 gene, one leucine-rich protein kinase, one MYB48 transcription factor, one bHLH61 transcription factor, one BTB/POZ domain-containing protein, and one NRT1/PTR FAMILY 5 protein. These genes are associated to plant responses to biotic and abiotic stress, either directly or indirectly. Among the 75 genes associated with salt damage traits, multiple cysteine-rich receptor-like protein kinases are located on chromosome 8, and sodium/hydrogen exchanger 2 is located on chromosome 5. A total of 13 genes were associated with cold tolerance traits, of which nine genes are located on chr2, P4H1, PPR, PRN, SAD2, and Ku70, and four genes have unknown functions, meanwhile, four genes are located on chr10, CAT2, PLGG1, PLAT, and BTB/POZ-MATH. A total of 44 SNP loci were mapped to the waterlogging tolerance traits of rootstocks, among them chr3 and chr10 genes hosted the largest number of them. We mapped WRKY20, WRKY72, and WRKY32, two MADS-box genes including MADS-box 3 and MADS-box X1, three F-box genes involved in abiotic stress, one ERF gene and one EIN gene involved in the ethylene regulation, and six endo-1,3 (Reisch et al., 2012)-beta-glucanase genes also involved in plant resistance to disease.

The P-values of the GWAS analysis of the traits of grape resistance to phylloxera and cold differed from the ideal P-values. We located six WRKY transcription factors related to phylloxera and three genes related to root growth, including ROOT INITIATION DEFECTIVE 3, ROOT PRIMORDIUM DEFECTIVE 1, and SODIUM POTASSIUM ROOT DEFECTIVE 2, which may be related to phylloxera infestation during root development, thus endangering grape growth. C-repeat binding factor 1 was associated with cold-tolerance traits on chr11, which is the master switch of plant cold-resistant crop mechanisms. In addition, we identified 22 MYB transcription factors and found PIF4 and cold shock proteins (CSP2) on chr13. CSP can enhance the ability of cells to resist cold shock stress.

In addition, the association analysis revealed that several attributes shared a large number of the same genes. 19 genes were associated with three traits, of which six genes were also associated with the resistance to phylloxera, salt, and cold traits, including pyruvate kinase isozyme A, armadillo repeat-containing protein LFR, and Very-long-chain 3-ketoacyl-CoA synthase. Six genes were instantaneously and associated with the resistance to phylloxera, cold, and root-knot nematode traits, including three serine/threonine-protein kinases, Protein LURP-one-related 5, and ethanolamine-phosphate cytidylyltransferase. Seven genes relating to resistance to waterlogging, phylloxera, cold, and waterlogging traits were linked, including aluminum-activated malate transporter 12, protein ABCI7, and Ribulose bisphosphate carboxylase/oxygenase activase 1, and two protein phosphatase 2C genes.