1 Introduction
==============

Experimental evolution has been used to investigate drug resistance, pathogenesis, disease emergence, and other important topics in virology, as well as aid vaccine development.  For instance, selection experiments in cell cultures can help to predict the emergence of drug-resistance mutations and their fitness costs in medically relevant viruses such as HIV-1, hepatitis C virus, and influenza virus ([@vev008-B25]; [@vev008-B7]; [@vev008-B39]; [@vev008-B42]; [@vev008-B46]).  In the context of viral emergence, experimental evolution has been used to study the selective constraints acting on host radiation ([@vev008-B43]), the relationship between ecological specialization and host-shift risk ([@vev008-B44]), and the effects of host alternation on the adaptability and genetic diversity of vector-borne viruses ([@vev008-B10]; [@vev008-B9]).  Viral evolution in the laboratory has been particularly useful for the development of live attenuated vaccines.  The pioneer procedures used for creating the first oral polio vaccines in the 1950s involved serial transfers in non-human hosts, plaque-to-plaque transfers, and viral propagation at low temperatures ([@vev008-B26]).  Similar techniques were used for the development of vaccines against measles, mumps, and yellow fever ([@vev008-B27]).  In these experiments, several evolutionary processes were probably involved in viral attenuation.  First, serial transfers in alternate host cells often lead to fitness losses in the original host, because mutations that are selectively advantageous in a given environment tend to be costly in alternate environments ([@vev008-B43]; [@vev008-B18]; [@vev008-B6]).  Second, plaque-to-plaque transfers typically promote the accumulation of deleterious mutations through random genetic drift ([@vev008-B8]; [@vev008-B17]; [@vev008-B16]).  Finally, relaxed selection in permissive environments should favor the accumulation of conditionally deleterious mutations. 

However, most of the above studies were designed with little knowledge of the molecular and cellular pathways involved.  Here, we sought to take a more gene-based approach to the experimental evolution of viral attenuation.  To explore this, we adapted a model virus to a permissive cellular environment by knocking-out the double-stranded RNA (dsRNA)-dependent protein kinase PKR, a major component of the innate immune response against viruses.  PKR is induced by type-I interferons and is directly activated by viral dsRNA present in the cytoplasm.  The best-known anti-viral effect of PKR is exerted through phosphorylation of the eukaryotic translation initiation factor 2, which induces translational shutoff, but PKR mediates additional antiviral responses through its effects on protein phosphorylation, mRNA stability, and apoptosis ([@vev008-B20]; [@vev008-B35]; [@vev008-B15]; [@vev008-B29]).  The ability to suppress PKR activity and other innate immune responses is thus a major determinant of viral fitness and has promoted the evolution of multiple immune evasion mechanisms including inhibition of upstream pathogen sensors such as RIG-I, signal transducers, and effector antiviral proteins such as PKR ([@vev008-B5]; [@vev008-B45]; [@vev008-B24]).  Therefore, we reasoned that relaxation of the selective pressure exerted by PKR should create a permissive environment that may promote the fixation of effectively neutral mutations usually not tolerated in cells with intact innate immunity.  In addition to this 'drift hypothesis', other processes may contribute to promoting viral attenuation.  For instance, some mutations may be selectively advantageous in PKR negative but neutral or deleterious in normal cells (selection hypothesis). 

To test these ideas, we chose vesicular stomatitis virus (VSV) as model system, a prototypic non-segmented negative-stranded RNA virus of the family Rhabdoviridae.  VSV can be easily propagated in most mammalian cells and is highly sensitive to PKR, as illustrated by the fact that the viral lethal dose is reduced by orders of magnitude in PKR-deficient mice, although this effect is more modest in cell cultures ([@vev008-B2]; [@vev008-B41]).  We found that VSV readily adapted to PKR-deficient mouse embryonic fibroblasts (MEFs), and that this adaptation was PKR-dependent, since the evolved viruses exhibited no fitness gains in PKR-positive MEFs.  We also sought to investigate the genetic basis of this adaptation by sequencing the full viral genome.  This revealed several parallel substitutions in the evolved lines, and we demonstrate the beneficial fitness effects of these substitutions by site-directed mutagenesis.  However, most parallel substitutions were also beneficial in PKR-positive MEFs and other cell types and thus were not responsible for PKR-dependent changes in fitness.  Therefore, some of the non-parallel substitutions may be deleterious in PKR-positive cells and counteract the generally beneficial effects of parallel substitutions.  Consequently, in the absence of PKR, relaxed selection may have allowed for the fixation of mutations with conditionally deleterious effects in cells with normal antiviral immunity, but these mutations did not reach fixation in a deterministic or repeatable manner.  Hence, the genetic data lend support to the drift hypothesis for viral attenuation.  We also found that one of the parallel substitutions appearing in multiple evolution lines was strongly deleterious when assayed as a single mutation, underscoring the importance of epistasis for adaptation.  Collectively, our results suggest that experimental evolution in cells with gene-specific defects in innate immunity may drive viral attenuation, but that this process can be obscured by the appearance of generally beneficial mutations, genetic drift, and by gene-by-gene interactions. 

2 Methods
=========

2.1 Virus and cell culturing
----------------------------

VSV was recovered from a cDNA clone, originally created by Whelan et al. ([ @vev008-B47]), by transfecting BHK-21 cells as previously described ([@vev008-B40]).  This virus was used as the founder of the evolution experiment (T0).  Primary MEFs derived from *pkr* −/− or wild-type C57BL6 mice were isolated as previously described ([@vev008-B34]) and were provided by Dr Carmen Rivas (Universidad de Santiago de Compostela, Spain).  Cells were cultured in Dulbeco's modified Eagle's Medium (Invitrogen) supplemented with 10 per cent fetal bovine serum (Invitrogen) at 37°C under 95 per cent humidity and 5 per cent CO~2~, and immortalized by serial passaging.  The absence of the *pkr* gene was verified by polymerase chain reaction (PCR). 

2.2 Experimental evolution
--------------------------

MEF monolayers containing approximately 10^5^ cells were inoculated with 10^4^ particle forming units (PFU) and incubated for 24 h, the time required to reach a titer plateau of approximately 5 × 10^7^ PFU/ml, as determined by growth curves.  The supernatants were titrated after each passage by the standard plaque assay in baby hamster kidney (BHK-21) cells (ATCC) and diluted as required to infect fresh cells with a constant number of PFUs throughout the course of the experiment.  Assuming a yield of approximately 1,000 PFU/cell, the estimated number of viral generations (i.e., infection cycles) per passage is ln(5 × 10^7^/10^4^)/ln 10^3 ^= 1.23 ([@vev008-B28]). 

2.3 Growth curves and model fitting
-----------------------------------

MEF monolayers containing approximately 10^5^ cells were inoculated with 10^4^ PFU, and small volumes were sampled at the indicated times and titrated by the standard plaque assay.  A logistic growth model with viral decay of the form $N_{t} = \frac{N_{\text{max}}}{1 + \text{e}^{c - rt}}\text{e}^{- dt}$ was fit to the data, where *N*~t~ is the observed titer, *r* the exponential growth rate, *N*~max~ the maximal viral yield in the absence of degradation, *c* sets the initial conditions, and *d* is the viral decay (degradation) rate.  The model was fit to the data by the least-squares method for each experimental replicate separately.  We obtained *r*^2^ values higher than 0.95 in all cases and higher than 0.99 in most cases. 

2.4 Fitness assays
------------------

MEF monolayers containing approximately 10^5^ cells were inoculated with 10^4^ PFU of the assayed virus (evolved or T0) and a common competitor which was isogenic to the T0 virus except for a single-nucleotide replacement (C3853A) in the surface glycoprotein gene G which controls sensitivity to a monoclonal antibody ([@vev008-B40]).  The proportion of each competitor at inoculation and after a single 24 h transfer was determined by titrating the mix in the presence and in the absence of antibody in BHK-21 cells.  Fitness was calculated as *W* = *R*~24~/*R*~0~, where *R* denotes the titer ratio of the assayed virus relative to the common competitor at inoculation (*T*0) and after 24 h.

2.5 Sequencing
--------------

Viral RNA was extracted from virus supernatants using the NucleoSpin RNA Virus kit (Macherey-Nagel), reverse transcribed with AccuScript (Agilent Technologies), and used for PCR amplification.  The viral genome was amplified using eight different, overlapping, primer pairs (sequence available upon request).  PCRs were column purified and used for Sanger sequencing.  Chromatograms were analyzed with the Staden package (<http://staden.sourceforge.net>). 

2.6 Site-directed mutagenesis
-----------------------------

Mutations were introduced using a pair of self-complementary primers carrying the desired change.  The cDNA clone was amplified for eighteen cycles using the high-fidelity Phusion DNA polymerase (New England Biolabs).  To remove template DNA, amplification products were digested with DpnI (New England Biolabs), which selectively cuts methylated GATC sites.  Products were used for transforming competent *Escherichiacoli* cells by the rubidium chloride heat-shock method.  Plasmid DNA was then purified using the Nucleospin Plasmid purification kit (Macherey-Nagel) and used for transfecting BHK-21 cells as previously described ([@vev008-B47] [@vev008-B40]).  Briefly, young 90 per cent confluent BHK-21 cells were infected with a recombinant vaccinia virus expressing T7 RNA polymerase and then co-transfected with the full-length VSV cDNA clone and three helper plasmids encoding the P, L, and N proteins.  Transfections were done using Lipofectamine LTX (Life Technologies), following manufacturer's instructions.  After 6 h, 25 µg/ml 1-β-D-arabinofuranosylcytosine was added to inhibit vaccinia replication.  After 3--4 days, supernatants were tested for the presence of infectious VSV particles by the standard plaque assay, vaccinia virus was removed by filtration, and one additional VSV infection cycle was performed to reach sufficient titer. 

2.7 RNA structure prediction
----------------------------

The predicted minimum free energy of the trailer sequence was obtained using the RNAfold server available online (rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi) with default parameter options. 

3 Results
=========

VSV was recovered from an infectious cDNA clone and used for propagating five independent lines (L1--L5) for forty transfers in *pkr*−/− MEFs, which corresponds to approximately fifty generations (infection cycles) of viral evolution.  To assess adaptation, we performed growth curves of the founder virus (T0) and each evolved line ([Fig.  1](#vev008-F1){ref-type="fig"}A).  Fit of a logistic growth model with viral decay showed that the exponential growth rate increased significantly in each of the lines (from 0.38 h^--1^ in T0 to 0.57--0.66 h^--1^; [Table 1](#vev008-T1){ref-type="table"}).  This change represents a strong fitness gain, since the amount of viral progeny produced in 10 h of growth (∼one infection cycle) ([@vev008-B12]) was approximately 10-fold higher for the evolved lines than for the founder virus.  The maximal viral yield also tended to increase (from two- to fourfold; [Table 1](#vev008-T1){ref-type="table"}) except in line L1, although the data showed higher variance, and differences were thus less significant.  Finally, we found no evidence for changes in particle stability, as inferred by viral decay rates.  We therefore conclude that VSV adapted to PKR-negative cells mainly by increasing its short-term growth rate.  Figure 1.VSV adaptation to pkr-/- MEFs. ( A) Growth curve of the five evolved lines (filled circles) and the founder (T0) virus (open circles).  Each line was assayed in triplicate.  Error bars indicate the standard error of the mean. ( B) Competition assays of the T0 and evolved viruses against a common, phenotypically marked competitor in *pkr* -/- MEFs.  Each line represents the change in the ratio of the assayed virus to the common competitor after a single 24 h transfer and corresponds to one replicate assay.  Four replicates were performed for each competition. ( C) Competition assays in *pkr* +/+ MEFs were performed in the same way as in (B).  Table 1.Growth parameters of the five lines adapted to PKR-negative MEFs. LineGrowth rate (h^--1^)Maximal yield (PFU/ml × 10^--7^)Decay rate (h^--1^)T00.38 ± 0.021.27 ± 0.210.02 ± 0.01L10.66 ± 0.00\*\*\*1.03 ± 0.080.03 ± 0.01L20.64 ± 0.00\*\*\*4.48 ± 1.250.04 ± 0.02L30.63 ± 0.02\*\*\*3.33 ± 0.91[\*](#vev008-TF1){ref-type="table-fn"}0.03 ± 0.01L40.57 ± 0.02\*\*\*2.61 ± 0.52[\*](#vev008-TF1){ref-type="table-fn"}0.04 ± 0.01L50.63 ± 0.03\*\*\*5.09 ± 0.93[\*](#vev008-TF1){ref-type="table-fn"}0.05 ± 0.02[^1]

To further test for fitness changes, we also performed head-to-head competition assays against a phenotypically marked competitor isogenic to the T0 virus.  Fitness (*W*), measured as the change in titer ratio relative to the common competitor after a single transfer, was \>10 times higher for the evolved lines (*W* = 18.5, *W* = 38.2, *W* = 17.5, *W* = 13.0, *W* = 15.2 for L1--L5, respectively; *t*-test of log*W* against zero: *P* \< 0.001 for all lines) than for the T0 virus (*W* = 1.15; *P* = 0.132; [Fig.  1](#vev008-F1){ref-type="fig"}B).  To ascertain whether the observed fitness changes were specific to PKR-negative cells, we performed identical competition assays in *pkr* +/+ MEFs.  Lines L1, L2, L4, and L5 showed marginally significant increases in fitness (*t*-test of log*W* against zero: 0.134 \> *P* \> 0.057), whereas L3 showed a more marked increase (*W* = 6.22; *P* = 0.002) but still less pronounced than in *pkr* −/− cells ([Fig.  1](#vev008-F1){ref-type="fig"}C). 

Therefore, viral adaptation was PKR dependent.  We also noticed that, throughout the course of the evolution experiment, there was a strong reduction in plaque size in BHK-21 cells.  Although *pkr*−/− MEF monolayers did not reach enough density to perform plaque assays, we verified the reproducibility of the small-plaque phenotype in *pkr* +/+ MEFs ([Fig.  2](#vev008-F2){ref-type="fig"}A).  The population frequency of this plaque phenotype gradually increased until reaching fixation (\>98%) in all lines except L1, in which the population was still polymorphic by passage 40 although small plaques were highly abundant (∼90%; [Fig.  2](#vev008-F2){ref-type="fig"}B).  Therefore, adaptation to *pkr* −/− cells was accompanied by a clear reduction in plaque size in PKR-positive cells.  Figure 2.Serial transfers in *pkr* -/- MEFs select for a small-plaque phenotype in PKR-positive cells. ( A) Plaque assays of the T0 virus and line L2 (transfer 40) in BHK-21 cells and in *pkr* +/+ MEFs. ( B) Fixation of the small-plaque phenotype (assayed in BHK-21 cells) during evolution.  Approximately 100 plaques from a single plaque assay were counted to determine the small/total plaque ratio.  Small plaques were not observed before passage 22 and reached quasi-fixation in all lines except L1 by transfer 35. 

To investigate the genetic basis of the observed changes in fitness and plaque size, we sequenced each of the evolved lines.  In total, we found 22 nucleotide substitutions compared with the T0 virus, of which seven (G2686A, G3351A, A3460C, G3837A, A8679C, C10224U, and A11090G) appeared in two or more lines and thus represent cases of parallel evolution ([Table 2](#vev008-T2){ref-type="table"}).  Previous work with VSV ([@vev008-B11]; [@vev008-B31]; [@vev008-B31], [@vev008-B33]; [@vev008-B36]; [@vev008-B13]) and several other viruses ([@vev008-B14]; [@vev008-B37]; [@vev008-B1]; [@vev008-B3]; [@vev008-B48]; [@vev008-B4]) has suggested that parallel substitutions tend to be selectively advantageous mutations, although their effects have not usually been assayed directly.  Six of these changes were non-synonymous, and most produced non-conservative amino acid replacements, thus further suggesting a significant fitness effect.  We noticed that the G3837A substitution occurred in lines L2--L5, all of which showed a small-plaque phenotype in *pkr* +/+ cells but that it was absent from L1 in which this phenotype was not fixed, thus suggesting that this mutation was responsible for small plaque formation.  Table 2.Genetic changes detected in the five lines adapted to PKR-negative MEFs. GenomeGeneResidueL1L2L3L4L5C1282UNNone−−−+−C2649UMH134Y−−+−−G2686AMG146E−+−−+U3030CInterNone+−−−−G3351AGE92K−++−+A3460CGQ128P+−−+−G3837AGE254K−++++C4180AGT368K+−−−−G4182AGE369K+−−−−U4397CGNone−−+−−A4730GInterNone−−−+−G5557ALNone+−−−−A5700ULD323V+−−−−A6931GLNone−−+−−U7369CLNone−−+−−U7522CLNone−−+−−G7591ALNone+−−−−A8320GLNone+−−−−A8679CLD1316A−+−−+C10224ULT1831M++−+−C10282ALNone−−+−−A11090GTrailerNone−+−−+

To dissect the effects of the parallel substitutions on fitness and plaque size, we introduced them individually into the parental cDNA clone by site-directed mutagenesis and transfected BHK-21 cells to recover mutant viruses.  We successfully recovered all single mutants except C10224U. As predicted, the G3837A mutation produced a small-plaque phenotype but, immediately after recovery, normal-plaque revertants were observed at 3, 3, 4, 40, and 90% population frequencies in five independent transfections, indicating that the G3837A mutation was severely deleterious in BHK-21 cells.  This was confirmed by competition assays in *pkr* +/+ MEFs, which yielded an extremely low fitness value of *W* = 0.081 (*t*-test of log *W* against zero: *P* = 0.003).  We then performed competition assays for each of the six single mutants in *pkr* −/− cells.  Beneficial fitness effects were confirmed for mutants G2686A (*W* = 2.18; *P* = 0.047), G3351A (*W* = 22.9; *P* = 0.008), and A3460C (*W* = 8.46; *P* = 0.025), whereas A8679C and A11090G did not deviate significantly from neutrality (*P* ≥ 0.08; [Fig.  3](#vev008-F3){ref-type="fig"}A).  Interestingly, the G3837A mutation also had a negative fitness impact on *pkr* −/− MEFs (*W* = 0.254; *P* = 0.031).  The fact that this mutation was deleterious regardless of cell type (BHK-21, *pkr* + /+ MEFS, or *pkr* −/− MEFs) but appeared in four out of five evolved lines indicates strong sign epistasis.  However, these genetic interactions were probably complex, since the G3837A change did not occur in systematic association with any other substitution ([Table 2](#vev008-T2){ref-type="table"}).  Figure 3.Individual fitness effects of parallel substitutions. ( A) Competition assays of the T0 and single mutants against a common, phenotypically marked competitor in *pkr* -/- MEFs.  Each line represents the change in the ratio of the assayed virus to the common competitor after a single 24 h transfer.  Three replicates were performed for each competition. ( B) Competition assays in *pkr* +/+ MEFs. ( C) Correlation between fitness in *pkr* −/− and *pkr* +/+ cells.  Error bars represent the standard error of the mean

We then performed competition assays of the single mutants in *pkr* +/+ MEFs to determine whether their effects were PKR specific.  However, we again found increased fitness for mutants G2686A (*W* = 6.00; *P* = 0.006), G3351A (*W* = 5.06; *P* = 0.050), and A3460C (*W* = 3.11; *P* = 0.041), thus indicating that these were generally beneficial mutations in MEFs regardless of PKR status.  Also reproducing the results obtained in *pkr* −/− cells, A8679C and A11090G were apparently neutral (*P* ≥ 0.19; [Fig.  3](#vev008-F3){ref-type="fig"}B).  Therefore, the hypothesis of an adaptive trade-off whereby the effects of these mutations in *pkr* +/+ and *pkr* −/− cells would be negatively correlated was not supported (Pearson *r* = 0.737; *P* = 0.059; [Fig.  3](#vev008-F3){ref-type="fig"}C).  We, therefore, conclude that PKR-specific adaptation was determined by other, non-parallel, mutations or a combination of mutations showing epistasis. 

4 Discussion
============

We found a large number of parallel sequence changes in the five evolved lines, similar to previous work ([@vev008-B11]; [@vev008-B31], [@vev008-B33]; [@vev008-B36]; [@vev008-B13]).  To further assess whether these changes were cell type -specific, host species-specific, or generally beneficial mutations under laboratory conditions, we reviewed sequences from previous experiments in which VSV was adapted to MEFs, BHK-21 cells, human cervix cancer HeLa cells, and dog kidney MDCK cells.  We found that the G2686A substitution appeared in two of four lines evolved in HeLa cells ([@vev008-B36]), whereas the A3460C mutation appeared in one line evolved in BHK-21 cells ([@vev008-B32]) and in one line evolved in p53-deficient MEFs ([@vev008-B21]).  Strikingly, the highly deleterious mutation G3837A appeared in five lines evolved in BHK-21 cells ([@vev008-B33]) and in three out of five lines evolved in MDCK cells ([@vev008-B13]).  The C10224U mutation, which we failed to rescue, was found in one line evolved in p53 knock-out MEFs ([@vev008-B21]).  Additionally, the C4180A mutation, which appeared only once in *pkr* −/− MEFs was found in more than ten lines evolved in BHK-21 cells ([@vev008-B33]) and in one of four lines evolved in HeLa cells ([@vev008-B36]).  Finally, although we did not find the G3351A mutation reported in previous works, it appeared under our laboratory conditions after serial transfers in BHK-21 cells (unpublished data).  Therefore, most parallel substitutions probably correspond to generally beneficial mutations conferring adaptation to cell culture conditions.  This is particularly likely considering that, in our experiments, the founder virus was recovered from an infectious cDNA clone with no previous laboratory passaging history.  A candidate factor that is common to all the above cell lines is immortalization, which often tends to constitute a permissive environment for viruses.  The various parallel mutations found in the glycoprotein G may also improve viral attachment/entry in mammalian fibroblasts such as BHK-21, MEF, or MDCK cells, which are not the natural target cell type of VSV. 

All parallel changes were non-synonymous except for A11090G, which mapped to the trailer sequence.  To explore possible effects of this mutation, we obtained the predicted minimum free energy secondary structure of the trailer sequence using the RNAfold program ([@vev008-B23]).  This gave a highly folded structure, with 66 out of 100 sites forming intramolecular base pairs ([Fig.  4](#vev008-F4){ref-type="fig"}).  The A11090G mutation mapped to position 28 of the trailer and introduced a new base pair which displaced the next 12 bases, thus producing a major change in the predicted structure.  The leader and trailer are the only VSV transcripts which are not 5'-capped and have a 5´ppp terminus ([@vev008-B22]).  This type of RNA can be a RIG-I agonist provided that it contains dsRNA regions formed by intramolecular base pairs, triggering interferon secretion and an innate immune response.  Furthermore, partially double-stranded 5´ppp RNA can directly activate PKR in an interferon-dependent manner ([@vev008-B30]), leading to EI2F alpha phosphorylation and translation arrest.  This suggests that slight changes in RIG-I induction may be tolerated in the absence of a major antiviral effector such as PKR but not in normal cells.  However, as shown above, the A11090G single mutant did not experience fitness defects in *pkr* −/− or *pkr* +/+ MEFs.  It is possible though that PKR-specific fitness effects may arise from epistatic interactions between A11090G and other mutations.  Interestingly, the A11090G, A8679C, and G2686A changes co-occurred in lines L2 and L5, suggesting that they may be an epistatic set.  This speculation is substantiated by the fact that G2686A produced a non-conservative (G146E) amino acid replacement in protein M, which is the only VSV protein known to block innate immunity ([@vev008-B19]; [@vev008-B38]).  Figure 4.Effect of the A11090G mutation on the predicted secondary structure of the trailer.  The minimum free energy (MFE) structure predicted by the RNAfold program is shown for the T0 (top) and A11090G (bottom) trailers.  Numbers indicate genome positions.  The mutated base is circled.  In the predicted MFE structure, the single-base substitution modified 14 neighboring base pairs, as indicated

In conclusion, our results show that serial transfers in cells with a specific genetic defect in innate immunity do not produce an easily identifiable, repeatable change in viral sequence which could be used for engineering an attenuated virus.  Parallel substitutions occurred, but these were generally beneficial and thus were not responsible for *pkr*-specific changes in fitness.  These changes should thus be driven by singleton substitutions, which may sometimes act epistatically.  Although the majority of these substitutions were synonymous, they may still have fitness effects, which remain to be characterized.  In contrast to directed evolution experiments showing deterministic changes at the molecular level, we suggest that the fixation of mutations under relaxed selection is governed by random genetic drift.  These types of mutations would be tolerated under the permissive environment but not in cells with intact antiviral responses, thus driving viral specialization and attenuation.  However, the stochastic nature of molecular evolution and epistasis make it inherently difficult to use a relaxed selection approach for the rational, molecularly informed attenuation of viruses, a desirable goal in the clinical setting. 

We thank Silvia Torres for technical assistance and Dr Carmen Rivas for the *pkr*−/− MEFs.  This work was supported by grants from the Spanish MINECO (BFU2013-41329) and the European Research Council (ERC-2011-StG-281191-VIRMUT) to R.S.

Data availability
=================

All data are available upon request. 

*Conflict of interest:* None declared. 

[^1]: \**P* \< 0.05; \*\*\**P* \< 0.001 (*t*-test against T0).