Homozygous TYLCD susceptible (ty-1/ty-1) and heterozygous TYLCD resistant (Ty-1/ty-1) tomato quasi-isogenic lines (Monci et al., 2002) (kindly provided by M. J. Díez, Universidad Politécnica de Valencia, Spain), S. nigrum (IHSM “La Mayora” germplasm bank) and common bean cv Donna (Nunhems, Haelen, The Netherlands) were used as experimental hosts. Plants at the three-leaf growth stage were inoculated with Agrobacterium tumefaciens LBA4404 containing infectious clones, and therefore single sequence variants, of isolates [ES:72:97] of the Mld strain of TYLCV (GenBank accession number AF071228) (García-Andrés et al., 2009), [ES:Mur1:92] of the ES strain of TYLCSV (Z25751) (Navas-Castillo et al., 1999), and [ES:421:99] of TYLCMaV (AF271234) (Noris et al., 1994). These isolates were derived from infected tomato plants except for TYLCMaV, which was derived from an infected common bean plant. TYLCMaV is the result of a genetic exchange between TYLCV-Mld and TYLCSV-ES, with recombinant genome sequences that share 99% nucleotide identity with that of the parental virus sequences. As depicted in Figure 1 in TYLCMaV, the 3′ half of the IR and the virus sense V2 and the CP ORFs derive from TYLCSV, whereas the remainder 5′ half of the IR and the Rep, C4, TrAP and most of REn ORFs derive from TYLCV, with the 3′ terminal 40 nucleotides of the REn ORF derived from TYLCSV (Monci et al., 2002). Viral infectious clones of TYLCSV, TYLCV or TYLCMaV were used to inoculate three plants each of a TYLCD-susceptible tomato, a quasi-isogenic resistant tomato containing the Ty-1 allele in heterozygosis, common bean and S. nigrum (Figure 1). Liquid cultures of A. tumefaciens were adjusted to an OD₆₀₀ of 1.00 and 20 μl were inoculated to each plant by the stem-puncture method as described elsewhere (Monci et al., 2002). After agroinoculation, plants were maintained for 45 days in a growth chamber (26°C during the day and 18°C at night, 70% relative humidity, with a 16 h photoperiod at 250 μmol s⁻¹m⁻² photosynthetically active radiation). Systemic virus infection of agroinoculated plants was confirmed on young non-inoculated leaves by tissue-print hybridization using a mixture of digoxigenin (DIG)-labeled probes specific to TYLCV and TYLCSV (Monci et al., 2002) that also recognizes TYLCMaV.

To monitor quasispecies in each infected plant, young apical leaves were harvested at 15, 30, and 45 dpi (Figure 1) and lyophilized. DNA was extracted following the Edwards' procedure (Edwards et al., 1991) with some modifications. From each sample, 8 mg of dried tissue were ground in a 1.5 ml Eppendorf tube with liquid nitrogen, and 400 μl of extraction buffer [200 mM Tris-HCl, 200 mM NaCl, 25 mM EDTA, SDS 0.5% (w/v), pH 7.5, 0.1% β-mercaptoethanol] added. After vortexing at 37°C for 2 min extracts were centrifuged at 12,000 × g at room temperature for 5 min. The supernatants were transferred to a new tube and mixed with 1 volume of phenol:chloroform (1:1). The aqueous phase was transferred to a new tube and mixed with 1 volume of isopropanol. After 2 min at room temperature, nucleic acids were recovered by centrifugation at 12,000 × g at room temperature for 5 min. The pellet was washed with 70% ethanol, dried and resuspended in 200 μl of RTE buffer (50 mM Tris-HCl pH 8.0, 10 mM EDTA pH 8.0, 0.1% RNase A) by incubation at 65°C for 10 min. Viral DNA from total DNA extracts was amplified by rolling-circle amplification (RCA) with ϕ 29 DNA polymerase according to Inoue-Nagata et al. (Inoue-Nagata et al., 2004) using the TempliPhi DNA Amplification Kit (GE Healthcare). For all samples, 400 to 500 ng/μl of the amplification product measured by UV-spectrophotometry and gel electrophoresis were obtained. Full-lenght genomic consensus sequences were obtained by direct sequencing (Macrogen Inc., South Korea) of RCA products amplified from total DNA extractions. Primers used to sequence consensus sequences of TYLCSV, TYLCV, and TYLCMaV are indicated in Table 1. Sanger sequencing was preferred over NGS due to the difficulty of the latter to discriminate between the error of the technology and the error rate of the viruses, which is paramount for quasispecies analysis.

RCA products were digested with single cutter (in viral genome) restriction enzyme BamHI to obtain the linear full-length DNA (c. 2.7 kbp) that was subsequently cloned into pBluescript II SK (+) (Stratagene). The product obtained was transformed into Escherichia coli DH5α by electroporation. Transformed colonies were analyzed by PCR using the universal M13 primers (Table 1). Only those clones containing inserts representing the 2.7 kbp full-length genome size were selected for the study. Amplification of molecular clones was also performed using RCA before sequencing. A total of 18–22 clones were sequenced per mutant spectrum from a single plant. Primers MA863 and M13 UNIV (-) were used for all virus isolates (Table 1) to sequence the 5′-end of Rep and C4 ORFs, and the intergenic region (IR) on the C strand. On the virus sense strand the sequenced region comprises the 5′-end of the V2 ORF and the complete CP ORF. Primers MA1197, MA1198, and MA1199 were used for sequencing molecular clones of TYLCMaV, TYLCSV, and, TYLCV, respectively, whereas primer MA1365 was used for all of them.

EditSeq, SeqBuilder and SeqMan software (DNAStar Inc., USA) were used for sequence assembly and analysis. Mutations were scored relative to the consensus sequence of each mutant spectrum. To analyze virus quasispecies genetic complexity (mutational composition of the ensemble) mutation frequencies were calculated by dividing the number of different mutations (repeated mutations in each quasispecies were computed only once) by the total number of nucleotides sequenced (Domingo et al., 2006). Virus quasispecies genetic complexity was also assessed by analyzing the average genetic distance (d), defined as to the average number of mutations per site between any pair of sequences chosen at random from the population. Genetic distances were estimated for each genomic region by Kimura's two-parameter method (Kimura, 1980) in MEGA version 4 (Tamura et al., 2007) in the different host-begomovirus combinations at 15, 30, and 45 dpi. Standard error was calculated by the bootstrap method with 1,000 repeats (Nei and Kumar, 2000).

Mutant spectra heterogeneity was determined using the normalized Shannon entropy (SE) calculation according to the formula − [Σi (pi x lnpi)]/lnN], where pi is the frequency of each sequence in the mutant spectrum and N is the total number of sequences compared (Volkenstein, 1994). SE values range from 0 (all sequences are identical) to 1 (all sequences are different).

To analyze the effect of selection on mutant spectra the rate of substitutions per synonymous (d_S) and nonsynonymous (d_NS) site was estimated using the Pamilo-Bianchi-Li method (Pamilo and Bianchi, 1993) implemented in MEGA. The structure-genetic (SG) matrix of Feng et al. (1984) was used to obtain the acceptability values of amino acid changes found in each genomic coding region. The values range from 0 (drastic amino acid changes) to 6 (synonymous replacements) and take into account structural similarities and probabilities of amino acid changes. The probability of occurrence of amino acid replacements found in quasispecies was calculated according to the PAM-250 substitution matrix (Feng and Doolittle, 1996).

Mutations in overlapping genes were only considered for the computation of mutation types or the estimation of d, mutation frequency or Shannon entropy of each ORF but not of the entire sequenced zone. Also, since the Rep and C4 ORFs are not in frame, neither V2 and CP, the resulting amino acid changes in the overlapping sequence were considered for the d_NS/d_S calculation.

The statistical significance of differences in ORF mutation frequencies between samples was evaluated by a generalized linear model with logit link function assuming a binary distribution, and followed by post hoc analysis of pairwise comparisons with least square (LSD) correction to test for differences between treatments, using IBM SPSS Statistics v. 22 software (P < 0.05 was considered statistically significant). To check for a random distribution of mutations in the overlapping Rep and C4 a Wald–Wolfowitz runs test was conducted.

To determine absolute fitness the accumulation ability of a virus genotype in a given host, that is, the number of viral genomes produced after a given time (per unit of total DNA extracted from the plant tissue) in apical leaves was measured. Absolute fitness was approximated as the Malthusian growth rate per day, m, calculated according to the formula m=1t logQ, where Q is the number of pg of begomovirus ssDNA per 100 ng of total plant DNA (Lalić et al., 2011), which is estimated from the qPCR determination of the accumulated viral load in the apical part of each plant being monitored during the 15 days that have elapsed between each collection time (15, 30, and 45 dpi).

Quantification of viral DNA in test plants was carried out on total DNA extracts previously analyzed by agarose gel electrophoresis and UV-spectrophotometry (ND-1,000 spectrophotometer, NanoDrop, Fisher Thermo). Quantification was done in triplicate with a final volume of 20 μl with SYBR Premix Ex Taq Perfect Real Time (Takara), 0.125 μM of specific primers and 10 ng of total DNA. An amplicon of 169 bp within the Rep ORF of TYLCV and TYLCMaV was amplified with primers MA1508 and MA1509 (Table 1). In the case of TYLCSV, a fragment of 240 bp also in Rep was amplified with primers MA1506 and MA1507. PCR amplification was performed in an iCycler PCR System (BioRad) and consisted of an initial denaturation step at 95°C for 30 s, followed by 40 cycles at 95°C for 5 s and 60°C for 34 s. Specificity of the amplified products was assessed by melting curve analysis. Data from total DNA plant extracts were normalized to the 25S ribosomal RNA gene amplified using primers 25S rRNA UNIV (+ and –) (Mason et al., 2008; Table 1). To obtain the total number of ssDNA molecules, values were multiplied by two as each standard dsDNA molecule (a phagemid harbors the infectious clone) contains two ssDNAs (virion sense and complementary sense).

For the absolute quantification of TYLCV-like viruses, standard curves were prepared using pBluescript II SK (+) containing full-length viral genomes. Phagemids were serially diluted tenfold in 5 ng/μl of DNA extracted from either non-inoculated tomato, common bean or S. nigrum plants, to obtain 10² to 10⁹ viral genome copies per μl. For each PCR system (primer pair and background host DNA), standard curves were obtained by linear regression analysis of the threshold cycle (Cq) values of three standard dilution replicates over the Log of the number of genome copies present in each sample.

A generalized linear model was fitted to the Malthusian growth rate and the log value of viral load data. Post hoc analysis of pairwise comparisons via Bonferroni correction were used to test for differences between hosts at 45 dpi (P < 5.56 × 10⁻³) or along the time course of viral infections (P < 0.010 or 6.25 × 10⁻³).

To evaluate how the host affects the behavior of the different begomoviruses first we measured the viral load of the three begomoviruses in the four hosts. Under the conditions studied, no effective amplification was possible from young apical non-inoculated leaves for TYLCSV in common bean or for TYLCV in S. nigrum. As shown in Figure 2, it was clear that the presence of the Ty-1 resistance allele repressed accumulation of TYLCSV, and TYLCV to some extent, whereas no such repression was observed for TYLCMaV. Significantly higher TYLCV accumulation was observed in common bean at 45 dpi than in susceptible (510-fold; P = 0.005) or resistant (1 × 10³-fold; P = 0.000) tomato plants. Also, significantly more TYLCSV molecules accumulated in the wild reservoir S. nigrum at 45 dpi than in susceptible (60-fold, P = 0.005) or resistant (3.6 × 10⁶-fold, P = 0.0004) tomatoes (Figure 2). Thus, TYLCV and TYLCSV infections in common bean and S. nigrum, respectively, resulted in the generation of larger virus populations suggesting higher viral fitness compared to tomato. Our results also suggest better adaptation of TYLCMaV to common bean (with an accumulation level similar to TYLCV) than to S. nigrum (much lower accumulation than TYLCSV).

To better understand the effect of the different hosts on the three begomoviruses absolute fitness was assessed. Since we used single-sequence variants as inoculum the Malthusian growth rate (m) during viral infection time-courses was calculated as an estimate for absolute fitness (Lalić et al., 2011). The results shown in Figure 3 indicate differences in fitness across hosts during time-course infections. The presence of the Ty-1 allele significantly restricted the growth rate of TYLCSV (P = 0.004) in apical leaves at 45 dpi, and of both TYLCSV and TYLCV during the time-course infection (P = 0.000). However, as previously shown (Monci et al., 2002), the presence of the Ty-1 resistance allele in tomato had no evident effect on TYLCMaV fitness at 45 dpi (P > 5.56 × 10⁻³) or during the infection (P > 6.25 × 10⁻³). Despite being a recombinant of the former two viruses, TYLCMaV maintained similar fitness throughout the whole experiment in all four hosts (P > 6.25 × 10⁻³ or 0.010). Moreover, viral fitness in common bean and in S. nigrum was sustained along the experiment (P >0.010). In summary, these results suggest better performance of TYLCV and TYLCSV in common bean and S. nigrum, respectively, whereas TYLCMaV accumulates similarly well in all hosts, including the resistant tomato.

Next we analyzed the quasispecies complexity and heterogeneity of the three begomoviruses during the time-course infections. Specifically, viral quasispecies complexity refers to the composition of the mutant spectrum, which can be described both by mutation frequency and average genetic distance (d), whereas viral quasispecies heterogeneity refers to the proportion of different viral genomes in the mutant spectrum, as measured by Shannon entropy (Domingo et al., 2006). DNA samples from one infected plant per virus/host combination at each time point, plus a replicate sample of each combination at 45 dpi, were subjected to rolling circle amplification (RCA) to specifically amplify circular DNA using the high-fidelity φ29 DNA polymerase with a reported error rate of 3–6 × 10⁻⁶ (Esteban et al., 1993; Nelson et al., 2002). As a measure to reduce amplification bias we amplified 3–7 × 10⁸ input viral molecules. Sequences obtained for each begomovirus region were compared to the consensus sequence to determine mutation frequencies. Despite large variation among viruses, the highest mutation frequency values calculated for the whole sequenced region were detected in common bean and in S. nigrum (Figure 4A). When analyzed by genomic region, the mutation frequency of the quasispecies varied from <0.616 × 10⁻⁴ to 8.91 × 10⁻⁴ mutations/nt (mut/nt) (Table 2). The most complex region was V2 (e.g., 8.91 × 10⁻⁴ for TYLCV in resistant tomato, 8.6 × 10⁻⁴ and 8.0 × 10⁻⁴ mut/nt in common bean for TYLCV and TYLCMaV, respectively, or 8.4 × 10⁻⁴ mut/nt in S. nigrum for TYLCSV). In contrast, the least complex mutant spectra were found in the Rep and C4 genomic regions with values generally below 3 × 10⁻⁴ mut/nt (Table 2). Higher complexity values were found in common bean or S. nigrum than in tomato for CP or IR in all viral species. Thus, CP mutation frequency of TYLCV in common bean was higher than in resistant tomato at 15 dpi (P = 0.025), and also higher at 30 dpi than in susceptible (P = 0.034) and resistant tomato (P = 0.034). Only in TYLCV, CP (P = 0.030) and V2 (P = 0.045) regions showed more complexity in resistant tomato than in common bean at 45 dpi. Additionally, our results suggest that CP and IR were differentially enriched in mutations respect to other regions within a particular quasispecies. In S. nigrum at 45 dpi the complexity of CP (P = 0.001) and IR (P = 0.045) of TYLCSV was higher than that of both Rep and C4 regions. Also in S. nigrum CP region of TYLCMaV was more complex than Rep, C4, and V2 at 15 (P = 0.045) and 30 dpi (P = 0.005). In susceptible tomato at 45 dpi the mutation frequency was higher in CP than in C4 (P = 0.025) for TYLCSV and in resistant tomato CP was more complex than Rep (P = 0.001) and C4 (P = 0.014) for TYLCV.

Quasispecies genetic complexity was also assessed by analyzing average genetic distance (d) values. As shown in Table 3, d values varied from 0.00012 to 0.00173. The maximum genetic distance values for the four coding regions and the IR were found in the V2 region (especially for common bean and S. nigrum). High d values were also observed for C4 in some cases (e.g., in resistant tomato for TYLCSV at 15 dpi or in common bean for TYLCV at 30 dpi).

Begomovirus mutant spectra heterogeneity was evaluated by calculating the normalized Shannon entropy of viral quasispecies. Shannon entropy values calculated for the whole sequenced region were highest in common bean at 15 dpi and no difference between hosts were detected a later times post infection (Figure 4B). Susceptible tomato begomovirus quasispecies were the least heterogeneous, followed by those of resistant tomato and S. nigrum. Analyzed by genomic region, as shown in Table 2, Shannon entropy values ranged from 0 to 0.340. The CP region showed the highest heterogeneity in all four hosts and viruses, although this was especially pronounced in common bean and S. nigrum. Conversely, the lowest Shannon entropy values in the four hosts were detected in the C4 genomic region, except for TYLCSV at 15 dpi in resistant tomato. Therefore, heterogeneity of quasispecies was high in all host/begomovirus combinations, with the most heterogeneous populations exhibited by common bean and S. nigrum. Notably, the wild host S. nigrum contained the most heterogeneous populations at later infection time points for the viruses that were able to infect it. In general, mutation frequency, d and Shannon entropy values at 45 dpi were similar for replicates 1 and 2 for each virus/host combination.

Altogether, our results support the idea that begomovirus infections are characterized by the generation of highly complex and heterogeneous ssDNA mutant spectra in the host species tested, especially in common bean and the wild reservoir S. nigrum. V2 region supported the most complex spectra, CP was the most heterogeneous, while C4 exhibited both the least complexity and least heterogeneity. Interestingly, despite the high genetic complexity and heterogeneity of TYLCD-associated begomovirus quasispecies, no changes were detected in their whole genomic consensus sequences in any of the four hosts at any time point during infection.

Mutation data, comprising 128 base substitutions and 16 indels (complete list shown in Table 4), were pooled together irrespective of sampling time, virus or host. Mutation types and their frequency (%) relative to the number of changes found in a given genomic region or the whole sequenced zone are indicated in Table 5. Independent analysis of Rep, C4, IR, V2, and CP regions, shows that the most abundant mutations were base substitutions (86.7–100%), with indels only absent in C4. The IR contained the highest fraction of insertions, representing 11.1% of all mutations in the region. On the other hand, the highest fraction of deletions was detected in Rep (5.3% of total mutations). Of the base substitutions, the most common transition was C → T, which was highly enriched in C4 and V2 regions. G → A was the next most common transition, although it was less common in the IR and absent from C4. With regards to transversions, G → T was the most abundant. Conversely, A → T and A → C and G → T transversions were exceptionally frequent in the IR compared to the other four regions. Interestingly, the transition:transversion ratio was close to 1 in all regions. Transition:transversion ratios in DNA-based organisms (Begun et al., 2007; Hodgkinson and Eyre-Walker, 2010) as well as in RNA and dsDNA viruses (Duchêne et al., 2015) are over 2. While it has not been shown for ssDNA viruses, lower ratios might suggest that in all the begomovirus regions analyzed, especially in the IR and CP (with ratios of 0.9), either negative selection is removing transitions or transversions are favored. Nonsynonymous changes were more frequent than synonymous mutations in all regions, ranging from 65.4 to 79.2% of total mutations in Rep and V2, respectively.

Next, we studied the frequency and acceptability of the amino acid changes found in the different gene products under study. Firstly, amino acid changes were analyzed using the SG matrix as described previously (Feng et al., 1984) generating values from 0 to 6, with 0 indicating the most drastic amino acid changes and 6 the least, and their relative frequencies shown in Figure 5A. The results show that amino acid changes in Rep were poorly tolerated, since 83.9% of changes have acceptability values above 3, with 38.7% having a value of 6 (synonymous changes). Although the overall degree of amino acid change was similar in all regions, the higher relative frequency of more pronounced changes (with values of 2 or less) found in the CP might indicate that this protein is more tolerant to drastic amino acid changes. We analyzed the probability of amino acid substitution occurrence using a PAM-250 substitution matrix (Feng and Doolittle, 1996). The results shown in Figure 5B indicate that the highest frequency of probable amino acid replacements occurred in Rep and that unlikely amino acid changes were most prevalent in the CP protein. Altogether the results indicate that the Rep protein tends to be more conserved and that the CP protein is more amenable to change.

Sequence variations explored during begomovirus infections seemed to be affected by specific virus-host interactions. Sequence space distribution of the mutations detected in our time-course assays is summarized in Figure 6, which shows the location of the mutations found in begomovirus quasispecies by assay time point, host and begomovirus (TYLCSV, TYLCV, or TYLCMaV). Considering the genomic regions explored at each sampling time, we observed that the tomato Ty-1 resistance allele did not significantly influence mutation frequency for any ORF, since no significant differences between susceptible tomato and resistant tomato were found at any time point (generalized linear model assuming binary distribution and using logit link function, P ≥ 0.05). However, differences between hosts were observed with respect to where mutations occurred over time. Thus, for example, mutation frequency of the CP increased for TYLCSV in susceptible tomato between 15 and 45 dpi (P = 0.048) and between 30 and 45 dpi (P = 0.003), and for TYLCV in resistant tomato between 15 and 45 dpi (P = 0.001). This increase in CP mutation frequency might be related to the decrease in fitness in tomato. In S. nigrum, TYLCMaV (replicates 1 and 2 at 45 dpi) showed a significant increase in the number of mutations in Rep between 15 and 45 dpi (P = 0.044) and between 30 and 45 dpi (P = 0.005) but not in the overlapping C4 ORF. A random distribution of mutations was found both in Rep and C4 (runs test, P > 0.05) so mutational hotspots were ruled out. This result suggests an active exploration of sequence space in Rep but not in C4 for the recombinant virus. No parallel mutations were found in replicate samples at 45 dpi. Our results suggest that begomoviruses explore different sequence space depending on host species, although further studies are needed to confirm this including additional replicates (at 15 and 30 dpi). They also suggest that the presence of the Ty-1 gene in tomato does not affect how begomoviruses move in sequence space at each time point. Furthermore, it is interesting to note that begomovirus consensus sequences were completely invariant throughout the infection time-courses regardless of mutant spectra or host, and despite the large number of mutations found.

To assess begomovirus mutant spectra adaptation to hosts, we studied the direction of selective forces on their coding regions. To this end the average rate of synonymous substitutions per synonymous site (d_S) and of nonsynonymous substitutions per nonsynonymous site (d_NS) for each mutant spectrum were estimated as reported previously (Pamilo and Bianchi, 1993). Where possible, the d_NS/d_S ratio was also calculated. The results in Table 6 show that d_NS/d_S ratios of begomovirus quasispecies evolving in tomato ranged from 0.149 to 0.616 in Rep, V2, and CP. Ratios calculated for TYLCSV and TYLCMaV were below 1, suggesting that all three regions were under negative selection, both in susceptible and resistant tomato. However, d_NS/d_S ratios ranged from 0.357 to 1.842 in common bean and from 0.464 to 1.778 in S. nigrum, showing that, in these hosts, selective forces acted differentially depending on the coding region. In fact, we observed positive selection (d_NS/d_S > 1) in the CP region of both common bean, at 15 and 45 dpi for TYLCV and TYLCMaV infections, respectively, and S. nigrum, at 45 and 15 dpi for TYLCSV and TYLCMaV infections, respectively. The results on d_NS/d_S ratios of replicates 1 and 2 at 45 dpi of each virus/host combination were similar, although additional replicates would add more robustness to our conclusions. Altogether, these results suggest differing selective constraints of begomovirus quasispecies in different hosts. Also, it seems that Ty-1 resistance does not alter the way begomoviruses are selected in tomato, as purifying selection dominates in all the coding regions analyzed. Our data also support positive selection of virus variants in the CP region in common bean and S. nigrum, whereas Rep, C4, and V2 coding regions are under purifying selection. Positive selection of CP region variants in common bean and S. nigrum together with higher virus fitness and diversity, and differential mutation distribution occur without alterations to begomovirus consensus sequences. These results suggest that the mutant spectrum plays an important role in begomovirus infections in each host.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.