Evolutionary interactions between haemagglutinin and neuraminidase in avian influenza
© Ward et al.; licensee BioMed Central Ltd. 2013
Received: 28 June 2013
Accepted: 16 September 2013
Published: 9 October 2013
Skip to main content
© Ward et al.; licensee BioMed Central Ltd. 2013
Received: 28 June 2013
Accepted: 16 September 2013
Published: 9 October 2013
Reassortment between the RNA segments encoding haemagglutinin (HA) and neuraminidase (NA), the major antigenic influenza proteins, produces viruses with novel HA and NA subtype combinations and has preceded the emergence of pandemic strains. It has been suggested that productive viral infection requires a balance in the level of functional activity of HA and NA, arising from their closely interacting roles in the viral life cycle, and that this functional balance could be mediated by genetic changes in the HA and NA. Here, we investigate how the selective pressure varies for H7 avian influenza HA on different NA subtype backgrounds.
By extending Bayesian stochastic mutational mapping methods to calculate the ratio of the rate of non-synonymous change to the rate of synonymous change (d N /d S ), we found the average d N /d S across the avian influenza H7 HA1 region to be significantly greater on an N2 NA subtype background than on an N1, N3 or N7 background. Observed differences in evolutionary rates of H7 HA on different NA subtype backgrounds could not be attributed to underlying differences between avian host species or virus pathogenicity. Examination of d N /d S values for each subtype on a site-by-site basis indicated that the elevated d N /d S on the N2 NA background was a result of increased selection, rather than a relaxation of selective constraint.
Our results are consistent with the hypothesis that reassortment exposes influenza HA to significant changes in selective pressure through genetic interactions with NA. Such epistatic effects might be explicitly accounted for in future models of influenza evolution.
The influenza A virus has its natural reservoir in wild waterfowl, who transmit it sporadically to other avian species along migratory flyways . The main antigenic influenza proteins - the surface proteins haemagglutinin (HA) and neuraminidase (NA) - are each encoded by a separate RNA segment and are classified into subtypes which do not cross-react serologically. Reassortment – the exchange of genetic segments between co-infecting parental viruses during replication – leads to novel combinations of HA and NA subtypes. There are currently 16 known HA subtypes (H1-H16) and 9 known subtypes of NA (N1-N9) circulating in birds . Whilst all of subtypes H1-H16 and N1-N9 can be found amongst wild waterfowl , viruses with certain HA/NA combinations occur frequently in nature whereas others are rarely observed [4–6]. This, combined with the failure of laboratory studies to produce viable reassortant viruses of particular subtype combinations, has led to the suggestion that there is a requirement for a functional match between the influenza HA and NA .
The HA and NA proteins play complementary roles in the life cycle of the influenza virus. Both HA and NA bind to host cell receptors containing sialic acid residues: HA to initiate viral entry into the host cell, and NA to permit the release of viral progeny from infected cells. Experimental studies have suggested that a fine balance between HA and NA activity must be achieved for productive viral infection . Such a balance may, in fact, be more important for viral fitness than high levels of activity per se. For example,  showed that when artificially generated reassortant viruses of the N1 NA subtype were cultured, several (e.g. H3N1) only gave low yields. However, when the low-yield H3N1 culture was passaged, a number of changes occurred in the HA which reduced its receptor binding affinity, apparently to match that of the NA in the reassortant rather than to return to the high levels of HA activity found in the H3N8 parent virus.
Both the HA and NA proteins are thought to determine sensitivity of naturally-occurring influenza viruses to neuraminidase-inhibiting drugs (NAIs) . In vitro studies have investigated genetic interactions between HA and NA in terms of NAI resistance. Evidence suggests that mutations in the HA which decrease receptor binding activity may compensate for a decrease in NA activity resulting from treatment with NAIs, thus restoring the balance between HA and NA function [7, 11–13]. In addition, HA and NA mutations which individually confer low-level resistance to NAIs have been found to combine synergistically to confer resistance at a higher level . Interdependence between the length of the NA stalk section and the number of HA glycosylation sites has been identified in laboratory strains [8, 15] and may also have direct consequences for the transmission of influenza viruses to other host species. For example, influenza A viruses which have become established in terrestrial poultry may possess additional HA glycosylation sites, accompanied by deletions in the stalk section of their NA [16, 17].
Reassortment has been implicated in the emergence of pandemic influenza viruses, including those of avian origin which were responsible for significant human mortality in the twentieth century [18, 19] and the 2009 H1N1 pandemic strain . Naturally-occurring reassortment events could affect the functional balance between the HA and NA proteins  and this could in turn affect their evolution. Whilst previous studies have investigated evolutionary rates of influenza (e.g. [21, 22]), few have focused on how rates of evolution are affected by genetic interactions between segments .
Evolution of protein coding sequences can be quantified in terms of rates of synonymous (d S ) and non-synonymous substitution (d N ) and their ratio, d N /d S , following the counting-based methods of  and . Departures from selective neutrality can be detected by a d N/d S ratio which differs from 1. Positive selection is inferred when d N/d S > 1. When d N/d S < 1, it is inferred that purifying selection is acting. However, gene-wide estimates of d N/d S which show overall purifying selection may mask a small number of sites experiencing positive selection. For example, while the overall rate of non-synonymous substitution across the influenza HA has been found to be lower than the synonymous substitution rate in birds and humans (e.g. [22, 26]), evidence has been provided for positive selection at certain amino acid sites, particularly those of antigenic significance (e.g. [27–30]).
Avian influenza viruses of the H7 HA subtype present an epidemiological and economic threat on a global scale. Along with H5, H7 is the only subtype associated with the highly pathogenic form of avian influenza and has been known to cause outbreaks in domestic poultry (e.g. [17, 31–33]), human infection [34–36] and even human mortality . The danger posed by H7 viruses is exemplified by recent human infections with H7N9 avian influenza, which had claimed at least 37 lives in China as of 28 May, 2013, and has been associated with an estimated 36% fatality rate amongst cases admitted to hospital . In particular, reassortment events between H7, N9 and H9N2 viruses have been suggested to have been important in the emergence of the outbreak-causing H7N9 lineage .
In this study, we adopted a Bayesian stochastic mutational mapping approach [39, 40] to investigate how the association with different NA subtypes influences the evolution of the HA-encoding segment of avian influenza. Specifically, d N /d S ratios of avian influenza H7 HA1 were evaluated for clades associated with different NA subtype backgrounds. We extended the mutational mapping approach of Nielsen [39, 40] by rescaling the inferred numbers of synonymous and non-synonymous changes to calculate d N /d S . Ancestral trait mapping was used to construct a clade-model that inferred background NA subtypes for branches across the tree, and d N /d S was averaged across all parts of the tree corresponding to a particular subtype. The ancestral trait mapping accounts for a lack of monophyly across the tree with respect to NA subtype background, which arises through repeated exposure of H7 HA to different NA backgrounds via reassortment. We find substantial differences between gene-wide d N /d S for avian influenza H7 HA on different NA subtype backgrounds, consistent with the hypothesis that the selective pressure experienced by HA can be affected by its genetic context.
Composition of avian H7 HA sequence dataset (background NA subtypes N1, N2, N3 and N7)
All subtypes (253)
N. America (99)
N. America (2)
N. America (66)
N. America (29)
N. America (2)
S. America (8)
S. America (0)
S. America (0)
S. America (8)
S. America (0)
For each background NA subtype, the H7 HA sequences covered a time-span of at least 25 years. There were roughly equal numbers of sequences from Eurasia and America (132 and 107 respectively), and sequences from Europe, Asia and North America were present for all four subtypes considered. The geographic spread of H7 avian influenza viruses of different background NA subtypes appeared to differ between continents. For example, 85% of the H7N1 sequences and 74% of the H7N7 sequences were from Europe, whilst 88% of the H7N2 isolates were from North America. H7N3 appeared to be the most ubiquitously sampled subtype, in terms of location, host order and pathogenicity. Overall, geographic and temporal diversity appeared to be captured in all subtypes.
Within the Eurasian clade, the Australian isolates formed a clade with 100% bootstrap support. The maintenance of a distinct Australasian lineage of H7 avian influenza within the Eurasian clade, with continued reassortment of different NA subtypes onto the H7 HA, has recently been reported . The phylogenetic position of early European fowl plague viruses (FPV) as a sister lineage to the Australian clade has been observed in other studies [43, 44, 48] and was observed in our ML and MrBayes phylogenies, although both methods appeared to have difficulty in placing this clade (which could account for the relatively low posterior probability observed for the Eurasian clade in the MrBayes consensus tree). Following other evolutionary studies , we excluded the FPV sequences from our mutational mapping analysis of evolutionary rates, since they have been highly cultured and may show artificially high rates of molecular change.
On a smaller geographic scale, H7 HA sequences from within avian influenza outbreaks, such as the Italian H7N1 outbreak of 1999–2000, clustered together. The observation that H7 HA sequences from viruses with different NA subtype backgrounds were distributed across the tree, rather than forming distinct clades, is indicative of repeated reassortment between H7 HA and NA of different subtypes. Avian H7 HA sequences did not cluster into distinct lineages corresponding to HP or LP viruses, or viruses from avian hosts of orders Anseriformes or Galliformes.
For all 1000 MrBayes phylogeny samples, the average d N estimate across all HA1 sites for a given NA background was plotted against the d S value for that tree sample (Additional file 1: Figure S2). This indicated that a phylogeny sample with a higher rate of synonymous substitution would also have a higher rate of non-synonymous substitution, although the rate of synonymous substitution was not an exact predictor of the corresponding non-synonymous substitution rate. It may be observed that, whilst the same d S value would lead to a similar expected d N for background NA subtypes N1, N3 and N7, there was little overlap between the d N values on the N2 background and on backgrounds N1, N2 and N3, with the d N values for N2 appearing to be higher than for the other NA background subtypes.
Average d N / d S across the H7 avian influenza HA1 on different NA backgrounds
Mean d N /d S
Lower 90% HPD limit for d N /d S
Upper 90% HPD limit for d N /d S
Comparing evolutionary rates for H7 avian influenza HA1 on different NA subtype backgrounds
d N /d S
For all NA subtype comparisons, the distributions of paired differences for d S were roughly centred on zero (i.e. approximately 50% of the paired differences were greater than zero, and 50% less than zero), indicating no substantial differences between the distributions, as suggested by the HPD interval plot. However, the pairwise difference comparisons indicated an elevated rate of non-synonymous change in H7N2, leading to a substantially higher d N /d S for H7N2 than for the other subtypes (split of p = 0.979/0.021 against H7N1; p = 0.991/0.009 against H7N3; p = 0.962/0.038 against H7N7).
Our results for the ordering of d N /d S values across H7 HA1 on different NA subtype backgrounds are consistent with the point estimates obtained by a previous study  which was based upon the single likelihood ancestor counting (SLAC) method . The results from  could not be statistically compared between subtypes and did not account for uncertainty in the phylogenetic or mutational history. Furthermore, estimating d N /d S separately for H7 HA datasets corresponding to different background NA subtypes, as was carried out in , implicitly assumes that the tree of all H7 HA sequences should split into distinct clades according to background NA subtype. Our phylogenetic analysis, along with previous studies (e.g. ), has shown that H7 HA sequences are not monophyletic with respect to viral NA subtype. It is therefore possible that error might be introduced into d N/d S estimates from datasets corresponding to individual NA subtype backgrounds, by incorrectly assuming that ancestral lineages were associated with a particular NA subtype.
Average d N / d S across H7 avian influenza HA1 for lineages corresponding to different viral pathogenicities
Mean d N /d S
Lower 90% HPD limit for d N /d S
Upper 90% HPD limit for d N /d S
Comparing H7 avian influenza HA1 evolutionary rates along lineages classified by viral pathogenicity
d N /d S
Average d N / d S across H7 avian influenza HA1 for lineages corresponding to different avian host orders
Avian host order
Mean d N /d S
Lower 90% HPD limit for d N /d S
Upper 90% HPD limit for d N /d S
Comparing H7 avian influenza HA1 evolutionary rates along lineages classified by avian host order
d N /d S
Ans. - Gall.
Ans. - other
Gall. - other
The largest number of sites under putative positive selection was observed on the N2 NA background (23 sites under putative positive selection, out of the 329 sites considered). This was approximately twice the number of sites with a mean d N /d S > 1 on N1, N3 or N7 backgrounds (13, 9 and 8 sites respectively). When the largest 50 mean d N /d S values across the HA1 codon sites were ordered by magnitude for each NA background subtype, the d N /d S value on the N2 background was higher than the d N /d S value of that rank on all other NA subtype backgrounds (Additional file 1: Figure S4a). The large d N /d S values observed at individual codon sites for H7 HA1 on the N2 NA background would have led to the elevated HA1-wide d N /d S observed on the N2 NA background; however, H7N2 also had many of the smallest d N /d S values out of the different subtypes at individual amino acid sites (Figure 6, Additional file 1: Figure S4b and Figure S5). For all NA subtype backgrounds, sites with mean d N /d S > 1 were observed in each of the fusion, vestigial esterase and receptor binding domains.
Although high d N /d S values were observed at two sites in the signal peptide region of H7 HA on NA backgrounds N2, N3 and N7, no sites with mean d N /d S > 1 were observed for the H7 HA signal peptide region on the N1 NA background. The signal peptide region appears to have been considered in previous gene-wide or HA1-wide calculations of d N /d S (e.g. [22, 28]), and the values we have reported across the alignment encompass the signal peptide and HA1. Note that we observed the same general pattern of average d N /d S across sites for H7 avian influenza on different NA backgrounds (i.e. a higher d N /d S when H7 HA was on an N2 NA background than on an N1, N3 or N7 NA background) when averaging across just the HA1 coding region, i.e. excluding the signal peptide region (data not shown).
Some commonality was observed between the H7 HA1 sites with mean d N /d S > 1 on different NA subtype backgrounds. One site (site 218 in H3 numbering) had mean d N /d S > 1 on all four NA subtype backgrounds; 3 amino acid sites had mean d N /d S > 1 on 3 out of the four NA subtype backgrounds and 10 sites had mean d N /d S > 1 in two out of the four background NA subtypes (Additional file 1: Table S1). Site 218 has been linked with receptor-binding specificity [52–54] and thus high levels of non-synonymous change at this site could signify a move towards viruses which are capable of infecting other host species.
Of the 75 H7N2 HA1 sequences studied, 66 were from viruses circulating in the North American live bird markets between 1994 and 2006, or from the many avian influenza outbreaks they seeded in commercial poultry in the Northeast United States during this period [41, 55]. It may also be noted that 88% of the North American H7N2 sequences possessed a deletion of 8 amino acids at the HA receptor binding site, and a recent study has put forward the idea that non-synonymous changes might have occurred in the HA to maintain functionality . This would be compatible with our observation that a large number of sites with mean d N /d S > 1 were found in the receptor binding domain for H7 HA on the N2 NA background (Figure 5 and Additional file 1: Table S1). If the elevated level of non-synonymous change only applied to H7N2 HA1 lineages associated with the receptor binding site deletion then our results could also be compatible with this hypothesis. It is possible that molecular changes at, or adjacent to, other sites in the receptor binding region (for example, the elevated d N /d S that we observed in H7N2 at sites 216 and 218 – H3 numbering) could be compensating for the HA deletion. Although this could indicate co-evolution at sites within the HA, again this could be to restore HA activity levels to match those of the NA.
H7N2 was the most common avian influenza subtype isolated from the North American live bird markets between 1994 and 2006 [57, 58], garnering attention as a potential source for a human pandemic virus [35, 59] after it proved capable of causing limited human infection [60, 61]. North American H7N2 viruses isolated between 2002 and 2003 were found to exhibit increased affinity towards human-like α-2,6-linked sialic acid receptors  which has also been associated with adaptation to certain terrestrial birds, such as chickens and quails [63–65]. While (like other known H7N2 avian influenza lineages) North American H7N2 only presented in a low pathogenic form, molecular evidence suggested a step-wise accumulation of basic amino acids at the North American H7N2 HA cleavage site towards those observed in highly pathogenic viruses . An elevated level of non-synonymous change amongst circulating avian influenza viruses could signify a heightened risk of molecular changes occurring which would increase the pathogenicity of the virus, or its ability to infect new species and become transmissible amongst humans. Although H7N2 avian influenza appeared to have been eradicated from domestic poultry in North America by mid-2006 , such findings might be particularly pertinent if the strain re-emerges.
Our stochastic mutational mapping method for calculating the d N /d S ratio provides many advantages for investigating selective pressure in influenza HA on different NA subtype backgrounds in the presence of reassortment. By using the rescalings described in Methods, we are able to estimate rates of synonymous substitution (d S ) and non-synonymous substitution (d N ), rather than merely counting the number of synonymous or non-synonymous changes along branches [39, 49]. Also, estimating d N and d S separately allowed us to attribute differences in the d N /d S ratio to underlying differences in the non-synonymous or synonymous rate. Our method also enabled us to estimate d N and d S along parts of the HA tree corresponding to different NA subtype backgrounds, despite sequences from viruses with different NA subtypes being distributed across the tree; this does not require the introduction of additional model parameters, but merely summarizes the relevant lineages. Finally, our rescalings allowed d N and d S to be compared between clades of different sizes and divergence.
Bayesian methods for phylogenetic inference and mutational mapping provide an advantage over parsimony and maximum-likelihood methods since they naturally accommodate uncertainty in the phylogenetic reconstruction (by considering multiple tree and model samples) and the mutational history (by sampling multiple histories for each site in each phylogeny sample). Failing to account for phylogenetic uncertainty can lead to artificially narrow confidence intervals for estimating substitution rates . We note that, whilst the topologies and relative branch lengths are consistent between our maximum likelihood and Bayesian phylogenies, the MrBayes trees had longer branch lengths. This is likely to be due to a known artefact of MrBayes ; however, our d N /d S estimates for H7 HA are consistent with those from a previous study  which used different phylogenetic inference methods.
Another advantage over parsimony is that non-parsimonious maps are not automatically excluded. Using parsimony to minimise the number of mutations required to produce the observed pattern in the data can lead to an underestimate in substitution rates, perhaps by a factor of over 20%, and can also bias d N /d S estimates by underestimating the number of synonymous changes in scenarios where synonymous mutations occur more frequently than non-synonymous mutations .
In addition to the ability to use a collection of trees and sample multiple mutational histories, our mutational mapping method also possessed advantages over the PAML maximum likelihood software [68, 69]. Although PAML can be used to estimate d N /d S along the branches of a phylogeny [70, 71], its use in our study would have led to an over-parameterised model with very little power for statistical testing using likelihood ratio tests, since parameters would be estimated for each branch in the tree. Furthermore, with stochastic mutational mapping we did not have to pre-specify branches with potentially positively-selected sites, which is a requirement of the branch-site models in PAML. In addition, PAML assigns d N /d S values for branches to a pre-determined number of rate classes (bins), which would lead to a loss of precision compared to the stochastic mutational mapping approach. Mutational mapping also records the timings of mutations across the tree, which we have used in calculating evolutionary rates, whereas existing maximum likelihood methods do not.
Assuming that all synonymous changes are essentially neutral, d S is independent of the effective size (N e ) of the population and is simply the mutation rate , although synonymous rates in RNA viruses can be affected by the virus’ secondary structure . Our finding that d S for H7 influenza HA1 did not vary across different NA subtype backgrounds therefore suggested that the mutation rate was constant for H7 HA1 across NA subtype backgrounds.
Under non-neutral models of evolution, differences in selective pressure could lead to differences between substitution rates . Since non-synonymous changes in the HA1 coding region are likely to be non-neutral, the elevated d N observed for avian influenza H7 HA1 on an N2 NA subtype background might be explained by a number of scenarios. Firstly, selection could be acting to fine-tune the functional HA-NA balance of H7 HA on an N2 NA background following reassortment. Secondly, a burst of positive selection could have occurred in the H7N2 lineages, which is not a consequence of the N2 NA background, but instead a consequence of an unrelated, co-varying factor such as avian host, demographic scenario, or an interaction with another gene segment. Thirdly, a relaxation of selective constraint could have taken place when H7 HA was exposed to the N2 NA background. The results of this study do not definitively distinguish between such scenarios and causality cannot be inferred. However, whilst d N /d S > 1 was observed in a larger number of HA1 sites on the N2 NA background than on N1, N3 or N7 backgrounds, at many sites the N2 viruses also had the lowest d N /d S values out of all NA subtype backgrounds (Figure 6 and Additional file 1: Figure S4b) and this is not indicative of an overall relaxation of selective constraint. One explanation for the observed pattern of site-by-site d N /d S values could be a larger effective population size in HA for the H7N2 viruses, which would allow selection to act more effectively in removing deleterious mutations, leading to a reduction of variation at some sites.
The results presented in this study are consistent with the hypothesis that reassortment exposes HA to significant changes in selective forces via association with different NA subtypes. However, establishing a causal relationship between background NA subtype and differences in evolutionary rates of HA is not straightforward. Mutational mapping analyses excluded underlying differences in evolutionary rates between viruses of different pathogenicity, or between different avian host orders, as causative factors in the elevated d N /d S observed in H7N2 avian influenza HA1. Nonetheless, other differences between the environments from which sequences were isolated may have influenced the selective pressure experienced. For example, it has been suggested that long term evolution in commercial poultry, which are not the natural reservoir of avian influenza, could lead to accelerated rates of evolution and the accumulation of point mutations in viruses in the live bird markets [74, 75].
Although we cannot exclude prolonged circulation of avian influenza viruses in non-natural avian hosts as a factor in observing an elevated d N /d S for H7 HA on an N2 NA background, it can be noted that 66% of the H7N1 sequences we analysed were sampled during an outbreak of LP and HP H7N1 avian influenza in domestic poultry in Italy, and that the elevated d N /d S did not appear to extend to this subtype background. However, Italian H7N1 sequences were sampled over a period of less than two years, compared to over 12 years for H7N2 in the North American live bird markets. The effect of continuous circulation amongst non-natural avian hosts on selective pressure could be investigated in H5N1 avian influenza, which is endemic in the live bird markets of East Asia . Given detailed information about the origin of the avian hosts from which viruses were collected, d N /d S could also be compared along lineages corresponding to wild or domestic avian hosts.
Future studies could investigate rate variation along individual branches of the H7 HA1 phylogeny to determine whether the elevated d N /d S extends to all lineages on the N2 NA subtype background (for example in both Eurasia and North America), or whether it is localised to particular parts of the tree (for example, to a particular geographical location such as the North American live bird markets, or specifically after transmission to a new avian species e.g. ). Further analyses could also consider whether the elevated d N /d S observed for H7N2 HA1 also extends to other segments, for example whether the NA for these viruses showed higher levels of non-synonymous change than the NA sequences for the H7N1, H7N3 or H7N7 viruses. Other investigations could consider interactions with other influenza proteins, such as the matrix protein, with which the HA and NA both interact closely. The precise nature of the genetic changes which take place when HA is placed in a novel NA background (or vice versa) could also be explored in the laboratory using reverse genetics experiments, to provide an insight into how the balance between HA and NA activity is regulated.
Future influenza modelling studies could explicitly incorporate genetic interactions between segments, rather than assuming that their evolution is independent. Such effects might be included in extensions to frameworks such as that of Zhang et al. , who model the impact of reassortment on the dynamics of novel human influenza strains. Although much modelling work has focused on human influenza rather than avian influenza, a recent study suggested that evolutionary changes mediating the HA-NA functional balance were an important determinant of the transmissibility of the 2009 H1N1 pandemic influenza strain , thus our result might find application in models of the emergence and spread of zoonotic influenza strains in human populations.
Reassortment of avian influenza segments creates novel combinations of influenza genes and repeatedly exposes segments to different genetic backgrounds. Our study has shown that the selective pressure experienced by the influenza HA can vary depending upon the genetic context in which a segment finds itself. In this case, the average d N /d S across avian influenza HA1 of subtype H7 differed according to the background NA subtype of the virus. Observed differences in selective pressure could not be accounted for by differences in the pathogenicity of the virus, or the taxonomic order of the avian host from which it was sampled. We believe that future influenza modelling studies could incorporate epistatic interactions between gene segments, for example when considering the impact of reassortment on the emergence dynamics of novel strains.
All available complete H7 avian influenza nucleotide sequences for the HA protein-coding region were downloaded from the NCBI database ( http://www.ncbi.nlm.nih.gov)  and labelled according to the corresponding NA subtype of the virus. Sequences were screened for identity and, in the case of identical sequences, only one such isolate was included. Only NA subtypes for which there were more than 20 sequences were analysed – these subtypes were N1 (62 sequences), N2 (75 sequences), N3 (69 sequences) and N7 (47 sequences) (Table 1). Sequences were also labelled according to the taxonomic order of the avian host from which the virus was isolated (Additional file 1: Table S2). Where possible, classification of the sequences into highly pathogenic (HP) or low pathogenic (LP) was made by searching the literature for studies confirming the pathogenic status of the strain using laboratory testing. Where no record of the pathogenicity of an isolate could be found, sequences were classified as HP if they possessed a motif at the HA1/HA2 cleavage region which was the same as that of a previously confirmed HP strain, in accordance with . Sequences with a novel cleavage site motif which had not been previously documented as either HP or LP were not labelled by pathogenicity.
Sequence alignment was performed manually, using BioEdit . The alignment of H7 HA sequences was split at the HA1/HA2 cleavage site  and just the HA1 coding region, which encompasses approximately two thirds of the length of the whole HA and has the major antigenic role for the virus , and the signal peptide region (17 amino acids immediately preceding the start of the HA1), were analysed in this study. A single breakpoint analysis  in the HyPhy software [86, 87] found no evidence of recombination in the alignment. Investigations using the method of Xia et al. (2003)  and plots of transitions and transversions against genetic distance in the DAMBE software  found no evidence of saturation at codon positions 1 and 2; whilst there was some evidence of saturation at the third codon position, this was not severe (Additional file 1: Figure S6).
A bootstrapped phylogenetic tree (with 1000 bootstrap replicates) was constructed for the avian influenza H7 HA1 coding region using maximum likelihood inference in the PhyML software . A GTR + Γ model of nucleotide substitution  was used, which allowed for gamma-distributed rate variation across sites. MrBayes version 3.1.2 [92, 93] was used to obtain posterior samples of topologies, branch lengths and substitution model parameters for the H7 HA1 alignment. A GTR + Γ model of nucleotide substitution was again selected. An outgroup sequence, A/Australian_shelduck/Western Australia/1756/1983(H15N2) [GenBank accession number: ABB90704], was used to root the trees. H15 been shown to be the closest HA subtype phylogenetically to H7 .
Three independent MrBayes runs were conducted, each with Markov Chain Monte Carlo (MCMC) searching over 2,000,000 generations. Trees and parameters were sampled every 1000 generations. The Tracer software  was used to inspect the chain traces, which indicated that a burnin period of 1,000,000 generations was sufficient to exclude samples taken before the chains had converged. Chain traces were compared across the three runs, with similar post-burnin values in all runs. A post-burnin sample of 1000 posterior trees and sets of parameter estimates was used for the analysis of selection.
Stochastic mutational mapping [39, 40, 95] was used to infer mutational histories (maps) using posterior phylogeny samples taken from MrBayes runs. Mutational histories describe the nature and location of molecular changes along the branches of a phylogeny (Additional file 1: Figure S7). Stochastic mutational mapping is a Bayesian approach in which mutational histories are sampled from the posterior distribution of mappings, given the observed nucleotide data.
Thus, mappings are sampled in proportion to their posterior probability. For a more detailed description see .
For each of the 1000 post-burnin MrBayes phylogeny and substitution model samples, 10 mutational mappings were simulated from the posterior distribution for each nucleotide site in the H7 HA1 alignment. Within each phylogeny sample and mutational mapping replicate, the mutational history of each amino acid site in the alignment was reconstructed by combining the mutational maps for the first, second and third codon positions. Branch lengths from the maps for codon positions 1 and 2 were rescaled to the branch lengths of position 3. This allowed us to identify codon substitutions and count the number of synonymous and non-synonymous changes (C s and C n respectively) along different parts of the tree, as well as to record their timing along the branches (Additional file 1: Figure S8).
c = number of codon intervals (distinct codon states) along a part of the tree. A new interval occurs every time there is a nucleotide change, even if it is silent, since this alters the codon state
j = position of nucleotide site in the codon (1, 2 or 3)
s ij = proportion of changes at the j th codon position of the codon at interval i which are synonymous
n ij = proportion of changes at the j th codon position of the codon at interval i which are non-synonymous
v ij = "mutational time interval" or "dwell time". This is obtained by multiplying the substitution rate r j with the length along the branch spent in each codon state. The parameter r j is drawn from a gamma distribution, whose parameters were sampled during the MrBayes analysis. A value of r j is sampled for each codon position (j = 1, 2, or 3) at the root from its respective posterior distribution and the stochastic mutational map is then sampled under this rate
V T = sum across all codon positions and over all codon intervals of the v ijs, i.e.
Together with the v ijs, this gives a time-weighted average which assigns more weight to codons with longer dwell times.
Note that, for a single codon interval, if the dwell time information is not used then our calculation of the number of synonymous and non-synonymous sites is the same as that of Nei and Gojobori , since our s ij is equivalent to their f i . However, unlike the Nei and Gojobori approach, by using the dwell time weighting we accommodate variation in branch lengths which may affect the counting procedure. Note also that Nei and Gojobori used the evolutionary distance formula of Jukes and Cantor (1969)  to estimate the expected number of synonymous changes per synonymous site (or non-synonymous changes per non-synonymous site) from the proportions of synonymous and non-synonymous differences between pairs of sequences. However, our method samples the full nucleotide state history across the phylogeny for each nucleotide in the alignment, thus d N and d S may be estimated directly by counting synonymous and non-synonymous changes along branches and rescaling by numbers of synonymous and non-synonymous sites, and dwell times, as described above. In addition, we account for uncertainty in the tree and model parameters by performing our analysis across 1000 MrBayes samples.
Values of C s , C n , S s and S n were used in calculating synonymous and non-synonymous evolutionary rates (d S and d N respectively) along different parts of the phylogeny, corresponding to background NA subtypes N1, N2, N3 and N7. In order to calculate d N and d S for H7 HA1 on different NA subtype backgrounds, parsimony mapping was used to assign ancestral NA subtypes at internal nodes along the MrBayes phylogeny samples, based on assignments at the tips of the phylogeny (i.e., the NA subtypes corresponding to the H7 HA sequences in our dataset). This allowed branches to be classified by NA subtype: N1, N2, N3 or N7 (Additional file 1: Figure S9). Branches where a subtype could not be unambiguously assigned from a single pass of the parsimony algorithm from the tips of the tree to the root were not used in the analysis. The use of parsimony avoids the possible confounding factor of incorrect lineage classification which could arise from methods which force ancestral states to be inferred for every branch, although the exclusion of ambiguous lineages potentially results in a loss of information. S s and S n were calculated as described above across all branches to which a particular NA subtype had been assigned, and numbers of synonymous and non-synonymous changes were counted along those parts of the tree.
Here, T is obtained by summing the branch lengths at all nucleotide positions in the amino acid site, with branch lengths for the first and second codon positions rescaled to the third codon position lengths (i.e. 3* sum of the third position branch lengths), for all branches in the phylogeny to which a particular NA subtype has been assigned. Rescaling by the length of the portion of the tree corresponding to each background NA subtype allowed for a comparison of evolutionary rates between clades of different sizes. This differs from the previous mutational mapping approaches of Nielsen and others [39, 40, 95], including those implemented in the SIMMAP software . By performing these calculations upon each of the 1000 MrBayes posterior phylogeny samples, we obtained approximations to the posterior distributions for d N and d S for each background NA subtype, at each codon site in the H7 HA1 alignment.
Estimates of d N and d S , obtained at each codon site for each background NA subtype (see Additional file 1: Table S3 for a list of sequences used in the mutational mapping analysis), were averaged over the 10 mutational mapping replicates for each phylogeny sample. Average values of d N across the sites in the HA1 alignment were obtained for each NA subtype by calculating the mean of the d N values across all codon sites in the alignment (and similarly for d S ). For all 1000 MrBayes phylogeny samples, we divided the HA1-wide d N estimate for a given NA subtype by the corresponding HA1-wide d S value for that subtype to obtain an approximation to the posterior distribution for the HA1-wide d N /d S for that subtype.
Estimates of d N /d S at individual codon sites in the H7 HA1 alignment were also calculated for each NA background subtype. For each site, d N and d S values were averaged over the 10 mutational mapping replicates for each tree, and then averaged over the 1000 MrBayes tree samples. To calculate the d N /d S ratio on a site-by-site basis, d N for each site was divided by the average d S value across the genome for that subtype. The gene-wide d S was used to avoid inflation of d N /d S values as a result of unobserved synonymous change at individual sites, and ensured that we were conservative in identifying sites under putative positive selection. Sites with a mean value of d N /(gene-wide d S ) greater than one were identified as being under putative positive selection. Sites in the H7 HA alignment were converted to H3 numbering prior to being reported, as is the convention for influenza, and numbering was based upon the alignment of Nobusawa et al.  (sites numbered -17 to -1 for the signal peptide region and 1 to 329 for HA1). The HA1 domain in which putatively positively selected sites were found was reported, using the alignment of Yang et al.  in which portions of the influenza HA corresponding to the fusion domain, vestigial esterase domain and receptor binding domain were identified.
Posterior distributions of d N /d S and rates of synonymous and non-synonymous substitution for avian H7 HA on different background NA subtypes could be visualised by plotting highest posterior density (HPD) intervals. A 100*(1-α)% credible interval for a posterior distribution for a parameter θ is any interval [a, b] in the domain of the distribution such that the posterior probability of θ lying between a and b is 1 – α. The highest posterior density (HPD) interval is the narrowest such credible interval. After checking the distributions for unimodality, 90% HPD intervals were calculated using the Chen and Shao algorithm  in the boa R package for the analysis of Bayesian output  and plotted using a custom R script (available on request). The overlap of the HPD intervals can be used as an indicator of whether the means of the distributions are significantly different.
In order to assess the overlap between posterior distributions of evolutionary rates for different background NA subtypes, the following comparison was implemented using 'distributions of differences’. For rate distributions corresponding to arbitrary NA background subtypes A and B, a comparison method was implemented as follows. Multiple pairings of evolutionary rate estimates were drawn randomly from across the 1000 posterior samples, with one observation from subtype A and one from subtype B in each pair. The proportion of pairings for which the observed rate from subtype A was greater than the observed rate from B (and vice versa) was recorded. For a null hypothesis that there is no difference between the distributions, the point of interest is where zero lies in the distribution of paired differences. If the distributions for A and B were identical then the corresponding distribution of paired differences should be centred on zero, as one would expect A > B for half of the paired samples and A<B for the other half. However, if the proportion of samples for which A>B is extremely skewed (e.g. less than 0.05 or greater than 0.95) then zero lies in the tail of the distribution of paired differences, providing evidence that the location of the distributions is different (Additional file 1: Figure S10). A total of 106 random pairings were sampled for each comparison of evolutionary rate distributions; this gave similar values to systematically comparing each of the 1000 observations for one subtype with each of the 1000 observations for the other subtype. Here we report the values from the randomized pairing approach.
In this study, avian H7 HA sequences were labelled according to the NA subtype of the virus and rates of evolution were calculated for lineages corresponding to different NA subtypes. In order to test whether a non-uniform distribution of host species or pathogenic viruses across different NA backgrounds could be confounding the ability to infer differences in d N /d S between subtypes, we performed two further analyses in an analogous manner to the NA subtype analysis. These analyses involved labelling sequences and performing stochastic mutational mapping to calculate and compare d N /d S between (a) HP and LP viruses and (b) viruses from different avian host orders. Bird orders compared were Galliformes (turkeys, chickens etc.) and Anseriformes (ducks, geese, etc.) (Additional file 1: Table S2), with all other avian host orders combined (classified as "other") due to a paucity of sequences. To further investigate the potential effect of uneven sampling of NA subtype backgrounds with respect to avian hosts, we also performed a Spearman’s rank correlation test between the proportion of sequences from terrestrial poultry and our mean d N /d S estimate for each background NA subtype.
A list of GenBank accession numbers is provided (Additional file 1: Table S3) for the sequence dataset analysed in this study.
Haemagglutinin subunit 1
We would like to thank Andrea Betancourt for discussion of re-scaling evolutionary rates and Anne Kupczok for helpful comments on a draft version of this manuscript.
This work was supported by the Biotechnology and Biological Sciences Research Council, the Government of the Republic of Panama, the Interdisciplinary Centre for Human and Avian Influenza Research ( http://www.ichair-flu.org) funded by the Scottish Funding Council, and the Institute for Science and Technology Austria.
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License ( http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.