The mutations that arise from adaptive passage map to receptor binding sites in a cell-type specific manner

Although the HeLa and A549 cell passage series presented similar mean genetic diversity at passage 40 (Fig ​(Fig1E1E–1G, P = 0.314), we mined the deep sequence data for all minority variants above 1% frequency that might explain adaptation in each condition (Table 1). In both cases, several mutations recurring in multiple replicates mapped to residues known to be part of the CAR receptor and DAF co-receptor footprints [16]; however, the distribution of these mutations was different for the HeLa- and A549-adapted viruses. In HeLa cells, the most abundant variants involved mutations strictly in the CAR footprint (VP1-259, 2–7% frequency), in the CAR/DAF shared footprint (VP2-138, 20–80%), and in the DAF footprint (VP3-234, 2–73% frequency; VP3-63, 3–14% frequency). In contrast, viruses passaged in A549 cells presented no variants strictly involved in the CAR footprint. Instead, in addition to CAR/DAF residue VP2-138 (2–81% total), a larger cluster of exclusively DAF-footprint variants were identified (VP3-63, 3–35% frequency; VP3-234, 3–82% frequency and VP1-271, 2–13% frequency). Strikingly, the most abundant mutation in every replicate passaged in A549 cells, and not observed in HeLa cell passage, was VP3-76 (61–95% total). This residue is not known to participate in either footprint.

Permissive and less permissive cell environments present different expression profiles of the primary and secondary receptors

Considering these CAR/DAF-specific differences, we hypothesized that expression levels of these molecules must vary between these cells. There were no significant differences in CAR expression in either cell by both flow cytometry (Fig 2A) or Western blot (Fig 2C). On the other hand, while DAF was very highly expressed in HeLa cells (Fig ​(Fig2B2B and ​and2D),2D), expression was one order of magnitude lower in A549 cells by flow cytometry (Fig 2B) and barely detectable by Western blot (Fig 2D). Normalized quantification of Western blot signals revealed that DAF expression was 4-fold higher in HeLa than in A549 cells (Fig 2E). To characterize the localization of CAR and DAF in these cell types, confocal microscopy was performed. In HeLa cells, CAR (Fig 2F) and DAF (S1 Fig) were expressed throughout the surface of the cell. Interestingly, in A549 cells CAR predominantly localized at cell-to-cell contacts (Fig 2G), while DAF was diffused throughout the surface of the cell (S1B Fig), albeit at considerably lower levels than in HeLa cells (see also, S1 and S2 Movies).

Open in a separate window

Fig 2

Differential expression of CAR and DAF in HeLa and A549 cells.

(A-B) Expression of CAR (A) and DAF (B) by flow cytometry in HeLa cells (solid line) or A549 cells (dashed line). Fluorescence of cells in absence of antibody (black) or in presence of CAR- or DAF-specific antibody (red) is shown. Figures are representative of three independent experiments. (C-D) Expression of CAR (C) and DAF (D) by Western blot analysis of HeLa or A549 cell extracts, revealed by CAR-, DAF- and GAPDH-specific antibodies. Molecular weight markers are indicated on the left, super-exposed section of each figure is included to reveal the presence of CAR and the two isoforms of DAF in each sample. Figures are representative of two independent experiments. (E) The relative expression of CAR and DAF with respect to GAPDH was quantified by imaging quantification using ImageJ. Localization of CAR, in red, in HeLa (F) or A549 (G) by confocal microsocopy. Nuclear (blue) staining was done with DAPI.

Population deep sequencing of longitudinal samples identifies multi-step adaption to novel host environments

It has been shown that adaptation is often a multi-step process, especially in asexual populations like viruses. Usually, until the mutation with the largest beneficial effect becomes fixed, secondary beneficial mutations with less effect are unable to compete, thereby rendering the adaptation process sequential [17–19]. Given the different frequencies of minority variants observed at passage 40, we examined the patterns of emergence and sequential adaptation to novel environments by deep sequencing every second passage in the A549 cell series (Fig 3 and S1 Dataset). A number of step-wise trends were observed across the replicates. In all six replicates, the VP3-E76G mutation was the first mutation to emerge (already above 0.1% at passage 1) and peak by passage 11 to become the most abundant mutation in the population (between 65 and 95%). Furthermore, the increase in frequency of E76G over the first 11 passages (Fig 3) correlated with the increases in fitness over the same period (Fig 4A), underscoring a considerable contribution of this single mutation to population fitness. Indeed, the E76G variant was generated by reverse genetics and found to confer significantly enhanced fitness in A549 cells (Fig 4B). Although the VP3-76 residue was not previously identified as a contact between virus and receptor, this residue maps to the icosahedral three fold region of the capsid (Fig 4C). In the three dimensional virus-receptor structure, there is an interaction between the C-terminal 6-His tags of three symmetry related DAF molecules [20]. This interference by the His-tag may have restricted the normal interactions of DAF SCR3/4 with the virus surface, masking a potential role for residue 76. In a biolayer interferometry assay using BLItz technology [21] (Fig 4D), E76G bound DAF with the same affinity as wildtype virus, but no increased binding could be detected.

Open in a separate window

Fig 3

Population sequencing reveals multi-step adaptation of residues in the DAF-specific footprint during passage in A549 cells.

(A-F) The virus populations were purified and deep sequenced for each replicate every two passages (p1, 3, 5, 7, etc...) across the P1 region coding for the four capsid proteins. The variant frequencies were determined by the ViVAN pipeline and all statistically significant mutations that showed a trend of increasing frequency were plotted in linear scale (left panels) and log scale (right panels). Each set of linear and log graphs represents the same data from one replicate, all six replicate passage series presented in this work are shown.

Open in a separate window

Fig 4

Mutation E76G confers fitness advantage in A549 cells.

(A) Increasing relative fitness (y-axis) of replicates a-f, over the initial stage of passage series (histograms indicate values for passage 5, 7, 9, 11), mean values with SEM, n = 3. (B) Relative fitness of the CVB3-parental or E76G mutant generated from an infectious cDNA. The relative fitness (y-axis) is the ratio of the viral RNA quantification at 24h and 0h, mean values with SEM, n = 3. (C) The surface rendered map of three CVB3 protomers shown at the icosahedral three-fold axis of the virus (grey). Three symmetry related bound molecules of DAF (blue) are shown with a transparent surface rendering to allow visualization of the location of VP3 76 (red) on the virus surface directly below the DAF. (D) Binding to DAF was measured by biolayer interferometry by dipping a DAF loaded sensor into aliquots of CVB3-parental (black) or CVB3-E76G (red).

Despite the clear fitness advantage conferred by the E76G mutation alone (approximately 10-fold increase with respect to wildtype), it could not account for the more significant fitness increases observed in the passaged populations, particularly at late passage stages (nearly 100-fold increases). Following the fixation of this first mutation, the frequency of VP3-E76G plateaued until a second increase after passage 27, in each of the six replicates, coinciding with the emergence of different combinations of mutations, all mapping to the DAF footprint (Fig 3). Closer examination of the patterns of emerging mutations revealed potential synergistic and antagonistic epistasis among variants. Despite the advantage afforded by extreme depth in coverage, the Illumina technology used here set a limit of 69 nucleotides for read length. In principle, this limitation forfeits the possibility of linking more distant mutations and identifying haplotypes. However, the longitudinal deep sequence data revealed that several mutations increase with parallel kinetics and frequency, suggesting that: i) each mutation appeared in individual genomes and the variants were selected as a group and/or ii) the mutations are accumulating in the same genome, resulting in the selection of single haplotypes.

To distinguish between the above cases, phylogenetic trees describing the mix of haplotypes at each passage were inferred from the longitudinal data using maximum likelihood estimation in a Bayesian model. The best-fit, predicted haplotypes were then generated by reverse genetics and their relative fitness values were measured. As expected, this analysis illustrated the quick rise of the E76G genotype that continues to coexist with, although generally dominating over, the original WT genotype (Fig ​(Fig5A5A–5F), with a significant increase in relative fitness (Fig 5G). For residue VP3-234, mutations were predicted to arise on either the WT background (Fig ​(Fig5C5C and ​and5F)5F) or the E76G genotype (Fig ​(Fig5A5A and ​and5E).5E). Mutations in residue VP2-138 were also predicted to arise on either WT (Fig 5B) or E76G (Fig ​(Fig5C,5C, ​,5D5D and ​and5F)5F) genotypes at approximately the same time in the passage series as residue 234. However, their frequencies on the WT background tended to remain lower than on the E76G background where they seemed to be better tolerated. To confirm these modelled predictions, we generated each mutation on both backgrounds and measured the relative fitness. Indeed, the D138G mutation alone on the WT background conferred nearly a ten-fold drop in fitness and the Q234K mutation conferred up to 100-fold drop in fitness (Fig 5G), while the fitness costs of these mutations on the E76G background resulted in neutral fitness relative to WT. The data suggest that epistasis between these mutations and the E76G mutation rescues the fitness of these double variants and permits their positive selection in the viral population. Interestingly, the D138G and Q234K variants seemed to entirely exclude one another from the population (Fig ​(Fig5A,5A, ​,5B,5B, ​,5D5D and ​and5E),5E), and the MLE analysis suggested that if both are present in the population, they must occur on separate genomes. We thus generated the E76G-D138G-Q234K haplotype and confirmed that the triple mutant bears a significant fitness cost, rendering this haplotype less fit than the original WT genotype (Fig 5G).

Open in a separate window

Fig 5

Maximum likelihood estimates of haplotypes present within the passaged virus populations and contribution of minority variants to population fitness.

MLE was performed using the frequency values of each mutation at each passage in the deep sequence data for replicates a-f (A-F). The trees indicate the presence of wildtype genotype (lowest line on tree) along with the predicted haplotypes as they emerge during the passage series (x-axis). Solid lines indicate haplotypes with highest MLE scores, dashed lines indicate alternative haplotypes that could exist, with lower scores. The thickness of grey area indicates the expected frequency of each haplotype in the population, according to deep sequence data. (G) Relative fitness of the wildtype (WT), single, double and triple variants, followed by combinations of single variants, followed by passaged population samples (p40 replicates A-F) and reconstituted quasipecies. Reconstituted quasispecies is composed of 50:30:10:10 of E76G, E76G+N63Y, E76G+D138G, and E76G+Q234K. Vertical dashed lines are placed to separate the aforementioned groups to facilitate the reader. Horizontal line is placed to facilitate reading of neutral fitness (value 0). The relative fitness (y-axis) is the ratio of the viral RNA quantification at 24h and 0h, mean values with SEM, n = 3. (H) A549 cells were transfected with in vitro transcribed infectious RNA corresponding to dawildtype (WT), single and double variants. At 8 hours post-transfection, the progeny virus was quantified by TCID₅₀ assay. No significant differences were observed between WT and variants (p = 0.203, 0.400, 0.365, 0.504, respectively, n = 5–6, two-tailed Mann Whitney test).

Finally, shortly after the appearance of residue 138 and 234 variants, a third mutation appeared during the passage series at residue 63 (Fig ​(Fig5A,5A, ​,5B5B and ​and5E).5E). Once again, computational modelling predicted that it too exists on the E76G background, and on separate haplotypes than the position 138 and 234 mutations. Furthermore, data from two replicates (Fig ​(Fig5A5A and ​and5E)5E) suggested that the E76G-N63Y double variant bears higher fitness than E76G-Q234K, as its frequency increased towards the end of the passage series as the other's decreased. Indeed, the residue 63 variants presented higher relative fitness than both the 138 and 234 variants on the E76G background (Fig 5G).

Although we could not determine whether residue 76 plays a role in receptor interaction directly, and the role of accompanying mutations are inferred from previously published studies, to rule out that these mutations impact virus fitness on activities downstream of receptor binding and entry, we transfected cells with in vitro transcribed RNAs corresponding to some of these variants and assayed virus production at 8 hours (before a new round of infection can occur). The data revealed that the fitness advantages (E76G, E76G-N63Y) and disadvantages (WT-Q234K, E76G-Q234K) observed during infection of cells (Fig 5G) do not appear to exist when the binding and entry step is bypassed (Fig 5H).

Viral titers by tissue culture infectious doses (TCID₅₀)

Ten-fold serial dilutions of virus were prepared in 96-well flat-bottom plates in DMEM. Dilutions were performed in octuplate and 100 μl of dilution were transferred to 10⁴ Vero cells plated in 100 μl of DMEM-10% NCS. After 5 days living cell monolayers were colored by crystal violet. TCID₅₀ values were determined by the Reed and Muensch method. No significant differences were observed between TCID₅₀ values when using Vero, HeLa or A549 cells as the cellular substrate.

Sequencing

5x10⁸ virion from passaged samples were RNA extracted and RT-PCR amplified by RT (Superscript III) and PCR (Phusion) using primers sets that covered the whole genome, in 3–4 kb fragments. For consensus sequencing, the resulting PCR products were purified, sequenced and analyzed using Lasergene software (DNAStar Inc). For deep sequencing, PCR fragments were purified via the Nucleospin Gel and PCR Clean-up kit (Macherey-Nagel) and total DNA was quantified by Nano-drop. PCR products were then fragmented (Fragmentase), linked to Illumina multiplex adapters, clusterized and sequenced with Illumina cBot and GAIIX technology. Sequences were demultiplexed by CASAVA with no mismatches permitted. Clipping was performed using the fastq-mcf tool, removing common adapter contaminants and trimming low quality bases (Phred<30). Clipped reads were aligned to the Coxsackie virus B3 Nancy sequence as reference with a maximum 2 mismatches per read, and no gaps, using BWA v0.5.9. Alignments were processed using SAMTools to obtain a pileup of the called bases at each position. An in-house pipeline, termed ViVAN (Viral Variant ANalysis) [28] was used to identify statistically significant variants above the background noise due to sequencing error, in every sufficiently covered site (>100x). Briefly, for each position throughout the viral genome, base identity and their quality scores were gathered. Each variant was determined to be true using a generalized likelihood-ratio test (used to determine the total number of minority variants) and its allele rate was modified according to its covering read qualities based on a maximum likelihood estimation. Additionally, a confidence interval was calculated for each allele rate. In order to correct for multiple testing, Benjamini-Hochberg false-discovery rate of 5% was set. The total allele rates passing these criteria, across the whole genome, were used to calculate the mean variation rates (diversity) at different passages. The variation rate at position i is defined as the proportion (F) of significant non-reference alleles (k) and is denoted V_i:

The region-wide variation rate is the averaged variation rate across all covered positions in the genome (denoted n):

Confocal microscopy

Hela and A549 cells were plated onto coverslips, fixed with 2% paraformaldehyde for 20 minutes at room temperature, and then washed with PBS. Before staining, non-specific staining was blocked by incubating the cells with 5% FCS and 0.05% saponin in PBS during 10 minutes. Staining antibodies were diluted also in this buffer, and it was used for washes between antibodies. Cells were incubated with either CAR (Santacruz) or DAF (Abcam) primary antibody, washed, stained with a secondary antibody coupled to the appropriate fluorophore and washed. Cells were analyzed using a Zeiss LSM-700 confocal microscope. 3D reconstruction of the images was performed using the Imaris software.

Western blot

One million cells of HeLa and A549 cells were lysed using a buffer containing 1% Triton-X and 1% sodium deoxycholate with protease inhibitors (Sigma). A fraction of the lysate was run for 1 hour on a 4–15% gradient gel (Biorad) on denaturalizing conditions. After the run, we performed the transfer to a nitrocellulose membrane. We washed the membrane with PBS-T (PBS 1X and 0.1% Tween-20) and blocked for 1 hour in PBS-T plus 5% milk. After the blocking we washed again with PBS-T and left overnight with each one of the antibodies, anti-CAR and anti-DAF (Santacruz Biotechnology). We washed the membrane and we added the secondary fluorescent antibodies (DyLight 680 and DyLight800 conjugated, Thermo Scientific) for 1 hour. We washed the membrane one last time and measure the fluorescence in the Odyssey system (Li-Cor).

Flow cytometry

In order to prepare cells for cell cytometry analysis, cells were washed twice with PBS 1x and trysinized, washed again twice in PBS. Cells were stained for 30 minutes with either CAR-PE or DAF-FITC antibodies (Millipore and Abcam) on ice. Unbound antibody was discarded and cells were washed again with PBS 1x. Cells were resuspended in 1% Parafolmaldehyde (PFA, Electron Microscopy Sciences) and kept in the dark for 15 minutes. After the incubation PFA was discarded and cells were resuspended in 200 μl of PBS 1x. Cells were kept at 4C until analysed. For each cell type used (HeLa and A549) specific instrument settings were set according to the size and complexity of the cell type, as well as antibodies fluorescence. Samples were analysed using the MACSquant flow cytometer (Miltenyi Biotec) using 96 well plates and obtaining 10,000 events per sample. Mock samples were also used in each plate to setup the baseline. Results were analysed using Flowjo software v10.

Fitness assays

Relative fitness values were obtained by competing each virus population with a marked reference virus that contains four adjacent silent mutations in the polymerase region introduced by direct mutagenesis. Co-infections were performed in triplicate at an MOI of 0.01 using a 1:1 mixture of each variant with the reference virus for 24 hours. The proportion of each virus was determined by real time RT-PCR on extracted RNA from the infection supernanant, using a mixture of Taqman probes labelled with two different fluorescent reporter dyes. MGB_CVB3_WT detects WT virus (including the fidelity variants) with the sequence CGCATCGTACCCATGG and labelled at the 5’ end with a 6FAM dye (6-carboxyfluorescein) and MGB_CVB3_Ref containing the four silent mutations: CGCTAGCTACCCATGG was labelled with a 5’ VIC dye. Each 25 μl-reaction contained 5ul of RNA, 900nM of each primer (forward primer, 5’-GATCGCATATGGTGATGATGTGA-3’; reverse primer, 5’-AGCTTCAGCGAGTAAAGATGCA-3’) and 150 nM of each probe. The relative fitness was determined by the method described by Carrasco et al. [15]. Briefly, the formula W = [R(t)/(R (0))] ^(1/t), represents the fitness, W, of each mutant genotype relative to the common competitor reference sequence, where R(0) and R(t) represent the ratio of mutant to reference virus densities in the inoculation mixture and t days post-inoculation, respectively. The fitness of the normal wildtype to reference virus was 1.019, indicating no significant differences in fitness due to the silent mutations engineered in the reference virus.

BLItz assay

CVB3 parental strain, and CVB3-E76G were propagated in HeLa cells and purified as described previously [14]. The membranes of infected cells were broken by three freeze-thaw cycles. Virus was concentrated by pelleting through sucrose and purified by tartrate step gradient ultracentrifugation. The virus bands were collected and exchanged into PBS and virus concentration and quality was estimated by measuring the absorbance at 260, 280, and 310 nm. Binding assays were performed by biolayer interferometry (BLI) measured by the BLItz from Fortebio [21]. Briefly, BLI measures binding by sending white light down a glass fiber-based biosensor, which is reflected back up to the instrument from two interfaces: 1) the interface between the glass fiber and the biosensor, and 2) the interface between the surface chemistry and solution. Since the two reflections come from the same white light source in the instrument, they both contain the same wavelengths. When molecules bind to the surface of the biosensor, the path length of the reflection (the one reflecting from the interface between surface and solution) increases while that of the other reflection remains the same.0.4% BSA and 0.08% Tween was added to virus and DAF to prevent non-specific binding to the sensor during the assay. Purified, His-tagged DAF was diluted to 0.1 mg/ml and attached to a Ni-NTA sensor. The DAF loaded sensor was then dipped into virus diluted to 0.5 mg/ml.

Mathematical modelling of emerging minority variants in deep sequence data

Let X be a phylogenetic tree, taken to include relative population sizes of the branches, and let X _t = (s _t p _t)^T be the state of the tree at time t = 1, 2,…,T. Here s _t denotes the structure of the tree, i.e. the set of branches existing at time t and p _t is a vector with the relative population sizes of all existing branches. Furthermore, let Y be the set of measurements, where Y _{t +} ^64xN is the number of reads supporting each of the 64 codons at each of the N different positions in the genome. By Bayes' theorem,

The output probability P(Y | X) follows a multinomial distribution for each position. The dynamics of the phylogenetic tree are naturally assumed to be Markovian, i.e.

where the transition probability can be expressed using the structural and population dynamics parts separately,

A simple random walk model is used for the population dynamics P(p _t | s _t, p _t-1), the population change from t−1 to t is taken to be lognormally distributed with standard deviation σh where h is the length of the time step. The structural part, P(s _t | s _t−1, p _t−1), depends on the mutation frequency of the virus and the population size of the haplotype in which the mutation occurs. The posterior probability is thus, assuming that s ₁ and p ₁ are known,

We would like to thank Dr. Isabel Puigdomenech for her help with protein expression analysis and for helpful discussions. We would like to dedicate this paper in her memory, a promising young scientist and friend, who left us abruptly and far too soon.