To obtain information on the CMV population present in Serbia, samples from various host plants showing symptoms resembling CMV-like infection (Supplementary Figure S1) were collected and tested using the DAS-ELISA kit with a commercial antiserum specific for the detection of CMV (Bioreba AG, Switzerland). Out of 189 samples tested, 51 samples (26.98%) were positive for CMV and in total 17 samples were selected for further analysis based on host plant symptoms and origin: six isolates from cucurbits, one from tobacco, two from pepper, three from tomato, one from bean, three from ornamentals, and one from weeds. In addition, two isolates (674–11 and 137-13), previously characterized at the molecular level (Zečević et al., 2023), were also included and used as reference isolates. Metadata regarding the CMV isolates, including their location of origin and the original host plant, are shown in Table 1.

Total RNAs were extracted from 100 mg of symptomatic plant materials using the RNeasy Plant Mini Kit (Qiagen, Germany). To obtain complete or partial fragments of different CMV genomic regions, reverse transcription (RT)-PCR was performed, using the One-Step RT-PCR kit (Qiagen) and previously described primers given in the Supplementary Table S1 (Zečević et al., 2023). Assays were performed in a reaction mixture of 25 µL total volume, containing 5 µL 5× Qiagen OneStep RT-PCR buffer (containing 12.5 mM MgCl2), 400 µM each dNTPs, 1 µL Enzyme Mix (Omniscript Reverse Transcriptase, Sensiscript Reverse Transcriptase and HotStarTaq DNA Polymerase), and 2 µL extracted RNA. The final primer concentrations were specific for each primer: 1.5 µM RNA1a-fwd, RV11 and CMVCPrev, 3 µM 2brev and CMV3a-rev, 4.5 µM RNA1a-rev and CMVCPfwd, 6 µM RW8, 2bfwd and CMVMP3. Reverse transcription was performed at 50 °C for 30 min, followed by a PCR denaturation step at 95 °C for 15 min, and final extension at 72 °C for 10 min. Cycling conditions and number of cycles were specific for each primer pair (Zečević et al., 2023). Amplicons were analyzed on 1% agarose gel stained with Midori Green Advance DNA Stain (Nippon Genetics Europe, Germany), at 100 V for 45 min and visualized using a UV transilluminator. Additionally, the RT-PCR products of all ORFs from each Serbian isolate were analyzed in five replicates using the previously described RT-PCR-RFLP protocol (Zečević et al., 2023). No mixed reaction patterns were detected in any of the five analyzed ORFs, confirming the absence of mixed infections in the selected Serbian CMV isolates (data not shown).

Products of an expected size for each genomic region/isolate combination, obtained in RT-PCR assays, were purified using a QIAquick PCR Purification Kit (Qiagen) and Sanger sequenced in both directions on an automated sequencer (Macrogen-Europe BV, the Netherlands), using the same primers as for amplifications. For each isolate, nucleotide sequences from both strands of each ORFs were aligned using ClustalW program (Thompson et al., 1994) implemented in MEGA12 software (Kumar et al., 2024), visually inspected and manually edited. Each consensus sequence was verified by BLASTn (Altschul et al., 1990) and deposited in GenBank under the accession numbers listed in Table 1.

For each of the five genomic loci, a multiple sequence alignment was generated using the ClustalW program implemented in MEGA12 software and manually trimmed to match the length of the shortest sequence. The sequences for each of the five ORFs (1a, 2a, 2b, MP and CP) were then compared pairwise to estimate genetic diversity by calculating nucleotide (nt) and amino acid (aa) identities using the p-distance model implemented in MEGA12 (Kumar et al., 2024). Single nucleotide polymorphisms (SNPs), defined as positions in the genome where a single nucleotide differs between sequences, were identified for each ORF using DnaSP 6.12.03 (Rozas et al., 2017). SNPs were used to assess genetic variation within and between CMV subgroups, providing insight into the diversity and evolutionary dynamics of the virus population.

DnaSP 6.12.03 was also used to: 1) identify the number of haplotypes (each represented by a set of SNPs inherited as a unit) and haplotype diversity (Hd); 2) analyze fixed differences and shared polymorphisms; and 3) estimate the total number of segregating sites (S), total number of mutations (Eta), nucleotide diversity (π) – average number of nucleotide differences per site between randomly chosen sequences, which accounts for the frequency of polymorphisms, nucleotide polymorphism (θ) – a measure of total diversity based on S, irrespective of the frequency of variants, as well as Tajima’s D value to assess deviations from neutrality and infer potential demographic or selective processes shaping genetic variation. Additionally, molecular evolutionary analyses based on the rates of synonymous/non-synonymous changes across the five coding regions were performed using MEGA12 to test the neutrality hypothesis.

To investigate phylogenetic relationships and evolutionary histories at individual loci, five different alignments previously described (1a, 2a, 2b, MP and CP) of Serbian sequences were further aligned with representative strains downloaded from the NCBI database (Supplementary Table S2) and phylogenetic trees were constructed using maximum likelihood (ML) method implemented in MEGA12. The best-fit substitution model was selected based on the Bayesian Information Criterion (BIC), and branch support was evaluated using 1,000 ultrafast bootstrap replicates and the SH-like approximate likelihood ratio test (aLRT) with 1,000 replicates. All branches with bootstrap value support <70% were collapsed. RNA sequences from peanut stunt virus strain P (NCBI GenBank accession numbers EU570236, EU570237, and EU570238 for RNA1, RNA2, and RNA3 segments, respectively) were included in the phylogenetic analyses as outgroup.

To assess the potential for genetic exchange through recombination within the Serbian CMV population, a phylogenetic network analysis was performed using SplitsTree 6.4.13 (Huson and Bryant, 2024). For all ORFs, CMV sequences used for phylogeny were aligned, trimmed to the length of the shortest sequenced fragment, and concatenated to obtain sequences of 735, 822, and 1,175 nt for the RNA 1, RNA 2 and RNA 3 segments, respectively. The Neighbor-Net method was applied on concatenated RNA 1, RNA 2 and RNA 3 datasets using default settings, and branch support was evaluated with 1,000 bootstrap replicates.

Recombination analyses were performed using RDP v4.56 software package, based on concatenated sequences of the RNA 1, RNA 2, and RNA 3 segments, derived from the corresponding individual loci as explained for splits tree analysis. Possible recombination breakpoints in all sequences were identified using RDP (Martin and Rybicki, 2000), GENECONV (Padidam et al., 1999), BOOTSCAN (Salminen et al., 1995), MAXCHI (Smith, 1992), CHIMAERA (Posada and Crandall, 2001), and SISCAN programs (Gibbs et al., 2000), all implemented in the RDP4 package (Martin et al., 2015), with default parameter settings and a Bonferroni corrected P-value of 0.01. Signals were confirmed as true if recognized by at least four detection methods with a high significance level (P-value of <1.0×10⁶) (Liu et al., 2009). These analyses also determined which sequences were the most similar to the recombinant ones, indicating them as ‘parental isolates’ - lineages from which the analyzed regions of the recombinant genomes likely originated.

Nucleotide (nt) and corresponding amino acid (aa) distances between 19 Serbian CMV isolates were examined at all five ORFs. Across all regions (1a, 2a, 2b, MP and CP), nt and aa identities ranged broadly indicating remarkable sequence diversity among Serbian isolates (Supplementary Table S3). Region 2b showed the highest nt and aa diversity (63.5% to 100% nt, 43.2% to 100% aa identities), followed by regions 2a, MP, CP, and 1a. Nucleotide distances allowed for the division of isolates into two groups, except in the 1a region, where nt identities suggested the existence of three groups. The classification of isolates, and thus the number of isolates per group, varied depending on the locus examined. Notably, when intragroup nt distances were considered, they remained above 94.9%, indicating much higher sequence identity among isolates within the same group (Supplementary Table S3).

To further characterize the genetic variation of the Serbian CMV population, several genetic diversity parameters were calculated (Table 2). Measures of genetic variation (such as nucleotide and haplotype diversity), both within subgroups and species-wide (Table 2), indicate a high overall level of genetic diversity in the tested Serbian CMV population. Gene 2b exhibited the highest level of genetic variation among all five analyzed regions, as illustrated by Hd and π values of 0.942 and 0.177, respectively (Table 2).

With a goal of identifying SNPs, we examined five coding regions representing 2712 nt, in 19 CMV isolates originating from 13 hosts in Serbia, yielding a dataset with over 51500 nt (Table 2). Based on nucleotide distances, as well as phylogenetic trees constructed for individual loci (see below), each isolate was classified into one of the two subgroups, IA and II (three subgroups for the 1a gene: IA, IB, and II), which were then treated as populations. Sequences at each locus were compared within and between subgroups to determine levels of intra- and inter-group variation in CMV populations. Analysis of the concatenated dataset revealed 738 SNPs, resulting in an overall SNP density of approximately one per 70 nt, with variation ranging from one per 54 nt in the 2b gene to one per 99 nt in the CP gene. However, the extent of polymorphism within subgroups varied distinctly. At the CP locus, for instance, one SNP per 293 nt and one per 385 nt were detected within subgroups I and II, respectively, illustrating lower level of genetic variation within each subgroup. The degree of polymorphism within subgroups also varied across loci: isolates belonging to subgroup II exhibited less genetic variation at the 1a, 2b, and CP loci, whereas the opposite trend was observed at the 2a and MP loci.

With the goal of determining the subgroup affiliation of Serbian CMV isolates originating from various hosts and geographic regions, as well as studying evolutionary histories of individual loci, phylogenetic analyses for each of the five ORFs were conducted, based on a comparison with 32 CMV isolates, retrieved from the GenBank and representing different genotypes. Phylogenetic trees constructed via the ML method, based on nucleotide sequences of the 2a, 2b, MP, and CP genes, designated all Serbian CMV isolates as either subgroup IA or subgroup II, with distinct isolate distributions between the two categories depending on the locus (Figures 1B-E). For instance, in case of the 2a gene, the vast majority of isolates analyzed in this study (17 out of 19) clustered into subgroup IA, while only two isolates clustered in subgroup II (Figure 1B). On the other hand, phylogenetic tree derived from the 1a gene showed that the Serbian CMV isolates separated into three subgroups (IA, IB, and II) (Figure 1A). Considering all five phylogenies, these findings are consistent with previously recognized CMV subgroups worldwide and indicate the presence of both subgroup I and II isolates within the population. It is evident that subgroup IA isolates are the most prevalent in the examined Serbian CMV population, while isolates belonging to the IB subgroup were detected only at the 1a gene. It is also important to note that isolates originating from Serbia clustered with international isolates, irrespective of their locations of origin, host, or year sampled.

In addition to the phylogenetic analyses, estimates of DNA divergence (fixed differences and shared mutations between CMV subgroups) were calculated to quantify genetic differentiation among CMV subgroups (Table 3). A fixed difference is defined as a nucleotide change that would be present in each sequence from one subgroup but would not be present in the sequences of the other subgroup. A high number of fixed differences and a virtual lack of shared mutations between subgroups I and II (Table 3), suggests that the majority of genetic variation observed in the dataset results from the divergence between these two major CMV lineages. Indeed, an examination of isolate-by-isolate fixed differences at individual loci revealed DNA divergence characterized by significantly fewer changes between isolates placed by phylogeny into the same phylogroup and significantly more between sequences originating from distinct phylogroups (Supplementary Figure S2). For example, at the 1a locus, 270-09 (IA) and 201-11 (IA) isolate comparison revealed four fixed differences, while 149 fixed polymorphisms were detected between isolates 270-09 (IA) and 1-12 (II) (Supplementary Figure S2). Additionally, it should be noted that, within subgroup 1A, isolates 674–11 and 286–12 formed a separate clade from the remaining nine Serbian isolates (Figure 1A), which is mirrored in larger number of fixed differences (up to a dozen differences between any pair of the nine isolates from one clade and approximately 30 mutations separating two clades within a subgroup) (Supplementary Figure S2). Isolate-by-isolate examination of shared mutations yielded zero shared polymorphisms at all five loci (data not shown). Altogether, this analysis enabled defining isolate- and subgroup- specific mutations in the CMV population.

The phylogenetic network analysis of concatenated sequences of all five genes (1a, 2a, 2b, MP, and CP) showed a reticulate topology, indicating possible recombination and/or reassortment events among the CMV isolates (Figure 2). The reticulation patterns confirmed the preliminary observations derived from the comparisons of individual gene phylogenetic trees. Notably, phylogenetic trees based on the 2a and 2b genes (RNA 2) suggested that isolates 101-08, 581-11, 1-12, and 267–13 may be recombinants. These isolates clustered within subgroup IA based on the 2a gene sequence but were placed in subgroup II according to the 2b gene. To further assess the putative recombination signals observed in phylogenetic analyses, we performed recombination detection analyses with the RDP4 program. The program confirmed recombination events in the examined isolates, supported by four to six different detection algorithms and highly significant P-values (Table 4). Isolates 101-08, 581-11, 1-12, and 267–13 were identified as recombinants between subgroups IA and II in RNA 2, with breakpoints in the 2b gene (positions 464/470–814).

A similar pattern was observed when comparing the phylogenetic trees of the MP and CP genes (RNA 3), indicating a recombination event in isolate 540-10, which grouped within subgroup IA based on the MP gene and within subgroup II based on the CP gene. This isolate was identified as a recombinant between subgroups IA and II within the CP gene (positions 634–1172) by five algorithms implemented in RDP4 (Table 4).

In addition to recombination, phylogenetic tree comparisons suggested at least three obvious reassortants among the Serbian isolates included in the current study. In particular, isolate 723–10 grouped within subgroup IB based on the 1a gene, but within subgroup IA for all remaining genes. The second reassortant, isolate 79-13, belonged to subgroup II based on the 1a, MP, and CP genes and to subgroup IA based on the 2a and 2b genes. The recombinant isolate 267–13 also exhibited reassortment event, as it was placed in subgroup IB based on the 1a gene, in subgroup IA based on the 2a, MP and CP genes, and in subgroup II based on the 2b gene.

To assess the genetic diversity of Serbian CMV isolates, each isolate was classified as a genetic variant (e.g., IA-IA, IA-IA, IA) based on the phylogenetic grouping of sequences obtained from three analyzed genomic regions. For the RNA 2 and RNA 3 segments, two sequences corresponding to the genes within each segment were analyzed, and the resulting genotypic assignments are presented sequentially. Overall, based on the phylogenetic grouping of each analyzed region, seven genetic variants were defined within Serbian CMV population: ten isolates (415-07, 650-07, 670-08, 270-09, 201-11, 570-11, 674-11, 58-12, 286–12 and 473-12) possess variant IA-IA,IA-IA,IA; two isolates (137–13 and 367-14) have variant II-II,II-II,II; three isolates (101-08, 581–11 and 1-12) have variant II-IA,II-II,II, while each of the remaining four isolates can be characterized as follows: IB-IA,IA-IA,IA (723-10), IB-IA,II-IA,IA (267-13), IA-IA,IA-IA,II (540-10) and II-IA,IA-II,II (79-13).

In efforts to learn about the mechanisms responsible for maintaining the levels of genetic diversity observed in the analyzed dataset, we estimated Tajima’s D, a statistical test of neutrality under infinite sites-constant population size model. It assesses the relationship between observed (π) and expected (θ) variation and its value can be instrumental in examining evolutionary history of the locus in question. Positive D values (a higher variation than expected) can be indicative of balancing selection, population structure or a population decreasing in size, negative D values (observed variation is lower than expected) are consistent with positive or negative selection, as well as a population expanding in size. Table 5 shows Tajima’s D results for the five analyzed loci within subgroups I and II, as well as species-wide. Calculated values were negative across all loci within subgroups, although not statistically significant. On the contrary, when Tajima’s D values were assessed at the level of an entire population (considering both subgroups together), they were positive in the majority of cases, and statistically significant in the case of the CP gene.

Finally, in attempts to assess selection pressure on CMV coding regions, we estimated rates of synonymous (dS) and nonsynonymous (dN) changes in coding regions and tested for violations of neutrality using a dN/dS (ω) ratio (Table 6). Nonsynonymous sites evolving faster than expected (ω>1) are indicative of diversifying (positive) selection, while nonsynonymous sites evolving slower than expected (ω<1) suggest purifying (negative) selection. At all five loci, dN/dS ratios were less than 1.0 (Table 6). Gene 1a exhibited the lowest ω values (along with MP and CP genes), approximately five times lower than the 2b locus, perhaps illustrating considerably greater stringency toward amino acid changes.