Role of host genetic diversity for susceptibility-to-infection in the evolution of virulence of a plant virus.

Predicting viral emergence is difficult due to the stochastic nature of the underlying processes and the many factors that govern pathogen evolution. Environmental factors affecting the host, the pathogen and the interaction between both are key in emergence. In particular, infectious disease dynamics are affected by spatiotemporal heterogeneity in their environments. A broad knowledge of these factors will allow better estimating where and when viral emergence is more likely to occur. Here, we investigate how the population structure for susceptibility-to-infection genes of the plant Arabidopsis thaliana shapes the evolution of Turnip mosaic virus (TuMV). For doing so we have evolved TuMV lineages in two radically different host population structures: (1) a metapopulation subdivided into six demes (subpopulations); each one being composed of individuals from only one of six possible A. thaliana ecotypes and (2) a well-mixed population constituted by equal number of plants from the same six A. thaliana ecotypes. These two populations were evolved for twelve serial passages. At the end of the experimental evolution, we found faster adaptation of TuMV to each ecotype in the metapopulation than in the well-mixed heterogeneous host populations. However, viruses evolved in well-mixed populations were more pathogenic and infectious than viruses evolved in the metapopulation. Furthermore, the viruses evolved in the demes showed stronger signatures of local specialization than viruses evolved in the well-mixed populations. These results illustrate how the genetic diversity of hosts in an experimental ecosystem favors the evolution of virulence of a pathogen.

1. Introduction

Since the term ‘emerging infectious disease’ was coined in the mid-1900s, its definition has evolved (Rosenthal et al. 2015). Woolhouse and Dye (2001) enunciated the most comprehensive definition of an emerging infectious disease as one ‘whose incidence is increasing following its first introduction into a new host population or whose incidence is increasing in an existing host population as a result of long-term changes in its underlying epidemiology’. Engering, Hogerwerf, and Slingenbergh (2013) rephrased this definition and suggested that emergence events can be classified into three groups: (1) viruses showing up in a novel host, (2) mutant viral strains displaying novel traits in the same host, and (3) already known viral diseases spreading out in a new geographic area. Following these definitions, viral emergence is usually associated with cross-species transmission but it can actually occur with or without species jump (Di Giallonardo and Holmes 2015).

The emergence and re-emergence of viruses cause serious threatening to public health (Morens and Fauci 2013) and are responsible for large yield losses in crops that compromise food security (Vurro, Bonciani, and Vannacci 2010). Gaining knowledge on the principles that govern viral emergence would allow to predict when and where such events are more likely to happen. Achieving an accurate prediction of the emergence of viral diseases is a challenging task, because emergence is governed by multiple and diverse factors that remain poorly understood or completely unknown. Predictions will become even more difficult under the ongoing climate change that will favor the conditions for development and dispersal of the virus vectors (Garrett et al. 2006; Baylis 2017).

The spectrum of disease severity can be attributed to heterogeneity in virus virulence or in host factors; the two are not necessarily independent explanations and they may actually complement and/or interact with each other. A problem faced by viruses is that host populations consist of individuals that had different degrees of susceptibility to infection (Schmid-Hempel and Koella 1994; Pfennig 2001). Therefore, adaptive changes improving viral fitness in one host may be selected against, or be neutral, in an alternative one. Genetic variability in susceptibility of hosts and infectiousness of viruses have been well studied in animals and plants (Schmid-Hempel and Koella 1994; Altizer 2006; Hughes and Boomsma 2006; Brown and Tellier 2011; Anttila et al. 2015; Parrat, Numminen, and Laine 2016). Variability within host populations can arise from nonrandom spatial distributions of genotypes: social groups of animals are more closely related to each other than to other members of the population, and plant crops are generally cultivated as genetically homogeneous plots. These situations facilitate virus transmission among genetically similar host genotypes. Spatial structure and local migration predict evolution of less aggressive exploitation in horizontally transmitted parasites (reviewed in Parrat, Numminen, and Laine 2016). For example, Boots and Mealor (2007) observed that Plodia interpunctella granulosis virus (PiGV) evolved in spatially structured lepidopteran host populations become less virulent than the one maintained in well-mixed host populations. More recently, Berngruber, Lion, and Gandon (2015) showed that a latent λ bacteriophage won competitions against a virulent one in a spatially structured host populations but lost in well-mixed populations. In a natural context, it has been observed that virulent Linum marginale fungi were more frequent in highly resistant Melampsora lini plant populations whereas avirulent pathogens dominated susceptible populations (Thrall and Burdon 2003; Thrall et al. 2012). Similar results were observed in laboratory evolution experiments with the pathosystem Arabidopsis thaliana–Tobacco etch virus (TEV): a negative association between plant natural accessions permissiveness to infection and TEV virulence was evolved (Hillung et al. 2014). Finally, host population structure also promotes coexistence of hosts and parasites by creating refugia due to limited dispersal of viruses; increased dispersal makes coexistence less stable (Brockhurst, Buckling, and Rainey 2006).

The role of host population heterogeneity has also received quite a lot of attention from theoreticians (Comins, Hassell, and May 1992; Gandon et al. 1996; Boots and Sasaki 1999, 2002; Gandon and Michalakis 2002; Tellier and Brown 2011), resulting in a number of interesting predictions that in many instances have not been properly tested experimentally. One of the most tantalizing predictions is that in the absence of host heterogeneity, parasites must evolve toward a host exploitation strategy that maximizes transmission with low virulence (Haraguchi and Sasaki 2000; Regoes, Nowak, and Bonhoeffer 2000; Rodríguez and Torres 2001; Ganusov, Bergstrom, and Antia 2002; Gandon 2004; Lively 2010; Moreno-Gámez, Stephan, and Tellier 2013). However, in the context of emerging diseases, right after the spill-over of the new pathogen into the heterogeneous host population and prior to adaptive evolution to take place, Yates, Antia, and Regoes (2006) showed that host heterogeneity has a small effect in the probability of establishing the disease. Very recently, Chabas et al. (2018) showed that evolutionary emergence is more likely to occur when the host population contains intermediate levels of resistant hosts, confirming this prediction using different phages and bacterial hosts with different alterations in the CRISPR/Cas immune systems that conferred the cells with different levels of resistance.

In this study, we use experimental evolution to explore the effect of host population structure for genes involved in susceptibility to infection in the evolution of infectiousness and virulence of a plant virus; in particular, we want to explore the extent to which host heterogeneity determines the rate of evolution of these two fitness-related traits. According to the above theoretical predictions and experimental observations in other pathosystems, we expect the virus infectiousness and virulence to evolve to lower levels in genetically homogeneous host subpopulations, but at a faster rate, than in genetically diverse well-mixed host populations. The pathosystem we have studied is composed of Turnip mosaic virus (TuMV; genus Potyvirus, family Potyviridae) as pathogen and six different natural accessions (hereafter ecotypes) of A. thaliana that differ in their susceptibility to infection as hosts.

2. Methods

2.1 Selection of A. thaliana ecotypes for evolution experiments

Before we could begin the evolution experiment, we sought to identify a set of A. thaliana ecotypes representative of the phenotypic variability in response to TuMV infection (Rubio et al. 2019) observed for a larger collection of ecotypes. Sixteen A. thaliana ecotypes were evaluated. Based on the similarity of their phenotypic responses to infection with TuMV (disease intensity over time measured as the area under the disease progression steps (AUDPS) curve; see Section 2.4 for a definition and explanation) during 14 days post-inoculation (dpi), ecotypes were clustered as shown in Fig. 1 (UPGMA). Six accessions were chosen as representative of the five clusters shown in Fig. 1: Col-0, Ga-0, Gy-0, Oy-0, Ta-0, and Wt-1. A progressive k-clustering algorithm confirmed that adding additional clusters did not result in significant improvement in model predictability (five vs six clusters partial-F test: F1,9 = 2.841, P = 0.058). Ta-0 showed the least intense disease progression (AUDPS = 0.490) while Ga-0 the most intense one (AUDPS = 6.500); Oy-0 and Wt-1 showed similar progression (AUDPS = 1.709).

Figure 1.

UPGMA clustering of A. thaliana ecotypes according to their response to TuMV infection. The six ecotypes selected for the evolution experiment are highlighted in gray. As a measure of virulence, the AUDPS of the ancestral TuMV isolate on each selected ecotype is indicated in the left scale.

Furthermore, the six ecotypes reached growth stage 3.5 in the Boyes’ scale (Boyes 2001) at the same time after germination (±1 day), and at that moment they were inoculated. This synchronization ensures that they all were at the same phenological state when inoculated.

In all experiments performed in this study, plants were maintained in a BSL-2 growing chamber at 8 h light:16 h dark cycles and temperature variation of 24 °C day:20 °C night.

2.2 TuMV original isolate and inocula

Infectious saps were obtained from TuMV-infected Nicotiana benthamiana plants inoculated with an infectious plasmid containing TuMV genome cDNA (GenBank accession AF530055.2) under the control of the Cauliflower mosaic virus 35S promoter. This TuMV sequence variant corresponds to YC5 isolate from calla lily (Zantedeschia sp.) (Chen et al. 2003). The same stock of plasmid was used to inoculate two batches of five N. benthamiana plants each, with a year of difference. After plants showed visible symptoms of infection, two independent infectious saps (viral stocks) were obtained by grinding infected tissues from N. benthamiana plants in a mortar with ten volumes of grinding buffer (50 mM KH2PO4 pH 7.0, 3% polyethylene glycol 6000).

2.3 Experimental evolution

TuMV lineages were evolved during twelve consecutive serial passages as schematized in Fig. 2 in two treatments that differ in the composition of the host population. The first treatment consists of a host metapopulation structured in six genetically homogeneous demes. TuMV was evolved in this genetically homogeneous host demes. Passages were made by harvesting the symptomatic plants at 14 dpi, preparing infectious sap as described above and inoculating the virus to the next demes (ten plants each) of the same ecotype. Three leaves from 21-day old plants were rub-inoculated each with 5 μl of infectious sap and 10 per cent Carborundum (100 mg/ml). One TuMV lineage was evolved on each one of the six ecotypes chosen (Fig. 1). In this treatment, each TuMV lineage in the metapopulation only experienced a particular host genotype along their evolution. In this treatment, we expect evolution to be dominated by rapid local adaptation.

Figure 2.

Schematic representation of a passage during the evolution experiments. A total of twelve such passages were done, in each of two replications conducted at different dates. See Section 2.3 for a detailed description. Arabidopsis thaliana draws were adapted from https://figshare.com/articles/Arabidopsis_Rosette_drawing_steps/4688839.

The second treatment consists of a host population without subpopulation structure composed of ten well-mixed plants from each of the six ecotypes. TuMV was thus evolved in a host population with maximal genetic heterogeneity. In this case, the passages were made harvesting all the symptomatic plants of all the ecotypes, pooling the tissues, preparing sap and transmitting it again to ten plants from each of the ecotypes (Fig. 2). In this treatment, the evolving population of symptomatic TuMV strains has an equal opportunity of infecting each plant ecotype at each passage. More susceptible ecotypes will contribute more to the next viral generation because more plants are infected and support more generations of replication, while more restrictive ones will contribute in a minor amount to the viral population. In this treatment, we expect a slower adaptation due to fluctuating selection and the evolution of generalist viruses.

This evolution experiment was replicated in two fully independent blocks, hereafter, referred as P and R, respectively. Each block was initiated with a different viral stock (as described in Section 2.2), different light sources (P with fluorescent tubes at PAR 100–150 μmol/m2/s and R with LED tubes at PAR 90–100 μmol/m2/s), inoculations were done by two different researchers, and started with more than year of difference in time (P started 2 December 2016 and R started 24 January 2018). In addition, the substrate mix used in passage 8 of block P was not the A. thaliana standard one in some of the pots but the one used for growing tomato plants, as a consequence A. thaliana plants in these pots grew slightly smaller.

2.4 Evaluation of infectiousness and virulence

Upon infection, plants were observed every day and the number of plants showing symptoms was recorded. Infectiousness, I, was evaluated as the frequency of plants showing symptoms at 14 dpi out of the ten plants inoculated with a standardized amount of infectious sap. Using the frequency time-series within the duration of a passage, the AUDPS (Simko and Piepho 2012) was evaluated as a proxy to virulence. AUDPS represents the intensity at which symptoms appear in a population of inoculated plants, and in our case, it is bounded between zero (no plant shows symptoms 14 dpi) and ten (all plants show symptoms at 1 dpi). In the case of the heterogenous host population treatment, I and AUDPS were evaluated in each of the ecotypes. I data were probit-transformed as hence the new variable f is now Gaussian distributed with mean zero and variance 1. For visualizing infectiousness data in a meaningful scale, I will be presented in figures while the f transformed data will be used in the statistical analyses described in Section 2.5.

At each passage, infection status of each plant was assessed by the presence of symptoms rather than by molecular detection methods. In our extensive experience with the TuMV/A. thaliana pathosystem, there is an almost one-to-one match between infection and the development of symptoms. Symptoms started with leaf curling and vein bleaching (∼5 to 6 dpi) that quickly developed to diverse grades of leaf chlorosis and/or necrosis (∼10 to 12 dpi). Plants also suffered a developmental arrest, with deformed new leaves, abortion of flowering button and abnormal growth of the caulinar apex.

2.5 Data analyses

AUDPS and f data were fitted to a fully factorial multivariate analysis of covariance model (MANCOVA), in which experimental block (B; P and R), population structure (D; homogeneous or heterogeneous), and plant ecotype (E) were treated as orthogonal factors and passage (t) as a covariable. Main effects and all the interactions between factors and the covariable were incorporated into the model. The full model equation thus reads where X = (AUDPS, f)T is the vector of phenotypic traits observed at time i, experimental block j, population structure k, and ecotype l, μ represents the vector of phenotypic grand mean values and ε stands for the vector of errors. In addition, univariate ANCOVA analyses were performed for AUDPS and f using the same model equation.

The magnitude of effects was evaluated using the statistic (proportion of total variability in the traits attributable to each factor in the model). Conventionally, values < 0.05 are considered as small, 0.05 ≤ < 0.15 as medium and ≥ 0.15 as large effects.

Unless otherwise mentioned, all statistical analyses described in this work were performed using SPSS version 25 software (IBM, Armonk, NY).

2.6 Evaluation of rates of phenotypic evolution

AUDPS and I data obtained for each lineage were fitted to the following first-order autoregressive integrated moving-average, ARIMA(1, 0, 0), model:

, where Y represents the variable being analyzed at passage k, measures the degree of self-similarity in the time-series data (correlation between values at passages t and t–1), represents the sampling error at passage t, and β represents the linear dependency of variable Y with passage number, that is, the rate of phenotypic evolution.

To explore the effect of factors D and E in the rates of AUDPS and I evolution, the following generalized linear model (GLM) was fitted to the β values: , were superscript Y again refers to the trait being analyzed, the grand mean value for the rate of evolution of trait Y and all other terms are defined in Section 2.5. A Normal distribution and identity link function were chosen based on the minimal BIC value among competing models. Notice that the error term is obtained from the differences in the estimates from both experimental blocks.

2.7 Infection matrices

To analyze the specificity of adaptation of each evolved TuMV lineage, we performed a full cross-infection experiment in which all the fourteen evolved lineages were inoculated into ten plants of all six ecotypes. In the case of the lineages evolved in the well-mixed host population, the virus isolated at passage 11 was inoculated into ten additional plants of each of the six ecotypes to separate it into sub-samples. Infection matrices were analyzed using tools borrowed from the field of community ecology to explore whether they show random associations between viral lineages and host genotypes, one-to-one associations, nestedness indicative of a gene-for-gene type of interaction, or modularity (Weitz et al. 2013). The statistical properties of the resulting infection matrices were evaluated using the bipartite version 2.11 package (Dormann, Gruber, and Fruend 2008) in R version 3.3.2 (R Core Team 2016). Three different summary statistics were evaluated: nestedness (Bascompte et al. 2003), modularity (Newman 2006), and overall specialization d′ index (Blüthgen, Menzel, and Blüthgen 2006). d′ is based in Kullback–Leibler relative entropy, that measures variation within networks and quantifies the degree of specialization of elements within the interaction network. Statistical significance of nestedness and modularity was evaluated using Bascompte et al. (2003) null model.

3. Results

3.1 Evolutionary dynamics of virulence and infectiousness

Figure 3 illustrates the evolutionary dynamics of virulence, measured as AUDPS, observed on each of the host ecotypes and for both population structures. Likewise, Fig. 4 shows the evolution of infectiousness, I, along the passages. In both figures black symbols and lines correspond to the results of experimental block P while red symbols and lines correspond to experimental block R. Both variables are significantly correlated (Pearson’s partial correlation coefficient controlling for B, D, E, and t: r = 0.640, 306 d.f., P < 0.001) and hence multivariate methods were used to analyze the data while increasing the power of the tests. Both datasets were fitted to the MANCOVA model equation described in Section 2.5. In all cases, an overall trend to increase virulence and infectivity along time can be observed in the time-series data (Figs 3 and 4). This trend is statistically supported by a significant net effect of the passage (covariable t) on both phenotypic traits (Table 1; P < 0.001).

Figure 3.

Evolution of virulence (AUDPS). Black circles and lines represent the results from experiment P and red triangles and lines from experiment R. Solid symbols represent evolution in the corresponding deme of host metapopulation while open symbols represent evolution in the well-mixed host populations.

Figure 4.

Evolution of infectiousness (I). Black circles and lines represent the results from experiment P and red triangles and lines from experiment R. Solid symbols represent evolution in demes within a metapopulation while open symbols represent evolution in the well-mixed host populations. Notice that infectiousness data were probit-transformed for statistical analyses.

Table 1.

MANCOVA analysis of the symptoms development (AUDPS) and infectiousness (f) data (Figs 3 and 4). See Section 2.5 for a description of the model equation and parameters. 1−β is the power of the corresponding test.

Source of variationa	Wilk’s Λ	F	d.f.	P	ηP2	1−β
Intersection (μ)	0.244	406.426	2,263	<0.001	0.756	1
t	0.353	241.176	2,263	<0.001	0.647	1
B	0.787	35.600	2,263	<0.001	0.213	1
B×t	0.911	12.805	2,263	<0.001	0.089	0.997
D	0.943	7.905	2,263	<0.001	0.057	0.952
D×t	0.972	3.816	2,263	0.023	0.028	0.690
E	0.713	9.678	10,526	<0.001	0.155	1
E×t	0.947	1.445	10,526	0.157	0.027	0.731
B×D	0.993	0.969	2,263	0.381	0.007	0.218
B×D×t	0.985	1.979	2,263	0.140	0.015	0.407
B×E	0.911	2.524	10,526	0.006	0.046	0.954
B×E×t	0.941	1.631	10,526	0.094	0.030	0.794
D×E	0.991	0.225	10,526	0.994	0.004	0.130
D×E×t	0.990	0.275	10,526	0.986	0.005	0.152
B×D×E	0.985	0.402	10,526	0.946	0.008	0.212
B×D×E×t	0.983	0.464	10,526	0.913	0.016	0.316

Model factors: passage ≡ t ∈ (0, …, 12), block ≡ B ∈ (P, R), population structure ≡ D ∈ (metapopulation, well-mixed population), host ecotype ≡ E ∈ (Col-0, Ga-0, Gy-0, Oy-0, Ta-0, Wt-1).

Next, given the differences in starting inocula, light conditions and experimenter responsible for performing each experimental block, we expect a priori to find significant block effects, either net or in combination with other factors in the model. Table 1 shows that experimental block has a net (B) effect on the phenotypic vector (P < 0.001). This effect changes along evolution passages (B × t) (P < 0.001) and depends on the particular ecotype (B × E) (P = 0.006). The observed net block effect is of large magnitude ( = 0.647), while the effect of its interaction with passage number ( = 0.089) and ecotype ( = 0.046) can be considered of small–medium size.

Most interestingly, the population structure in which the TuMV lineages had evolved has a highly significant net (D) effect on the phenotypic vector (Table 1; P < 0.001), though it is of small–medium size in the multivariate analysis ( = 0.057), and of similar effect for both AUDPS and f (univariate ANCOVAs shown in Supplementary Table S1; = 0.042 and = 0.051, respectively). As it can be seen in Figs 3 and 4, this effect comes from a pattern (consistent across both experiments) that AUDPS and I values are usually larger for the lineages evolving in the well-mixed population during the early passages of evolution than for those evolved in the metapopulation (solid symbols are above open ones). Indeed, on average AUDPS was 7.42 per cent larger for TuMV lineages evolved in well-mixed population than in metapopulation. Likewise, f was 12.25 per cent higher for TuMV lineages evolved in the well-mixed plant population than in the genetically homogeneous demes. Furthermore, the effect of population structure changed along passages (D × t) (Table 1; P = 0.023), though this effect is of small size in the multivariate analysis ( = 0.023) and also in the univariate ones (Supplementary Table S1; = 0.016 and = 0.028, respectively).

The host ecotype in which lineages evolved has a highly significant effect on the magnitude of the phenotypic vector (Table 1; P < 0.001), the size of this effect being large ( = 0.155). Indeed, the univariate analyses show that the effect is much larger for AUDPS than for f (Supplementary Table S1; = 0.233 and = 0.085, respectively). The TuMV evolved in Ga-0 was the most virulent (AUDPS = 9.492), followed by viruses evolved in Ta-0 and Gy-0, while the less virulent infection corresponds to a homogenous group formed by TuMV lineages evolved in Col-0, Wt-1, and Oy-0 (AUDPS in the range: 5.938–6.525) (sequential Bonferroni’s post hoc test, P ≤ 0.017). This ranking of ecotypes according to virulence slightly differs from the ranking observed for the ancestral TuMV isolate (Fig. 1). Likewise, the virus evolved in Ga-0 was also the most infectious (f = 1.264), followed by lineages evolved in Col-0 and Ta-0. The less infectious viruses were those evolved in Wt-1, Oy-0, and Gy-0 (f in the range: 0.905–0.972).

Therefore, we conclude in this first section that the presence of maximal host genetic diversity for genes involved in susceptibility to infection selects for more virulent and infectious viruses. Differences are mostly due to the dynamics of the initial passages. For both phenotypic traits studied, this effect depends on the presence of particular ecotypes.

3.2 Rates of phenotypic evolution

In the previous section, we found a net effect of passage on the virulence-related traits. More interestingly, this effect depended on the population structure (D) and was consistent across both experimental blocks (non-significant B × D × t effect; Table 1), which allows us to treat the estimates of rates of evolution from each block as replicates in a GLM analysis. Next, we sought to get a better understanding of the effect of population structure for susceptibility to infection on the rates of phenotypic evolution. To do so, we have estimated evolution rates as described in Section 2.6. Summary statistics for the ARIMA(1, 0, 0) model fitting are shown in Supplementary Table S2. Figure 5 compares the estimated rates of evolution for AUDPS and I for both population structures analyzed. Rates of evolution were fitted to the GLM described in Section 2.6. The results from these analyses are shown in Table 2. First, lets focus in the rate of AUDPS evolution. In all cases except for the lineages isolated from Wt-1, the rate of evolution for lineages evolved in the corresponding deme of the metapopulation is larger than for the corresponding lineages evolved in the well-mixed host population (Fig. 5A). This trend is statistically significant (Table 2; P = 0.015) and the size of effect must be considered as large ( = 0.219).

Figure 5.

Mean rates of phenotypic evolution for the two traits studied, AUDPS (A) and I (B). Rates of evolution were estimated from the ARIMA(1, 0, 0) model described in Section 2.6. Error bars represent ±1 SEM.

Table 2.

GLM analyses of the rates of evolution of AUDPS and I data (Fig. 5). See Section 2.6 for a description of the model equation and parameters. 1−β is the power of the corresponding test.

Source of variationa	Rate of evolution (βY)	χ 2	d.f.	P	ηP2	1−β
Intersection (β^Y)	AUDPS	80.518	1	<0.001	0.965	1
	I	25.879	1	<0.001	0.660	0.993
D	AUDPS	5.936	1	0.015	0.219	0.393
	I	4.979	1	0.026	0.187	0.334
E	AUDPS	8.835	5	0.116	0.308	0.259
	I	6.176	5	0.289	0.227	0.179
D×E	AUDPS	3.378	5	0.642	0.131	0.110
	I	1.071	5	0.975	0.044	0.067

Model factors: population structure ≡ D ∈ (metapopulation, well-mixed population), host ecotype ≡ E ∈ (Col-0, Ga-0, Gy-0, Oy-0, Ta-0, Wt-1).

Second, a similar result has been observed for the rates of I evolution (Fig. 5B): rates of infectiousness evolution are faster in each deme of the metapopulation than in the well-mixed population. Again, differences among population structures are significant (Table 2; P = 0.026) and of large effect ( = 0.187).

For both traits, no significant differences in rates of evolution exist between ecotypes (E) nor for the interaction between ecotypes and population structure (D × E) (Table 2), with the power of the tests being too low (1−β ≤ 0.259) as to fully discard we are not accepting a false null hypothesis of no effects. Low statistical power must clearly be due to the small sample size used (only two experimental blocks).

As a conclusion for this second section, evolution of virulence and infectiousness took place at a faster pace in genetically homogeneous demes of a metapopulations than in a well-mixed genetically heterogeneous population.

3.3 Comparison of infection matrices

Finally, we evaluated the degree of specificity of adaptation to the six ecotypes for all the evolved lineages. To build up infection matrices, TuMV lineages evolved in each of the ecotypes, or isolated from each ecotype in the case of the two blocks evolved in the well-mixed host population, were inoculated into each one of the six ecotypes and AUDPS was estimated as described above. Figure 6 shows a binary representation of the infection matrices estimated (averaging across experiments) for the two host population structures. Black squares represent host-virus combinations in which virulence was equal or greater than for the value observed for the viral lineage in its corresponding local host (host of isolation in the case of the well-mixed host population). White squares represent less virulent combinations. The first row in the matrices correspond to the most generalist TuMV lineage, in this case those evolved in Gy-0, whereas the last row represents the most specialist TuMV lineage, here those evolved in Oy-0. This observation is consistent for both matrices. Likewise, both matrices are also consistent regarding the rank order in susceptibility to infection of the different ecotypes: Oy-0 is the most susceptible ecotype, being infected by almost all viral lineages with similar virulence (exception being the virus isolated from Col-0 in the well-mixed host population) and Gy-0 being the most resistant ecotype, infected with high virulence only by the lineages evolved in or isolated from Gy-0. In other words, this suggests that more permissive ecotypes select for less virulent viruses while more restrictive ecotypes select for more virulent ones.

Figure 6.

Infection matrices obtained from the virulence data (mean AUDPS values of the two experiments). Black cells represent cases in which virulence was equal or greater than the value estimated for the corresponding viral lineage on its local host ecotype. In the case of the well-mixed population, it corresponds to the value estimated on the ecotype from which the virus was isolated in the last evolution passage. As described in Section 3.2, in each matrix evolved viruses are ordered from the most virulent (upper row) to the less virulent (bottom row) and host ecotypes from the most sensitive (most left column) to the most resistant (rightest column).

To further test this hypothesis, we evaluated the nestedness and modularity of the two matrices. The infection matrix estimated for the demes of the metapopulations shows a significant nestedness (Fig. 6, left; P = 0.035). This suggests that virus evolution in a single host genotype selects for a gene-for-gene interaction mechanism. This model of host-virus interaction is fully compatible with the above hypothesis. However, the matrix estimated for the well-mixed population did not show significant nestedness (Fig. 6, right; P = 0.149), suggesting that gene-for-gene interactions have not being selected under this ecological situation. Indeed, odds ratios indicate that the left matrix in Fig. 6 is 3.16 per cent more nested that the right one.

Both matrices show significant modularity (P = 0.022 for the metapopulation and P = 0.012 for the well-mixed population), though the modularity in the well-mixed host population matrix is slightly larger (odds ratio: 0.24%). A module is a group of viruses and hosts for which the viruses in the set are more likely to be virulent in these hosts that to any other host outside the group, and that hosts in the group are more likely to be infected with similar virulence by viruses from within the group. Such modules suggest common selective constraints imposed by the hosts and similar evolutionary solutions found by the viruses.

Finally, we sought to quantify the degree of specialization in the host ecotype—viral lineage, or partner diversity. To this end, we applied Blüthgen, Menzel, and Blüthgen (2006)d′ index. For the infection matrix estimated for viruses evolved in the metapopulation (Fig. 6 left), d′ = 0.00207, while for the matrix estimated for viruses evolved in the well-mixed population (Fig. 6 right), d′ = 0.00061. Therefore, 3.39 times greater specialization evolved when the virus was facing a single host genotype during the evolution experiment.

4. Discussion

4.1 Host population structure and the evolution of specialist and generalist viral strategies

Most plant viruses are generalists capable of infecting more than one host species and thus generalism versus specialism should be properly defined on the basis of variance in infectiousness across different hosts (Leggett et al. 2013). In this sense, a generalist virus will be characterized by a low variance in infectiousness over different host genotype and/or species; in contrast a specialist virus will have a high variance. It is logical to expect that virus evolution should proceed faster in homogenous than in heterogeneous host environments because viruses with a narrower host range (specialist) have higher probabilities of fixing beneficial alleles, taking less time to do so (Gavrilets and Gibson 2002; Whitlock 2002, 2003; Whitlock and Gomulkiewicz 2005; Papaïx et al. 2013). Consequently, viral species or genotypes with broader host ranges (generalists) must have slower rates of evolutionary response (Bennett, Lenski, and Mittler 1992; Fry 1996; Whitlock 1996; Kassen and Bell 1998; Kassen 2002). If evolution occurs in too many hosts, then the total selection for any particular host-specific trait would not be as effective, and viruses specializing into a single host would then evolve faster and outcompete their generalist counterparts. As a result of this faster evolutionary rate, specialist viruses may persist longer in time in a constant host landscape. In this study, we have directly challenged this hypothesis using experimental evolution. Lineages of TuMV, a prototypical RNA virus of the picorna-like superfamily (sensuKoonin et al. 2008), were evolved either in six alternative single host demes constituting a metapopulation or in a well-mixed genetically heterogeneous host environment composed by equal numbers of the same six host genotypes. Giving support to the above hypothesis, we observed that rates of phenotypic evolution were significantly faster for lineages evolved in the individual demes of the metapopulation than in the well-mixed population.

Pfennig (2001) made the distinction between polymorphism and polyphenism as causes of pathogens generalism. On one hand, polymorphism means that different strains of pathogens evolve specialized virulence strategies in different host genotypes. On the other hand, polyphenism means that pathogens facultatively express alternative virulence strategies depending on the host phenotypes. In case of viruses with compacted genomes and multifunctional proteins, polymorphisms seem a more plausible explanation, though we cannot rule out polyphenism. In this sense, highly polymorphic viral populations would behave as generalists owed to the diversity of specialists they contain.

4.2 Host population structure, selection of recognition mechanisms and evolution of virulence

We have also characterized the evolution of two-virulence-related traits, AUDPS and infectiousness. Before engaging ourselves in further discussion, it is worth mentioning that our definition of virulence, AUDPS, is to some extent different from the classical one used by evolutionary ecologists, namely the reduction in host fitness due to infection, which is the one generally used in most theoretical treatments of the evolution of virulence in infectious diseases. However, our definition is more akin with definitions used by plant pathologists: development of symptoms (AUDPS) or susceptibility to infection (infectiousness). Given that TuMV is a castrating pathogen, and hence infection has a strong impact in the plant fitness, AUDPS turns out to be a good predictor of the reduction in plant fitness.

According to theoretical predictions (Haraguchi and Sasaki 2000; Regoes, Nowak, and Bonhoeffer 2000; Gandon 2004; Lively 2010; Moreno-Gámez, Stephan, and Tellier 2013), one should expect to observe viruses to evolve toward mild infections in genetically homogeneous host populations, as a more efficient exploitation strategy that maximizes transmission. In contrast, host heterogeneity is postulated to select for increased virulence (Ganusov, Bergstrom, and Antia 2002). Indeed, in the absence of trade-offs between the virulence in the different host ecotypes, it is expected that virulence increases without bounds (Regoes, Nowak, and Bonhoeffer 2000). Our results also provide support to this hypothesis. TuMV lineages evolved in demes from a metapopulation are, on average, less virulent on their local hosts than those lineages evolved in the well-mixed host population. However, a note of caution must be added here: the aforementioned models assume a trade-off between virulence and transmission rates because more virulent viruses will reduce their host’s lifespan or production of viable progeny. In our experiments, this trade-off has been broken since transmission rate is controlled by us and, therefore, virulence can increase without paying a cost. Therefore, our results are consistent with the expectation that mild infections are selected in genetically homogenous host populations, though the precise mechanism does not rely in the transmission–virulence trade-off and needs to be explored in future work.

Genetic variability in host susceptibility to infection and the evolution of infectiousness and virulence of parasites has been well documented in animals and plants (see Schmid-Hempel and Koella 1994; Parrat, Numminen, and Laine 2016 for reviews). Resistance is not always achieved against all parasite genotypes but results from the specific interaction between the genotypes of the host and the pathogen in agreement with the so-called gene-for-gene relationship. This mechanism proposes that for each locus determining resistance in the plant host, there is a corresponding locus for virulence in the virus. This model predicts that infection matrices must show significant nestedness: the most susceptible host carrying few resistance genes will be successfully infected by every virus, while the host genotypes carrying large number of resistance loci will be infected by very few viral genotypes carrying more virulence loci (Weitz et al. 2013). To test whether the gene-for-gene relationship has evolved in our experiment, we evaluated the structure of the infection matrices. In agreement with the gene-for-gene model, we found that the infection matrix for the TuMV lineages evolved in homogeneous host demes was significantly nested, with viruses evolved in the most restrictive host Gy-0 deme being also the most generalist ones able of infecting all alternative ecotypes equally well. In contrast, viruses evolved in the most permissive host Oy-0 deme were the most specialized ones, infecting alternative hosts with less efficiency. Interestingly, the lineages evolved in the well-mixed host population did not generate a nested matrix, suggesting that fluctuating selection avoided the fixation of beneficial alleles in all required virulence loci. It is also noteworthy mentioning that, regardless whether or not nested, both infection matrices were significantly modular, suggesting that the mechanism underlying TuMV–A. thaliana interaction is far more complex than expected under a simple gene-for-gene model. Two possible mechanisms will contribute to create modularity: first, plant ecotypes sharing alleles in given sets of defense-response genes will likely impose similar selective constraints to the virus. Second, negative pleiotropic effect of mutations in small and compacted RNA genomes limits the number of alternative evolutionary pathways (Belshaw et al. 2008).

4.3 Are our findings in agreement with those from other previous evolution experiments?

We would like to compare our findings with those from studies with small RNA and large DNA viruses. First, Cuevas, Moya, and Elena (2003) simulated in vitro the migration of Vesicular stomatitis virus (VSV) between three alternative cell types, and observed that migration rate had a negative impact on the rate of fitness improvement on each cell type. With no migration, viruses evolved as specialists with high fitness in their local host, whilst at 50 per cent migration rate per generation VSV evolved as a lower fitness generalist metapopulation. This result is apparently at odds with our findings and with theoretical expectations. However, it is worth noticing that Cuevas, Moya, and Elena (2003) did not measure virulence but fitness in competition against a common reference strain. Since fitness and virulence are not perfectly correlated traits in cell cultures for VSV (Furió et al. 2012), this may account for the discrepancy.

Second, Hillung et al. (2014) evolved TEV, another member of the Potyvirus genus, in five different ecotypes of A. thaliana in a set up similar to our deme host metapopulations treatment. In full agreement with our results, the infection matrix of the evolved TEV genotypes on the five ecotypes was significantly nested and thus compatible with a gene-for-gene relationship.

Third, Boots and Mealor (2007) evolved lineage of the large DNA betabaculovirus PiGV in lines of its host lepidopteran at varying viscosity of the food media in which larvae grew. Increasing viscosity decreased the mobility of the larvae and thus created stronger spatial structure. Viruses evolved in spatially structured host populations evolved to be less virulent than those maintained in well-mixed host populations, matching our findings.

Fourth, Kubinak et al. (2014) and Cornwall et al. (2018) observed that serial passages of the mouse-specific Friend virus (FV) complex in mouse populations that differed in diversity in the major histocompatibility complex (MHC) resulted in different fitness and virulence. Both traits rapidly increased when FV was evolved only in inbred mouse strains of the same genotype. However, both viral fitness and virulence were significantly reduced when FV was evolved in mouse populations containing different MHC genotypes.

Finally, Berngruber, Lion, and Gandon (2015) showed that a latent λ phage won competitions against a virulent one in a spatially structured bacterial host populations but lost in well-mixed ones. This result agrees with our observation that evolution in isolated homogeneous host populations selected for less virulent TuMV while evolution in well-mixed heterogeneous population selected for more virulent virus.

4.4 Are our results compatible with field observations?

There is a number of field studies performed with plant viruses whose results are fully compatible with those from our experimental approach. We are just mentioning a few. Malpica et al. (2006) analyzed the prevalence of five plant viruses on twenty-one wild plant species and found that specialization was the most successful viral strategy, with highly selective viruses (specialist) being the most prevalent ones. In a recent follow up study, McLeish et al. (2017) have applied ecological networks theory to explore the association between plant host diversity and virus spatial distributions. Two important observations were made: (1) an association between host abundance and viral prevalence and (2) most prevalent viruses behave as facultative generalists, meaning that they have the widest host range breadth but can narrow it to the most permissive host without paying a fitness cost.

Rodelo-Urrego et al. (2013) explored the association between the population structure of Capsicum annuum glabriusculum and the prevalence of pepper golden mosaic and pepper huasteco yellow vein begomoviruses, founding that landscape heterogeneity affected the epidemiology and genetic structure of the begomoviruses. When host population diversity was removed by human intervention, the prevalence of the viruses increased. Higher levels of anthropization and loss of plant biodiversity (i.e. similar to our genetically homogenous host populations) resulted in an increase in the genetic diversity of the begomoviruses (Rodelo-Urrego, García-Arenal, and Pagán 2015).

5. Concluding remarks

Understanding of ecological factors, such as the spatial distribution of host genotypes, is critical for deciphering and predicting the evolutionary and epidemiological dynamics of infectious diseases. The role of spatial host structure has attracted quite a lot of attention from theoreticians that made a number of interesting testable predictions. Unfortunately, many of these still remain untested. Here, we have tried to provide experimental evidences for a plant-virus pathosystem to test the hypothesis that host population homogeneity would promote fast local adaptation and low virulence of the virus, whereas the presence of host genetic heterogeneity for susceptibility to infection will slow down the rate of viral evolution and favor more virulent viruses. Supporting these predictions, we found faster TuMV adaptation to homogeneous than to heterogeneous A. thaliana experimental populations. However, viruses evolved in heterogeneous host populations were more pathogenic and infectious than viruses evolved in the homogeneous population. Furthermore, the viruses evolved in homogeneous populations showed stronger signatures of local specialization than viruses evolved in heterogeneous populations. These results illustrate how the genetic diversity of hosts in an experimental ecosystem favors the evolution of virulence of a pathogen and may help agronomists to handle crops in ways that will minimize the rise and spread of virulent viral strains.

Data availability

Raw data, including numbers of symptomatic plants, computation of AUDPS and infectiousness values at each passage and experiment are available at LabArchives under doi: 10.25833/q8y2-f433.

Funding

We thank Francisca de la Iglesia for continuous excellent technical support. Work was supported by Spain’s Agencia Estatal de Investigación—FEDER grant BFU2015-65037-P and Generalitat Valenciana grant GRISOLIA/2018/005 to S.F.E. R.G. was supported by Spain’s Agencia Estatal de Investigación pre-doctoral contract BES-2016-077078.

Conflict of interest: None declared.

Click here for additional data file.