Samples of 50 wild grapevines were collected in three protected areas of the Black Sea region ( Supplementary Figure S1 ). There were no symptoms of viral disease. In the Crimean Nature Reserve located on the southern coast of the Crimean Peninsula, 18 grapevine samples were collected at two locations near the cities of Alushta and Yalta. In the Utrish State Nature Reserve located on the Abrau Peninsula and adjacent territories, 27 samples were collected at four locations: Lobanova Shchel, Vodopadnaya Shchel, Shirokaya Shchel, and Shkolny Khutor. The samples were transported to the laboratory at +4°C to 8°C and then stored at −20°C. An additional five samples were collected from the mountain forests of Abkhazia. These samples were stored in dried form. Detailed sample information is provided in Supplementary Table S1 .

Plant DNA was extracted from the vine and leaves using the CTAB method (Rogers and Bendich, 1985). DNA quality was verified using a Nano-500 spectrophotometer (Allsheng, Hangzhou, China).

Genotyping was performed at nine microsatellite markers: VVS2, VVMD5, VVMD7, VVMD25, VVMD27, VVMD28, VVMD32, VrZAG62, and VrZAG79, as recommended by the OIV for DNA profiling of grapevine (This et al., 2004). Pinot noir and Chardonnay were used as reference cultivars of V. vinifera, the DNA profiles of which at nine loci are available in the international database Vitis International Variety Catalogue (VIVC). Flower sex was determined using the VVIb23 marker associated with the sex locus (Merdinoglu et al., 2005). PCR products 284 bp in size were considered to correspond to the hermaphroditism allele (H), 288 bp corresponded to the female allele (F), and 304 bp corresponded to the male allele (M) (Battilana et al., 2013). Kishmish Vatkana was used as a reference cultivar for the VVIb23 marker. The primers used for the analysis are listed in Supplementary Table S2 . PCR was carried out using reagents from SibEnzyme (Moscow, Russia). The 20 µl reaction mixture included 1× Taq polymerase buffer, 0.1% BSA, 0.4 mM each dNTP, 1 U Taq polymerase, 0.3 µM each primer, and 10–40 ng of DNA template. PCR was carried out using a MasterCycler Nexus GX2 (Eppendorf, Hamburg, Germany). The amplification products were separated via capillary electrophoresis using a Nanofor 05 Genetic Analyzer and built-in software (Institute of Analytical Instrumentation of RAS, St. Petersburg, Russia). The sizes of the target alleles of the SSR loci were determined via comparison with the molecular weight marker CD-600 (Syntol, Moscow, Russia).

A population genetic analysis of samples of wild grapevine from the Black Sea region was carried out using DNA profiles at nine SSR loci. The analysis also included 61 SSR profiles of previously genotyped samples from the Crimean and Utrish nature reserves ( Supplementary Table S3 ) (Ilnitskaya et al., 2023, 2024; Korniliev, 2023), 119 SSR profiles of cultivars from the VIVC database, and 540 SSR profiles of V. vinifera ssp. sylvestris from previous studies ( Supplementary Tables S4 and S5 ) (De Andrés et al., 2012; Riaz et al., 2018; D’Onofrio, 2020; Zdunić et al., 2020; Rahimi et al., 2021; Perko et al., 2024). All of them were normalized by SSR profiles of the Cabernet Sauvignon, Pinot Noir, Syrah, Chardonnay Blanc, and Riesling Weiss varieties. To obtain more even sampling, the final dataset was formed in such a way that it covered a comparable number of samples from each population and characterized its genetic diversity.

A statistical analysis of data on polymorphisms of microsatellite loci was carried out for the groups of grapevine using the average values of the number of alleles per locus (Na), the effective number of alleles per locus (Ne), and the Shannon diversity index (I). The parameters of expected heterozygosity (He), observed heterozygosity (Ho), and Wright’s fixation index (F) were calculated using the Genalex v. 6.5 software (Peakall and Smouse, 2012).

The degree of genetic differentiation (Fst) and level of gene flow (Nm) between wild grapevine populations from the Black Sea region and other groups of grapevine were calculated using Genalex v. 6.5 software. Principal Coordinates Analysis (PCoA) was performed with PAST v. 4.17 software using the Dice similarity index (Hammer et al., 2001). The PCoA projections were represented in a two-dimensional scatter plot. Genetic relatedness between grapevine groups was assessed by the UPGMA method using POPGENE 1.32 software (Yeh and Boyle, 1997).

A Bayesian analysis was performed using STRUCTURE 2.3.4 software with the following settings: burn-in period = 100,000, reps = 100,000, and number of iterations = 5 (Pritchard et al., 2000). The optimal number of clusters was calculated with Structure Harvester software using deltaK values, as defined in Evanno et al. (2005).

Total RNA was extracted from grapevine leaves and petioles using CTAB buffer (Morante-Carriel et al., 2014). Total RNA extraction was conducted on individual plants. The quality of the RNA was verified using an Eppendorf BioSpectrometer (Eppendorf, Hamburg, Germany) and electrophoresis in a 1.2% agarose gel. The RNA concentration was measured on a Qubit 3.0 fluorometer (Invitrogen, Waltham, MA, USA) using the Qubit RNA BR Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA).

Prior to library preparation, total RNA was subjected to DNase I treatment (Thermo Fisher Scientific, Waltham, MA, USA), and plant ribosomal RNA was depleted using the RiboMinus Plant Kit (Thermo Fisher Scientific, Waltham, MA, USA).

cDNA sequencing libraries were prepared using the QIAseq Stranded RNA Library Kit (Qiagen, Hilden, Germany) and NEBNext Ultra II DNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA). The DNA concentration was measured using a Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA) on the Qubit 3.0 fluorometer (Invitrogen, Waltham, MA, USA). The quality of the libraries was assessed using a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). The prepared libraries were paired-end sequenced (2 × 150 bp) on a NovaSeq 6000 System (Illumina, San Diego, CA, USA). The bioproject containing the sequenced libraries was submitted to GenBank under accession number PRJNA1148509.

Bioinformatics processing of the HTS data was performed using Geneious Prime v. 2023.2.1 software (Biomatters, Ltd., Auckland, New Zealand). Low-quality and adapter sequences were removed using the BBDuk plugin (minQ = 30, min length = 10 bp). The trimmed reads were merged, and duplicates were removed at kmer = 30.

Taxonomic classification of the preprocessed reads was performed using the Kraken2 tool (Wood et al., 2019) as described by Haegeman et al. (2023). The PlusPFP database (12/1/2024) with a size of 171 GB containing reference genomes of archaea, bacteria, viruses, plasmids, protozoa, fungi, plants, and humans was used (Index zone by Ben Langmead, 7 October 2024).

Preprocessed reads were mapped to reference genomes of grapevine viruses and viroids (latest update of the genome list: March 2024) with medium-low sensitivity (maximum mismatches per read = 20%, word length = 18). To prevent false-positive identification, a threshold of 10 reads with 15% reference coverage was set for viruses with a large number of reads: GRSPaV, grapevine pinot gris virus (GPGV; Trichovirus pinovitis), and grapevine virus T (GVT; Foveavirus tafvitis).

De-novo assembly was performed using two assemblers, SPAdes, and Geneious, with default settings. Contigs were compared to the NCBI reference viral genome database (accessed on 15 September 2022) using the TBLASTX algorithm. In further analysis, only contigs with E-value ≤ 1 × 10⁻⁴⁰ were considered.

Consensus sequences obtained after mapping reads to reference genomes or whole-genome contigs (if they were available) were used to assemble the genomes of viruses and viroid. If there were gaps in the consensus sequence, we performed a BLASTN search against the GenBank database and selected the closest genome. This genome was used as a reference for remapping. The consensus sequences contained only sample-specific data. Genomes were considered assembled when they covered more than 90% of the reference sequence and had ambiguities of less than 5% of the consensus length. The GenBank accession numbers of the assembled genomes are provided in Supplementary Table S6 .

The results of the bioinformatics analysis were visualized using the ggplot2 and ggnet2 packages in the R programming language in the RStudio integrated development environment.

Plant virus contigs that showed matches with unexpected grapevine viruses after TBLASTX analysis were analyzed to search for novel grapevine viruses. These contigs were searched against the GenBank database using BLASTN and BLASTX. For further analysis, we used contigs that did not give any matches to sequences from the database using BLASTN but were similar to viral sequences using BLASTX. These contigs were used to map reads from each library; upon successful mapping, all preprocessed reads from these libraries were combined into one data pool. The reads that mapped to the V. vinifera genome (GCF_000003745) were removed from this pool, followed by de-novo assembly using SPAdes and Geneious assemblers. The resulting contigs were analyzed using TBLASTX against the NCBI reference viral genome database as well as a local database containing nucleotide sequences of viruses of the family or the order to which the previously identified contigs belong.

All the selected contigs were extended by mapping reads from all the libraries onto them with several iterations and aligned to the genome of the closest virus. Gaps between contigs were filled by bidirectional Sanger sequencing using the BigDyeTM Terminator v3.1 Cycle Sequencing Kit (Thermo Fisher Scientific, Waltham, MA, USA) on an ABI PRISM 3730 automated sequencer (Applied Biosystems, Foster, CA, USA). The primers are listed in Supplementary Table S2 . The sequences were analyzed using Finch TV 1.4.0 (Geospiza Inc., Seattle, WA, USA).

The 5’ and 3’ ends of the novel virus were identified using the Step-Out rapid amplification of cDNA ends (Step-Out RACE) (Matz, 1999); cDNAs were synthesized using the Mint cDNA Synthesis Kit (Cat. #SK001, Evrogen, Moscow, Russia) and total RNA as a template. Step-by-step cDNA amplification was then performed using the external priming technique with universal adapter primers from the Mint RACE primer set (Cat. #SK004, Evrogen, Moscow, Russia) and 5’ and 3’ ends specific primers ( Supplementary Table S2 ). The primers were checked for secondary structures using the Beacon Designer free tool (Premier Biosoft International, San Francisco, CA, USA). The polymerase from the Encyclo Plus PCR Kit was used for amplification (Cat. #PK001, Evrogen, Moscow, Russia). The amplification products were separated by electrophoresis in a 1.2% agarose gel, excised using the Cleanup Standard Kit (Evrogen, Moscow, Russia), and cloned into the pAL2-T vector using the Quick-TA Kit (Cat. #TAK02, Evrogen, Moscow, Russia). After chemical transformation of competent Escherichia coli XL1-Blue cells (Cat. #CC001, Evrogen, Moscow, Russia), the cells were cultured on LB medium supplemented with ampicillin, XGal, and IPTG (Casali et al., 2003). Cloned DNA fragments were sequenced using the primers M13F and M13R from the Quick-TA Kit.

The full genome sequence of a novel virus was assembled based on the nucleotide sequences from the same library to eliminate polymorphisms of different isolates from the sequence. The assembled genomes were deposited in GenBank ( Supplementary Table S6 ).

The genome annotation of novel viruses was carried out using the NCBI ORFfinder tool. The prediction of conserved protein domains was performed using the NCBI Conserved Domain Search and InterPro database (Paysan-Lafosse et al., 2023). A pairwise comparison of the sequences of novel viruses with the closest known viruses was performed using the Clustal W alignment algorithm in the Sequence Demarcation Tool (SDT) (Muhire et al., 2014).

The presence of each bioinformatically identified virus and viroid was verified in all samples using RT-PCR. cDNAs were synthesized using 1 μg of total RNA as a template, random hexamers, and RevertAid H Minus Reverse Transcriptase (Thermo Fisher Scientific, Waltham, MA, USA). As a control for successful cDNA synthesis, a fragment of the 18S rRNA gene was amplified. The reaction mixture for RT-PCR comprised 0.2 mM of each dNTP, 1× Taq Buffer with (NH₄)₂SO₄, 2.5 mM MgCl2, 1 µM of each primer, and 0.375 U of Taq polymerase (Thermo Fisher Scientific, Waltham, MA, USA). All primers are listed in Supplementary Table S2 . The PCR products were visualized in a 1.2% agarose gel with the 100+ bp DNA Ladder (Evrogen, Moscow, Russia).

In cases where bioinformatics identification coincided with PCR validation, the sample was considered virus-positive.

Nucleotide sequences for phylogenetic analysis of the viruses identified in this study were downloaded from GenBank. Multiple sequence alignments were performed using the Clustal Omega 1.2.2 algorithm. The optimal model of DNA sequence evolution was determined in MEGA11 software (Tamura et al., 2021). Dendrograms were constructed in MEGA11 using the maximum likelihood method with 1,000 bootstrap replicates. The tree was rooted with an outgroup that was selected individually for each virus among related taxa. Detailed information about the settings for each dendrogram is provided in Supplementary Table S7 .

Population genetic parameters were estimated for the most widely represented viruses in the studied samples (VCV, GPGV, GVT, and GRSPaV). For the analysis of VCV, GPGV, and GVT, we used sequences prepared for phylogenetic analysis. Representatives of small populations and isolates with unknown countries of origin were removed. For GRSPaV analysis, we used >95% assembled genomes from GenBank. Sequences aligned by the Clustal Omega algorithm were cut off at both ends.

Probable recombination events were determined in RDP v. 4.101 software (Martin et al., 2015) using RDP, GENECONV, Chimaera, MaxChi, BootScan, SiScan, and 3Seq algorithms. In further analysis, only recombination events confirmed by five or more algorithms were considered.

The assessment of nucleotide diversity in populations was performed using MEGA11. The diversity index π was calculated as the average number of nucleotide differences per site between any two sequences randomly selected from the population (Nei, 1987). To test the hypothesis of whether populations are affected only by neutral evolution or they are under selection pressure, we performed Tajima’s D analysis (Tajima, 1989) in Arlequin v. 3.5.2.2 software (Excoffier and Lischer, 2010). At D > 0, the population was considered to be under balancing selection, which maintains genetic polymorphism; at D < 0, the population was assumed to have recently passed through a bottleneck or was under purifying selection; at D = 0, the population was considered to evolve under the influence of random processes in accordance with neutral evolution.

Analysis of molecular variance and calculation of the proportion of genetic diversity due to differences among populations (Φst) were performed in Arlequin v. 3.5.2.2 software with the «Compute distance matrix» setting. The Φst index is an analog of Fst for haplotype data and can have values from 0 to 1 (Holsinger and Weir, 2009). At Φst → 0, populations were considered to be connected by frequent gene flow; at Φst → 1, the populations were considered to be genetically separated.

An analysis of wild grapevine samples from the Black Sea region collected in this study using nine microsatellite markers revealed that the total number of identified alleles ranged from 8 (for VVMD25, VVMD27) to 13 (VVMD28), with an average of 10 alleles per marker ( Supplementary Table S8 ). An analysis using the VVIb23 marker revealed that 37 samples carried the F allele (female), and 27 carried the M allele (male). The H allele (hermaphrodite) was found in one sample from the Utrish Nature Reserve ( Supplementary Table S9 ). Another six alleles of different lengths, for which a statistically significant association with flower sex had not previously been established (Battilana et al., 2013), were detected in 20 samples. As a result, dioecy was reliably determined for 30 plants.

The highest genetic diversity was established for the wild grapevine population from the Utrish Nature Reserve: Na = 8.333, Ne = 4.211, whereas for the population from the Crimean Nature Reserve, these parameters were 6.000 and 3.191, respectively ( Supplementary Table S8 ). The population from Abkhazia was the least diverse in terms of allele composition (Na = 3.333, Ne = 2.765). This is probably due to the uneven distribution of this population and the limited number of samples analyzed. The fixation index F for the three analyzed populations had negative values, indicating a high degree of heterozygosity.

An analysis of the genetic kinship of the wild grapevine genotypes from the Black Sea region collected in this study and genotyped previously ( Supplementary Table S3 ) by PCoA revealed that the population of grapevine from the Utrish Nature Reserve was the most genetically heterogeneous ( Figure 1A ). The distribution areas of samples from Crimean and Utrish nature reserves on the graph partially overlap, which may indicate the kinship of these populations.

Bayesian analysis of the population structure of wild grapevine from the Black Sea region generally revealed similar results. According to the data of the Structure Harvester, the population structure in this sample selection was best described at three and six clusters ( Supplementary Table S10 ). At K = 3, the samples from Abkhazia grouped together with some of the samples from Utrish, and the other two clusters generally corresponded to the geographical division into the Utrish and Crimean populations ( Supplementary Figure S2 ). At K = 6, the populations from Utrish and Crimea divided into three clusters each, and one of the clusters consisted of both Crimean and Utrish samples, which may indicate their genetic relationship ( Figure 1B ).

The populations of wild grapevine from the Black Sea region were used for pairwise comparisons of the Fst and Nm indices with the cultivars of V. vinifera ssp. vinifera and V. vinifera ssp. sylvestris from different geographic regions ( Supplementary Table S5 ). Populations from the Crimean and Utrish nature reserves were genetically most similar to each other (Fst = 0.050) and were characterized by a high value of gene flow (Nm = 4.785). The closest to them in values of Fst and Nm were populations of wild grapevine from Italy and the Balkan Peninsula countries, along with a group of European cultivars ( Table 1 ). At the same time, native varieties and wild grapevine from Crimea were genetically isolated from each other (Fst = 0.201, Nm = 0.992). Using the present dataset, the population of Abkhazian wild grapevine was genetically close to Abkhazian native varieties (Fst = 0.049, Nm = 4.897) and isolated from other populations of wild grapevine from the Black Sea region (Fst = 0.200–0.206, Nm = 0.967–0.997).

As a result of PCoA of wild and cultivated grapevines ( Supplementary Table S5 ), the samples were divided into three groups ( Figure 2A ). The first group included genotypes of V. vinifera ssp. sylvestris from Transcaucasia (Azerbaijan, Armenia, Georgia). The second group comprised cultivars (of both European and Asian origin), wild grapevines from Israel, and, partly, Spain. The largest, third group included V. vinifera ssp. sylvestris from Europe, some European cultivars, and wild grapevine from the Black Sea region.

The dendrogram reflecting the hierarchical division of the sample groups in general clustered similarly to the PCoA plot ( Figure 2B ). Wild grapevine from the Crimean and Utrish nature reserves clustered together with European V. vinifera ssp. sylvestris, whereas wild grapevine from Abkhazia clustered together with native varieties from Abkhazia and Georgia.

A detailed study of the population structure of wild and cultivated grapevine was made possible by Bayesian analysis in STRUCTURE software. The optimal number of clusters for describing the population structure in the available sample selection ( Supplementary Table S5 ) was three, four, eight, and ten ( Supplementary Table S10 ). The analysis at K = 3 in general made it possible to divide grapevines into cultivated, wild European, and wild Asian ( Figure 3 ). However, wild grapevines from Israel and southern Spain were assigned to the same cluster as cultivated grapevines. Almost every geographic population had samples distributed among different clusters. At higher K values, wild grapevines split into subgroups based on geography. The most complex structure was observed in the wild grapevine samples from Croatia and Utrish.

At K = 10, most of the wild grapevines from Crimea formed a single cluster ( Figure 3 ). It also included a few samples of V. vinifera ssp. sylvestris from Slovenia, Croatia, and Utrish. The genotypes of the wild grapevine from Abkhazia belonged to a cluster that included some samples of V. vinifera ssp. sylvestris from Georgia, Armenia, Azerbaijan, and Utrish, as well as some Caucasian varieties. The samples of wild grapevine from Utrish did not form a single subgroup and belonged to clusters that included wild grapevine from both Europe and Asia.

The results of the analyses indicate that wild grapevine populations from the Crimean and Utrish nature reserves are genetically close to V. vinifera ssp. sylvestris from Europe. This observation, together with the absence of the hermaphroditism allele in most of the studied samples, allows us to attribute wild grapevine from Utrish and Crimea to V. vinifera ssp. sylvestris. The population of wild grapevine from Abkhazia was related to both the native varieties from Abkhazia and Georgia and to a number of samples of V. vinifera ssp. sylvestris from this region.

The sequencing of 50 total RNA libraries of grapevine from the Black Sea region yielded a total of approximately 793 million reads, with an average of 16 million reads per library. Owing to preprocessing, the number of reads decreased by more than three times and, on average, comprised 5 million reads per library. The detailed sequencing and preprocessing results are provided in Supplementary Table S11 .

A smaller number of reads were obtained from the Abkhazian samples than from the other samples; on average, there were 7 million raw reads per library and only 2.5 million reads after preprocessing. Sequencing of the Crimean samples made it possible to obtain an average of 17 million reads per library and 5.7 million after preprocessing. The data obtained from sequencing samples from the Utrish Nature Reserve were comparable to the results of sequencing the Crimean Nature Reserve samples. On average, 17 million raw reads per library were obtained; after preprocessing, their number was reduced to ~4.4 million per library. Most of the preprocessed reads in all the libraries were assigned by Kraken2 as plant reads ( Supplementary Figure S3 and Supplementary Table S12 ). Viral reads ranged from 0% to 0.14%. Some of the preprocessed reads were assigned by Kraken2 to phytopathogenic bacteria, such as Agrobacterium tumefaciens, Allorhizobium ampelinum, Agrobacterium vitis, Pseudomonas syringae, and several species of phytoplasmas ( Supplementary Figure S4 ).

Preprocessed reads from each library were used for de novo assembly. Using Geneious assembler, between 20 thousand and 92 thousand contigs were assembled. Using the SPAdes assembler, from 200 to 8,000 contigs were assembled ( Supplementary Table S11 ).

The summarized results of read mapping, TBLASTX contig analysis, Kraken2 read classification, and PCR validation are provided in Supplementary Table S13 . In 50 libraries, a total of ten known grapevine viruses and one viroid were identified.

Depending on the area, the composition of the identified viruses and viroids varied greatly ( Supplementary Table S14 ). Only VCV was detected in all three regions, and it was the only virus found in the Abkhazian samples. Wild grapevines from all four locations of the Utrish Nature Reserve were infected with viruses, among which there were no economically significant ones: VCV, GRSPaV, GPGV, and GVT. The most diverse composition of viruses and viroids was observed in samples collected in Alushta (Crimea); along with the viruses found in Utrish, it included the economically significant GLRaV-1 and GVA, as well as grapevine-associated tymo-like virus (GaTLV), grapevine satellite virus (GV-Sat), grapevine rupestris vein feathering virus (GRVFV), grapevine red globe virus (GRGV), and grapevine yellow speckle viroid 1 (GYSVd-1; Apscaviroid alphaflavivitis).

The most common virus in the wild grapevine samples from the Black Sea region was VCV, which infected 76% of the plants ( Figure 4 ). Approximately 45% of the plants were infected with GVT, GRSPaV, and GPGV, with the latter two viruses being more common in Utrish than in Crimea.

Samples of wild grapevine were found to be infected with up to nine viruses per plant ( Supplementary Figure S5 ). No viruses were found in five plants. The smallest number of preprocessed reads was found in the Abkhazian samples, and the largest was found in the Crimean samples. There was no relationship between the number of viral reads and the number of viruses detected. For example, the same number of viruses–four–was found both in a library with 2 million reads, of which fewer than 5,000 were viral, and in a library with 5.5 million reads, of which about 18 thousand were viral. Thus, the available number of reads in the libraries was sufficient for accurate virus identification.

As a result of this study, 91 complete or nearly complete virus genomes and one viroid genome were obtained, all of which were used in phylogenetic analysis ( Supplementary Table S7 ). The assembly of viruses such as GRSPaV, GRVFV, GVA, GLRaV-1, and GPGV required special settings for mapping and/or selection of the optimal reference isolate ( Supplementary Table S15 ). The genomes of the other viruses were assembled using read mapping to a reference genome with default settings.

Phylogenetic analysis of GLRaV-1 revealed that the isolates from libraries A14, A25, and A26 formed a phylogenetic group with two Russian isolates, whereas the isolate from library A15 clustered closer to an isolate from France ( Supplementary Figure S6 ). A similar situation was observed with GV-Sat isolates: the isolates from libraries A14, A19, and A25 formed one phylogroup, whereas the isolate from library A15 was located at a distance from them ( Supplementary Figure S7 ). The GYSVd-1 isolate belongs to type 3 ( Supplementary Figure S8 ). The only GRVFV isolate was phylogenetically very distant from the other Russian isolates ( Supplementary Figure S9 ). The GVA isolate was assigned to phylogroup I, which, to date, has not been associated with any symptoms in grapevine ( Supplementary Figure S10 ) (Wu et al., 2023).

The results of the population genetic analysis of grapevines from different geographical regions obtained in this study indicate a high level of genetic differentiation between Asian and European populations of V. vinifera ssp. sylvestris. Similar results were reported in a recent study of wild grapevine in Slovenia (Perko et al., 2024). Many European, Abkhazian, and Georgian varieties were genetically close to the wild grapevine of their regions. This confirms the assumption of the influence of wild forest grapevine on the formation of the gene pool of cultivated grapevine in Europe (Magris et al., 2021) and is consistent with the presence of a center of grapevine domestication in the Caucasus along with the West Asian center (Dong et al., 2023).

Populations of wild grapevine from the Crimean and Utrish nature reserves are genetically isolated from most V. vinifera ssp. vinifera cultivars and are genetically related to V. vinifera ssp. sylvestris from southeastern Europe. The samples from the Utrish Nature Reserve were distinguished by a complex population structure and were grouped together with wild grapevine from both Crimea and Abkhazia. This is probably because the Utrish Nature Reserve is located on the Black Sea coast between Crimea and Abkhazia. Moreover, wild grapevines from Crimea are not genetically related to the native Crimean varieties, which means that they represent geneplasm introduced from outside. This finding is consistent with previously obtained data on the origin of native Crimean varieties through introduction from other regions (Sekridova et al., 2022).

A different situation was observed in the analysis of wild grapevine samples from Abkhazia. A high level of gene flow between them and local native varieties, as well as genetic isolation from other groups of cultivated and wild grapevine, indicates that this population is represented by forms related to local varieties and a few samples of Caucasian V. vinifera ssp. sylvestris. An earlier study of grapevines in this geographical region also revealed a genetic relationship between wild forms and native varieties (Imazio et al., 2013; Ekhvaia et al., 2014; Riaz et al., 2018).

Thus, the population genetic analysis of wild grapevine from different areas of the Black Sea region that we performed made it possible for the first time to determine its place in the world gene pool in comparison with cultivated and wild grapevine from different geographical regions.

Recently, Haegeman et al. demonstrated that the Kraken2 metagenomic classifier can be applied to RNA-seq data for rapid analysis of preprocessed reads and detection of cellular organisms. We applied this method to our data and tested Kraken2 for virus detection (Haegeman et al., 2023). For some viruses (GPGV, GaTLV, and GV-Sat), the number of reads classified by Kraken2 was comparable to the number of reads mapped to the reference genome in Geneious Prime software ( Supplementary Table S13 ). Approximately two times fewer GRGV, GRVFV, and GLRaV-1 reads were classified by Kraken2 than were mapped to the reference genomes. For the detection of GRSPaV and GVA, the number of reads mapped to RefSeq was several orders of magnitude greater than the number of reads classified by Kraken2. Nonetheless, when the genome of the detected virus was present in the database used by Kraken2, its reads were identified in all virus-positive samples.

Compared with commercial and collection vineyards in southern Russia (Navrotskaya et al., 2021; Shvets et al., 2022a, 2022b; Belkina et al., 2023), wild grapevines from the Black Sea region were generally infected with a small number of viruses. The number of viruses per plant averaged three ( Supplementary Figure S5 ), whereas in cultivated grapevine from southern Russia, an average of 5–10 viruses per plant were detected. Moreover, in this study, no virus was detected in five plants, four of which grew in Abkhazia.

In the wild grapevine from Abkhazia, which is genetically distinct from the grapevine from the Crimean and Utrish nature reserves, only one virus (VCV) was detected, which limits our ability to conduct population genetic analysis of grapevine viruses from all three study areas. The almost complete absence of viral infections in the wild grapevine from Abkhazia is interesting and unusual; further studies will help to shed light on the virome of V. vinifera ssp. sylvestris in this area.

Grapevines from the other two areas were predominantly infected with VCV, GVT, GRSPaV, and GPGV ( Supplementary Table S14 ). The distribution and population genetics of these viruses are discussed in detail later in this article. In the wild grapevine of the Crimean Nature Reserve, we detected, along with the four viruses mentioned above, GRVFV (22%), GRGV (17%), GV-Sat (22%), GaTLV (33%), GYSVd-1 (6%), and two economically significant viruses: GLRaV-1 (33%) and GVA (6%) ( Figure 13A ). These viruses were previously identified on V. vinifera ssp. sylvestris in Tuscany (Italy): GLRaV-1 infected 11% of the analyzed grapevines, and GVA infected 7% (Sabella et al., 2018). Phylogenetically, GLRaV-1 isolates from the Crimean Nature Reserve grapevines clustered with Russian isolates from the Anapa ampelographic collection and, with a high bootstrap, formed a separate clade ( Supplementary Figure S6 ). GLRaV-1 and GVA are known to be transmitted by mealybugs ( Figure 13B ), and most other viruses found in Crimean grapevine also have known vectors (Fuchs, 2020). Notably, the viticulture industry is well developed on the Crimean Peninsula, and the planting material is usually imported from European countries. Therefore, it is likely that the listed viruses spread to the Crimean Nature Reserve grapevine from nearby commercial vineyards. Phylogenetic analysis of the isolates of these viruses found in this study demonstrated the absence of unique phylogroups specific to wild grapevine from the Black Sea region.

Only some of the viruses found in wild grapevine from the Black Sea region are associated with grapevine diseases ( Figure 13A ). The most dangerous pathogens are GLRaV-1, a causative agent of grapevine leafroll disease, and GVA, a causative agent of rugose wood disease and possibly Shiraz disease. The manifestation of these and other viral diseases occurs only when a number of factors are combined, such as a virulent genetic variant of the virus, a susceptible grapevine variety, and specific environmental conditions (Meng et al., 2017; Fuchs, 2020). Viruses such as GRSPaV, GPGV, and GRVFV are associated mainly with latent infections, but under certain conditions, they can cause rupestris stem pitting (RSP), grapevine leaf mottling and deformation (GLMD), and vein feathering disease, respectively. The other viruses detected in this study (VCV, GVT, GaTLV, GV-Sat, and GRGV) are currently not associated with the manifestation of any symptoms in grapevine. This may be due to insufficient studies and difficulty in identifying symptoms in mixed infections. However, the possibility that these viruses do not have any negative effects on grapevines cannot be excluded.

The geographic region of this study is part of the natural range of distribution of wild grapevine (Grassi and De Lorenzis, 2021). Moreover, most of the samples were collected in specially protected natural areas; therefore, these ecosystems can be expected to be minimally affected by human activity. Viruses that are widespread in wild grapevine from the Black Sea region undoubtedly play a role in the ecosystem. Apparently, they do not interfere with the growth of the studied grapevines, which have been growing in this area for decades. It is possible that certain viruses are widespread in wild grapevines precisely because they provide the host plant with a competitive advantage for survival. GRSPaV-infected grapevines under greenhouse conditions were previously shown to resist drought more effectively than uninfected plants did (Pantaleo et al., 2016). In addition, viruses from the family Partitivitidae, which includes VCV, are known as harmless persistent plant viruses, and some of them are capable of repelling herbivorous insects (Safari et al., 2019).

Thus, in some cases, viral infections can increase the ability of the host plant to withstand abiotic and biotic stresses. This advantage can be very beneficial in intraspecies competition and can lead to the survival of grapevines that are infected with mutualist viruses. This may be the reason why some viruses have been evolutionarily preserved in the Black Sea grapevine population. However, this issue requires in-depth study, since different viruses can adhere to different strategies of interaction with the host plant (Lefeuvre et al., 2019; Rodríguez-Nevado et al., 2020), and their widespread distribution in the wild grapevine population does not in itself imply a positive effect.

This study revealed that VCV, GPGV, GVT, and GRSPaV are very common in the wild grapevine from the Black Sea region. To determine how typical this situation is for this region and to trace the routes of virus migration between cultivated and wild grapevine, we used the results of metavirome studies of commercial vineyards (Navrotskaya et al., 2021) and ampelographic collections (Shvets et al., 2022a, 2022b; Belkina et al., 2023) from southern regions of Russia ( Figure 14 ). Notably, the Anapa ampelographic collection is located at a distance of approximately 25 km from the Utrish State Nature Reserve, where we collected some of the samples for this study, whereas the Magarach ampelographic collection is located on the Crimean Peninsula at a distance of approximately 50 km from the Crimean Nature Reserve.