Biased mutational pattern and quasispecies hypothesis in H5N1 virus.

Like other RNA viruses, influenza viruses are subject to high mutation rates. Carrying segmented RNA genomes, their genetic variability is even higher. We aimed at analyzing the mutational events occurring during the infection of chickens by the Highly Pathogenic Avian Influenza (HPAI) H5N1 virus. We therefore studied the different sequences of two surface proteins, hemagglutinin (HA) and neuraminidase (NA), as well as two internal proteins, PB2 and NS. Three organs (lung, spleen, brain) were obtained from a chicken, experimentally infected with a lethal dose of HPAI H5N1 virus. Cloning these PCR fragments enabled us to investigate the mutations undergone by the virus after several replicative cycles. The first outcome is the presence of a strong mutational bias, resembling host-driven ADAR1 adenosine deamination, which is responsible for 81% of all mutations. Whereas the frequency of RNA dependent RNA polymerase-related mutations is compatible with the survival of the virus, the ADAR1-like activity usually strongly increases the mutation frequency into a level of "error catastrophe" in theory incompatible with virus survival. Nevertheless, the virus was successfully infective. HPAI H5N1 virus displayed traits in agreement with the quasispecies theory. The role of this quasispecies structure in successful infection and the superposition with the ADAR1-like response is discussed.

Introduction

A feature common to all RNA viruses is their high genetic variability due to high mutation rates of around 10−4–10−5 misincorporations per nucleotide per replication cycle (Drake and Holland, 1999), facilitated by the lack of proof-reading activity. This high genetic variability has been associated, among other things, with the concept of quasispecies for the last two decades (Eigen, 1993, Eigen and McCaskill, 1988). The quasispecies approach is based on the concept of quasi-infinite error-prone populations displaying a high mutation rate (Lauring and Andino, 2010, Wilke, 2005). In this situation a virus will generate a large progeny with genomes highly similar to the parental molecule, each one differing from this original genome by one or few mutated sites. Other rounds of replication will generate additional variants also bearing limited mutations which with time will occupy the “sequence space”, a geometrical representation of all possible sequences. The outcome is a “cloud” of variants behaving almost (quasi) as a unique entity (species), or quasispecies (Lauring and Andino, 2010, Wilke, 2005). However, this is the outcome in the mathematical modeling of the growth of a quasi-infinite population. In a true biological context, this dynamics will be strongly counterbalanced by negative, purifying, selection to generate smaller clouds of closely related sequences, each one differing from the master sequences by very few mutations. Quasispecies theory is therefore a specific case of population genetics and mutation-selection concept in which a strong negative selection generates a group of closely related sequences described as “coupled” in the sequence space (Wilke, 2005). Several studies have investigated the applicability of this model to various RNA viruses, mostly to human immunodeficiency virus (HIV) (Bordería et al., 2010, Briones et al., 2006, Briones et al., 2003), foot-and-mouth disease virus (FMDV) (Domingo et al., 2002, Ojosnegros et al., 2010), and hepatitis viruses (Domingo and Gomez, 2007, Figlerowicz et al., 2010, Fishman and Branch, 2009, Sulbaran et al., 2010). The quasispecies structure of influenza virus was also investigated through different approaches such as PCR amplicons (Chen et al., 2010, Kongchanagul et al., 2008, Iqbal et al., 2009, Watanabe et al., 2011), next generation sequencing (Kuroda et al., 2010, Ramakrishnan et al., 2009, Höper et al., 2009, Höper et al., 2011), realtime RT-PCR (Van Borm et al., 2010) and pyrosequencing (Lackenby et al., 2008, Duwe and Schweiger, 2008). Other studies were conducted on viruses such as the severe acute respiratory syndrome-associated coronavirus (SARS-CoV) (Seligman, 2007), flaviviruses (Brackney et al., 2010, Jerzak et al., 2005) and plant viruses (Hasiów-Jaroszewska et al., 2010).

Like all influenza viruses, the Highly Pathogenic Avian Influenza (HPAI) H5N1 virus is highly variable. Its segmented RNA genome allows it to evolve rapidly via mutation, reassortment, recombination and mixed infection. The consequences of this high genetic variability are numerous: resistance to drugs, modification of virulence, antigenic drift, escape from the immune response, ability to adapt to new host species, etc. Given the animal health and public health hazard that the HPAI H5N1 virus presents, these features have been investigated and documented for over 10 years. Nonetheless, nothing is known about their behavior as defined in the quasispecies framework. Furthermore, the HPAI H5N1 virus being responsible for acute fatal infections, as opposed to chronic infections caused by most viruses previously studied in the quasispecies model, data obtained from the latter do not necessarily apply to the HPAI H5N1 virus.

Here, we studied the viral subpopulations of an HPAI H5N1 inoculum and their in-host mutations when infecting an avian host, specifically in three organs. We described the mutational pattern observed, and confronted it to the conditions of the quasispecies concept. We also reported the presence of a very strong mutational bias, present both in the inoculum and in infected organs, potentially associated to a host ADAR1-like activity.

Materials and methods

Biosafety and ethics statement

All experiments using H5N1 virus were performed within the biosafety level 3 laboratory (BSL-3) of Institut Pasteur in Cambodia. All animal experiments were carried out in the animal biosafety level 3+ laboratory (ABSL-3+) of the Institut Pasteur in Cambodia, complying with the Animal Committee regulations of Institut Pasteur in Paris, France, in accordance with the EC 86/609/CEE directive, and approved by the Animal Ethics Committee of Institut Pasteur in Cambodia.

Virus

The HPAI H5N1 virus used in these studies was the virus A/Chicken/Cambodia/LC1AL/2007 (GenBank accession numbers HQ200574 to HQ200581). The virus stock was obtained by two passages in Specific Pathogen Free (SPF) 9–11-days old embryonated hen eggs. Incubation time was of 48 h at 37 °C. The eggs were then chilled overnight at 4 °C, and the amnio-allantoic fluid was then harvested and stored at −80 °C until further use. Virus titer was determined by calculating the 50% egg infectious dose (EID50) per milliliter of virus stock. Titration endpoints were calculated using the Reed and Muench method (Reed and Muench, 1938).

Animals

SPF chickens, of a local Cambodian breed, 4–6-weeks old, were kindly provided by the National Veterinary Research Institute of Cambodia (NaVRI). During the experiments, all birds were housed in self-contained isolation units (Cap Engineering B.V. HM1500 Isolator, Montair, Sevenum, The Netherlands) ventilated under negative pressure with High-Efficiency Particulate Air (HEPA)-filtered air, with 12 daily hours of artificial lighting. Commercial food and water were provided ad libitum.

Experimental setting

Six SPF chickens were inoculated via oral, ocular and cloacal routes with 106 EID50 of virus in 500 μL of sterile water on day 0 (D0). All birds were checked for expression of clinical signs of illness on a daily basis. Post-mortem examination and collection of samples were conducted on all animals after their death. The samples collected post-mortem included organs such as spleen, lung, and brain. All samples were stored in phosphate buffer saline (PBS) at −80 °C prior to testing.

RNA extraction and retro-transcription (RT)

Four different specimens were used for this study: the inoculum (initial viral solution inoculated to the chickens), and spleen, lung and brain samples collected from an infected chicken at the time of necropsy. After homogenization of the organs, viral RNA was extracted using the QIAamp® Viral RNA Mini Kit (QIAGEN GmbH, Hilden, Germany) according to the manufacturer’s recommendations. Complementary DNA (cDNA) fragments targeting the hemagglutinin (HA), neuraminidase (NA), PB2 and NS genes were synthesized using the following primers: H5R-1773 for the HA gene, NA- for the NA gene, PB2R-2341 for the PB2 gene, NSR-890 for the NS gene. All primers used in this study are listed with their corresponding sequences in Supplementary Table S1. The reverse transcription reaction was performed using AMV Reverse Transcriptase (Promega, Madison, WI, USA) in a reactional mix containing 6 μL of AMV Reverse Transcriptase 5× Reaction Buffer (Promega), 1 μL of dithiothreitol (DTT, 0.1 M), 1 μL of primer (working stock concentration: 10 μM), 1 μL of mixed dNTPs (10 mM), 1 μL of AMV Reverse Transcriptase (Promega), 5 μL of RNAse-free water, and 5 μL of sample’s RNA. The reaction was carried out in a MyCycler Thermal Cycler (Bio-Rad, Applied Biosystems, Foster City, CA, USA): 5 min at 25 °C, 60 min at 50 °C, 15 min at 70 °C and 10 min at 4 °C.

Polymerase chain reaction (PCR)

The HA, NA, and PB2 genes were amplified in overlapping halves named HAa and HAb, NAa and NAb, and PB2a and PB2b. The NS gene was amplified in a single segment. Primers used for this PCR step included: H5F-1 and H5R-1111 for the HAa segment; H5F-751 and H5R-1773 for the HAb segment; N1F-1 and N1R-835 for the NAa segment; N1F-459 and NA- for the NAb segment; PB2F-1 and PB2R-1250 for the PB2a segment; PB2F-1105 and PB2R-2341 for the PB2b segment; NSF-1 and NSR-890 for the NS segment. Each reaction was conducted using the GoTaq Flexi DNA Polymerase (Promega) in a reactional mix containing 11 μL of 2× Reaction Mix for SuperScript III Platinum One-Step Quantitative RT-PCR system (Invitrogen, Carlsbad, CA, USA), 1.5 μL of 5′ primer (H5F-1, H5F-751, N1F-1, N1F-459, NSF-1, PB2F-1; working stock concentration: 10 μM), 1.5 μL of 3′ primer (H5R-1111, H5R-1773, N1R-835, NA-, NSR-890, PB2R-1250, PB2R-2341; 10 μM), 0.25 μL of GoTaq Flexi DNA Polymerase (Promega), 33.75 μL of water, and 2 μL of cDNA. The PCR was then carried out in a MyCycler Thermal Cycler (Bio-Rad): 1 cycle of 2 min at 95 °C, followed by 40 cycles of 30 s at 95 °C, 30 s at 55 °C and 1 min at 72 °C, and a final cycle of 5 min at 72 °C, before reaching a final storage temperature of 4 °C. Amplified PCR products were then purified using QIAquick Gel Extraction kit (QIAGEN GmbH), according to the manufacturer’s recommendations.

Cloning and sequencing

All amplicons were ligated into pCR 2.1-TOPO plasmid vector (Invitrogen) and used for transformation using the TA TOPO Cloning Kit with One Shot MAX Efficiency DH5α-T1 Escherichia coli (Invitrogen), following the manufacturer’s instructions. Seventeen to 20 clones were selected per PCR amplicon (locus) and per sample, and sent for sequencing to a contract sequencing facility (Macrogen, Seoul, Korea). Sequencing reactions were performed in a PTC-225 Peltier Thermal Cycler (MJ Research, Waltham, MA, USA) using ABI PRISM BigDye Terminator Cycle Sequencing Kit with AmpliTaq DNA Polymerase (Applied Biosystems, Foster City, CA, USA), following the protocols supplied by the manufacturer. Single-pass sequencing was performed on each template using the following primers: H5F-1, H5F-219, H5R-515 and H5R-1111 for the HAa segment; H5F-751, H5F-1150 and H5R-1773 for the HAb segment; N1F-1, N1F-459, N1R-560, and N1R-835 for the NAa segment; N1F-459, N1F-814, N1R-835, and NA- for the NAb segment; PB2F-1, PB2R-370, PB2F-940, PB2R-960, and PB2R-1250 for the PB2a segment; PB2F-1105, PB2R-1753, PB2F-1860 and PB2R-2341 for the PB2b segment; NSF-1, NSF-467, NSR-580, and NSR-890 for the NS segment. Final electrophoresis was performed on an ABI PRISM 3730XL Analyzer (Applied Biosystems).

Sequence alignments and phylogenetic trees

For each clone, sequences on both strands were assembled into contigs to obtain a consensus sequence using CLC Main Workbench 5.7.1 (CLC bio, Aarhus, Denmark). The clones corresponding to each locus, i.e. amplicon, were aligned with CLUSTAL W (Thompson et al., 1994) in BioEdit Sequence Alignement Editor 7.0.9.0 (Ibis Therapeutics, Carslbad, CA, USA). Each locus was treated independently. Because of the very low number of mutations per clone (typically 1 or 2) with respect to the consensus, phylogenetic trees from nucleotide sequences were generated with the Neighbour-Joining method with observed differences and 100 bootstrap replicates, using MEGA 5.05 (Tamura et al., 2011). Trees were edited using MEGA 5.05 and Inkscape 0.48 (Free Software Foundation, Boston, MA, USA). Sequences used were all nucleotidic variants, truncated when necessary so as to be of equal lengths. All positions described on the clones sequenced here refer to H5N1 numbering as defined by alignment with the sequence of A/Goose/Guangdong/1/1996 virus (GenBank accession number AAD51922 to AAD51931).

Data analysis and mutation models

Sequences showing nonsense and/or frameshift mutations were excluded from all subsequent analyses. The number of mutations per clone was modeled in the following way: for a given clone of size L (typically L  ∼ 800 bp), the number of observed mutations between the clone sequence and the consensus sequence was called n dif. Because of the low value of n dif (typically 1 mutation), one may assume that there is no saturation and that one observed difference corresponds exactly to one mutation. Under the hypothesis of independence of the mutation sites, n dif can be simply modeled by a Poisson process with parameter λ  = f ·  L where f is the mutation frequency

Considering now N independent clones (of the same size L), the number N(k) of clones with exactly k mutations follows a binomial distribution with parameter n  =  N and . The mean and variance of N(k) are therefore given by: and . P-values of observed data can be readily calculated from the binomial distribution.

The same analysis can be performed on the number of mutations per site on a multiple alignment of N clones of same size L. This is simply done by permuting the roles of N and L in the previous equations. The number of mutations per site follows a Poisson distribution with parameter λ′ =  f  ·  N where N is typically 20 (for one organ or inoculum) and f is the same mutation frequency as before. Considering a set of L independent sites, the number M(k) of sites with exactly k mutations follows a binomial distribution with parameter n  =  L and .

A more detailed derivation of the distribution of N(k) and M(k) is given in the Supplementary Annex in the case of two superposed Poisson processes.

Results

Animal infection and sample collection

All six chickens, inoculated via oral, ocular and cloacal routes, died of H5N1 infection, within a mean death time (MDT) of two days, without exhibiting any clinical symptoms. No gross lesions were observed on the organs collected post-mortem. The samples chosen for the study were all collected from one single chicken, randomly selected. Although the passage through eggs was reported to change the quasispecies composition (Ramakrishnan et al., 2009), it did not interfere in this work since the objective was not to study the evolution of a given sequence through successive infections but to compare the mixture present in the inoculum with the sequences found in different organs. This approach is quite similar to what could be found in the wild with interspecies transmission.

Phylogeny of individual loci and divergence of amino acid sequences

The phylogenetic trees are represented on Supplementary Fig. S1. All were left unrooted. Bootstraps values were generally low (⩽65) although few exceptions were observed, e.g. on the NAb tree with values of 85, 87, 88 and 99. Phylogenetic analysis of sequenced clones yielded similar results for each individual locus. No clear clusters appeared and all trees displayed a similar unresolved star structure when shown under a radial shape (Supplementary Fig. S1). This topology was observed for all loci, all organs and inoculum.

Isolates variation according to locus and organ

Mutation matrix and mutation frequencies

When considering, for each locus, the distribution of the number of mutations observed per sequence with respect to the consensus, or master sequence, a global mutation frequency of 1.58 × 10−3substitutions per nucleotide (nt) was obtained (Table 1 ). Owing to the low number of mutations observed, there was no saturation and each variation corresponded to a mutation. However, the mutation matrix was not symmetrical and was strongly biased towards the presence of A > G (259/556) and T > C (190/556) mutations which represented 81% of all mutations observed (Table 2 ). This bias was observed in all organs as well as in the inoculum (data not shown). Despite this strong mutational bias towards G and C, the average G + C content of sequences analyzed was 44%. These two classes of mutations were complementary, suggesting that only one class of mutation actually occurred but both on the genome and anti-genome. Such strongly biased A > G mutational patterns cannot be explained by RNA dependent RNA polymerase (RdRP) misincorporations but are usually associated with Adenosine Deaminases that act on RNA (ADAR) (Bass, 2002). However, the landmark of ADAR1 mutation is 5′ neighbor position preference (usually U = A > C > G) (Bass, 2002), a feature which was not observed here with the exception of 3 data out of 16 (data not shown), suggesting that the activity observed might be an ADAR1-related deamination process rather than strict ADAR1. Therefore, three different mutation frequencies were calculated (Table 1). The first one was simply the overall mutation frequency: f tot  = (number of mutations)/(number of nucleotides) = #mutations/#nt = 1.58 × 10−3 mutations/site. The second one tried to estimate the fraction of mutations due to the deamination process solely. We assumed that (i) deamination only produced A > G and T > C mutations; (ii) the total number of observed A > G + T > C mutations was the sum of mutations due to deamination on one hand and to RdRP misincorporations on the other hand, and (iii) the mutation matrix of RdRP misincorporations was symmetrical, i.e. that the actual number of A > G + T > C due to RdRP equalled the observed number of G > A + C > T transitions. Therefore the mutation frequency due to the deamination process solely was: fdeam = (A > G + T > C−G > A−C > T)/(A + T) = 1.88 × 10−3 mutations/A-or-T-sites. Finally we estimated the mutation frequency due to RdRP solely by substracting the deamination component from A > G and T > C, using the same assumption as before: fpol = (#mutation−A > G−T > C + G > A + C > T)/#nt = 5.19 × 10−4 mutations/nt. This corrected mutation frequency was similar to the values ranging from 3.6 × 10−4 to 5.0 × 10−4 mutations/nt reported by Iqbal et al. for H5N1 virus in chickens and in very similar experimental conditions (Iqbal et al., 2009).

Table 1

Mutation frequencies by locus.

Locus	#seqa	#ntb	A + Tc	#mutd	A > G + T > Ce	G > A + C > Tf	ftotg	fdeamh	fpoli
HAa	71	63,048	36,920	99	87	9	1.57 × 10−3	2.11 × 10−3	3.33 × 10−4
HAb	73	54,093	31,463	101	85	13	1.87 × 10−3	2.29 × 10−3	5.36 × 10−4
NAa	75	51,525	29,250	71	55	10	1.38 × 10−3	1.54 × 10−3	5.05 × 10−4
NAb	65	48,945	27,755	71	56	12	1.45 × 10−3	1.59 × 10−3	5.52 × 10−4
NS	71	44,943	23,643	49	39	6	1.09 × 10−3	1.40 × 10−3	3.56 × 10−4
PB2a	54	45,360	24,678	94	68	17	2.07 × 10−3	2.07 × 10−3	9.48 × 10−4
PB2b	49	44,688	25,186	71	59	9	1.59 × 10−3	1.99 × 10−3	4.70 × 10−4
Total	458	352,602	198,895	556	449	76	1.58 × 10−3	1.88 × 10−3	5.19 × 10−4

Total number of sequences analyzed.

Total number of nucleotides.

Total number of A and T bases.

Total number of mutations observed.

Total number of A > G and T > C mutations observed.

Total number of G > A and C > T mutations observed.

Overall frequency: #mut /#nt.

Estimated frequency of deamination process: (A > G + T > C−C > T−G > A)/A + T in mutations/A-or-T.

Estimated frequency of mutations caused by the RNA dependent RNA polymerase only : (#mut - A > G−T > C + G > A + C > T)) #nt in mutations/nt.

Table 2

Mutation matrix (all loci).

The mutation matrix represents mutations of nucleotides of the left column into nucleotides of the upper row.

It should be noted that the overall mutational process was actually modeled by two independent Poisson processes, one due to RdRP (with frequency f pol) and one due to deamination (with frequency f deam). As for the number of mutations per clone, the resulting process was therefore a Poisson process with parameter (see Supplementary Annex for details).

From these frequencies and from the total genome size of ∼13.5 kb the fraction of genomes, i.e. sequenced clones, displaying at most k mutations (cumulative Poisson) was calculated. These clones displaying no or very few mutations could correspond to the viable infective progeny compatible with the negative selective pressure associated to quasispecies. These values are given in Table 3 for k  = 0, 1, 2.

Table 3

Fraction of genomes displaying at most k mutations, with and without ADAR1-like deamination.

k	All mutationsa	RdRp onlyb
0	5.45 × 10−10	9.06 × 10−4
1	1.22 × 10−8	7.25 × 10−3
2	1.36 × 10−7	2.95 × 10−2

Considering all mutations, with a mutation frequency of 1.58 × 10−3 mutations/nucleotide and a total genome size of 13.5 Kb.

Considering mutations caused by the RNA dependent RNA polymerase only, with a mutation frequency of 5.19 × 10−4 mutations/nucleotide and a total genome size of 13.5 Kb.

Distribution of mutations per clone

The distribution of the number of mutations per clone was compared to a Poisson distribution in order to determine whether the occurrence was stochastic or adaptive. Data are presented in Supplementary Table S2. For each locus, the distribution was considered for each organ (lung, spleen, brain) and inoculum separately as well as for the set of all sequences (“total”). The comparison of the observed mutations per locus for all sets fitted quite well with a Poisson distribution. Few variations from Poisson distribution were observed, most of them occurring on PB2a. With the exception of PB2 there was no significant difference between loci with respect to the number of mutations. Similarly, there was no significant difference between organs nor between organs and inoculum. This suggested that these mutations were purely stochastic and therefore not fixed by selective pressure. This was also supported by two observations: (i) the distribution of the number of mutations was uniform across the three reading frames (Supplementary Table S6), whereas one would expect more mutations on the third frame if selective fixation had occurred; (ii) in the same way, the observed ratio of non-synonymous over synonymous mutations was 367/189 = 1.94, a value very close to the expected ratio of 2 for purely random mutations (the actual expected value for random mutations on these particular sequences was 2.09). Finally, we performed the same analysis by removing all the A > G and T > C mutations (Supplementary Table S3). Again, the number of mutations per clone did not show any discrepancy from a purely stochastic Poisson process.

Distribution of mutations by site

In order to try to detect a potential fixation process, the previous analysis was complemented by the analysis of the number of mutations per site, for a given locus. The comparison of this distribution with Poisson is presented in Supplementary Table S4 (for all mutations) and C2 (for all mutations except A > G and T > C). Although some few variations were observed, the distribution per site displayed a good fit with the Poisson distribution supporting further the stochastic occurrence of the mutations and the absence of fixation with time. No significant difference could be observed between organs, nor between organs and inoculum, suggesting that the occurrence of mutations followed the same stochastic trend in all organs and in eggs. This also confirmed that there was no evolution from the inoculum by fixation of mutations through time.

Discussion

The main objective of this work was to analyze the mutational events occurring during the infection of chickens by the highly variable and heterogeneous Highly Pathogenic Avian Influenza (HPAI) H5N1 virus. However, this work is based on the use of amplicons in order to cover a large portion of the viral genome. A limitation is thus that only a small proportion of the overall population of variants was analyzed. This work should therefore be completed by an in-depth high-throughput approach which allows deeper coverage of the whole population.

The study of this virus in multiple organs was made possible because of the fatal systemic infection it causes in chicken. Although the mean death time of chickens was short (2 days), the average duration of the replication cycle of the influenza A virus being of 5–6 h (Baccam et al., 2006, Zambon and Potter, 2009), it allowed at least 6–8 cycles to be accomplished during the span time of the experiment. This makes the experimental data reported in this work comparable to what would be expected in a situation of quasispecies evolution. The quasispecies nature of influenza has been investigated on a variety of serotypes, including the pandemic 2009 influenzaA virus H1N1 but also H3N1 and H5N1. These investigations were performed by using a wide range of complementary approaches including in depth next generation sequencing (Kuroda et al., 2010, Ramakrishnan et al., 2009, Höper et al., 2009, Höper et al., 2011), pyrosequencing (Lackenby et al., 2008, Duwe and Schweiger, 2008), realtime RT-PCR (Van Borm et al., 2010) and PCR amplicons sequencing (Chen et al., 2010, Kongchanagul et al., 2008, Iqbal et al., 2009, Watanabe et al., 2011). This work relates to this last approach.

The first expected feature is the very limited evolution of variants from the master sequence in order to occupy the sequence space. Furthermore, the phylogenetic analysis of the different loci on a per-organ basis yielded trees with similar unresolved topologies, i.e. non-structured or star-shaped trees, indicative of parallel and independent evolution of each clone from the same original master sequence (Biebricher and Eigen, 2006, Lauring and Andino, 2010). The data reported in this work indicate that the number of observed mutations is of the order of one mutation per clone, corresponding to an overall mutation frequency of 1.58 × 10−3 mutations/site. Owing to a total genome size of ∼13.5 kb, this would lead to an overall mutation frequency of 21.3 mutations/genome. This value might be too high to be compatible with the creation of viable variants (Crotty et al., 2001). However, the purely replicative part of the mutations was calculated to be about 5.19 × 10−4 mutations/site which would lead to a genomic frequency of 7.0 mutations/genome. Although this value is still rather high, it would lead to a fraction of e−7.0  ∼ 1 ‰ of genomes displaying no mutation and 3% of genomes with at most 2 mutations (Table 3) during the life time of infected chickens. This indicates that the mutation frequency attributable to RdRp is compatible with the survival of viable viruses in the population. Furthermore, both the overall mutations and the restricted mutations follow a random (Poisson) process with no shared mutations or specific patterns.

The main proportion of the observed mutations is attributable to a host-driven deamination process. HPAI H5N1 virus displayed a very strong A > G mutational bias. 81% of the observed mutations were A to G and T to C transitions whereas reverse G to A and C to T transitions as well as transversions were more limited (13.6%). A very similar bias was recently reported on SRAS, although to a lower extent with A > G and T > C substitutions making about 60% of all mutations observed (Eckerle et al., 2010). This bias was never reported for in vitro activity of AMV reverse transcriptase despite extensive misincorporation survey (Operario et al., 2005, Skinner and Eperon, 1986, Taube et al., 1998). This strong mutational bias is suggestive of the activity of a family of host cell proteins known as Adenosine Deaminases that Act on RNA or ADAR (Bass, 2002). ADAR-mediated RNA editing occurs through the conversion of adenosine to inosine in both coding and noncoding regions of double-stranded RNA. This finally leads to the editing into guanosine. The presence of T to C transitions on the sequence may correspond in fact to A to G mutations on the complementary strand indicating thus that 81% of the observed mutations would be ADAR-like mediated RNA editing. However, ADAR editing is dependent on two factors, the first one being specificity which is directly dependent on the stability of the double-stranded RNA (Bass, 2002). ADAR editing usually occurs on a limited number of adenosine residues but can extend up to 50% under some situations of high stability of the double-stranded RNA editing, a process known as hypermutation (Bass, 2002, Suspène et al., 2011). In this work, the overall editing is only 0.2 ‰ of all adenosines available and an average of two residues per sequence. The second main trait of ADAR is the neighbor position preference. ADAR-mediated editing displays a slight preference for certain adenosines over others depending upon their immediate neighbor base. ADAR1 displays no 3′ preference but a clear 5′ preference where U ∼ A > C = G (Bass, 2002, Suspène et al., 2011). If the data reported here underline the overwhelming presence of A > G editing among the HPAI H5N1 virus variants isolated from chicken, this does not correspond to true ADAR1 activity owing to the lack of neighbor preference and should be referred to as ADAR1-like deamination. However, the inoculum, which was coming from egg passage, displayed the same bias indicating thus that this ADAR1-like process is also active in hen eggs. The data reported in this work most likely result from the superposition of two mechanisms: the RdRP intrinsic errors during the replication of the virus on one hand and an antiviral ADAR-like activity from the host on the other hand.

ADAR-1 proteins also display antiviral activity and the p150 isoform was shown to inhibit measles virus but also to protect against several orthomyxoviruses and other paramyxoviruses (Samuel, 2011, Ward et al., 2011). This is in particular true for influenza A. ADAR1-depleted mice were more susceptible to influenza virus (Ngamurulert et al., 2009) while ADAR1 proteins were produced in higher amounts in influenzaA infected cells (Emmott et al., 2010). Furthermore, H5N1 NS1 protein was reported to specifically interact with ADAR1 proteins, to inhibit interferon production and to manipulate host proteins to facilitate infection and confer an unusual severity (Ngamurulert et al., 2009, Randall and Goodbourn, 2008, Samuel, 2011). Antiviral ADAR1 activity seems to be based at least on two complementary actions: interferon activity and accumulation of deleterious mutations in virus protein genes. This brings the mutation frequency to a level which might be too high to maintain viral genome infectivity (1 for 10 million) pushing thus viral population towards a model of “error catastrophe” (Crotty et al., 2001).

Despite the presence of an ADAR1-like antiviral response and an overall mutation frequency beyond the viability limit, the H5N1 virus studied in this work was able to infect the chickens and kill them within two days. This demonstrates its high virulence in a susceptible host but also its ability to efficiently counteract the host innate immunity. The answer might be found in some additional key concepts in the quasispecies theory. The first one is mutational robustness (Van Nimwegen et al., 1999), in which a population minimizes the detrimental impact of mutations by gathering sequences with low probability of bearing deleterious mutations (Bornberg-Bauer and Chan, 1999, Lauring and Andino, 2010, Van Nimwegen et al., 1999, Wilke, 2001, Wilke and Adami, 2003). The second concept to consider is that of functional cooperation (Lauring and Andino, 2010, Más et al., 2010) in which a global advantage is given, here pathogenicity versus host defense, to the cloud of variants by the association of complementary functions borne by separate variants. Although intuitively sound and in line with the reality, which is that the virus was successful in evading host defenses and kill the chickens, this is difficult to demonstrate and must thus remain at a somehow speculative level. This is nevertheless an expected outcome of the quasispecies concept. Unlike Darwinian evolution which relies on the preexistence in the population of rare but suitable mutations which can be fixed by selection, and thus can be traced, there is no such preexisting fixed mutation in quasispecies. There is instead a continuous stochastic occurrence of mutations in the overall population and the cooperation of a limited number of variants to provide a temporary functional advantage in a given situation (e.g. circumventing host defense, adapting to a new host, etc.). These mutations are temporary but through cooperation and mutational robustness they may confer a functional advantage in a novel environment or host. The occurrence of such “preferred” sites was recently reported in a comparative analysis of H5N1 virus mutation frequency in different hosts (Iqbal et al., 2009). In the work presented here, several of the mutated sites observed were previously described in the literature as bearing mutations of special functional importance (Table 4 ). These regions fell roughly into two categories: those involved in adaptation to human host, and those involved in virulence and pathogenicity. Many substitutions were observed on well-known sites in the hemagglutinin protein, such as antigenic sites, receptor-binding sites, the cleavage site, or the fusion peptide. Substitution of the glutamine for an arginine on position 192 in the H5 HA1 protein was shown to be responsible for the shift from sialic acid α2,3 galactose (SAα2,3Gal) to SAα2,6Gal recognition (Yamada et al., 2006) and thus adaptation to human hosts (Matrosovich et al., 1999, Matrosovich et al., 1997). Mutations described as potentially affecting membrane fusion (Chen et al., 1999, Langley et al., 2009, Skehel and Wiley, 2000) were also detected in the HA2 protein. The three substitutions observed on relevant sites of the NA protein were located on amino acids involved in the neuraminidase active site (Collins et al., 2008, Russell et al., 2006, Yen et al., 2005, Yen et al., 2007). As for the NS protein, the main element to report was the mutation Asp92Gly. Better known at this position is the Asp92Glu substitution, responsible for increased virulence and resistance to antiviral cytokines (Seo et al., 2002, Zhou et al., 2010). Position 92 is also located within the translation initiation factor eIF4GI binding domain (positions 81–113) of NS1. Truncations of this domain were shown to have an impact on virulence and pathogenicity of H5N1 viruses (Zhou et al., 2010). It would be out of scope to speculate on a possible effect but their occurrence on the limited number of clones analyzed in this work out of the numerous possibilities offered by the sequence space make them worth to be at least mentioned.

Table 4

Mutations observed in the region of functional importance.

Protein	Sample	Mutation observed	Mutation previously reported	Altered function or relevant site affected	References
HA1	Lung	S128P	None	Antigenic site (position 133 in H3-numbering)	Kaverin et al., 2007, Kaverin et al., 2002, Philpott et al., 1990, Stevens et al., 2006
	Inoculum	K153R	None	Antigenic site (position 157 in H3-numbering)
	Spleen	T159Ta	T159A (T160A in reference H5 numbering)	Receptor specificity (avian vs mammal)	Gao et al., 2009, Stevens et al., 2006
	Inoculum	N182Na	N182 K	Receptor specificity (avian vs human)	Yamada et al. (2006)
	Spleen	N182T
	Inoculum	Q192R	Q192R
	Brain	S223G	S223N (S227N in H3-numbering)	Receptor specificity (avian vs human) and/or antibody-antigen binding	Gambaryan et al., 2006, Hoffmann et al., 2005, Shinya et al., 2005
	Brain	P321Pa	None	Cleavage site	Bosch et al., 1981, OIE World Organisation for Animal Health, 2009, Perdue et al., 1997, Skehel and Wiley, 2000, Steinhauer, 1999, Walker and Kawaoka, 1993
	Lung	R326Ra
HA2	Brain	F9S	None	Fusion peptide	Langley et al., 2009, Skehel and Wiley, 2000
	Lung	Y34H	None	Membrane fusion intermediate	Chen et al. (1999)
	Spleen	A35Aa
NA	Lung	R118G	None	NA active site	Collins et al., 2008, Russell et al., 2006, Yen et al., 2005, Yen et al., 2007
	Inoculum	D151N
	Inoculum	R293G	R292K
PB2	Brain	A44Aa	A44S	Adaptation to human host	Miotto et al. (2010)
	Lung	A588T	A588I
	Spleen	A661Aa	A661T
NS	Lung brain	V86A	None	Virulence	Zhou et al. (2010)
	Spleen	L90La M93T W102R
	Lung	D92G	D92E	Virulence and resistance to host’s immune response	Seo et al., 2002, Zhou et al., 2010
	Brain	M98L K108E	None	Virulence	Zhou et al. (2010)
	Inoculum	S99L

Synonymous mutations.