resolved conflict
[synmut.git] / synmut.tex
1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 % ARTICLE ABOUT FATE OF SYNONYMOUS MUTATIONS IN HIV
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 \documentclass[rmp, twocolumn]{revtex4}
5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 \newcommand{\Author}{Fabio~Zanini and Richard~A.~Neher}
8 \newcommand{\env}{\textit{env}}
9
10 \newcommand{\Title}{Deleterious synonymous mutations hitchhike to high frequency in HIV \env~evolution}
11 \newcommand{\Keywords}{{HIV}, {synonymous}, {population genetics}}
12 \usepackage[english]{babel}
13 \usepackage[utf8x]{inputenc}
14 \usepackage{amsmath,amsfonts,amssymb,eucal,eurosym,textcomp}
15 \usepackage{color}
16 \usepackage{graphicx}
17 \usepackage[caption=false]{subfig}
18 \usepackage{natbib}
19 \usepackage{pslatex}
20 \usepackage[colorlinks,linkcolor=red,citecolor=red]{hyperref}
21 \hypersetup{pdfauthor={\Author}, pdftitle={\Title}, pdfkeywords={\Keywords}}
22 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
23 \graphicspath{{./figures/}}
24 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
25 %\DeclareMathOperator\de{d\!}
26 \newcommand{\comment}[1]{\textit{\textcolor{red}{#1}}}
27 \newcommand{\mut}{\mu}
28 \newcommand{\mfit}{\langle F\rangle}
29 \newcommand{\mexpfit}{\langle e^{F}\rangle}
30 \newcommand{\ox}{r}
31 \newcommand{\co}{\rho}
32 \newcommand{\gt}{g}
33 \newcommand{\locus}{s}
34 \newcommand{\locuspm}{t}
35 \newcommand{\OO}{\mathcal{O}}
36 \newcommand{\rev}{\textit{rev}}
37 \newcommand{\FIG}[1]{Fig.~\ref{fig:#1}}
38 \newcommand{\FIGS}[2]{Figs.~\ref{fig:#1} and~\ref{fig:#2}}
39 \newcommand{\shankaregion{C2-C5}
40
41 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
43 \begin{document}
44 \title{\Title}
45 \author{\Author}
46 \date{\today}
47 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48
49 \begin{abstract}
50 \noindent
51 Intrapatient HIV evolution is governed by selection on the protein level in the
52 arms race with the immune system (killer T-cells and antibodies). Synonymous
53 mutations do not have an immunity-related phenotype and are often assumed to be
54 neutral. In this paper, we show that synonymous changes in epitope-rich regions
55 are often deleterious but still reach high frequencies. We analyze time
56 series of viral sequences from the \shankaregion~part of {\it env} within individual
57 hosts and observe that synonymous derived alleles rarely fix in the
58 viral population. Simulations suggest that such synonymous mutations
59 have a (Malthusian) selection coefficient of the order of $-0.001$, and that
60 they are brought up to high frequency by hitchhiking with neighboring beneficial
61 nonsynonymous alleles. We detect a negative correlation between fixation of an allele and
62 its involvement in RNA stem-loop structures.
63 Deleterious synonymous mutations are observed in other parts of the HIV genome, in which
64 selective sweeps are less dense and hitch-hiking not as strong.
65 \end{abstract}
66
67 \end{abstract}
68 \maketitle
69
70 \section{Introduction}
71
72 HIV evolves rapidly within a single host during the course of the infection.
73 This evolution is driven by strong selection imposed by the host immune system
74 via killer T cells (CTLs) and neutralizing antibodies
75 (ABs)~\citep{pantaleo_immunopathogenesis_1996} and facilitated by the high
76 mutation rate of HIV~\citep{mansky_lower_1995}. When the host develops a CTL or
77 AB response against a particular HIV epitope, mutations in the viral genome that
78 reduce or prevent recognition of the epitope frequently emerge. Escape mutations
79 in epitopes targeted by CTLs typically evolve during early infection and spread
80 rapidly through the population~\citep{mcmichael_immune_2009}. During chronic
81 infection, the most rapidly evolving part of the HIV genome are the so called
82 variable loops of the envelope protein gp120, which need to avoid recognition by
83 neutralizing ABs. Mutations in \env~, the gene encoding for gp120, spread
84 through the population within a few months (see \figurename~\ref{fig:aft}, solid
85 lines). The (Malthusian) effect size of these beneficial mutations is of the
86 order of $s_a \sim 0.01$~\citep{neher_recombination_2010}.
87
88 These escape mutations are strongly selected for their effect on the amino acid
89 sequence of the viral proteins. Conversely, synonymous mutations are commonly
90 used as approximate neutral markers in studies of viral evolution. Neutral
91 markers are very useful in practice, because they can be used to make inferences
92 about the stochastic forces driving evolution~\citep{yang_statistical_2000}.
93 The viral genome, however, needs to satisfy further constraints in addition to
94 immune escape, such as efficient processing and translation, nuclear export, and
95 packaging into the viral capsid: all these processes operate at the RNA level
96 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 skewed codon usage bias (CUB) patterns 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, while 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. Their ability to hitchhike for extended times, which is a
123 core requirement for our analysis, is rooted in the small recombination rate of
124 HIV~\citep{neher_recombination_2010, batorsky_estimate_2011}. Extending the
125 analysis of fixation probabilities to the nonsynonymous mutations, we show that
126 time dependent selection or strong competition of escape mutations inside the
127 same epitope are necessary to explain the observed patterns of fixation and
128 loss.
129
130 \section{Results}
131
132 The central quantity we investigate is the probability of fixation of a
133 mutation, conditional on its population frequency. A neutral mutation
134 segregating at frequency $\nu$ has a probability $P_\text{fix}(\nu) = \nu$ to
135 spread through the population and fix; in the rest of the cases, i.e. with
136 probability $1-\nu$, it goes extinct. This is a simple consequence of the fact
137 that (a) exactly one of the $N$ individuals in the current population will be
138 the common ancestor of the entire future population at a particular locus and
139 (b) this ancestor has a probability $\nu$ of carrying the mutation, see
140 illustration in \FIG{fixp}. Deleterious or beneficial mutations fix less or
141 more often than neutral ones, respectively. Time series sequence data enable a
142 direct observation of both the current frequency $\nu$ of any particular
143 mutation and its future fate (fixation or extinction). They therefore represent
144 a simple way to investigate average properties of different classes of
145 mutations.
146
147 \subsection{Synonymous polymorphisms in \env, C2-V5 are mostly deleterious}
148
149 \FIG{aft} shows time series data of the frequencies of all mutations observed
150 \env~, C2-V5, in patient p8~\citep{shankarappa_consistent_1999}. Despite many
151 synonymous mutations reaching high frequency (dashed lines), very few fix. This
152 observation is quantified in panels \FIG{fixp1} and \ref{fig:fixp2}, which
153 stratify the data of 7 (resp. 10) patients according to the frequency at which
154 different mutations are observed (see methods). Considering all mutations in a
155 frequency interval around $\nu_0$ at some time $t_i$, we calculate the fraction
156 that is found at frequency 1, at frequency 0, or at intermediate frequency at
157 later times $t_f$. Plotting these fixed, lost, and polymorphic fraction against
158 the time interval $t_f-t_i$, we see that most synonymous mutations segregate for
159 roughly two years and are lost much more frequently than expected. The long-time
160 probability of fixation versus extinction is shown as a function of the initial
161 frequency $\nu_0$ in panel~\ref{fig:fixp2} (red line). In contrast to synonymous
162 mutations, the nonsynonymous seem to follow more a less the neutral expectation
163 (blue line) -- a point to which we will come back below.
164
165 \begin{figure}
166 \begin{center}
167 \includegraphics[width=\linewidth]{Shankarappa_allele_freqs_trajectories_syn_nonsynp8}
168 \caption{Synonymous mutations rarely fix in \env, C2-V5: mutation frequency
169 trajectories observed in patient 8~\cite{shankarappa_consistent_1999};
170 Nonsynonymous and synonymous mutations are shown as solid and dashed lines,
171 respectively. Colors indicate the position of the site along the C2-V5 region
172 (red to blue) MAYBE MAKE FIGURE WITH SYNONYMOUS AND NONSYN
173 SEPARATELY. While nonsynonymous mutations frequently fix, very few synonymous
174 mutations do even though they are frequently observed at intermediate
175 frequencies.}
176 \label{fig:aft}
177 \end{center}
178 \end{figure}
179
180 \citet{bunnik_autologous_2008} present a longitudinal dataset on the entire
181 \env~gene of 3 patients at $\sim 5$ time points with approximately 5-20
182 sequences each (see methods). Repeating the above analysis separately on the
183 C2-V5 region studied above and the remainder of \env~ reveal striking
184 differences (see \FIG{fixp}). Within C2-V5, this data fully confirms the
185 observations made in the data set by \citet{shankarappa_consistent_1999} (red
186 line). In the remainder of \env, however, observed synonymous mutations behave
187 as if they were neutral (orange line).
188
189 %ARE OBSERVED SYNONYMOUS MUTATIONS OUTSIDE C2-V5 NEUTRAL? (?? SOME!)
190 %DOES LOSS/FIX CORRELATE WITH CONSERVATION? YES.
191 %MAYBE WE COULD HAVE ONE -- COMPLETELY CIRCULAR -- FIGURE SHOWING LOSS/FIX VS CONSERVATION: SUPPLEMENTARY?
192 %CAN WE LOOK AT THE AVERAGE LEVEL OF CONSERVATION STRATIFIED BY MAX FREQ? TRICKY: the maximal freq is achieved by hitchhiking...
193
194 These observations suggest that many of the synonymous polymorphisms in the part
195 of \env~that includes the hypervariable regions are deleterious, while outside
196 this regions polymorphisms are mostly roughly neutral. Note that this does not imply that all
197 synonymous mutations are neutral -- only those mutations observed at high frequencies tend
198 to be neutral.
199
200 \begin{figure}
201 \begin{center}
202 \subfloat{\includegraphics[width=0.9\linewidth]{Shankarappa_fix_loss_dt_times}
203 \label{fig:fixp1}}\\
204 \subfloat{\includegraphics[width=0.9\linewidth]{Bunnik2008_fixmid_syn_ShankanonShanka}
205 \label{fig:fixp2}}
206 \caption{Left panel: time course of loss and fixation of synonymous mutations
207 observed in a frequency interval $\nu_0$. The ultimate fraction of synonymous
208 mutations that fix as a function of intermediate frequency $\nu_0$ is the
209 fixation probability. Right panel: fixation probability of derived synonymous
210 alleles is strongly suppressed in C2-V5 versus other parts of the {\it env}
211 gene, and of nonsynonymous ones. Data from
212 Refs.~\cite{shankarappa_consistent_1999, bunnik_autologous_2008}.}
213 \label{fig:fixp}
214 \end{center}
215 \end{figure}
216
217 \subsection{Synonymous mutations in C2-V5 tend to disrupt conserved RNA stems}
218
219 One possible {\it a priori} explanation for lack of fixation of synonymous
220 mutations in C2-V5 are secondary structures in the viral RNA. If any RNA
221 secondary structures are relevant for HIV replication, mutations in nucleotides
222 involved in those base pairs are expected to be deleterious and to revert
223 preferentially. Many functionally important secondary structure elements have
224 been characterized, including the RRE~\citep{fernandes_hiv-1_2012} and the 5'
225 UTR pseudoknot interacting with the host
226 tRNA$^\text{Lys3}$~\citep{barat_interaction_1991, paillart_vitro_2002}. It has
227 been suggested early on that parts of the viral genome that has the potential to
228 form stems is better conserved than the
229 remainder~\citep{forsdyke_reciprocal_1995}.
230
231 Recently, the propensity of nucleotides of the HIV genome to form base pairs has
232 been measured using the SHAPE assay (a biochemical reaction preferentially
233 altering unpaired bases)~\citep{watts_architecture_2009}. The SHAPE assay has
234 shown that the variable regions V1 to V5 tend to be unpaired, while the
235 conserved regions between those variable regions form stems. We partition all
236 synonymous alleles observed at intermediate frequencies above 10-15\% depending
237 on their final destiny (fixation or extinction). Subsequently, we align our
238 sequences to the reference NL4-3 strain used in
239 ref.~\citep{watts_architecture_2009} and assign them SHAPE reactivities. As
240 shown in \FIG{SHAPEA} in a cumulative histogram, the reactivities of fixed
241 alleles (red line) are systematically larger than of alleles that are doomed to
242 extinction (blue line) (Kolmogorov-Smirnov test, $P\approx
243 2~\text{\textperthousand}$). In other words, alleles that are likely to be
244 breaking RNA helices are also more likely to revert and finally be lost from the
245 population. As a control, the average over non-observed but potentially
246 available polymorphisms lies between the two curves (green line), as expected
247 (because only some of them will be helix breakers). Furthermore, as a
248 complementary analysis, we split the synonymous mutations in the extended V1-V5
249 region further into conserved and variable regions and found that the biggest
250 depression in fixation probability is observed in the conserved stems, while the
251 variable loops show little deviation from the neutral signature, see
252 \FIG{SHAPEB}.
253
254 In addition to RNA secondary structure, we have considered other possible
255 explanations for a fitness effect of synonymous mutations, in particular codon
256 usage bias (CUB). HIV is known to prefer A-rich codons over highly expressed
257 human housekeeping genes~\citep{jenkins_extent_2003}. Moreover, codon-optimized
258 and -pessimized viruses have recently been generated and shown to replicate
259 better or worse than wild type strains,
260 respectively~\citep{li_codon-usage-based_2012, ngumbela_quantitative_2008,
261 coleman_virus_2008}. We do not find, however, evidence for any contribution of
262 CUB to the ultimate fate of synonymous alleles. Several lines of thought support
263 this result. First of all, although codon-optimized HIV seems to perform better
264 {\it in vitro}, the distance in CUB between HIV and human genes is not shrinking
265 at the macroevolutionary level. Second, within a single patient, we do not
266 observe any bias towards more human-like CUB in the synonymous mutations that
267 reach fixation rather than extinction. Third, it is a common phenomenon for
268 retroviruses to use variously different codons from their hosts, and CUB effects
269 on fitness are thought to be so small that divergent nucleotide composition has
270 been suggested as a possible mechanism for viral
271 speciation~\citep{bronson_nucleotide_1994}. Fourth, CUB in the V1-C5 region is
272 not very different from other parts of the HIV genome, whereas the reduced
273 fixation probability is only observed there. In conclusion, although we cannot
274 exclude an effect of CUB on fitness as a general rule, we expect it to be a
275 minor effect in our context.
276 \begin{figure}
277 \begin{center}
278 \subfloat{\includegraphics[width=0.9\linewidth]{mixed_Shankarappa_Bunnik2008_Liu_fixation_reactivity_Vandflanking_fromSHAPE}
279 \label{fig:SHAPEA}}\\
280 \subfloat{\includegraphics[width=0.9\linewidth]{Shankarappa_fixmid_syn_V_regions.pdf}\label{fig:SHAPEB}}
281 \caption{Watts et al. have measured the reactivity of HIV nucleotides to {\it
282 in vitro} chemical attack and shown that some nucleotides are more likely to
283 be involved in RNA secondary folds. C1-C5 regions, in particular, show
284 conserved stem-loop structures~\citep{watts_architecture_2009}. We show that
285 among all derived alleles in those regions reaching frequencies of order one,
286 there is a negative correlation between fixation and involvement in a base
287 pairing in a RNA stem (left panel). The rest of the genome does not show any
288 correlation (right panel). There might be too few silent polymorphisms in the
289 first place, or the signal might be masked by non-functional RNA
290 structures. Data from Refs.~\cite{shankarappa_consistent_1999,
291 bunnik_autologous_2008, liu_selection_2006}.}
292 \label{fig:SHAPE}
293 \end{center}
294 \end{figure}
295
296
297 \subsection{Deleterious mutations are brought to high frequency by hitch-hiking}
298
299 While the observation that some fraction of synonymous mutations is deleterious
300 is not unexpected, it seems odd that we observe them at high population
301 frequency -- at least in some regions of the genome. The region of \env~ in
302 which we observe deleterious mutations at high frequency, however, is special in
303 that it undergoes frequent adaptive changes to evade recognition by neutralizing
304 antibodies~\cite{williamson_adaptation_2003}. Due to the limited amount of
305 recombination in HIV~\cite{neher_recombination_2010,batorsky_estimate_2011},
306 deleterious mutations that are linked to adaptive variants can reach high
307 frequency~\citep{smith_hitch-hiking_1974}.
308
309 The potential for hitchhiking is already apparent from the allele frequency
310 trajectories in \FIG{aft}, where many mutations appear to change rapidly in
311 frequency as a flock. Deleterious synonymous mutations can be amplified
312 exponentially by selection on linked nonsynonymous sites, a process known as
313 {\it genetic draft}~\citep{gillespie_genetic_2000, neher_genetic_2011}. In order
314 to be advected to high frequency by a linked adaptive mutation, the deleterious
315 effect of the mutation has to be substantially smaller than the adaptive effect.
316 The latter was estimated to be on the order of $s_a = 0.01$ per day~\citep{neher_recombination_2010}.
317 The approximate magnitude of the deleterious effects can be estimated from
318 \FIG{fixp1}, that shows the distribution of times for synonymous
319 alleles to reach the fix or get lost starting from intermediate frequencies. The
320 typical time to loss is of the order of 500 days. If this loss is driven by the
321 deleterious effect of the mutation, this corresponds to deleterious effects of
322 roughly $s_d \sim - 0.002$ per day.
323
324 To get a better idea of the range of parameters that are compatible with the
325 observations and our interpretation, we perform computer simulations of
326 evolving viral populations under selection and rare recombination. For this
327 purpose, we use the recently published package FFPopSim, which includes a module
328 dedicated to intrapatient HIV evolution~\citep{zanini_ffpopsim:_2012}. We
329 analyze many combinations of parameters such as population size, recombination
330 rate, selection coefficient and density of escape mutations, and deleterious effect
331 of synonymous mutation.
332
333 The main result of the simulations is that genetic draft can indeed bring weakly
334 deleterious mutations to high frequencies and result in a dependence of the
335 fixation probability on initial frequency that is compatible with observations.
336 First of all, since neutral mutations are much more likely to rise to high
337 frequency than deleterious ones, the majority of the synonymous mutations needs
338 to be slightly deleterious observe a significant reduction of $P_\text{fix}$.
339 In order to further quantify the reduction in fixation probability, we look at
340 the difference between the neutral curve ($P_\text{fix}(\nu) = \nu$) and the
341 measured fixation probability and calculate its area (see inset of
342 \FIG{simfixpvar}). The minimal and maximal values for this area are zero
343 (neutral-like curve) and 0.5 (no fixation at all), respectively. The HIV data
344 correspond to an area under the diagonal of approximately 0.2 for synonymous
345 changes, and a very small area over the diagonal for nonsynonymous changes.
346 Various simulation curves are shown in \FIG{simfixpvar}. Then, in
347 \FIGS{simheat1}{simheat2}, we explore the parameter space: the combinations that
348 yield areas close to the experimental result are roughly indicated by ellipses.
349 The two crucial parameters that control the fixation probability are the
350 following: (a) the deleterious effects of hitchhikers compared to the beneficial
351 effects of escape mutants, and (b) the density of escape mutations. Intuitively,
352 a higher density of escape mutations (i.e., epitopes) enables a larger degree of
353 genetic draft, because escape mutations from different epitopes start to combine
354 and their effects add up. In \FIG{simheat1}, we show that this is indeed the
355 case in simulations.
356
357 \begin{figure}
358 \begin{center}
359 \subfloat{\includegraphics[width=0.9\linewidth]{fixation_loss_shortgenome_distance_ada_frac_del_eff_coi_various.pdf}
360 \label{fig:simfixpvar}}\\
361 \subfloat{\includegraphics[width=0.9\linewidth]{fixation_loss_shortgenome_area_ada_frac_del_eff_coi_0_01_nescepi_6_heat.pdf}
362 \label{fig:simheat1}}\\
363 \subfloat{\includegraphics[width=0.9\linewidth]{fixation_loss_shortgenome_area_ada_frac_del_eff_coi_0_01_nescepi_6_nonsyn_heat.pdf}
364 \label{fig:simheat2}}
365 \caption{\comment{USE A GOOD FIG WITH THE CORRECT $s_a$!} The depression in $P_\text{fix}$ depends on the deleterious effect size
366 of the synonymous alleles (panel A). Simulations on the escape competition
367 scenario show that the density of selective sweeps and the size of the
368 deleterious effects of synonymous mutations are the main driving forces of the
369 phenomenon. A convex fixation probability is recovered, as seen in the data,
370 along the diagonal (panel B): more dense sweeps can support more deleterious
371 linked mutations. The density of sweeps is limited, however, by the
372 nonsynonymous fixation probability, which is quite close to neutrality (panel
373 C). Moreover, strong competition between escape mutants is required, so that
374 several escape mutants are ``found'' by HIV within a few months of antibody
375 production.}
376 \label{fig:simheat}
377 \end{center}
378 \end{figure}
379
380 However, if hitchhiking is driven by nonsynonymous mutations that are
381 unconditionally beneficial, we should find that nonsynonymous mutations almost
382 always fix once they reach high frequencies -- in contrast with \FIG{fixp} that
383 shows that nonsynonymous mutations fix as if they were neutral. We know,
384 however, that nonsynonymous variation in the variable regions is driven by
385 positive selection. Inspecting the trajectories of nonsynonymous mutations
386 suggest the rapid rise and fall of many alleles. We test two possible
387 mechanisms that are biologically plausible and could explain the transient rise
388 of nonsynonymous mutations: time-dependent selection and within-epitope
389 competition. If the immune system recognizes the escape mutant before its
390 fixation, the mutant might cease to be beneficial and disappear despite its
391 quick initial rise in frequency. In support of this idea,
392 \citet{richman_rapid_2003, bunnik_autologous_2008} report antibody responses to
393 escape mutants. These respones are delayed by a few months, roughly matching the
394 average sweep time of an escape mutant. Alternatively, several different escape
395 mutations in the same epitope can arise almost simultaneously and start to
396 spread. Their fitness benefits are not additive, because each of them is
397 essentially sufficient to escape. As a consequence, several escape mutations rise to
398 high frequency, while the escape with the smallest cost in terms of replication,
399 packaging, etc. is most likely to
400 eventually fix. In simulations, this kind of epistatic interactions within
401 epitopes reduces the fixation probability. The emergence of
402 multiple sweeping nonsynonymous mutations in real HIV infections has been shown
403 previously~\citep{moore_limited_2009, bar_early_2012}.
404 See the supplementary material for examples of successful simulations in both scenarios.
405
406 \section{Discussion}
407 Despite several known functional roles for RNA secondary structure in the HIV
408 genome, synonymous mutations are often used as approximately neutral markers in
409 evolutionary studies of viruses. We have shown that the majority of synonymous
410 mutations in the conserved regions C2-C5 of the \env~gene are deleterious.
411 Comparison with recent biochemical studies of binding propensity of bases in RNA
412 genome suggest that these mutations are deleterious, at least in part, because they disrupt
413 stems in RNA secondary structures. Furthermore, we provide evidence that these
414 mutations are brought to high frequency through linkage to adaptive mutations.
415 The latter mutations are only transiently adaptive, either through a
416 coevolution with the immune system or redundant escape within an epitope.
417
418 Our observations and conclusion rely heavily on longitudinal data in which the
419 dynamics of mutations can be explicitly observed. The fact that deleterious
420 mutations can be brought to high frequencies through hitchhiking underscores
421 the intensity of the coevolution with the immune system. The fact that
422 multiple escape mutations in the same epitope -- as is indeed observed in
423 studies of antibody escape~\citep{moore_limited_2009, bar_early_2012} -- are
424 necessary to explain the patterns of fixation of nonsynonymous mutations points
425 towards a large populations size that rapidly discovers adaptive mutations. A
426 similar point has been made recently by Boltz {\it et al.} in the context of
427 preexisting drug resistance mutations~\citep{boltz_ultrasensitive_2012}.
428
429 The observed hitchhiking highlights the importance of linkage due to infrequent
430 recombination for the evolution of HIV
431 \citep{neher_recombination_2010,batorsky_estimate_2011,
432 josefsson_majority_2011}. The recombination rate has been estimated to be on the
433 order of $\rho = 10^{-5}$ per base and day. It takes roughly $t_{sw} = s_a^{-1}
434 \log \nu_0$ generations for an adaptive mutation with growth rate $s_a$ to rise
435 from an initially low frequency $\nu_0\sim \mu$ to frequency one. This implies
436 that a region of length $l = (\rho t_{sw})^{-1} = s_a / \rho \log \nu_0$ remains
437 linked to the adaptive mutation. With $s_a=0.01$, $l\approx 100$ bases which is
438 consistent with strong linkage between the variable loops and the stems in
439 between. Furthermore, we do not expect hitchhiking to extend far beyond
440 the variable regions consistent with the lack of signal out side of C1-V5. In
441 case of much stronger selection -- such as observed during early CTL escape or
442 drug resistance evolution -- the linked region is of course much larger.
443
444 The functional significance of the insulating RNA structure stems between the
445 hyper variable loops has been proposed
446 previously~\citep{watts_architecture_2009, sanjuan_interplay_2011}.
447 \citet{sanjuan_interplay_2011} have shown that insulating stems are relevant for
448 viral fitness {\it in vivo}. Our analysis is limited by the availability of
449 longitudinal data which requires a focus on the the variable regions of \env.
450 Conserved RNA structures most likely exist in different parts
451 of the HIV genome (several are known). In absence of repeated adaptive substitutions in the vicinity
452 that cause hitchhiking, the deleterious synonymous mutations remain at low
453 frequencies and can only be observed by deep sequencing methods.
454
455 As far as population genetics models are concerned, our study uncovers the
456 subtle balance of evolutionary forces governing intrapatient HIV evolution. The
457 fixation and extinction times and probabilities represent a rich and simple
458 summary statistics to test sequencing data and computer simulation upon. A
459 similar method has been recently used in a longitudinal study of
460 influenza~\citep{strelkowa_clonal_2012}. The propagators suggested in that
461 paper, however, represent ratios between (certain kinds of) nonsynonymous
462 mutations and synonymous ones, hence they are inadequate to investigate
463 synonymous changes themselves. Those authors also conclude that several
464 beneficial mutations segregate simultaneously in influenza, a scenario
465 remarkably similar to our within-epitope competition picture. These results
466 jointly suggest that viral evolution proceeds by multiple concurrent sweeps
467 rather then by successive fixation~\citep{desai_beneficial_2007, neher_rate_2010}.
468
469 Finally, our results emphasize the inadequacy of independent site
470 models of HIV evolution, especially in the light of transient effects on
471 sweeping sites, such as time-dependent selection and within-epitope negative
472 epistasis. Although a final word about which mechanism is more
473 widespread is yet to be spoken, both intuition and biological evidence from the
474 literature support a mixed scenario~\citep{richman_rapid_2003,
475 moore_limited_2009, bar_early_2012}. Note also that, unlike influenza, HIV does
476 recombine if rarely, hence clonal interference as studied in
477 ref.~\citep{strelkowa_clonal_2012} is only a short-term effect.
478
479 \section{Methods}
480 \subsection{Sequence data collection}
481 Longitudinal intrapatient viral RNA sequences were collected for published
482 studies~\citep{shankarappa_consistent_1999,
483 liu_selection_2006, bunnik_autologous_2008} and downloaded from the Los Alamos
484 National Laboratory (LANL) HIV sequence database~\citep{LANL2012}. The sequences from
485 some patients showed signs of HIV compartimentalization into subpopulations and
486 were discarded; a grand total of 11
487 patients with approximately 6 time points each and 10 sequences per time point
488 were analyzed. The time interval or resolution between two ocnsecutive sequences
489 was approximately 6 to 18 months.
490
491 \subsection{Sequence analisys}
492 The good sequences were aligned within each patient
493 via the translated amino acid sequence, using
494 Muscle~\citep{edgar_muscle:_2004}, and to the NL4-3 reference sequence probed
495 by \citet{watts_architecture_2009}. Within each patient, a consensus RNA
496 sequence at the first time point was used to classify alleles as ancestral or
497 derived at all sites. Problematic sites that included large frequencies of gaps
498 were excluded from the analysis, because variable regions are known to be
499 subject to frequent indels, while our analysis is limited to nucleotide
500 substitutions. Time series of allele frequencies were extracted from the
501 sequences.
502
503 The synonymity of a mutation was assigned if the rest of the codon was
504 in the ancestral state and using the standard genetic code. Cases where more
505 than one mutation within the codon was observed were discarded. Slightly
506 different criteria for synonymous/nonsynonymous discrimination yielded similar
507 results.
508
509 \subsection{Fixation probability and secondary structure}
510 For the estimate of times to fixation/extinction, polymorphisms were
511 binned by frequency and the time to reaching the first boundary (fixation or
512 extinction) was stored. For the fixation probability, the long-time limit of the
513 resulting curves was used, excluding polymorphisms that arose late in the
514 clinical history (and would have had no time to reach either boundary).
515
516 For the correlation analysis with RNA secondary structure, the SHAPE scores were
517 downloaded from the journal website~\citep{watts_architecture_2009}. By virtue
518 of the alignment of the longitudinal sequences with the reference used by
519 \citet{watts_architecture_2009}, SHAPE reactivities were assigned to most sites.
520 Problematic assignments in indel-rich regions were excluded from the analysis.
521 In order to restrict the analysis to synonymous polymorphisms, a lower frequency
522 threshold of 0.15 was used (other thresholds yielded the same results). Since
523 very few polymorphisms hitchhike beyond, say, a frequency of 0.5, this pool is
524 enriched for to-be-lost mutations; hence the "lost" curve in \FIG{SHAPEA}
525 contains much more points than the "fixed" one.
526
527 The V loops and flanking regions were identified manually starting from the
528 annotated reference HXB2 sequence from the LANL HIV database~\citep{LANL2012}. A
529 similar approach was used to label the C2-V5 region sequenced in
530 ref.~\citep{shankarappa_consistent_1999}.
531
532 \subsection{Computer simulations}
533 Simulations were performed using the recently published software
534 FFPopSim~\citep{zanini_ffpopsim:_2012}. Both full-length HIV genomes and
535 \env{}-only simulations were performed and yielded comparable results. For each
536 set of parameters, approximately 100 simulation runs were averaged over. In each
537 run, a random fitness landscape with specified statistical properties (e.g.
538 density of beneficial sites, average deleterious effect of synonymous changes) was generated.
539 Although the curves shown in \FIG{simfixpvar} are not very smooth, small
540 parameter changes resulted in overall consistent trends across many repetitions.
541
542 For the discussion of simulation parameters, the areas below or above the neutral
543 diagonal were estimated from the binned fixation probabilities using the linear
544 interpolation between the bin centers. This measure is sufficiently precise for
545 our purposes, because the HIV data are quite scarse themselves.
546
547 \section*{Acknowledgements}
548 \comment{to be written\dots}
549
550
551 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
552 \bibliographystyle{natbib}
553 \bibliography{bib}
554 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
555 \end{document}
556 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
557