We used the backcross 1 (BC1) segregating population of 135 individuals (Table S1) from the cross between the hybrid VRH8771 [(‘Cabernet-Sauvignon’ x ‘Alicante Bouschet’) x M. rotundifolia] and V. vinifera cv. ‘Cabernet-Sauvignon’ (VRH8771 x CS). This cross has been initiated by Alain Bouquet at INRA Montpellier since 2005 and continued at INRA UMR EGFV (Bordeaux, France) after 2008. VRH8771, the female parent, is resistant to D. vitifoliae and X. index whereas V. vinifera ‘Cabernet-Sauvignon’ (CS), the male parent, is susceptible (Bouquet et al. 2000). Parental genotypes and BC1 individuals are conserved at the INRA germplasm repository at INRA UMR EGFV (Bordeaux, France).
Insect and nematode material
For D. vitifoliae, all experiments were conducted with the isofemale clone ‘Pcf7’. The original population was sampled in 2010 at Pineuilh (France) on V. vinifera cv. ‘Cabernet franc’ scions grafted on SO4 rootstock (V. berlandieri x V. riparia) and maintained at INRA UMR SAVE (Bordeaux, France) on leaves of the American variety ‘Harmony’, a complex hybrid between Dog-Ridge (V. champinii) and ‘1613C’ (V. labrusca x V. riparia x V. vinifera) and on root pieces of V. vinifera cv. ‘Cabernet Sauvignon’.
For X. index, all experiments were conducted with the isofemale line ‘Frejus’. The original population has been sampled in a GFLV-infected grapevine field in Frejus (Provence, France). Using the population grown on grapevine in the greenhouse, the line had been created from a single female inoculated on a fig plant previously grown from in vitro.
Experimental designs and phenotyping
An in planta assay was performed with 89 BC1 individuals. Two control genotypes were also tested: V. vinifera cv. ‘Pinot Noir’ and the rootstock ‘Börner’ (V. riparia x V. cinerea) as susceptible and resistant genotypes respectively. This experiment was organized according to a randomized complete block design with three experimental blocks, each block being an independent randomization of one replicate of each BC1 individual. Plants were grown in individual pots of 1L in a soil substrate composed of at least 50 % clay with pebbles at the bottom to improve drainage. The soil was sterilized by autoclaving in order to prevent cross contamination with other phylloxera population sources. Each pot was covered with an insect-proof transparent plastic bell and an automatic watering system was adapted (Fig. S6 A). One hundred phylloxera eggs of ‘Pcf7’ clone, previously grown on root pieces of V. vinifera cv. ‘Cabernet Sauvignon’ for two generations, were deposited on a moistened and sterilized filter paper near the root (~3 cm depth) of each pot. Three months after inoculation, the plants were uprooted and the nodosity number was counted (Fig. S6 B). The maximal number of nodosities scored among the three replicates of each BC1 individual was considered as the most reliable indicator of the quantitative resistance phenotype.
An in vitro assay was also performed according to Pouget (1975) to assess the larval development of phylloxera on 37 BC1 individuals. Five woody root pieces per individual, 6 to 7 cm long, were arranged in small bundles (a contact between all the roots is required) on a disk of dampened blotting paper in a Petri dish sealed with parafilm. Three replicates were realized for each individual. A fungicide treatment (Ridomil Gold, Syngenta ® - 2.3g/L) was carried out on the roots to prevent Botrytis cinerea infection. In each Petri dish, 50 phylloxera eggs were deposited on the roots pieces. Then Petri dishes were incubated at 25°C in the dark. The number of larvae that have developed on the five root pieces was counted one month after inoculation. The maximal number of larvae scored among the three replicates was considered as the best indicator of the quantitative resistance phenotype.
X. index assays
A total of 60 BC1 individuals were evaluated during five successive experiments conducted between 2010 and 2017 (Table S4). Three reference genotypes were also tested: the susceptible V. rupestris cv. ‘du Lot’, the resistant BC1 rootstock ‘Nemadex Alain Bouquet’ from the cross ‘VRH8773 (V. vinifera x M. rotundifolia) x 140 Ru (V. berlandieri x V. rupestris)’ and the susceptible to intermediate genotype ‘VRH8624’. VRH8773 is a brother clone of VRH8771, and VRH8624 in an F1 hybrid (V. vinifera x M. rotundifolia) whose muscadine parent is the accession ‘Trayshed’ (Esmenjaud et al. 2010).
For each evaluation, homogenous hardwood cuttings of the individuals were rooted annually in alveolated plates in the nursery at INRA UMR EGFV (Bordeaux, France) in February. In May, plants were delivered to INRA UMR ISA (Sophia-Antipolis, France), planted individually into 2-L pots with six replicates per individual and grown in a greenhouse. At the end of June, each pot was inoculated with a fixed number of nematodes that ranged from 300 to 900 depending on the year of experiment (Table S5). The plants were grown for ten to twelve months, which is the time that allows approximately three to four nematode developmental cycles over two successive calendar years.
At harvest, the aerial part of the plant was cut at the collar level and removed and each pot was hermitically placed into a plastic bag and stored in a cold chamber at 6°C. This stopped plant and nematode development simultaneously in all individual replicates until plant and nematode ratings. Ratings were done sequentially, i.e. replicate after replicate. Total soil of each 2-L pot was recovered in a 10-L bucket, over which plant roots were washed individually with caution under tap water. The entire root system of each plant was rated for its root development (RD) based on a 0-5 scale and its fresh root weight (RW) was also measured. Root galling was rated for each plant using a 0-5 gall index (GI) scale derived from studies for resistance to the root-knot nematode Meloidogyne spp. : 0 = no gall; 1 = 1-10 %; 2 = 11-30 %; 3 = 31-70 %; 4 = 71-90 %; 5 > 90 % of root system galled (Barker 1985). Nematodes of each plant were extracted from the total soil suspended in the bucket using an adapted Oostenbrink method (De Goede and Verschoor 2000). In the two first experiments (2010-2011 and 2011-2012), final nematode numbers were counted under a binocular microscope. Then the ratio between nematode final and initial numbers was calculated to evaluate the mean nematode reproduction factor (RF) for each BC1 individual and the parental and reference genotypes. Individuals were classified as resistant (R) when their RF value was lower than 1 and susceptible (S) when their RF value was equal or above 1. As the two first experiments showed that RF ratings were significantly correlated with GI ratings (see Results section), no nematode extraction was performed for the three last experiments (2012-2013, 2015-2016 and 2016-2017) and individuals were classified directly as R or S from their visual symptoms.
Treatment of phenotypic data sets
The number of nodosities and the number of larvae were explored using the following generalized mixed model:
Pij = µ + individual + blockj + εij
where Pij is the observed phenotype, µ is the overall mean of the phenotypic data, ‘individual’ corresponds to the genetic differences among the BC1 individuals, ‘block’ accounts for the differences in microenvironmental conditions among the three blocks, and εij is the residual term (R LME4 package). The factor ‘individual’ was treated as a random factor, whereas the factor ‘block’ was treated as a fixed factor. The ‘Best linear unbiased predictions’ (BLUP) of random effects were extracted from the selected generalized mixed model. The BLUP values were noted by the name of the trait preceded by the word BLUP.
For each quantitative phenotypic trait, broad sense heritability was estimated using the formula h²=σ²g/[ σ²g+ σ²e/n], where σ²g is the genetic variance, σ²e is the environmental variance and n is the number of plants per accession.
All principal component analyses (PCA) and parametric and non-parametric statistical tests were performed using R version 3.4.0. Statistical significance was set at P < 0.05.
Mapping of the resistance traits
Simple Sequence Repeats (SSR) markers
DNA was extracted from 50 to 100 mg of young leaves grown in a greenhouse from 90 BC1 individuals (Table S1). Leaves were ground in 5 mL of a first buffer extraction containing sodium metabisulfite [sodium metabisulfite 20 mM, Tris-HCl pH8 0.2 M, EDTA pH8 70 mM, NaCl 2 M]. 500 µL of this homogenate was incubated with 450 µL of a second buffer with CTAB [CTAB 2 %, NaCl 1.4 M, EDTA pH8 20 mM, Tris-HCl pH8 0.1 M] during one hour at 65°C. Then the solution was centrifuged for 30 minutes at 13,000 rpm at 4°C. 500 µL of supernatant was sampled and cleaned with the same volume of a solvent chloroforme:octanol (24:1). After 20 min of spin at 13,000 rpm at 4°C, we recovered 300 µL of the aqueous phase in 450 µL of a solution of isopropanol: ammonium acetate (2:1). We left one hour at 4°C and then we centrifuged at 13,000 rpm for 20 min at 4°C. We discarded the solution and washed the pellet with 70% ethanol v/v (two times). After centrifugation (13,000 rpm for 20 min at 4°C), we discarded the ethanol, and precipitated the pellet in 100 to 200 µL of TE 0.1 X [Tris HCl 10mM, EDTA 1mM]. DNAs were quantified by NanodropTM.
A total of 217 primers pairs were used: 61 VMC (Vitis Microsatellite Consortium, managed through AGROGENE, Moissy Cramayel, France), 7 VVMD (Bowers et al. 1999), 1 VVS (Thomas and Scott 1993), 6 SC08 (Cipriani et al. 2011), 1 VrZAG (Sefc et al. 1999), 7 VVC (Thomas and Scott 1993), 59 VVI (Merdinoglu et al. 2005), 11 MRBX (Zah-Bi et al. 2010), 21 UDV (Di Gaspero et al. 2005) and marker Gf13-9 (Zhang et al. 2009). Two new series of SSRs, VVBX and VVBX-A, were designed from the genome 12X of V. vinifera cv. ‘Pinot Noir’ (PN40024) using Primer3 software (Untergasser et al. 2012). Primers characteristics of VVBX and VVBX-A markers are reported in Table S6.
PCRs were performed by a single reaction with M13-tailed forward primer (Oetting et al. 1995) conjugated with four different dyes (6FAM™, VIC®, NED™ and PET®) in 15 µl reaction volume containing : 5 ng of DNA template, 1.5 µL of 10xPCR reaction buffer, 2 mM of MgCl2, 0.2 mM of each dNTP, 0.05 µM of M13 tailed SSR forward primer, 0.2 µM of reverse primer, 0.2 µM of dye conjugated with M13 primer, and 0.2 U of JumpStart™ Taq DNA Polymerase (Sigma-Aldrich). Amplification conditions were as follows: 5 min initial denaturation step at 95°C followed by 2 cycles (30 s denaturation at 95°C, 1.5 min annealing at 60°C or 56°C and 1 min extension at 72°C) followed by 35 cycles (30 s denaturation at 95°C, 30 s annealing at 60°C or 56°C and 1 min extension at 72°C) then followed by 10 min final extension at 72°C. Visualization was performed by a 3730 DNA analyzer (Applied BiosystemTM). Eight to sixteen PCR products were pooled, according to the size of SSRs and the dyes, and analyzed in a single run. Electropherograms were analyzed using free software STRand. Any ambiguous genotypes were re-run, re-amplified or left as unknown.
Single Nucleotide Polymorphisms (SNP) markers
Two foliar disks (1.5 cm of diameter) of 128 BC1 individuals grown in a greenhouse were sampled (Table S1). DNA extractions were realized on dried leaf tissues using the same protocol as described by Cormier et al. (2019). Genotyping-by-sequencing (GBS) was performed as described by Elshire et al. (2011) (Keygene N.V. owns patents and patent applications protecting its Sequence Based Genotyping technologies) integrating two 96-well plates across 96 barcodes for library preparation. The genomic library was prepared using ApeKI restriction enzyme. Paired-end sequencing of 150 bp reads was performed on an Illumina HiSeq3000 system (at the GeT-PlaGe platform in Toulouse, France).
Raw reads were checked with FastQC (Andrews 2016), demultiplexed with a custom script (https://github.com/timflutre) and cleaned with CutAdapt (Martin 2011). Cleaned reads were then mapped to the V. vinifera cv ‘Pinot Noir’ (PN40024) genome assemblies for SNP calling. Alignment on this genome was performed using Burrows-Wheeler Aligner maximal exact match (BWA-MEM) with default parameters (Li and Durbin 2009), SAMtools and Picard (http://broadinstitute.github.io/picard/). SNP calling was performed with GATK using the hardfilter parameters (McKenna et al. 2010; DePristo et al. 2011; Van der Auwera et al. 2013).
VCFtools was used to remove SNPs with a quality score < 200 and with depth values < 10 (Danecek et al. 2011). Quality-filtered SNPs were analyzed with the major_minor and get_pseudo_test_cross scripts from HetMappS pipeline to identify pseudo-testcross markers (Hyma et al. 2015). In the variant call format (VCF) output file only sites with less than 10 % missing data were retained. Individuals with more than 50% missing data and those with genotype frequencies different from expected 1:1 marker segregation were discarded. Two sets of markers were obtained corresponding to the two parental genotypes.
Individual maps were constructed for each parental genotype following a double pseudo-testcross strategy (Grattapaglia and Sederoff 1994). Marker segregation was analysed with regard to goodness-of-fit to the expected Mendelian ratio using the Chi-square test (P < 0.05). Marker types of lm x ll and nn x np were retained for construction of maternal and paternal maps, respectively. Genetic maps and marker order were determined using the maximum likelihood (ML) algorithm with Haldane function and default parameters of JoinMap®4.1 software (Van Ooijen 2006). Linkage groups (LGs) were constructed with a minimum threshold logarithm of odds (LOD) score of 6.0. LGs were grouped and numbered based on their corresponding physical chromosome numbers (The French–Italian Public Consortium for Grapevine Genome Characterization 2007).
Methods for resistance mapping
The maximum number of nodosities (in planta assay), the maximum number of larvae (in vitro assay) and the BLUP values associated were used as the quantitative scores of susceptibility/resistance response.
Detection of quantitative trait loci (QTLs) was performed with the one-dimension scan function, scanone, of R/qtl software using a normal model but also a non-parametric analysis and the expectation-maximization (EM) algorithm method depending on the normality of the data (Dempster et al. 1977; Broman et al. 2003). Multipoint genotype probabilities were calculated beforehand using calc.genoprob with step = 1 and default parameters. Logarithm of odd score (LOD) significance threshold was estimated with 1000 permutations and for a significant level α of 0.05. An interval estimate of the location of each QTL was calculated using the 1.5-LOD intervals method of Rqtl (Broman and Sen 2009). The percentage of the phenotypic variation explained by a QTL corresponds to the regression value R² taken at its peak LOD score.
Based on SSR markers, an in silico method derived from the Bulked Segregant Analysis (BSA) was performed (Michelmore et al. 1991). In the VRH8771 x CS progeny, polymorphic markers were distributed into either 2 (ab), 3 (abc) or 4 (abcd) allelic forms. From each of the markers screened, markers retained were those for which an allele was detected in the resistant parent (VRH8771) and in part of the resistant BC1 individuals but in none of all susceptible BC1 individuals.
The resistant and susceptible phenotypes were converted to values equal to 0 and 1, respectively. The one-dimension scan function, scanone, of R/qtl software was performed with the argument model = ‘binary’ (Broman and Sen 2009). LOD significance threshold, the QTL interval and the R² value were obtained with the procedure described in the previous paragraph. A number of 47 phenotyped individuals were used for the BSA analysis with SSR markers whereas the 60 total individuals phenotyped were used for QTL analysis with SSR and SNP markers.