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