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