Merge branch 'master' of inn.eb.local:/ebio/ag-neher/home/fzanini/phd/papers/synmut
[synmut.git] / synmut.tex
1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 % ARTICLE ABOUT FATE OF SYNONYMOUS MUTATIONS IN HIV
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 \documentclass[rmp, twocolumn]{revtex4}
5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 \usepackage[english]{babel}
8 \usepackage[utf8x]{inputenc}
9 \usepackage{amsmath,amsfonts,amssymb,eucal,eurosym,textcomp}
10 \usepackage{color}
11 \usepackage{graphicx}
12 \usepackage[caption=false]{subfig}
13 \usepackage{natbib}
14 \usepackage{pslatex}
15 \usepackage[colorlinks,linkcolor=red,citecolor=red]{hyperref}
16 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
17 \graphicspath{{./figures/}}
18 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
19 \newcommand{\comment}[1]{\textit{\textcolor{red}{#1}}}
20 \newcommand{\mut}{\mu}
21 \newcommand{\mfit}{\langle F\rangle}
22 \newcommand{\mexpfit}{\langle e^{F}\rangle}
23 \newcommand{\ox}{r}
24 \newcommand{\co}{\rho}
25 \newcommand{\gt}{g}
26 \newcommand{\locus}{s}
27 \newcommand{\locuspm}{t}
28 \newcommand{\OO}{\mathcal{O}}
29 \newcommand{\rev}{\textit{rev}}
30 \newcommand{\FIG}[1]{Fig.~\ref{fig:#1}}
31 \newcommand{\FIGS}[2]{Figs.~\ref{fig:#1} and~\ref{fig:#2}}
32 \newcommand{\env}{\textit{env}}
33 \newcommand{\shankaregion}{C2-C5}
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35 \renewcommand{\thesubfigure}{\Alph{subfigure}}
36 \newcommand{\Author}{Fabio~Zanini and Richard~A.~Neher}
37 \newcommand{\Title}{Deleterious synonymous mutations hitchhike to high frequency in HIV \env~evolution}
38 \newcommand{\Keywords}{{HIV}, {synonymous}, {population genetics}}
39 \hypersetup{pdfauthor={\Author}, pdftitle={\Title}, pdfkeywords={\Keywords}}
40 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
41 \begin{document}
42 \title{\Title}
43 \author{\Author}
44 \date{\today}
45 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
46
47 \begin{abstract}
48 \noindent
49 Intrapatient HIV evolution is governed by selection on the protein level in the
50 arms race with the immune system (killer T-cells and antibodies). Synonymous
51 mutations do not have an immunity-related phenotype and are often assumed to be
52 neutral. In this paper, we show that synonymous changes in epitope-rich regions
53 are often deleterious but still reach high frequencies in the viral population. We analyze time
54 series of viral sequences from the \shankaregion~part of {\it env} within individual
55 hosts and observe that synonymous derived alleles rarely reach fixation.
56 Simulations suggest that such synonymous mutations
57 have a (Malthusian) selection coefficient of the order of $-0.001$, and that
58 they are brought up to high frequency by hitchhiking on neighboring beneficial
59 nonsynonymous alleles. We detect a negative correlation between the fixation of an allele and
60 its involvement in RNA stem-loop structures.
61 Deleterious synonymous mutations are not observed as abundantly in other parts of the HIV genome, in which
62 selective sweeps are less dense and hitchhiking not as strong; this behaviour is
63 confirmed by extensive computer simulations.
64
65 \end{abstract}
66 \maketitle
67
68 \section{Introduction}
69
70 HIV evolves rapidly within a single host during the course of the infection.
71 This evolution is driven by strong selection imposed by the host immune system
72 via killer T cells (CTLs) and neutralizing antibodies
73 (ABs)~\citep{rambaut_causes_2004} and facilitated by the high
74 mutation rate of HIV~\citep{mansky_lower_1995}. When the host develops a CTL or
75 AB response against a particular HIV epitope, mutations in the viral genome that
76 reduce or prevent recognition of the epitope frequently emerge. Escape mutations
77 in epitopes targeted by CTLs typically evolve during early infection and spread
78 rapidly through the population~\citep{mcmichael_immune_2009}. During chronic
79 infection, the most rapidly evolving part of the HIV genome are the so called
80 variable loops of the envelope protein gp120, which need to avoid recognition by
81 neutralizing ABs. Mutations in \env, the gene encoding for gp120, spread
82 through the population within a few months (see \figurename~\ref{fig:aft}, solid
83 lines). Consistent with this time scale, it is found that serum from a
84 particular time typically neutralizes virus extracted more than 3-6 month
85 earlier \citep{richman_rapid_2003}.
86
87 These escape mutations are strongly selected for their effect on the amino acid
88 sequence of the viral proteins. Conversely, synonymous mutations are commonly
89 used as approximately neutral markers in studies of viral evolution. Neutral
90 markers are very useful in practice, because there evolution can be compared to
91 that of putatively functional sites to detect purifying or directional selection
92 \citep{Bhatt:2011p43255,Hurst:2002p32608,Chen:2004p22606}. In most species,
93 however, it is known that synonymous codons don't evolve completely neutrally
94 but that some codons are favored over others \citep{plotkin_synonymous_2011}.
95 In addition to maintaining protein function and avoiding immune
96 recognition, the HIV genome has to ensure efficient processing and translation,
97 nuclear export, and packaging into the viral capsid: all these processes operate at the RNA
98 level and are sensitive to synonymous changes. A few functionally important RNA
99 elements are well characterized. For example, a certain RNA sequence in the HIV
100 genome, called \rev{} response element (RRE), enhances nuclear export of viral
101 transcripts~\citep{fernandes_hiv-1_2012}. Another well studied case is the
102 interaction between viral reverse transcriptase, viral ssRNA, and the host
103 tRNA$^\text{Lys3}$: the latter is required for priming reverse transcription
104 (RT) and bound by a specifical pseudoknotted RNA structure in the viral 5'
105 untranslated region~\citep{barat_interaction_1991, paillart_vitro_2002}.
106 Nucleotide-level fitness effects have been observed beyond RNA structure as
107 well. Recent studies have shown that genetically engineered HIV strains with
108 a skewed codon usage bias towards more or less abundant tRNAs
109 replicate better or worse \comment{replicate better or produce more protein?},
110 respectively~\citep{ngumbela_quantitative_2008, li_codon-usage-based_2012}. A similar conclusion has been reached about
111 influenza, and codon-pessimized influenza strains have been shown to be good
112 live attenuated vaccines in mice~\citep{mueller_live_2010}. Purifying selection
113 beyond the protein sequence is therefore expected
114 \citep{forsdyke_reciprocal_1995,snoeck_mapping_2011}, but it seems reasonable
115 that the bulk of positive selection through the immune system be restricted to amino acid sequences.
116
117 %SYNONYMOUS CONSERVATION. DO WE HAVE A PLOT OF GENOME WIDE CONSERVATION, MAYBE
118 %FOR SUPPLEMENT? YES
119
120 In this paper, we characterize the dynamics of synonymous mutations in \env{}
121 and show that a substantial fraction of these mutations is deleterious. We
122 further show that, although such synonymous mutations cannot be used as neutral
123 markers, the degree to which they hitchhike with nearby nonsynonymous mutations
124 is very informative both about their on deleterious effect, as well as the
125 strength and nature of selection for nearby positively selected sites.
126 Their ability to hitchhike for extended times, which is a core requirement for
127 our analysis, is rooted in the small recombination rate of
128 HIV~\citep{neher_recombination_2010, batorsky_estimate_2011}. Extending the
129 analysis of fixation probabilities to the nonsynonymous mutations, we show that
130 time dependent selection or strong competition of escape mutations inside the
131 same epitope are necessary to explain the observed patterns of fixation and
132 loss.
133
134 \section{Results}
135
136 The central quantity we investigate is the probability of fixation of a
137 mutation, conditional on its population frequency. A neutral mutation
138 segregating at frequency $\nu$ has a probability $P_\text{fix}(\nu) = \nu$ to
139 spread through the population and fix; in the rest of the cases, i.e.~with
140 probability $1-\nu$, it goes extinct. As illustrated in the inset of \FIG{aftsyn},
141 this is a simple consequence of the fact that
142 (a) exactly one of the $N$ individuals in the current population will be
143 the common ancestor of the entire future population at a particular locus and
144 (b) this ancestor has a probability $\nu$ of carrying the mutation (assuming
145 the neutral mutation is not preferentially associated with genomes of high or
146 low fitness).
147 Deleterious or beneficial mutations fix less or
148 more often than neutral ones, respectively. Time series sequence data enable a
149 direct observation of both the current frequency $\nu$ of any particular
150 mutation and its future fate (fixation or extinction).
151
152
153 \subsection{Synonymous polymorphisms in \env, C2-V5, are mostly deleterious}
154
155 \FIG{aft} shows time series data of the frequencies of all mutations observed
156 \env~, C2-V5, in patient p8~\citep{shankarappa_consistent_1999}. Despite many
157 synonymous mutations reaching high frequency, very few fix
158 (panel~\ref{fig:aftsyn}); however, many nonsynonymous mutations fix
159 (panel~\ref{fig:aftnonsyn}).
160 These observations are quantified in \FIG{fixp}, which stratifies the data
161 of 7 (resp.~10) patients according to the frequency at which
162 different mutations are observed (see methods). Considering all mutations in a
163 frequency interval around $\nu_0$ at some time $t_i$, we calculate the fraction
164 that is found at frequency 1, at frequency 0, or at intermediate frequency at
165 later times $t_f$. Plotting these fixed, lost, and polymorphic fraction against
166 the time interval $t_f-t_i$, we see that most synonymous mutations segregate for
167 roughly one year and are lost much more frequently than expected (panel
168 \ref{fig:fixp1}). The long-time probability of fixation versus extinction of
169 synonymous mutations is shown as a function of the initial frequency $\nu_0$
170 in panel~\ref{fig:fixp2} (red line). The nonsynonymous seem to follow more or
171 less the neutral expectation (blue line) -- a point to which we will come back
172 below.
173
174 \begin{figure}
175 \begin{center}
176 \subfloat{\includegraphics[width=\linewidth]{Shankarappa_allele_freqs_trajectories_syn_p8.pdf}
177 \label{fig:aftsyn}}\\
178 \subfloat{\includegraphics[width=\linewidth]{Shankarappa_allele_freqs_trajectories_nonsyn_p8.pdf}
179 \label{fig:aftnonsyn}}
180 \caption{Time series of allele frequencies in \env, C2-V5, from
181 patient 8~\cite{shankarappa_consistent_1999},
182 divided by synonymous (panel A) and nonsynonymous (panel B) derived alleles.
183 While nonsynonymous mutations frequently fix, very few synonymous
184 mutations do even though they are frequently observed at intermediate
185 frequencies. Colors indicate the position of the site along the C2-V5 region
186 (blue to red). Inset: the fixation probability $P_\text{fix}(\nu)$ of a neutral mutation
187 is simply the likelihood that the future common ancestor is currently carrying
188 it, i.e. its frequency $\nu$.}
189 \label{fig:aft}
190 \end{center}
191 \end{figure}
192
193 \citet{bunnik_autologous_2008} present a longitudinal dataset on the entire
194 \env~gene of 3 patients at $\sim 5$ time points with approximately 5-20
195 sequences each (see methods). Repeating the above analysis separately on the
196 C2-V5 region studied above and the remainder of \env~ reveal striking
197 differences (see \FIG{fixp2}). Within C2-V5, this data fully confirms the
198 observations made in the data set by \citet{shankarappa_consistent_1999} (red
199 line). In the remainder of \env, however, observed synonymous mutations behave
200 neutrally (green line).
201
202 These observations suggest that many of the synonymous polymorphisms at
203 intermediate frequencies in the part of \env~that includes the hypervariable
204 regions are deleterious, while outside this regions polymorphisms are mostly
205 roughly neutral. Note that this does not imply that all synonymous mutations are
206 neutral -- only those mutations observed at high frequencies tend to be neutral.
207
208 \begin{figure}
209 \begin{center}
210 \subfloat{\includegraphics[width=0.9\linewidth]{Shankarappa_fix_loss_dt_times}
211 \label{fig:fixp1}}\\
212 \subfloat{\includegraphics[width=0.9\linewidth]{Bunnik2008_fixmid_syn_ShankanonShanka}
213 \label{fig:fixp2}}
214 \caption{(Panel A) Time course of loss and fixation of synonymous mutations
215 observed in a frequency interval $\nu_0$. The ultimate fraction of synonymous
216 mutations that fix as a function of intermediate frequency $\nu_0$ is the
217 fixation probability. (Panel B) Fixation probability of derived synonymous
218 alleles is strongly suppressed in C2-V5 versus other parts of the {\it env}
219 gene. The horizontal error bars on the abscissa are bin sizes, the vertical ones
220 the standard deviation after 100 patient bootstraps of the data. Data from
221 Refs.~\cite{shankarappa_consistent_1999, bunnik_autologous_2008}.}
222 \label{fig:fixp}
223 \end{center}
224 \end{figure}
225
226 \subsection{Synonymous mutations in C2-V5 tend to disrupt conserved RNA stems}
227
228 One possible {\it a priori} explanation for lack of fixation of synonymous
229 mutations in C2-V5 are secondary structures in the viral RNA. It has
230 been suggested early on that parts of the viral genome that has the potential to
231 form stems is better conserved than the
232 remainder~\citep{forsdyke_reciprocal_1995,snoeck_mapping_2011}.
233
234 Recently, the propensity of nucleotides of the HIV genome to form base pairs has
235 been measured using the SHAPE assay (a biochemical reaction preferentially
236 altering unpaired bases)~\citep{watts_architecture_2009}. The SHAPE assay has
237 shown that the variable regions V1 to V5 tend to be unpaired, while the
238 conserved regions between those variable regions form stems. We partition all
239 synonymous alleles observed at intermediate frequencies above 10-15\% depending
240 on their final destiny (fixation or extinction). Subsequently, we align our
241 sequences to the reference NL4-3 strain used in
242 ref.~\citep{watts_architecture_2009} and assign them SHAPE reactivities. As
243 shown in \FIG{SHAPEA} in a cumulative histogram, the reactivities of fixed
244 alleles (red line) are systematically larger than of alleles that are doomed to
245 extinction (blue line) (Kolmogorov-Smirnov test, $P\approx 0.002$).
246 In other words, alleles that are likely to be
247 breaking RNA helices are also more likely to revert and finally be lost from the
248 population. As a control, the average over non-observed but potentially
249 available polymorphisms lies between the two curves (green line), as expected
250 (because only some of these mutations will interfere with stem loop formation).
251 To test the hypothesis that mutations in C2-V5 are lost since break stems
252 in the conserved between the variable loops, we split the synonymous mutations
253 in the extended V1-V5 region further into conserved and variable regions and
254 found that the biggest depression in fixation probability is observed in the
255 conserved stems, while the variable loops show little deviation from the
256 neutral signature, see \FIG{SHAPEB}. This is consistent with important stem
257 structures in conservered regions between loops.
258
259 In addition to RNA secondary structure, we have considered other possible
260 explanations for a fitness defect of synonymous mutations, in particular codon
261 usage bias (CUB). HIV is known to prefer A-rich codons over highly expressed
262 human codons~\citep{jenkins_extent_2003}. Moreover, codon-optimized
263 and -pessimized viruses have recently been generated and shown to replicate
264 better or worse than wild type strains,
265 respectively~\citep{li_codon-usage-based_2012, ngumbela_quantitative_2008,
266 coleman_virus_2008}. We do not find, however, any evidence for a contribution of
267 CUB to the ultimate fate of synonymous alleles. Several lines of thought support
268 this result. First of all, within a single patient, we do not observe any bias
269 towards more human-like CUB in the synonymous mutations that reach fixation
270 rather than extinction. Second, although codon-optimized HIV seems to perform
271 better {\it in vitro}, the distance in CUB between HIV and human genes is not
272 shrinking at the macroevolutionary level. Third, it is a common phenomenon for
273 retroviruses to use variously different codons from their hosts, and CUB effects
274 on fitness are thought to be so small that divergent nucleotide composition has
275 been suggested as a possible mechanism for viral
276 speciation~\citep{bronson_nucleotide_1994}. Fourth, CUB in the V1-C5 region is
277 not very different from other parts of the HIV genome, whereas the reduced
278 fixation probability is only observed there. In conclusion, although we cannot
279 exclude an effect of CUB on fitness as a general rule, we expect it to be a
280 minor effect in our context. \comment{Need to read literature}
281
282 \begin{figure}
283 \begin{center}
284 \subfloat{\includegraphics[width=0.9\linewidth]{mixed_Shankarappa_Bunnik2008_Liu_fixation_reactivity_Vandflanking_fromSHAPE}
285 \label{fig:SHAPEA}}\\
286 \subfloat{\includegraphics[width=0.9\linewidth]{Shankarappa_fixmid_syn_V_regions.pdf}\label{fig:SHAPEB}}
287 \caption{Permissible synonymous mutations have
288 high SHAPE reactivities. \comment{Is this C2-V5 only?}
289 Panel A) Synonymous substitions in the C2-V5 region have substantially higher
290 SHAPE reactivities than mutations that reach frequencies above 15\%, but are
291 susequently lost. This difference is quantified by a cumulative histogram
292 ($p=0.002$, Kolmogorov-Smirnov two sample test). Mutations that are never
293 observered show an intermediate distribution of SHAPE values.
294 Panel B shows the fixation probablility of synonymous mutations in C2-V5
295 separately for variable loops and the connecting conserved regions, that
296 harbor conserved RNA stems. As expected, fixation probability is lower inside
297 the conserved regions. Data from Refs.~\cite{shankarappa_consistent_1999,
298 bunnik_autologous_2008, liu_selection_2006}.}
299 \label{fig:SHAPE}
300 \end{center}
301 \end{figure}
302
303
304 \subsection{Deleterious mutations are brought to high frequency by hitch-hiking}
305
306 While the observation that some fraction of synonymous mutations is deleterious
307 is not unexpected, it seems odd that we observe them at high population
308 frequency -- at least in some regions of the genome. The region of \env~ in
309 which we observe deleterious mutations at high frequency, however, is special in
310 that it undergoes frequent adaptive changes to evade recognition by neutralizing
311 antibodies~\cite{williamson_adaptation_2003,richman_rapid_2003}. Due to the limited amount of
312 recombination in HIV~\cite{neher_recombination_2010,batorsky_estimate_2011},
313 deleterious mutations that are linked to adaptive variants can reach high
314 frequency~\citep[genetic hitchhiking,][]{smith_hitch-hiking_1974}.
315
316 The potential for hitchhiking is already apparent from the allele frequency
317 trajectories in \FIG{aft}, where many mutations appear to change rapidly in
318 frequency as a flock. Deleterious synonymous mutations can be brought to high
319 frequency by selection on linked nonsynonymous sites, a process
320 known as {\it genetic draft}~\citep{gillespie_genetic_2000, neher_genetic_2011}. In order
321 to be advected to high frequency by a linked adaptive mutation, the deleterious
322 effect of the mutation has to be substantially smaller than the adaptive effect.
323 The latter was estimated to be on the order of $s_a = 0.01$ per day~\citep{neher_recombination_2010}.
324 The approximate magnitude of the deleterious effects can be estimated from
325 \FIG{fixp1}, that shows the distribution of times for synonymous
326 alleles to reach the fix or get lost starting from intermediate frequencies. The
327 typical time to loss is of the order of 500 days. If this loss is driven by the
328 deleterious effect of the mutation, this corresponds to deleterious effects of
329 roughly $s_d \sim - 0.002$ per day.
330
331 \begin{figure}
332 \begin{center}
333 \subfloat{\includegraphics[width=\linewidth]{simulations_graduallydel.pdf}
334 \label{fig:simfixpvar}}\\
335 \subfloat{\includegraphics[width=\linewidth]{simulations_syn.pdf}
336 \label{fig:simsfig}}
337 \caption{Evolutionary parameters compatible with observations.
338 Panel A) The depression in $P_\text{fix}$ depends on the deleterious
339 effect size of the synonymous alleles (in HIV data is about $-0.2$).
340 Increasing this parameter also reduces synonymous
341 diversity, measured by the number of mutations at frequencies $0.25 < \nu < 0.75$,
342 normalized by the number of possible mutations (first inset).
343 Panel B) To assess the parameter space compatible with observations, we run many
344 simulations with random parameters for deleterious effect size, fraction of
345 deleterious synonymous sites, average size of escape mutations (color, blue ($10^{-2.5}$) to red ($10^{-1.5}$)), and rate of
346 introduction of new epitopes (marker size, from $10^{-3}$ to $10^{-2}$ per
347 generation). Mostly simulations with small deleterious effects and relatively
348 high deleterious fractions of synonymous mutations reproduce synonymous HIV-like
349 fixation probability and diversity. Parameters are chosen uniformly on a
350 logarithmic scale from the indicated windows. \comment{should we extend the
351 parameters space maybe to $10^{-4}$?}}
352 \label{fig:simheat}
353 \end{center}
354 \end{figure}
355
356 To get a better idea of the range of parameters that are compatible with the
357 observations and our interpretation, we perform computer simulations of
358 evolving viral populations assuming a mix of positive and purifying selection
359 and rare recombination.
360 For this purpose, we use the simulation package FFPopSim, which includes a
361 module dedicated to intrapatient HIV evolution~\citep{zanini_ffpopsim:_2012}.
362
363 The first result of the simulations is that genetic draft can indeed bring weakly
364 deleterious mutations to high frequencies and results in a dependence of the
365 fixation probability on initial frequency that is compatible with observations
366 (see \FIG{simfixpvar}).
367 The fixation probability of synonymous alleles decreases from the neutral
368 expectation to zero as their deleterious effect increases
369 (\FIG{simfixpvar}, from $-0.0003$ in blue to $-0.003$ in red). Along with the
370 reduction in fixation probability, the number of alleles that reach high
371 frequencies, which is a proxy for synonymous diversity, is reduced as well (\FIG{simfixpvar}, top inset).
372
373 To assess the reduction in fixation probability quantitatively, we calculate
374 the area between the measured fixation probability and the diagonal, which is
375 the neutral expectation (\FIG{simfixpvar}, lower inset). This area spans a range
376 between $-0.5$ (no fixation at all), zero (neutral-like curve), and $+0.5$ (always
377 fixed). The HIV data correspond to an area of approximately $-0.2$ for synonymous
378 changes, and essentially zero for nonsynonymous changes.
379
380 We pick many random combinations of the following four parameters: effect size
381 and rate of appearance of beneficial escape mutations, and effect size and
382 fraction of deleterious synonymous mutations.
383 Among all simulations, we select the one that show a large depression in
384 fixation probability (>-0.??) of synonymous mutations but, simulteneously, a
385 high synonymous diversity (>0.0? per site), as observed in HIV data. This
386 selection imposes a strong constraint on the parameters, in particular on the
387 ones related to synonymous changes, i.e.~their deleterious fraction and effect
388 size (see \FIG{simsfig}).
389 A high fraction ($\gtrsim 0.8$) of only slightly deleterious synonymous sites,
390 with effect size $|s_d| \lesssim 0.002$, is required. This result fits well the
391 expectation based on the fixation/extinction times above (see \FIG{fixp1}).
392 Moreover, it corresponds to {\it a priori} theoretical expectations:
393 (a) since high frequency bins are enriched in neutral mutations (it is easier
394 for them to hitchhike), the hallmark of deleterious mutations will be visible
395 there only if they are pervasive; (b) genetic draft stops working if the
396 deleterious effect size becomes comparable to the adaptive effect size of escape
397 mutations, because the double mutant has little or no fitness advantage.
398
399 As long as only features of synonymous mutations are filtered, as performed
400 above, no strong correlation of the two nonsynonymous parameters, effect size
401 and rate of appearance of epitopes, is observed; in \FIG{simsfig}, the points
402 shown possess a variety of shapes and colors. In most of these simulations,
403 however, nonsynonymous mutations almost always fix once they reach high frequencies
404 -- their nonsynonymous fixation area is well above zero. This is incompatible
405 with the blue line in \FIG{fixp}: in an HIV infection, nonsynonymous mutations fix as if they
406 were neutral. Although there are no doubts that nonsynonymous variation in the
407 variable regions is driven by positive selection and not by genetic drift, the
408 biological reason for the neutral-like curve in \FIG{fixp} is unclear.
409 Inspecting the trajectories of nonsynonymous mutations
410 suggests the rapid rise and fall of many alleles. We test two possible
411 mechanisms that are biologically plausible and could explain the transient rise
412 of nonsynonymous mutations: time-dependent selection and within-epitope
413 competition.
414
415 The former explanation can be formulated as follows: if the immune system
416 recognizes the escape mutant before its fixation, the mutant might cease to be
417 beneficial and disappear soon, despite its quick initial rise in frequency. In support of this idea,
418 \citet{richman_rapid_2003, bunnik_autologous_2008} report antibody responses to
419 escape mutants. These respones are delayed by a few months, roughly matching the
420 average time needed by an escape mutant to cross frequencies of order one.
421 In the alternative explanation, several different escape
422 mutations within the same epitope might arise almost simultaneously and start to
423 spread. Their benefits are not additive, because each of them is
424 essentially sufficient to escape. As a consequence, several escape mutations rise to
425 high frequency rapidly, while the one with the smallest cost in terms of replication,
426 packaging, etc.~is most likely to eventually fix. The emergence of
427 multiple sweeping nonsynonymous mutations in real HIV infections has been shown
428 previously~\citep{moore_limited_2009, bar_early_2012}.
429 We tested both hypotheses in our simulation framework and they both seem to
430 achieve a reduced, neutral-like fixation probability while maintaining a high
431 rate of substitutions, compatibly with HIV measurements. In real HIV infections,
432 both mechanisms are likely to be playing a role.
433
434 \section{Discussion}
435 Despite several known functional roles for RNA secondary structure in the HIV
436 genome, synonymous mutations are often used as approximately neutral markers in
437 evolutionary studies of viruses. We have shown that the majority of synonymous
438 mutations in the conserved regions C2-C5 of the \env~gene are deleterious.
439 Comparison with recent biochemical studies of binding propensity of bases in RNA
440 genome suggest that these mutations are deleterious, at least in part, because they disrupt
441 stems in RNA secondary structures. Furthermore, we provide evidence that these
442 mutations are brought to high frequency through linkage to adaptive mutations.
443 The latter mutations are only transiently adaptive, either through a
444 coevolution with the immune system or redundant escape within an epitope.
445
446 Our observations and conclusion rely heavily on longitudinal data in which the
447 dynamics of mutations can be explicitly observed. The fact that deleterious
448 mutations can be brought to high frequencies through hitchhiking underscores
449 the intensity of the coevolution with the immune system. The fact that
450 multiple escape mutations in the same epitope -- as is indeed observed in
451 studies of antibody escape~\citep{moore_limited_2009, bar_early_2012} -- are
452 necessary to explain the patterns of fixation of nonsynonymous mutations points
453 towards a large populations size that rapidly discovers adaptive mutations. A
454 similar point has been made recently by Boltz {\it et al.} in the context of
455 preexisting drug resistance mutations~\citep{boltz_ultrasensitive_2012}.
456
457 The observed hitchhiking highlights the importance of linkage due to infrequent
458 recombination for the evolution of HIV
459 \citep{neher_recombination_2010,batorsky_estimate_2011,
460 josefsson_majority_2011}. The recombination rate has been estimated to be on the
461 order of $\rho = 10^{-5}$ per base and day. It takes roughly $t_{sw} = s_a^{-1}
462 \log \nu_0$ generations for an adaptive mutation with growth rate $s_a$ to rise
463 from an initially low frequency $\nu_0\sim \mu$ to frequency one. This implies
464 that a region of length $l = (\rho t_{sw})^{-1} = s_a / \rho \log \nu_0$ remains
465 linked to the adaptive mutation. With $s_a=0.01$, $l\approx 100$ bases, hence we
466 expect strong linkage between the variable loops and their surrounding stems,
467 but none far beyond the variable regions, consistent with the lack of signal
468 outside of C1-V5. In case of much stronger selection -- such as observed during
469 early CTL escape or drug resistance evolution -- the linked region is of course
470 much larger.
471
472 The functional significance of the insulating RNA structure stems between the
473 hyper variable loops has been proposed
474 previously~\citep{watts_architecture_2009, sanjuan_interplay_2011}.
475 \citet{sanjuan_interplay_2011} have shown that insulating stems are relevant for
476 viral fitness {\it in vivo}. Our analysis is limited by the availability of
477 longitudinal data which requires a focus on the the variable regions of \env.
478 Conserved RNA structures exist in different parts of the HIV genome (several are
479 known). In absence of repeated adaptive substitutions in the vicinity that cause
480 hitchhiking, the deleterious synonymous mutations remain at low frequencies and
481 can only be observed by deep sequencing methods.
482
483 Our study uncovers the
484 subtle balance of evolutionary forces governing intrapatient HIV evolution. The
485 fixation and extinction times and probabilities represent a rich and simple
486 summary statistics useful to characterize sequencing data and compare to
487 models via computer simulations.
488 A similar method has been recently used in a longitudinal study of
489 influenza~\citep{strelkowa_clonal_2012}. The propagators suggested in that
490 paper, however, represent ratios between nonsynonymous
491 mutations and synonymous ones, hence they are inadequate to investigate
492 synonymous changes themselves. The authors also conclude that several
493 beneficial mutations segregate simultaneously in influenza, a scenario
494 remarkably similar to our within-epitope competition picture. While
495 recombination hardly affects evolution within an epitope, it might facilitate
496 immune escape by easing competition between different epitopes
497 \citep{neher_rate_2010,Rouzine:2005p17398}.
498
499 Our results emphasize the inadequacy of independent site
500 models of HIV evolution and the common assumption that selection is time
501 independent. If genetic variation is only transiently beneficial, existing
502 methods to quantify selection will yield substantial underestimates
503 \citep{williamson_adaptation_2003,neher_rate_2010,OTHER}. To explain the
504 observations regarding the fixation probabilities of non-synonymous mutations,
505 either transient selection, or substational within-epitope competition are
506 necessary. Which mechanism is more widespread is not clear as of now,
507 there is evidence for both~\citep{richman_rapid_2003, moore_limited_2009,
508 bar_early_2012}.
509
510 \section{Methods}
511 \subsection{Sequence data collection}
512 Longitudinal intrapatient viral RNA sequences were collected for published
513 studies~\citep{shankarappa_consistent_1999,
514 liu_selection_2006, bunnik_autologous_2008} and downloaded from the Los Alamos
515 National Laboratory (LANL) HIV sequence database~\citep{LANL2012}. The sequences from
516 some patients showed signs of HIV compartimentalization into subpopulations and
517 were discarded; a grand total of 11
518 patients with approximately 6 time points each and 10 sequences per time point
519 were analyzed. The time interval or resolution between two consecutive sequences
520 was approximately 6 to 18 months.
521
522 \comment{How did you determine who to discard. Maybe we should include the PCA
523 plots into the supplement.}
524
525
526 \subsection{Sequence analysis}
527 The good sequences were aligned within each patient
528 via the translated amino acid sequence, using
529 Muscle~\citep{edgar_muscle:_2004}, and to the NL4-3 reference sequence used
530 by \citet{watts_architecture_2009} in the SHAPE assay. Within each patient, a
531 consensus RNA sequence at the first time point was used to classify alleles as ancestral or
532 derived at all sites. Problematic sites that included large frequencies of gaps
533 were excluded from the analysis to avoid artefactual substitutions due to
534 alignment errors. Time series of allele frequencies were extracted from the
535 sequences.
536
537 The synonymity of a mutation was assigned if the rest of the codon was
538 in the ancestral state and using the standard genetic code. Cases where more
539 than one mutation within the codon was observed were discarded. Slightly
540 different criteria for synonymous/nonsynonymous discrimination yielded similar
541 results.
542
543 \subsection{Fixation probability and secondary structure}
544 For the estimate of times to fixation/extinction, polymorphisms were
545 binned by frequency and the time to reaching the first boundary (fixation or
546 extinction) was stored. For the fixation probability, the long-time limit of the
547 resulting curves was used, excluding polymorphisms that arose late in the
548 clinical history (and would have had no time to reach either boundary).
549
550 For the correlation analysis with RNA secondary structure, the SHAPE scores were
551 downloaded from the journal website~\citep{watts_architecture_2009}. By virtue
552 of the alignment of the longitudinal sequences with the reference used by
553 \citet{watts_architecture_2009}, SHAPE reactivities were assigned to most sites.
554 Problematic assignments in indel-rich regions were excluded from the analysis.
555 In order to restrict the analysis to synonymous polymorphisms, a lower frequency
556 threshold of 0.15 was used (other thresholds yielded the same results). Since
557 very few polymorphisms hitchhike beyond, say, a frequency of 0.5, this pool is
558 enriched for to-be-lost mutations; hence the ``lost" curve in \FIG{SHAPEA}
559 contains much more points than the ``fixed" one.
560
561 The V loops and flanking regions were identified manually starting from the
562 annotated reference HXB2 sequence from the LANL HIV database~\citep{LANL2012}. A
563 similar approach was used to label the C2-V5 region sequenced in
564 ref.~\citep{shankarappa_consistent_1999}.
565
566 \subsection{Computer simulations}
567 Simulations were performed using the recently published software
568 FFPopSim~\citep{zanini_ffpopsim:_2012}. Both full-length HIV genomes and
569 \env{}-only simulations were performed and yielded comparable results. For each
570 set of parameters, approximately 100 simulation runs were averaged over. In each
571 run, a random fitness landscape with specified statistical properties (e.g.
572 density of beneficial sites, average deleterious effect of synonymous changes) was generated.
573 Although the curves shown in \FIG{simfixpvar} are not very smooth, small
574 parameter changes resulted in overall consistent trends across many repetitions.
575
576 For the discussion of simulation parameters, the areas below or above the neutral
577 diagonal were estimated from the binned fixation probabilities using the linear
578 interpolation between the bin centers. This measure is sufficiently precise for
579 our purposes, because the HIV data are quite scarse themselves.
580
581 \section*{Acknowledgements}
582 \comment{to be written\dots}
583
584
585 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
586 \bibliographystyle{natbib}
587 \bibliography{bib}
588 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
589 \end{document}
590 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
591