Molecular and Functional Bases of Selection against a Mutation Bias in an RNA Virus

The selective pressures acting on viruses that replicate under enhanced mutation rates are largely unknown. Here, we describe resistance of foot-and-mouth disease virus to the mutagen 5-fluorouracil (FU) through a single polymerase substitution that prevents an excess of A to G and U to C transitions evoked by FU on the wild-type foot-and-mouth disease virus, while maintaining the same level of mutant spectrum complexity. The polymerase substitution inflicts upon the virus a fitness loss during replication in absence of FU but confers a fitness gain in presence of FU. The compensation of mutational bias was documented by in vitro nucleotide incorporation assays, and it was associated with structural modifications at the N-terminal region and motif B of the viral polymerase. Predictions of the effect of mutations that increase the frequency of G and C in the viral genome and encoded polymerase suggest multiple points in the virus life cycle where the mutational bias in favor of G and C may be detrimental. Application of predictive algorithms suggests adverse effects of the FU-directed mutational bias on protein stability. The results reinforce modulation of nucleotide incorporation as a lethal mutagenesis-escape mechanism (that permits eluding virus extinction despite replication in the presence of a mutagenic agent) and suggest that mutational bias can be a target of selection during virus replication.


Introduction
RNA viruses exploit high mutation rates to generate complex mutant spectra (ensemble of genetic variants) which allow finding mutational pathways towards adaptation (Andino and Domingo 2015;Domingo 2016;Domingo et al. 1978;Lauring et al. 2013). Increasing the mutation rate above the basal level determined by the replication machinery can generate altered viral RNA and proteins in detriment of infectious progeny production. This is the basis of an antiviral strategy termed lethal mutagenesis in which virus-specific mutagenic agents such as nucleotide analogues increase mutation rates above a survival threshold (Andino and Domingo 2015; Dapp et al. 2013;Domingo 2005;Domingo and Schuster 2016a;Domingo et al. 2012;Graci and Cameron 2008;Holland et al. 1990;Loeb et al. 1999). However, viruses display molecular mechanisms to overcome the effect of mutagenic activities mediated by nucleotide analogues. When dependent on the polymerase, resistance can be achieved through a general increase of polymerase fidelity, a specific restriction to incorporate the mutagenic nucleotide, or modulation of nucleotide incorporation to achieve a balance among mutation types (Agudo et al. 2010;Borderia et al. 2016;Ferrer-Orta et al. 2010;Pfeiffer and Kirkegaard 2005;Vignuzzi et al. 2006;Zeng et al. 2013). Interestingly, modulation of mutation types in the picornavirus foot-and-mouth disease virus (FMDV) to confer resistance to the purine analogue ribavirin can be achieved not only through amino acid replacements in the viral polymerase (termed 3D; Agudo et al. 2010) but also by a single amino acid substitution in non-structural protein 2C, without the need of mutations in the polymerase (Agudo et al. 2016). Thus, modulation of nucleotide incorporation is a well established mechanism to escape ribavirin-induced lethal mutagenesis of FMDV.
The pyrimidine analogue 5-fluorouracil (FU) is mutagenic for several RNA viruses (Agudo et al. 2009;Holland et al. 1990), and replication in the presence of FU can lead to virus extinction in cell culture and in vivo (Grande-Pé rez, G omez-Mariano, et al. 2005;Grande-Pé rez, L azaro, et al. 2005;Grande-Pé rez et al. 2002). FU is intracellularly converted into fluorouridine-triphosphate (FUTP), which is then incorporated as fluorouridine-monophosphate (FUMP) into RNA and can compete with incorporation of the standard nucleotides. FU exerts both an inhibitory and a mutagenic activity on FMDV (Agudo et al. 2008). As inhibitor, FU in its FUTP form blocks initiation of viral RNA synthesis. As mutagen, FUTP is incorporated during viral RNA elongation, and it evokes mainly A ! G and U ! C transitions. Here, we describe the selection, isolation and characterization of a FMDV mutant with substitution V173I in the F motif of 3D which is able to survive the mutagenic activity of FU. V173I does not produce direct resistance to the inhibition of FMDV replication by FU, but compensates for the increase of A ! G and U ! C transitions that the wild-type (wt) virus endures in the presence of FU. Compensation in the mutant virus entails an increase of G ! A and C ! U transitions in the presence of FU that approximates the mutational pattern to that of the wt virus replicating in the absence of FU. The same incorporation bias is observed with purified mutant polymerase using in vitro nucleotide incorporation assays. The structural modifications associated with the altered nucleotide incorporation have been examined by comparing the threedimensional structures of the wt and mutant polymerases by X-ray crystallography. Structural changes have been identified at the N-terminal region of the protein. The codon composition of FMDV RNA predicts a higher frequency of non-synonymous mutations as a result of an increase of the proportion of G and C in the genome, with negative consequences for protein 3D stability. Experimental results and predictive algorithms of RNA and protein stability support that in the course of FU mutagenesis, a polymerase that can limit the frequency of A ! G and U ! C transitions will display a selective advantage over its wt counterpart. The results show that mutant spectrum composition can be a target of selection, and evoke mechanisms to escape lethal mutagenesis that do not involve a significant modification of polymerase fidelity or mutant spectrum complexity.

Cells, Viruses, and Infections
Baby hamster kidney 21 (BHK-21) cells were grown in Dulbecco's modified Eagle's medium (DMEM) (Gibco) supplemented with non-essential amino acids and 5% fetal calf serum, as previously described (Domingo et al. 1980). The origin of the FMDVs used in the present study is described in table 1. FMDV-RP [named Bp35 in Sierra et al. 2007] was subjected to 15 additional passages at a MOI of 0.05-0.1 plaque-forming units (PFU) per cell in the absence or presence of 50 lg/ml or 200 lg/ml of FU, leading to populations RP-15, RP(FU50)-15, and RP(FU200)-15, respectively. All infections were cytolytic and were allowed to progress until overt cytopathology (>80% cell killing), generally at 30-50 h post infection (p.i.). The procedures for infection of BHK-21 cell monolayers in liquid medium or under an agar overlay for titration of infectivity have been previously described (Domingo et al. 1980;Sobrino et al. 1983).
Treatment with 5-FU FU was prepared in DMEM at a concentration of 2.5 mg/ml, and stored at À70 C. Cells were preincubated with the desired concentration of FU for 7-13 h prior to infection; infections were allowed to proceed in the presence of the same drug concentration (Sierra et al. 2000).
Preparation of FMDV with Substitution V173I pMT28 (table 1) was used to construct plasmid pMT28-3D(V173I), encoding FMDV with amino acid substitution V173I in the polymerase (3D)-coding region. To this aim, pMT28 DNA was amplified with PFU ultra polymerase and primers with the desired nucleotide changes (the oligonucleotide primers used in the present study are listed in supplementary table S1, Supplementary Material online). A first amplification with V173Iplus and Pol1XbaI was followed by a second amplification with V173Iminus and PolCKpnI; the two distinct mutations in the relevant codon served as marker for the engineered virus. The products of both amplifications were shuffled and amplified with primers 5 0 3D and AV2New. Digestion of the DNA with BamHI (that cleaves at position 7427) and ClaI (position 7004) allowed ligation to pMT28 linearized with the same enzymes. Ligation, transformation of Escherichia coli DH5a, colony screening, nucleotide sequencing, transcription, and RNA transfections were carried out as previously described (Sierra et al. 2007). The effect of substitution V173I in 3D on FMDV replication was studied by comparing FMDV-wt with FMDV-3D(V173I), that differs from FMDV-wt only in the presence of mutations G7126A and A7128C (the first mutation leads to 3D with substitution V173I, and the second mutation serves as a marker for the construct). In this study we refer to FMDV-wt to mean a FMDV whose genomic sequence is identical to that of C-S8c1 or its molecular clone pMT28 (table 1). Likewise, wild type 3D refers to a polymerase 3D whose amino acid sequence is identical to the amino acid sequence of 3D encoded in C-S8c1 or its molecular clone pMT28.

Fitness Assays
Relative fitness between FMDV-wt and FMDV-3D(V173I) was determined by growth-competition experiments in the absence or presence of FU, for a total of eight serial passages. The initial infection was carried out with a 1:1 ratio of infectivity of FMDV-wt:FMDV-3D(V173I). The proportion of the two viruses was determined by RT-PCR, by using discriminatory primers, as previously described (Agudo et al. 2010). The logarithm of the ratio of the two viruses was plotted against the passage number, and the fitness vector was adjusted to the exponential equation y ¼ a . e bx ; the antilogarithm of the vector slope gave the fitness value of FMDV-3D(V173I) relative to FMDV-wt. Control experiments showed that under the same passage conditions of the fitness assays, substitution V173I did not revert upon subjecting FMDV-3D(V173I) to 10 passages in absence of FU, and that FMDV-wt did not acquire substitution V173I after 10 passages in the presence of 400 lg/ml of FU.

Quasispecies and Clonal Analyses
Quasispecies is the ensemble of related mutants subjected to a continuous process of genetic variation, competition among variants and selection of those that constitute a mutant spectrum at any given time. To sample the composition of the mutant spectrum generated in the viral populations and to obtain sequences for mutant polymerase preparation, cDNAs from single genomes (obtained from viral RNA through RT-PCR amplification with primers PolCKpnI and Pol1XbaI, spanning most of the 3D-coding region) were cloned into pGEM-T or pGEM-T Easy (Promega). To avoid redundant cloning of the same genome, an excess of viral RNA was ensured by using samples from which a 1/100 dilution produced a visible DNA band as the product of RT-PCR amplification (Airaksinen et al. 2003). Escherichia coli DH5a was transformed with the ligation products and DNA from individual positive colonies was amplified with Templiphi (GE Healthcare) and sequenced (Macrogen, Inc.).

Computational Analyses
The Quickfold program on the DINAMelt server (http://una fold.rna.albany.edu/, Zuker 2005, 2008) was used to predict the secondary structure stability of the wild type RNA sequences encoding the 3D protein, and of the corresponding mutant RNA sequences observed in the cloned samples (supplementary tables S3 and S4, Supplementary Material online), or used to express the mutant 3D proteins analyzed by ThermoFluor. For each sequence, Quickfold predicted multiple structures with similar folding free energy, and we considered the lowest predicted free energy as an estimate of the RNA stability.
The impact of mutations on the folding free energy (DDG) of the 3D protein was computed using two different programs: PoPMuSiC and DeltaGREM. PoPMuSiC (http://bab ylone.ulb.ac.be/popmusic, Dehouck et al. 2011) relies on a set of statistical potentials derived from known protein structures, which describe the correlations between amino acid types, pairwise inter-residue distances, backbone torsion angles, and solvent accessibilities. It was shown to yield a correlation coefficient of 0.63 between predicted and experimental DDG values, in cross-validation on a set of 2,648 mutations in 131 different proteins (Dehouck et al. 2009). The changes in free energy resulting from substitutions in 3Dwt have been evaluated on the basis of seven different structures (PDB codes: 1U09, 1WNE, 2D7S, 2E9T, 2E9Z, 2EC0, 2F8E). The DDG predictions were very similar for all structures, and we report only the average values. In the case of 3D(V173I), the same structures were used (after introducing the V173I substitution), in addition to the two structures of 3D(V173I) presented here.
DeltaGREM (http://ub.cbm.uam.es/software/Delta_GREM. php) estimates the free energy of the native and the nonnative state adopting the contact matrix representation of protein structures with the contact interaction parameters determined in Bastolla et al. (2000). The non-native state is modeled as the combination of the ensemble of unfolded structures, which dominates when the hydrophobicity is low, and the ensemble of compact, wrongly folded structures, which dominates when the hydrophobicity is high. The misfolded state is modeled through a Random Energy Model (Derrida 1981) of a large set of protein-like compact contact matrices, as described in Minning et al. (2013) [for computational details, see Arenas et al. 2015]. In the present study, the unfolded state was modeled through the configuration entropy of the protein chain. The free energy of a mutant was computed assuming that the native contact matrix does not change. This method yielded a correlation coefficient of 0.72 between predicted and experimental DDG for a set of 195 mutants that fold with two-state thermodynamics (Bastolla 2014).
For the DeltaGREM analysis, we used the available structures of the 3D protein (Agudo et al. 2010;Ferrer-Orta et al. 2010, 2004, 2015 as a template, adopting for each sequence the structural model that minimizes the predicted free energy, which in all cases corresponded to the structure with PDB code 1U09 (Agudo et al. 2010;Ferrer-Orta et al. 2004). Because all sequences generated by the mutant polymerase included also substitution V173I, we repeated the calculations considering mutant sequences where the substitution V173I had been computationally reverted; the mutation did not affect predicted stabilities. We also computed changes in hydrophobicity of the coded protein sequences, using the hydrophobicity scale derived by (Bastolla et al. 2005).

Measurement of Polymerase Stability
The stability of wt and several mutant 3Ds that have been cloned, expressed and purified in the present study was compared using the ThermoFluor assay (Semisotnov et al. 1991), using a 96-well PCR plates (BIORAD). A total of 36 sample wells were prepared: three replicas of the 3Dwt, of the ten 3D mutants, and of the buffer control. Each well contained 42.5 ll of protein solution in 0.5M NaCl, 01.M Tris pH 8.5, 8% glycerol, 0.8 mM DTT and 0.8 mM EDTA buffer plus and 7.5 ll of diluted SYPRO Orange solution (ThermoFischer Scientific, S-6651) at 300Â in the same buffer to reach a final volume of 50 ll and a final protein concentration of 1 mg/ml. The plates were then sealed with an adhesive optical clear seal (MicroAmp Optical Adhesive Film). Data were collected using IQ5 Multicolor Real-Time PCR Detection System (Biorad) and analyzed with iQ TM 5 Optical System Software. The fluorescence in each well was measured at regular intervals with a gradient of 0.5 C per minute over a temperature range spanning from 30 to 95 C. The interaction of the dye with the protein leads to a sigmoidal curve. The temperature at which 50% of the protein is unfolded and bound to the fluorescent dye (apparent melting temperature, T m ) corresponds to the inflexion point of the slope. Raw data were plotted and analyzed using Excel-Microsoft.

Preparation of FMDV Polymerases
To prepare the FMDV polymerase with amino acid substitution V173I, the 3D gene was cloned and expressed from vector pET-28a 3Dpol [vector pET-28a (Novagen) with the 3D-coding region] using described procedures (Arias et al. 2005;Ferrer-Orta et al. 2004). 3D(V173I) was expressed from the corresponding mutant 3D constructed by sitedirected mutagenesis with oligonucleotides V173Iplus and V173Iminus, that included nucleotidic changes G7126A and A7128C (supplementary table S1, Supplementary Material online), using the QuickChange site-directed mutagenesis kit (Stratagene). Enzyme expression and purification by affinity chromatography were performed as previously described (Agudo et al. 2010;Arias et al. 2005;Ferrer-Orta et al. 2004). Enzymes were >95% pure, according to analytical SDS-PAGE and Coomassie brilliant blue staining. The enzyme used for the crystallographic studies was subjected to an additional step of purification by size exclusion chromatography on a Superdex 200 HR column (GE Healthcare).
The same procedure was used to express and purify 3Ds containing two or three amino acid substitutions for studies of thermal stability. The polymerases were chosen among those encoded in genomes from the mutant spectra of FMDV populations passaged in the presence of FU or among ribavirin-resistant and ribavirin-mutagenized FMDVs mutants previously characterized in our laboratory; the aim was to test polymerases that span a broad range of predicted stability effects, their origin and amino acid sequence are given in Results. Temperature required for efficient expression in E. coli, protein yield and aggregation state were recorded for each polymerase.

Polymerization Assays Using Heteropolymeric Symmetrical/Substrates (Sym/Sub) Template-Primers
Incorporation of standard nucleoside-5 0 -triphosphates or 5fluorouridine-5 0 -triphosphate by the wt (3D) and mutant [3D(V173I)] polymerases was quantified using sym/sub RNAs (Dharmacon Research; Arnold and Cameron 2000;Castro et al. 2005;Gohara et al. 2000). The sym/sub RNAs are named according to the two nucleotides nearest to the duplex region; their sequence is given for each experiment. The oligonucleotides were purified, 5 0 -32 P-end-labeled with [c-P 32 ] ATP by polynucleotide kinase (New England Biolabs), purified and annealed in 500 mM Tris-HCl (pH 7.8) and 1 M KCl, by heating at 95 C for 5 min and slow cooling to room temperature. Nucleotide incorporation by 3D was monitored by incubating sym/sub RNA with 1 mM 3D (concentration as active site equivalents, measured by rapid quench flow, described below) and increasing concentration of the required NTP substrates for different time periods (Agudo et al. 2010;Arias et al. 2008). Polymerization products were resolved by electrophoresis on a denaturing 23% polyacrylamide, 7 M urea gel in 90 mM Tris-base, 90 mM boric acid, 2 mM EDTA, scanned with a Phosphorimager (BAS-1500; Fuji), and quantitated with ImageJ 1.45s (NIH).

Rapid Quenched-Flow Assays
To measure incorporation of cognate nucleotides under presteady-state conditions, a model RQF-3 chemical quench-flow apparatus (KinTek Corp., Austin, TX) was used. Experiments were performed at 37 C in 50mM HEPES (pH 7.5), 5 mM MgCl 2 and 10 mM 2-mercaptoethanol. Protein 3D was preincubated with the RNA template/primer and the first incoming NTP (10 mM) to allow formation of 3D-RNA nþ1 product complexes. This reaction mixture was loaded in syringe A (sample loop volume 17.2 ml), and mixed with varying concentrations of NTP loaded in syringe B (sample loop volume 17.4 ml). Reactions were allowed to proceed and stopped with 0.5 M EDTA at different time periods (0.01-5 s). Reactions with each NTP concentration were tested for >6 time points. Products were resolved on a 23% polyacrylamide 7 M urea gel, and scanned with Typhoon FLA 9000 (GE Healthcare). The reaction products were quantitated with ImageJ 1.45s. Data were fitted by non-linear regression to the exponential equation P ¼ A(1e -kobst ) (where A is the amplitude of the burst, k obs is the observed burst rate constant for NTP incorporation, and t is the reaction time), using the program GraphPad Prism 4 (GraphPad Inc.).
To obtain the dissociation constant K d.NTP for NTP binding to the 3D-sym/sub complex, the observed burst rates (k obs ) were plotted as a function of NTP concentration and the data were fitted by non-linear regression to the hyperbolic equation where k pol is the maximal rate for nucleotide incorporation) using GraphPad Prism 4.
To determine the proportion of active sites in the polymerase preparations, incorporation rates in 0.25 s were determined for 1 mM 3D in the presence of varying concentrations of RNA (0.1, 0.2, 0.5, 1, and 2 mM) and saturating NTP concentration (2 mM). Values of total amount of RNA elongated were plotted and fitted to a quadratic equation to obtain the concentration of active enzymes.

Other Assays with 3D
The VPg uridylylation activity was measured employing two different protocols that include either Mn 2þ and poly(A), or Mg 2þ and cre as template, plus 3CD. In the assay in the presence of Mn 2þ the reaction mixture contained 30 mM MOPS, pH 7.0, 33mM NaCl, 0.6mM MnCl 2 , 40 ng/ml poly(A) (average length of 300 residues), 150 mM VPg, 8% glycerol, 0.4 mg/ml bovine serum albumin (New England Biolabs), and 50 mM [a-32 P UTP] (0.01 mCi/ml; 200 mCi/ nmol) in the presence of different concentrations of FUTP. The reaction was carried out for 30 min at 37 C, and it was stopped by adding EDTA to a final concentration of 83 mM. Relative activity was measured from the amount of VPg-pU(pU) determined electrophoretically as previously described (Agudo et al. 2008). For the assay in the presence of Mg 2þ the reaction mixture contained 30 mM MOPS, pH 7.0, 33 mM NaCl, 5 mM MgCl 2 , 20nM cre, 0.16mM 3CD, 150 mM VPg, 8% glycerol, 0.4 mg/ml bovine serum albumin, and 5 mM [a-32 P UTP] (0.01 mCi/ml; 200 mCi/nmol).
Binding of 3D to RNA was quantified by a gel mobility shift assay as previously described (Arias et al. 2005). The RNA binding mixture included 285 nM of the labeled sym/sub-UG, 100 mM MOPS (pH 7.0), 10 mM Mg(CH 3 COO) 2 , 5% polyethylene glycol, and increasing concentrations (0, 250, 500, 1,000, and 2,000 nM) of 3D. The mixtures were loaded onto a non-denaturing 4% polyacrylamide gel in TAE retarding buffer (120 mM Tris-acetate, 7 mM EDTA, pH 7.5) and 5% glycerol and electrophoresed at 200 V, 4 C for 1 h in TAE retarding buffer; then the gel was fixed, dried and the free and complexed RNA was measured in a Phosphorimager (Arias et al. 2005).

Crystallization and Structure Solution of 3D(V173I)
Purified 3D(V173I) was stored in a buffer containing Tris-HCl (50 mM pH 8.0), NaCl (500 mM), DTT (0.8 mM), EDTA (0.8 mM), and glycerol (8%), at a concentration of 5 mg/ml. Prior to the crystallization experiments, the polymerase was incubated overnight at 4 C with oligonucleotide 5 0 ACGFuGGGCCC 3 0 (NWG-Biotech) (oligonucleotide with FU-monophosphate at position 4), in an equimolar proportion, in presence of 2 mM MgCl 2 . Crystals were obtained at 20 C by the hanging drop vapor diffusion method in limbro plates, using a precipitant/well solution containing 36% PEG 4000, 0.2 M ammonium acetate, 0.1 M MES[2-(N-morpholino) ethanesulfonic acid] pH 6.0, and 4% c-butyrolactone. Two distinct crystal forms, which were suitable for X-ray analysis, were visible between 4 and 5 days after solution preparation. In order to obtain the ternary complexes, crystals were soaked for 6 h in a harvesting solution containing the crystallization buffer, 20 mM ATP and 5 mM MnCl 2 . Crystals were then transferred to a cryoprotecting solution, containing 20% glycerol in the crystallization buffer immediately before flash freezing in liquid nitrogen.
Two different data sets were collected at 100 K, using synchrotron radiation at ALBA-CELLS, beamline XALOC on a Pilatus 6M DECTRIS detector (k ¼ 0.979 Å ). Trigonal crystals, belonging to the space group P3 2 21 diffracted to 2.4 Å resolution and tetragonal crystals, space group P4 1 2 1 2, diffracted to 2.6 Å resolution. All data were processed with XDS (Kabsch 2010) and internally scaled with SCALA [CCP4i, Potterton et al. 2003]. Data collection statistics are given in supplementary table S2, Supplementary Material online.
The initial maps for the tetragonal crystal structures were obtained after rigid-body fitting of the coordinates of the isolated wt polymerase [crystallized in the tetragonal P4 1 2 1 2 (PDB id. 1U09)] to the new unit cells, using the program Refmac5 [CCP4i, Murshudov et al. 1997]. Initial maps for de la Higuera et al. the trigonal crystal form were obtained following the same procedure but using the P3 2 21 coordinates of 3Dwt (PDB id. 1WNE) as starting model (supplementary table S2, Supplementary Material online). In both structures, the weighted 2jFoj -jFcj and jFoj -jFcj difference maps clearly allowed the rebuilding of the mutated residue and other regions presenting conformational changes as a consequence of the mutation. Several cycles of automatic refinement, performed with Refmac5 (Murshudov et al. 1997), were alternated with manual model rebuilding using Coot (Emsley et al. 2010

Selection of a FMDV Mutant Displaying Increased Fitness in the Presence of 5-FU
FMDV-RP is a ribavirin-resistant population whose mutant spectrum included mutations that led to amino acid replacements M296I and P169S in the polymerase (3D) (at 100% and 50% dominance, respectively; Agudo et al. 2010; Sierra et al. 2007). To investigate a possible cross-resistance to FU, FMDV-RP was subjected to 15 passages in the absence or presence of FU (50 or 200 lg/ml) at a MOI of 0.05-0.1 PFU/ cell ( fig. 1A). Sequencing of the 3D-coding region indicated dominance (presence without alternative residues), at passage 15 in the presence of 200 lg/ml FU, of a new 3D substitution, V173I in the polymerase motif F (resulting from transition G7126A in the 3D-coding region), and loss of replacements M296I and P169S ( fig. 1B). To determine if V173I alone conferred FU resistance, FMDV-3D(V173I) was prepared as described in Materials and Methods, and its response to FU compared with the response of FMDV-wt. FU displays a dual mutagenic and inhibitory activity on FMDV (Agudo et al. 2008) but FMDV-3D(V173I) did not exhibit a significant resistance to the inhibitory effect of FU (supplementary fig. S1, Supplementary Material online). The lack of effect on the inhibition by FU establishes a distinction between polymerase residues that can affect inhibition by FU and those that can affect its mutagenic properties. Growth-competition fitness assays indicated a fitness cost inflicted by replacement V173I in 3D in the absence of FU but a FU concentration-dependent fitness gain in the presence of FU ( fig. 1C and supplementary  fig. S2, Supplementary Material online); the selective strength (fitness þdrug /fitness -drug ) was 2.23 for substitution V173I with the maximum concentration of FU tested (400 mg/ml). Thus, 3D replacement V173I in FMDV 3D confers a selective advantage to the virus when it replicates in the presence of FU.

Mutation Types and Mutant Spectrum Complexity of wt and Mutant Virus
We investigated how replacement V173I in 3D might affect quasispecies complexity in response to FU. Comparative analyses of the mutant spectra showed that passages in the presence of FU increased the mutation frequency for both FMDV-wt and FMDV-3D(V173I) (table 2; Pearson's chisquared test: P < 0.0001), resulting in a larger heterogeneity of the sequences in the mutant spectra. Comparison of mutation types and predicted amino acid substitutions indicated a predominance of transitions and non-synonymous mutations (supplementary tables S3 and S4, Supplementary Material online). The proportion of amino acid replacements with negative PAM 250 "(point accepted mutation matrix)" values ranged from 8% to 26% suggesting a good average acceptability of the amino acid substitutions in the polymerases of the viruses sampled from the mutant spectrum of the 3D-coding region; amino acid substitutions were found along the 3D-coding region (supplementary fig. S3, Supplementary Material online). The most salient difference between the two viruses regarding the mutation types was the proportion of A ! G plus U ! C relative to all mutations scored. For FMDV-wt such proportion was 0.67 in the mutant spectrum of the virus passaged 10 times in the absence of FU, and 0.84 in the virus passaged in the presence of FU (P ¼ 0.04; Fisher's test), a bias diagnostic of FU mutagenesis (Agudo et al. 2008;Grande-Pé rez et al. 2002). Instead of an increase, a decrease of the mutational bias was observed in the mutant spectrum of mutant FMDV-3D(V173I): 0.76 in the absence of FU versus 0.73 in its presence (P > 0.05; Fisher's test). The different response of the two viruses to FU can be expressed by the mutation type ratio normalized to the frequency of each nucleotide (named transition bias in table 2). The transition bias in the presence of FU increased 1.9-fold for FMDV-wt, and decreased 1.4-fold for FMDV-3D(V173I). Thus, the analyses suggest that the mutant polymerase tends to counteract the mutational bias toward A ! G and U ! C produced by FU mutagenesis.

Kinetics of Nucleotide Incorporation by wt and Mutant Polymerase
The difference of composition of mutant spectra generated during the replication of FMDV-wt and FMDV-3D(V173I), prompted us to compare nucleotide incorporation properties of the corresponding purified polymerases. To explore the enzymatic basis and to confirm the restriction of A ! G and U ! C transitions by FMDV-3D(V173I) in the presence of FU, the wt (3Dwt) and substituted 3D(V173I) polymerases were expressed, purified and studied biochemically and structurally. incorporation of FUMP (k pol /K d,NTP measures the catalytic efficiency of the reaction, taking into account the velocity of nucleotide incorporation and the affinity of the enzyme for the substrate). This gives a 3.2-fold higher selectivity (incorporation of UMP relative to FUMP) for the mutant than for the wt enzyme, contributed mainly by an increase in the dissociation of FU from the elongation complex resulting in a restriction of FUMP incorporation ( fig. 2B; table 3). This result is in agreement with the reduction of transitions A ! G and U ! C observed in the mutant virus population compared with the wt when replicating in the presence of FU. In contrast, the two enzymes did not differ significantly regarding the relative incorporation of FU and C opposite G in the template and showed a strong preference for C incorporation (table 3;  The incorporation of AMP and GMP opposite FU in the template indicated a 14-fold higher selectivity for A versus G as a result of substitution V173I, largely a consequence of increased dissociation of G relative to A ( fig. 2D; table 3). Thus, the parameters of presteady state kinetics suggest that substitution V173I attenuates the mutagenic activity of FU by forcing recognition of the cognate nucleotide. Transitions G ! A and C ! U (which are the mutation types opposite to those favored by FU) will occur when FUMP is incorporated instead of CMP, directed by a G in the template. The kinetics parameters indicate that this misincorporation should not differ significantly between the two polymerases (supplementary fig. S5C and D, Supplementary Material online). However, for a FU in the template (originated by misincorporation opposite G) to give rise to a mutation it should direct the incorporation of A. 3D (V173I) exhibits higher selectivity than 3Dwt for the incorporation of A instead of G ( fig.  2C and D; table 3). The end result is that 3D(V173I) will tend to generate a higher frequency of G ! A transitions than 3Dwt thus compensating the tendency of FU to induce A ! G and U ! C mutations when replication is catalyzed by 3D-wt. Discrimination of nucleotides is greater when FU has to be incorporated rather than copied as template residue (table 3), suggesting that FU is more likely to induce mutations when it serves as template. Thus, the kinetic parameters determined with purified polymerase are consistent with the mutation repertoire quantified in the mutant spectra of FMDV-wt and FMDV 3D(V173I) passaged in presence of FU (table 2).
Substitution V173I did not affect the VPg uridylylation activity (supplementary fig. S6, Supplementary Material online), thereby establishing a structural and functional distinction between the mutagenic and inhibitory activities of FU.

Structural Modifications of the Polymerase with Replacement V173I
To investigate the structural modifications associated with the changes in nucleotide affinity, 3D(V173I) was crystallized and solved in two different crystal forms at 2.6 Å (P4 1 2 1 2 crystals) The mutant spectrum analysis involved nucleotides 6,610-7,980 in the 3D-coding region; the number of individual clones analyzed is given in parenthesis. c Mutation frequency [min] is the number of different mutations found in the mutant spectrum (relative to the sequence of the parental, unpassaged FMDV clone), divided by the total number of nucleotides sequenced (given in the second column); values in parenthesis are the total number of unique mutations. The statistical significance of the differences between values is the following: FMDV-wt, no FU versus FU: v 2 ¼ 42.7, P ¼ 6.5 Â 10 À11 ; FMDV-3D(V173I), no FU versus FU: v 2 ¼ 47.6, P ¼ 5.0 Â 10 À12 ; no FU, FMDV-wt versus FMDV-3D(V173I): v 2 ¼ 0.79, P ¼ 0.37; FU, FMDV-wt versus FMDV-3D(V173I): v 2 ¼ 0.14, P ¼ 0.71. (In all cases, Pearson's chi-squared test, df ¼ 1). d Mutation frequency [max] is the number of total mutations found in the mutant spectrum (relative to the sequence of the parental, unpassaged FMDV clone), divided by the total number of nucleotides sequenced (given in the second column); values in parenthesis are the total number of mutations. The statistical significance of the differences between values is the following: FMDV-wt, no FU versus FU: v 2 ¼ 84.1, P ¼ 4.5 Â 10 À12 ; FMDV-3D(V173I), no FU versus FU: v 2 ¼ 62.2, P ¼ 2.9 Â 10 À15 ; no FU, FMDV-wt versus FMDV-3D(V173I): v 2 ¼ 0.95, P ¼ 0.33; FU, FMDV-wt versus FMDV-3D(V173I): v 2 ¼ 4.2, P ¼ 3.9 Â 10 À2 . (In all cases, Pearson's chi-squared test, df ¼ 1). e Shannon entropy (S) is calculated by the formula S ¼ À [S i (p i Â ln p i )]/ln N, in which p i is the frequency of each sequence in the quasispecies, and N is the total number of sequences compared. It is a measure of heterogeneity of the sequences sampled in the mutant spectrum (S ¼ 0, all sequences are identical; S ¼ 1, each sequence is different). f Transition bias is defined here as the ratio between the frequency of transitions A ! G plus U! C, relative to the frequency of transitions G-A plus C-U: and 2.4 Å (P3 2 21 crystals) resolution, respectively (supplementary table S2, Supplementary Material online). The analysis of electron density maps allowed the tracing of the mutated and surrounding residues that were omitted from the initial model to eliminate model bias. The structure of 3D(V173I) was similar to 3Dwt (fig. 3). The superimpositions of all 476 amino acids of the two enzymes showed root mean square deviation (rmsd) values of 0.33 Å and 0.62 Å for space groups P4 1 2 1 2 and P3 2 21, respectively. Despite overall structural similarity, some significant changes were noted in the polymerase Nterminus (residues 14-18), and loop b9-a11 (residues 293-302, at the N-terminal side of polymerase motif B). In both space groups, the electron density indicates high flexibility around residues H14-K18 and this region could not be traced in the trigonal P3 2 21 crystals. In the P4 1 2 1 2 crystals the electron density was compatible with the tracing of two alternative conformations of the M16-R17 dipeptide main and side chains, with $50% occupancy each one. The first orientation was coincident with the position found for this region in the wt enzyme (Ferrer-Orta et al. 2004), with the side chain of R17 oriented towards the template channel, in a position accessible to the template RNA. In the second position, the main and side chains of residues M16 and R17 appeared reoriented, with the R17 side chain pointing to the polymerase interior, and interacting with residues N41 and Y285. Interestingly, this second orientation is coincident with that previously observed in the M296I, P44S, P169S (SSI) polymerase mutant, resistant to high concentrations of ribavirin (Agudo et al. 2010).
Two types of changes were observed in loop b9-a11 depending on the crystal type. While in the trigonal crystals the G299-S301 stretch was disordered and not visible in the electron density, the tetragonal crystals showed a 2.5 Å movement of G299 to S301 from their standard position in the wt enzyme, packed against the fingers (Ferrer-Orta et al. 2004), to a configuration where the loop protrudes towards the catalytic cavity in an orientation that would disturb the binding of the RNA template. Cocrystallization with sym/sub RNAs indicated that the bound RNA was disordered in the crystals, precluding a comparison with complexes with 3Dwt. A soaking of the crystals in 20 mM ATP resulted in the presence of the nucleotide bound to the NTP binding channel. However, only the triphosphate moiety of this incoming nucleotide was clearly seen in the electron density, interacting with the motif F residues R168 and K172, neighbor of the substituted residue I173 (fig. 3). Thus, localized structural changes are associated with altered nucleotide incorporation tendencies of 3D(V173I).

Probing Targets of Selection with Predictive Algorithms and Experimental Determinations
The biochemical and structural results prompted us to explore the possible consequences of the mutational effects on fitness  Nucleotidic substrate used as a triphosphate for incorporation at the second position of the sym/sub; the exception is sym/sub-FuG for which incorporation was measured at the first position.  2.-Presteady-state kinetics of nucleotide incorporation into sym/sub-UA and sym/sub-FuG by 3Dwt and 3D(V173I). (A) Sequence of 5 0 -endlabeled, annealed sym/sub-UA. Arrows indicate the template residue at which nucleotide incorporation is measured. (B) 3D (0.5mM active sites) was preincubated at 37 C with sym/sub-UA (0.5 mM duplex) and ATP (10 mM) for 900 s to allow the formation of 3D-RNA product complex, with AMP incorporated at the first position. The 3D-RNA product was then mixed with the indicated concentration of UTP or FUTP using a rapid chemical quenchflow apparatus, and reactions were quenched by the addition of EDTA (0.3 M). Independent time points at 0, 0.01, 0.05, 0.1, 0.25, 0.5, 1, and 2 s were taken for each nucleotide concentration tested. Time courses at fixed nucleotide concentrations were fit to an exponential curve to obtain the observed rate constant for nucleotide incorporation at the second position, k obs . The observed rate constants were then plotted as a function of nucleotide concentration, and the data were fit to a hyperbola to obtain k pol and K d,app . Note that the scale at ordinate is different in the two panels. (C) Sequence of 5 0 -end-labeled, annealed sym/sub-FuG (Fu means FU). Arrows indicate the template residue at which nucleotide incorporation is measured (FU is written F). (D) Incorporation of AMP or GMP; no nucleotide was preincubated with 3D and RNA. Procedures are those described in (B). For the assays with GTP, the concentration of 3D used was 1 mM active sites, and no rapid chemical quench-flow apparatus was needed. Duplicate samples at 0, 10, 30, 60, 120, 300 and 900 s time points were taken for each nucleotide concentration tested. Note that the scale in ordenate is different in the two panels. Procedures are further detailed in Materials and Methods. The left and right panels show two views of 3D(V173I) protein rotated by 90 . The polymerase is depicted in blue ribbons with substituted amino acid I173 shown in sticks in red. The triphosphate moiety of the bound ATP and the motif F contacting side chains R168 and K172 are also shown as sticks and explicitly labeled. (B) Stereoview of rA-weighted jFoj À jFcj electron density map (contoured at 3r) around the mutated residue I173. The substituted residues and surrounding amino acids were omitted from the phasing model. The model is placed inside in ball and stick representation and colored in atom type code. (C) Stereoview of the structural changes around the mutated I173 residue. 3D(V173I) polymerase is shown in white (with I173 depicted in red) and the superimposed 3Dwt in cyan (TP refers to the ATP triphosphate moiety). (D) Superimposition of the de la Higuera et al. at the molecular level. This objective has as difficulties that in RNA viruses both the RNA and the proteins are part of the phenotype, the genomic base composition can affect several steps in the virus life cycle, proteins are multifunctional, and selection may act on complex mutant spectra (Martin and Domingo 2008), and not necessarily on one or several dominant sequences [these points have been reviewed in Domingo 2016]. Since selection acts on phenotypes we sought to explore viral traits that might have driven the selection of FMDV with V173I in 3D upon virus replication in the presence of FU. According to the codon composition of the entire FMDV genome and of the 3D-coding region (Toja et al. 1999), the fraction of non-synonymous mutations is significantly larger for U ! C than C ! U mutations (Pearson's chisquared test: df ¼ 1, v 2 ¼ 64.0, P < 10 À14 ), and for A ! G than G ! A mutations (Pearson's chi-squared test: df ¼ 1, v 2 ¼ 36.9, P < 10 À8 ; fig. 4A). Furthermore, the impact on 3D's folding free energy (DDG) of every possible nonsynonymous mutation was evaluated using PoPMuSiC and DeltaGREM (described in Materials and Methods). Both methods indicate that, on average, U ! C mutations have a much stronger destabilizing effect on the protein than C ! U mutations (Welch's two-tailed t-test, for the comparison of the mean DDG as predicted by PoPMuSiC: df ¼ 321, t ¼ 17.5, P < 10 À15 ). Similarly, A ! G mutations are predicted to be on average more destabilizing than G ! A mutations, although the differences are smaller (Welch's two-tailed t-test, df ¼ 521, t ¼ 0.88, P ¼ 0.38; fig. 4B). The relevance of the stability predictions is supported by evidence of stronger selection against mutations with a larger predicted destabilizing effect in the mutant spectra of FMDV and FMDV-3D(V173I; supplementary table S5, Supplementary Material online).
To further investigate the importance of protein stability as a target of selection, we expressed 10 mutants of the 3D polymerase and evaluated their melting temperatures (T m ) via ThermoFluor experiments (table 4). Most mutants show a T m decrease of at least 2 C, three of them have a T m lower than the temperature at which the passages were performed (37 C), and three more were found to aggregate. In addition, we applied PoPMuSiC and DeltaGREM to estimate, for each mutant, the change in folding free energy at room temperature (DDG). An important reduction of the stability is predicted for 6 out of 10 mutants, by either or both computational methods (table 4). The genomes isolated from the mutant spectra are a surviving subset that eluded negative selection. Despite this bias, both the ThermoFluor experiments and the DDG predictions indicate that clones produced under FU-enhanced mutagenesis may present important alterations of 3D protein stability, with likely consequences for viral fitness. This result is in agreement with previous reports of decreasing fitness in subsets of genomes that remain viable in the process of ribavirin mutagenesis (Arias et al. 2013).
The average predicted RNA secondary structure stability was increased in the mutants displaying a larger content of the strongly interacting bases G and C (table 4 and    The numbering of FMDV genomic residues is that described previously (Escarmis et al. 1999). The mutant 3Ds were expressed from plasmid pET-28a 3Dpol that included 3D with the amino acid substitutions given in the third column. Expression in E. coli was performed previously described (Ferrer-Orta et al. 2004) except that the temperature at which E. coli was grown had to be decreased to 20 C to obtain 3D in sufficient yield and monomeric (M) form.
i The yield is expressed as the total amount of 3D (in mg) obtained from an E. coli culture of 500 ml grown as previously described (Ferrer-Orta et al. 2004). j The forms observed are indicated. D in parenthesis in the triple P44S, P169I, M296I triple mutant means preponderance of the monomeric (M) form. de la Higuera et al. two-tailed t-test: df ¼ 59.2, t ¼ 2.2, P ¼ 0.03), while the opposite trend holds for mutations found in virus replicated in the presence of FU (supplementary fig. S7, Supplementary Material online). Besides the individual effect of each mutation, the increased mutation frequency resulting from FU mutagenesis markedly affects the average stability of 3D in each mutant spectrum. In consequence, according to DeltaGREM, the fraction of clones in which the stability of 3D is seriously altered (> 2.0 kcal/mol) increased from 4% for FMDV-wt replicated without FU to 32% with FU (Fisher's test: P ¼ 0.01), and from 12% to 23% for FMDV-3D(V173I) (Fisher's test: Interestingly, the average change of hydrophobicity is not significantly different from zero in the FMDV mutant spectrum passaged in the absence of FU (mean DH ¼ À0.007; One-sample t-test (two-tailed): df ¼ 27, t ¼ À0.2, P ¼ 0.85), despite the larger frequency of U ! C mutations (supplementary fig. S7, Supplementary Material online). This suggests that increasing hydrophobicity may be less deleterious than decreasing it, and/or that compensating hydrophobicity changes may be profitable. The fine balance between mutations increasing and decreasing hydrophobicity is compromised in presence of FU, resulting in a significant loss of hydrophobicity (mean DH ¼ À0.11; Welch's two-tailed t-test: Thus, experimental determinations with mutant viral polymerases and predictive algorithms of the deleterious effects of amino acid substitutions suggest decreased functionality of the viral polymerase encoded by components of a mutant spectrum enriched in G and C residues.

Discussion
The control of viral diseases using antiviral agents constitutes an important challenge. Lethal mutagenesis as an antiviral strategy arose as an application of error threshold concept included in quasispecies theory (Domingo and Schuster 2016b;Eigen and Schuster 1979;Perales and Domingo 2016;Schuster 2016;Tejero et al. 2016). The interest in applying lethal mutagenesis to combat viral disease was reinforced with the realization that some antiviral agents in use or in advanced clinical investigations such as ribavirin or favipiravir (T-705) may act as lethal mutagens against some viruses (Arias et al. 2014;Baranovich et al. 2013;Crotty et al. 2000;de Avila et al. 2016;Graci and Cameron 2006;Holland et al. 1990). Unfortunately the pertinacious adaptability of viral quasispecies has manifested also in the form of selection of viral mutants resistant to mutagenic nucleotide analogues (Agudo et al. 2010;Borderia et al. 2016;Ferrer-Orta et al. 2010;Pfeiffer and Kirkegaard 2005;Vignuzzi et al. 2006;Zeng et al. 2013).
The present study has addressed the selection of a FMDV mutant with amino acid substitution V173I in its polymerase that confers the virus a selective advantage during replication in the presence of FU. The selective advantage has been evidenced by the gradual increase in the frequency of amino acid substitution V173I in 3D when FMDV-RP was passaged in the presence of FU ( fig. 1B), and the consistent FU concentrationdependent fitness increase of FMDV-3D(V173I) relative to FMDV-wt (figs . 1C and supplementary fig. S2, Supplementary Material online). Under the passage conditions of the fitness assays, neither reversion of V173I upon passage of FMDV-3D(V173I) in absence of FU, nor acquisition of V173I by FMDV-wt in the presence of FU, were observed (see Materials and Methods). The selective advantage of FMDV-3D(V173I) over FMDV-wt was associated with attenuation of the mutational bias in favor of A ! G and U ! C inflicted by FU upon FMDV. The increase of A ! G and U ! C transitions induced by FU was observed in previous experiments with FMDV-wt (Pariente et al. 2001;Sierra et al. 2000). The altered mutational preferences observed in the mutant spectra of FMDV-wt and FMDV-3D(V173I) are in agreement with the kinetics of nucleotide incorporation by the corresponding mutant polymerases. The mutant polymerase exhibits a compensatory increase of the mutation frequency from C to U in the presence of FU while maintaining a broad mutant spectrum in progeny virus, a hallmark of virus adaptability (Pfeiffer and Kirkegaard 2005;Vignuzzi et al. 2006). At the structural level (supplementary table S2, Supplementary Material online), the localized structural changes in 3D are similar to those previously found with other FMDV 3D mutants, which display decreased or increased incorporation of ribavirin-triphosphate (Agudo et al. 2010;Ferrer-Orta et al. 2010, 2015. These comparisons have revealed that very subtle structural adjustments may lead to either an increase or a decrease of nucleotide analogue incorporation during viral RNA synthesis. The codon composition of an RNA viral genome is the result of a long evolutionary history, often in response to the host codon composition to ensure modulation of the rate of viral protein synthesis for virion production (Bosch et al. 2010;Mueller et al. 2010). Seeking to explore additional phenotypic traits that may have constituted a target for the selection of FMDV(V173I), we have applied some predictive and biochemical assays to document multiple effects at the RNA and protein level both in general for the FMDV genome, and for a sample of the mutants found in the relevant FUmutagenized viral populations. It should be stressed that mutants detected among the molecular clones from mutagenized FMDV correspond to replication competent viruses, while the clones with the most deleterious mutations cannot be detected. Nevertheless, measurements with 10 mutant proteins revealed that many of them experienced substantial reductions of folding stability (table 4). Our predictive computations suggest that the effects of mutations on protein stability can be modulated by the mutation bias, in accordance with previous studies (Bastolla et al. 2006;Mendez et al. 2010). In the genetic background of FMDV, the mutations from U to C and from A to G, which are induced by FU, have a more deleterious effect than the mutations from C to U and from G to A, respectively. This is due to both a smaller fraction of synonymous mutations, and a greater impact of non-synonymous mutations on protein stability (figs. 4 and supplementary fig. S7, Supplementary Material online). These two effects may thus contribute to the fitness advantage displayed by FMDV-3D(V173I) in presence of FU, since the 3D(V173I) enzyme tends to diminish the frequency of G and C introduced in the genome (tables 2 and 3). The assumption that synonymous mutations are less likely to be deleterious is supported by the significantly larger than expected frequency of such mutations in the mutant spectra of the viral populations replicating in presence of FU, providing evidence of a stronger selection against non-synonymous mutations (table 2). Similar evidence indicates that the selection is even stronger against non-synonymous mutations with a large destabilizing impact (supplementary table S5, Supplementary Material online). Additional support for the validity of protein stability predictions used in the present study comes from a previous report of a significant correlation between PoPMuSiC predictions and protein fitness on a set of $1,000 mutants of b-lactamase TEM-1 (Jacquier et al. 2013).
The relationship between mutation bias and protein stability may be, at least in part, explained by the hydrophobicity of the different amino acid types (Bastolla et al. 2006;supplementary fig. S7, Supplementary Material online). Indeed, codons containing U at second position correspond to hydrophobic residues; on average, U ! C mutations decrease hydrophobicity, while C ! U mutations increase it. At the RNA level, U to C mutations tend to increase the secondary structure stability of the RNA, with possible influence on the stability of the viral genome and on protein expression (Plotkin and Kudla 2011).
Although protein stability is not the sole determinant of viral fitness, the predicted fractions of clones in which the stability of 3D is altered (supplementary fig. S7, Supplementary Material online) are qualitatively consistent with the comparative fitness assays ( fig. 1C and supplementary fig. S2, Supplementary Material online); a quantitative investigation would require computational simulations based on more precise estimations of the mutational frequencies and biases. The complementary virological, enzymological, structural, and computational results presented here support the conclusion that the V173I substitution in 3D restores a required mutational balance, opening new perspectives on the role of mutational biases in viral and cellular evolution. For virology the results reinforce mutation modulation as a mechanism to escape lethal mutagenesis by purine and pyrimidine analogues, further illustrating the resources that RNA viruses have to confront hostile environments. For cell biology the results raise the intriguing possibility of external mutational activities affecting the replicative machinery over extended periods of evolution. These results support the view that mutation biases are under selective control, and they are consistent with the proposal, formulated in the context of microorganism evolution, that mutation biases can evolve through changes in the replicative machinery that are selected for their effects over extended periods of time (Mendez et al. 2010).

Supplementary Material
Supplementary data are available at Genome Biology and Evolution online.