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