The Assessment of Agrobiological and Disease Resistance Traits of Grapevine Hybrid Populations (Vitis vinifera L. × Muscadinia rotundifolia Michx.) in the Climatic Conditions of Crimea

The Crimean autochthonous grape varieties are unique by their origin and serve as a valuable source for breeding new cultivars with increased salt and frost resistance, as well as high-quality berries. However, they suffer from fungal pathogens, as the dry and hot summer months contribute to the epiphytotic course of diseases. An increase in the resistance of Crimean grape varieties is currently achieved through interspecific hybridization. In this study, we describe the genetic and agrobiological diversity of three hybrid populations obtained using the Vitis interspecific hybrid ‘Magarach 31-77-10′ as a female parent and Muscadinia rotundifolia × Vitis vinifera BC5 hybrid plants as male parents. The hybrid nature of the populations was assessed using RADseq high-throughput genotyping. We discovered 12,734 SNPs, which were common to all three hybrid populations. We also proved with the SSR markers that the strong powdery and downy mildew resistance of the paternal genotypes is determined by the dominant Run1/Rpv1 locus inherited from M. rotundifolia. As a result, the disease development score (R, %) for both mildew diseases in the female parent ‘Magarach 31-77-10’ was three times higher than in male parents 2000-305-143 and 2000-305-163 over two years of phytopathological assessment. The highest values of yield-contributing traits (average bunch weight ~197 g and 1.3 kg as yield per plant) were detected in the population 4-11 (♀M. No. 31-77-10 × 2000-305-163). Despite the epiphytotic development of PM, the spread of oidium to the vegetative organs of hybrids 4-11 did not exceed 20%. Some hybrid genotypes with high productivity and resistance to pathogens were selected for further assessment as promising candidates for new varieties.


Introduction
The Crimean Peninsula belongs to the northern part of the distribution area of Vitis vinifera L. subsp. sylvestris (Gmel.) Hegi, which is considered the wild forest ancestor of them, the famous hybrid NC6-15, was used as the resistant parent in a series of pseudobackcrosses with V. vinifera varieties [12]. The purpose of the performed backcrosses was to saturate the genome of interspecific hybrids with the genes of cultivated grapes, but with the preservation of the resistance loci inherited from Muscadinia rotundifolia. To this end, at each backcross step, a resistant individual was crossed with a different susceptible V. vinifera genotype to prevent inbreeding depression. As a result, BC5 and BC6 populations were obtained, in which resistance to powdery and downy mildew segregated in a 1:1 ratio [12,13]. The BC5 progeny Mtp 3294 (VRH 3082-1-42 × V. vinifera cv. Cabernet Sauvignon) was employed for the fine mapping, positional cloning, and functional characterization of the unique resistance locus from M. rotundifolia, containing two closely related genes MrRUN1 and MrRPV1, which confer strong resistance to powdery mildew and downy mildew, respectively [14].
The half-siblings of Mtp 3294, the genotypes 2000-305-143 and 2000-305-163, were obtained from crossing the same female parent VHR 3082-1-42 and V. vinifera cv. Regent. These genotypes were kindly provided by Prof. R. Eibach (Institut für Rebenzüchtung Geilweilerhof, Germany) to the Research Institute of Viticulture and Winemaking 'Magarach' in Yalta (Crimea). In 2011, the 2000-305-143 and 2000-305-163 genotypes were involved in the local breeding program of the 'Magarach' Institute, serving as donors of loci of resistance in the crosses with the elite parental female Magarach No. 31-77-10, which showed the highest breeding value when selecting for resistance to oidium in all crosses carried out earlier.
In this article, for the first time, we genetically assessed the obtained hybrid populations using SNP markers generated by restriction-site-associated DNA sequencing (RADseq) [15]. Discovering thousands of SNP loci, RADseq has become a valuable method in population studies, outperforming other molecular marker techniques (e.g., microsatellites) in applications that require individual-level genotype information, such as quantifying relatedness and individual-level heterozygosity [16]. Previously, we successfully employed the ddRADseq method for genotyping of non-model species [17,18] and found the technique as an efficient and low-cost approach to de novo SNP discovery.
The aim of our study was the assessment of the agrobiological and disease resistance traits of hybrid populations of grapevine, which presumably carry the introgressions of the species M. rotundifolia, in ecological conditions of the Crimea, as a specific area of viticulture.  Figure 1). Both resistant paternal genotypes are expected to carry MrRUN1 and MrRPV1 genes from M. rotundifolia, since their half-sibling Mtp 3294 was used to dissect those resistance loci [14]. Additionally, 2000-305-143 and 2000-305-163 may also inherit some additional resistance loci from their parent cv. Regent-a cultivar with quantitative resistance against downy and powdery mildew released in Germany in 1996 [19]. For cv. Regent, at least two resistance loci for Erysiphe necator (Ren3, Ren 9) [20,21] and one resistance locus for Plasmopara viticola (Rpv 3.1) [22] have been reported. Forty-three progenies were obtained from the cross

Assessment of the Hybrid Population Relationship Using Genotyping by Sequencing with IlluminaHiSeq2500
To assess the genetic identity and confirm the hybrid nature of the progenies obtained from three crosses, 139 progeny genotypes (66 + 43 + 30) and parental genotypes (female M. No. 31-77-10; male 2000-305-143, 2000-305-163) were genotyped using the double-digest ddRADseq approach [15]. The precise description of the procedure performed has been published previously [17,18]. We constructed 151 ddRAD libraries from the 139 progenies of three hybrid populations and three parents (each in 4 replicates). Sequencing was performed on Illumina HiSeq2500 with single end reads of 150 base pairs.
In total, 309,762,340 high-quality reads were obtained. For the female parent, which was common for all of the three populations, 6,717,571 reads were generated; for male parents 2000-305-143 and 2000-305-163, 4,414,422 and 9,010,505 reads were obtained, respectively. The progeny of population 2-11 provided 140,445,297 reads, the progeny of population 3-11 provided 97,219,515 reads, and the progeny of population 4-11 provided 51,955,030 reads.
Aligned reads were subjected to the SNP calling procedure using the Tassel v.5.2.40 GBSv2 plug-in [24] and Stacks v. 2.53 bioinformatic software [25]. Next, the filtration of raw SNPs using VCFtools was applied to determine high-quality SNPs. As shown in Table 1, Tassel provided a higher number of raw SNPs, while Stacks kept a higher number of filtered SNPs; therefore, subsequent analysis was performed on the Stacks-filtered SNP data set. According to the Venn diagram in Figure 3, the three hybrid populations shared a common set of 12,734 SNPs, which were further employed for principal component analysis (PCA). The SNP data set is available in Supplementary Materials, Table S1.
First, we examined whether the hybrids of the populations 3-11 and 4-11, for which both female and male parents were determined, show their intermediate genetic background when compared to the maternal and paternal genotypes (Figure 4a Figure S1 and Table S2).
First, we examined whether the hybrids of the populations 3-11 and 4-11, for which both female and male parents were determined, show their intermediate genetic background when compared to the maternal and paternal genotypes (Figure 4a The PCA plot for all three populations taken together shows that populations 3-11 and 4-11 are closely related, as expected, and are genetically separated from population 2-11 ( Figure 4d).

Transit to Generative Stage
All crosses to obtain hybrid populations 2-11, 3-11, and 4-11 were carried out in 2011. In 2012, the seeds were germinated, providing seedlings for the first year of life. In 2013, the seedlings were planted as self-rooted plants in the experimental field. Thus, the hybrid populations have been monitored since 2013.
Long-term observations showed that the seedlings from the three crosses varied greatly in many traits, but most significantly by the year of life, when the first bunches were recorded. As shown in Figure 5, about 40% of the progeny of population 2-11 still did not enter fruiting by 2020, while the progeny from populations 3-11 and 4-11 were more than 90% fertile in 2020. We assume that the remarkable difference can be explained by the complex genetic origin of the 2-11 progeny. In our practice, even in intraspecific crosses, up to 30% of seedlings do not begin to flower, and this is influenced by a number of factors, both genetic and environmental. Long-term observations showed that the seedlings from the three crosses varied greatly in many traits, but most significantly by the year of life, when the first bunches were recorded. As shown in Figure 5, about 40% of the progeny of population 2-11 still did not enter fruiting by 2020, while the progeny from populations 3-11 and 4-11 were more than 90% fertile in 2020. We assume that the remarkable difference can be explained by the complex genetic origin of the 2-11 progeny. In our practice, even in intraspecific crosses, up to 30% of seedlings do not begin to flower, and this is influenced by a number of factors, both genetic and environmental.

Productivity Trait Variation
Focusing only on the plants that were flowering by 2019 in each hybrid population, we compared the productivity-related traits between the populations based on the field evaluation data collected in 2019-2020. As shown in Table 2, all three populations were comparable in terms of the number of latent buds, slightly varied by the number of developed shoots, but significantly differed in the number of fertile shoots, the number of inflorescences, and, correspondingly, the number of bunches per plant. The highest values of these yield-contributing traits were shown by population 4-11 ( Table 2). As a result, for the progeny of the hybrid population 4-11, the maximum yield per plant was recorded both in 2019 and in 2020 ( Figure 6).

Productivity Trait Variation
Focusing only on the plants that were flowering by 2019 in each hybrid population, we compared the productivity-related traits between the populations based on the field evaluation data collected in 2019-2020. As shown in Table 2, all three populations were comparable in terms of the number of latent buds, slightly varied by the number of developed shoots, but significantly differed in the number of fertile shoots, the number of inflorescences, and, correspondingly, the number of bunches per plant. The highest values of these yield-contributing traits were shown by population 4-11 ( Table 2). As a result, for the progeny of the hybrid population 4-11, the maximum yield per plant was recorded both in 2019 and in 2020 ( Figure 6).     Here, we present the results of phytopathological analyses for 2017 and 2019, since 2018 was unusual in its climatic conditions, which resulted in diseases developing moderately across the entire sample. In general, the resistance score of all three populations to both PM and DM pathogens varied from 5 to 9 according to the OIV-455 and OIV-452 resistance scale [26], meaning intermediate or high resistance. To expose more clearly the interpopulation diversity in terms of resistance to fungal diseases, we also recorded for each population the disease development score (R%) reflecting the percentage of affected leaves per plant and the degree of disease development (see Section 4, "Materials and Methods").
In this regard, there was a significant difference between the parental genotypes for resistance to both downy and powdery mildew: the percentage of leaves affected by the pathogens in the female parent M. No. 31-77-10 was more than three times higher than in the male parents 2000-305-143 and 2000-305-163 in both 2017 and 2019 (Figures 7 and 8).    (Figure 7).
We noticed that in individual genotypes, resistance to each specific pathogen can be expressed to varying degrees. There were genotypes with high resistance to PM but medium resistance to DM, and vice versa. There were also some outstanding genotypes with high resistance to both mildew diseases.

Evaluation of Parental Genotypes with SSR Markers Linked to PM and DM Resistance Loci
The strong powdery mildew resistance of the parental genotypes 2000-305-143 and 2000-305-163 is presumably supported by a single dominant locus Resistance to Uncinula necator (Run1) inherited from M. rotundifolia. To confirm the presence of resistant alleles at this locus in the paternal genotypes, we employed the closest flanking SSR markers VMC4f3.1 and VMC8g9 reported for the Run1 locus [27] (Table 3). For the VMC8g9 marker, which completely co-segregates with Run1, two alleles were registered; among them a fragment of 159 bp has been described as a resistant allele, tracing back to the NC6-15 hybrid [28]. For the second flanking marker, VMC4f3.1, localized at a distance of 0.6 cM from the Run1 locus, just one allele size of 192 bp was recorded, suggesting the homozygous state of the marker for both paternal genotypes 2000-305-143 and 2000-305-163. The 192 bp allele has also been described as being associated with powdery mildew resistance [28].
The chromosomal interval between SSR markers VMC4f3.1 and VMC8g9, flanking Run1, is located in linkage group 12 of the V. vinifera consensus map [29]. The locus of

Evaluation of Parental Genotypes with SSR Markers Linked to PM and DM Resistance Loci
The strong powdery mildew resistance of the parental genotypes 2000-305-143 and 2000-305-163 is presumably supported by a single dominant locus Resistance to Uncinula necator (Run1) inherited from M. rotundifolia. To confirm the presence of resistant alleles at this locus in the paternal genotypes, we employed the closest flanking SSR markers VMC4f3.1 and VMC8g9 reported for the Run1 locus [27] (Table 3). For the VMC8g9 marker, which completely co-segregates with Run1, two alleles were registered; among them a fragment of 159 bp has been described as a resistant allele, tracing back to the NC6-15 hybrid [28]. For the second flanking marker, VMC4f3.1, localized at a distance of 0.6 cM from the Run1 locus, just one allele size of 192 bp was recorded, suggesting the homozygous state of the marker for both paternal genotypes 2000-305-143 and 2000-305-163. The 192 bp allele has also been described as being associated with powdery mildew resistance [28]. Table 3. Assessment of parental genotypes of hybrid populations with SSR markers associated with loci of resistance to powdery mildew and downy mildew. Ren3 ScORA7 -720 -720 -- 1 The allele size associated with resistance is shown in bold. 2 The resistance-associated haplotype in chromosome 15 (Ren9-Ren3) inherited from cv. 'Chambourcin' via cv. Regent is marked by light gray.

Loci of Resistance
The chromosomal interval between SSR markers VMC4f3.1 and VMC8g9, flanking Run1, is located in linkage group 12 of the V. vinifera consensus map [29]. The locus of resistance to downy mildew (Resistance to Plasmopara viticola (RPV1)) is tightly linked to the Run1 locus in M. rotundifolia × V. vinifera BC2 hybrid plants [30]. Merdinoglu et al. [30] were able to link the Rpv1 locus to the SSR marker VMC1g3.2, while Katula-Debreceni et al.
suggested 122 bp as an Rpv1-specific allele size [31]. However, Prazzoli et al. [32] signified 118 bp as the allele size linked to resistance at the Rpv1 locus. In our experiments, we documented the presence of a resistant allele of 118 bp in both paternal genotypes 2000-305-143 and 2000-305-163 but not in the susceptible female parent M. No. 31-77-10 (Table 3).
For a deeper assessment of Rpv1, we also employed SSR markers VVIb32 иVVIM11 mapped in the same chromosome interval as VMC1g3.2 [33].  [20,21], there are at least two resistant loci to E. necator loci on the chromosome 15: Ren9 (1.1-3.5 cM) and Ren3 (9.2-9.6 cM), spanning large physical distances between each other. The authors proposed to distinguish between Ren3 (interval ScORGF15-02-ScORA7*) and Ren9 (interval CenGen7-GF15-10). The SSR marker GenGen6 closely linked to CenGen7 was suggested for selection of the resistance-mediating haplotype of Ren9, while the SCAR marker ScORA7 can be used to detect the resistant allele of Ren3 inherited by Regent from its parent cv. 'Chambourcin'. We analyzed the parental genotypes of our hybrid populations with the markers developed for the Ren9 and Ren3 loci and revealed that male genotypes 200-35-143 and 200-35-163 inherit different haplotypes from Regent in the Ren9-Ren3 interval assessed by markers GenGen6, UDV116, and ScORA7. The male parent 200-35-143 carries the resistanceassociated haplotype identical to the 'Chambourcin' grandparent [20], while 200-35-163 has inherited the chromosome fragment from the susceptible predecessor 'Diana' (Table 3).

Discussion
The climate in the south coast of Crimea is close to dry subtropical (Mediterranean) conditions that have favored viticulture on the peninsula since ancient times. Currently, the area of vineyards in Crimea has reached the level of 31,000 ha, with an average yield of 4.2 t ha −1 , and 82.8% of the vineyard area is occupied by technical grape varieties and 17.2% by table varieties [34]. Fungal diseases have been and remain the main problem for viticulture in this area.
Following the classification of phytopathogens by harmfulness, powdery mildew should be attributed to the main diseases in the western foothill-coastal region of Crimean viticulture, since it causes yield losses from 10% to 50% [35]. Oidium, according to this classification, is definitely the dominant disease, because the yield loss in some years reaches 100%. The development of oidium is epiphytotic almost every year. For example, 8 of 10 years of observation (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) at the stationary experimental site 'Livadia' showed that the development of the disease on bunches of 'Muscat White' variety before harvesting was 85.9-100%. Only in 2005 and 2007 did the disease develop moderately (42.8% and 46.5% of bunches were affected, respectively) due to high daytime air temperatures and air drought in the first half of July. Thus, seven to eight sprayings of grapes with highly effective fungicides during the growing season are a common practice in vineyard care [35].
The discovery of the Run1 gene in American Muscadine grapes (M. rotundifolia), immune against PM, opened up new possibilities in viticulture [31]. The BC4 hybrid VRH3082-1-42, derived from the M. rotundifolia × V. vinifera cross [12], was employed in breeding resistant grapevines using a gene-pyramiding scheme e.g., [31,36]. It was reported that the defense mechanism by which Run1 confers resistance appears to be different from another dominant gene for PM resistance, Ren1, identified in the Central Asian V. vinifera variety Kishmish vatkana [31].
In another study, when evaluating F1 progeny derived from the VHR 3082-1-42 × 'Regent' cross, it was stated that the effect of the resistance genes originating from 'Regent' is higher for downy mildew, while for powdery mildew, the effect of the Run1 gene is much higher, and it even covers phenotypically the effect of the resistance genes originating from 'Regent' [36]. This report is consistent with our results obtained by phytopathological assessment and by the SSR marker assay of two hybrids of this F1 progeny, the genotypes 2000-305-143 and 2000-305-163. Only the latter carries resistant alleles of Ren3 and Ren9 genes inherited from Regent, while the Run1 resistant gene was detected in both of them. However, 2000-305-143 and 2000-305-163 showed no difference in the field resistance to both mildew diseases.
Nevertheless, the incorporation of other sources of resistance to PM into a breeding program makes possible the long-term support of the Run1 gene and decreases the danger of breakdown of resistance in the event that a new pathogen race appears [10,31]. Thus, about 40% of 102 grape accessions from the international, collaborative project VITISANA showed pyramiding of R-loci against powdery mildew [10].
In this study, we described for the first time, to the best of our knowledge, the genetic and agrobiological diversity of populations of remote hybrids obtained using Magarach 31-77-10 as a female parent by crossing with donors of Muscadinia rotundifolia introgressions. The genetic relatedness and distinctiveness of three hybrid populations were assessed using 12,734 SNP markers, generated by Illumina-based genotyping by sequencing (GBS) approach. For all SNP markers, the corresponding position in the reference V. vinifera genome was determined and its allelic state in the analyzed grape genotypes was described. The available SNP data set can further be employed as a tool for evaluation of the M. rotundifolia × V. vinifera BC populations that are involved in grapevine-breeding programs around the world.
Although our results of the resistance studies and field evaluation of the hybrid genotypes are preliminary (two years of study), they suggest that the three hybrid populations could be a valuable source for breeding new grape varieties well adapted to both abiotic and biotic stresses of the Crimean environment. First, genotypes with hermaphrodite flowers predominate in all three populations: 81% in population 2-11, 52% in population 3-11, and 65% in population 4-11. Thus, the stability of fertilization, which determines the yield, can be less influenced by unfavorable weather conditions. Additionally, long-term field monitoring of disease resistance against a natural infectious background showed that all hybrid populations demonstrate high resistance to powdery mildew. The infectious pressure is always present in the experimental field where other grapevine varieties grow together with the hybrid populations. Here, for some grape varieties, up to 91% of the spread of oidium on vegetative organs and up to 100% of damage to the generative organs were recorded, while the spread of oidium on vegetative organs of the hybrids, e.g., from the M. No. 31-77-10 × 2000-305-163 cross (population 4-11) never exceeded 20% (Figure 7).
In the present study, we experimentally proved that the progeny  Figure S2).
The current breeding programs at the 'Magarach' Institute focus on creation of new varieties with complex immunity to the insect pest Daktulosphaira vitifoliae and pathogens Plasmopara viticola and Erysiphe necator based on the introgression of resistance genes from the highly resistant (immune) species Muscadinia rotundifolia. Among hybrids from the crossing of the hybrid genotype Magarach 31-77-10 with the donors of Muscadinia rotundifolia genes, progeny have been obtained in which genotypes with high productivity and high resistance to pathogens are distinguished. These are recommended for further agrobiological and technological evaluation as promising candidates for new varieties. The area of the field is formed by valleys, ridges, and slopes with different exposure. Soils of the experimental fields are characterized as brown mountainous non-carbonate, with a heavy, loamy gravelly texture. The mild climate provides the possibility of uncovered grape culture: the average annual air temperature is 13.2 • C, the sum of active air temperatures (above 10 • C) is 3700-4055 • C, the number of days with temperature above 10 • C varies between 195 and 218, and the annual precipitation is 300-400 mm. In summer, rains are short, in the form of downpours. Plants often suffer from a lack of moisture because of the dry summer.

Plant Material and Field Evaluation
Assessment of agrobiological traits was conducted in 2019-2020 according to the method of Lazarevsky [37]. In short, for each genotype, the following traits were recorded: number of latent buds, number of developed shoots, number of fertile shoots, number of inflorescences, number of bunches, average bunch weight (g), and yield per plant (kg). The actual yield from a bush was recorded by counting all bunches on the bush and weighing them.
Phytopathological field evaluation was conducted by the examination of untreated plants against a natural infection pressure. In each season, two counts were carried out: the first in 15-20 days after flowering of grapes and the second at the beginning of grape ripening.
The nature and percentage of leaf damage were scored according to the recommended method [38]. Precisely, on each counting bush, up to 50 leaves were evaluated from both sides for signs of infestations. The percentage of affected leaves and the degree of disease development on each leaf were determined using a scale: 0-no signs of infestation 1-single, hardly visible spots on leaves (OIV resistance-9 points) 2-up to 10% of leaf surface affected (OIV resistance-7 points) 3-11-25% of leaf surface affected (OIV resistance-5 points) 4-26-50% of leaf surface affected (OIV resistance-3 points) 5-more than 50% of leaf surface affected (OIV resistance-1 point) Using the scale, the disease development score (R, %) in a particular genotype was calculated by the formula where a is a score of the scale, according to which the lesion was evaluated in the experiment; b is the number of affected leaves within the range of this score; N is the total number of leaves evaluated (pcs); K is the highest score of the scale; and 100 is the conversion factor. For example, of 48 leaves evaluated for a plant, 0 degree of disease development was determined for 27 leaves, 1 for 18 leaves, 2 for 3 leaves, and 3, 4, and 5 for 0 leaves. The disease development score (R, %) was calculated as follows: Fresh leaves were collected from parental genotypes and hybrid progenies and frozen in liquid nitrogen. DNA was isolated from the plant material using the protocol described by Rahimah et al. [39] with some modifications. Specifically, 160 mg of frozen leaf material was ground with Precellys 24 tissue homogenizer (Bertin Technologies, Montigny le Bretonneux, France) at 2500 rpm for 30 s. To each of the ground samples, 800 µL of modified CTAB buffer (2% CTAB w/v, 20 mM EDTA (pH 8.0), 1.4 M NaCl, 100 mM Tris-HCl (pH 8.0), 5 mM ascorbic acid, 4mM diethyldithiocarbamic acid sodium salt, and 2% polyvinylpyrrolidone-40), 3.2 µL of 2-mercaptoethanol, and 0.4 µL RNase A were added. The samples were incubated at 37 • C for 30 min, followed by another 30 min at 65 • C for enzyme inactivation. Then, 800 µL of chloroform:isoamyl alcohol (24:1) was added and mixed thoroughly. The mixture was centrifuged at 10,000 rpm for 15 min at room temperature, and the upper aqueous phase was transferred into a sterile 1.5 mL Eppendorf tube. DNA was precipitated by adding 0.7 volume of ice-cold isopropanol. The solution was stored at −80 • C for 15 min and centrifuged at 12,000 rpm at 4 • C for 15 min. The DNA pellet obtained was washed in 200 µL of wash buffer (76% ethanol and 10 mM ammonium acetate) and centrifuged at 14,500 rpm for 5 min. Then, the wash buffer was removed with an FTA-1 Aspirator (BioSan, Riga, Latvia). The pellet was dried at room temperature and suspended in 160 µL of TE buffer (10 mM tris-HCl (pH 8.0) and 1 mM EDTA (pH 8.0)). Extracted DNA was purified by adding to each DNA sample 0.5 volume of 7.5 M ammonium acetate (pH 7.7) and incubating in ice for 20 min. After centrifugation at 12,000 rpm for 15 min at 4 • C, the supernatant was transferred into a clean 1.5 mL Eppendorf tube. The DNA in the supernatant was re-precipitated by adding 2.5 volumes of ethanol, mixed thoroughly, and centrifuged at 12,000 rpm at 4 • C for 15 min. The supernatant was removed, and the DNA pellet was washed with 200 µL of 76% ethanol, dried at room temperature for 15 min, and suspended in 40 µL of TE buffer.
The quality of DNA samples was evaluated with a SPECTROstar Nano spectrophotometer (BMG LABTECH, Ortenberg, Germany), and the optical density was determined at A260, A280, and A350. The integrity of the DNA was further examined by electrophoresing approximately 5 µg of DNA in 1% agarose gel in 1× TAE buffer (0.04 M tris base, 20 mM acetic acid, and 2 mM EDTA). The concentration of each DNA sample was estimated using Qubit 4.0 Fluorometric Quantitation (Thermo Fisher Scientific, Waltham, MA, USA).

Construction and Sequencing of ddRAD Libraries
ddRADseq libraries were prepared, as described previously [21]. ddRAD libraries were constructed from 151 individual DNA samples. The concentration of DNA for each individual sample was normalized to 10 µg/µL, and then 10 µL of the obtained DNA solution was digested for 1 h at 37 • C with 0.7 units of HindIII (NEB, Ipswich, MA, USA) in a 20 µL reaction volume. Digestion was stopped by incubating the reaction mixture at 65 • C for 20 min. Next, the barcode adapter was ligated to the end generated by HindIII, allowing pooling of the samples. The ligation reaction was carried out at 22 • C for 1 h in a volume of 50 µL with 1.6 units of T4 ligase (NEB, Ipswich, MA, USA). T4 ligase inactivation was reached by heating at 65 • C for 10 min. Next, 10 µL of each sample was pooled and simultaneously purified with AMPure XP beads. The second restriction digestion was performed with 0.7 units of NlaIII (NEB, USA) in a 20 µL reaction volume at 37 • C for 1 h with subsequent inactivation at 65 • C for 20 min. The second common adapter was ligated to the overhanging end of NlaIII. The reaction mixture (50 µL) was dispersed into 14 aliquots, and PCR was performed on the 14 aliquots to amplify the DNA fragments. The PCR mix (50 µL) contained 1× HF buffer, 25 µL of each primer, 12.5 µM dNTPs, and 1 unit of Phusion ® High-Fidelity DNA Polymerase (NEB, USA). PCR conditions were as follows: 98 • C for 30 s, 14 cycles of 98 • C for 10 s, 65 • C for 30 s, 72 • C for 15 s, and then 72 • C for 2 min. All 14 PCR reactions were pooled and purified with AMPure XP beads. The concentration and size distribution of the prepared libraries were checked by automated electrophoresis with the 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Sequencing was performed on Illumina HiSeq2500 with single end reads of 150 base pairs at the CERBALAB Company (St. Petersburg, Russia, https://www.cerbalab.ru/, accessed on 18 December 2020).

SSR Marker Analysis
For SSR marker analysis, DNA was isolated from 100 mg of fresh young leaves using the CytoSorb kit (Syntol, Moscow Russia) according to the manufacturer's recommendations.
Amplification reactions were performed in a total volume of 25 µL with 20 ng of template DNA, 10 pmol of each primer, and 10 µL of 2.5× PCR reaction mixture containing 0.25 mM dNTPs, 3 mM MgCl 2 , and 3 units of Taq polymerase with antibodies inhibiting enzyme activity. The forward primer of each pair was 5 -end-labeled with a fluorescent dye (FAM, TAMRA, or R6G) (Syntol, Moscow, Russia).
PCR was carried out using a T100 thermal cycler (Bio-Rad Laboratories, Hercules, CA, USA). The cycling program consisted of the following steps: 5 min at 95 • C, followed by 35 cycles of 15 s at 95 • C, 25 s at 58-66 • C, and 25 s at 72 • C, and a final extension step of 7 min at 72 • C.