Evolution of the M gene of the influenza A virus in different host species: large-scale sequence analysis

Background Influenza A virus infects not only humans, but also other species including avian and swine. If a novel influenza A subtype acquires the ability to spread between humans efficiently, it could cause the next pandemic. Therefore it is necessary to understand the evolutionary processes of influenza A viruses in various hosts in order to gain better knowledge about the emergence of pandemic virus. The virus has segmented RNA genome and 7th segment, M gene, encodes 2 proteins. M1 is a matrix protein and M2 is a membrane protein. The M gene may be involved in determining host tropism. Besides, novel vaccines targeting M1 or M2 protein to confer cross subtype protection have been under development. We conducted the present study to investigate the evolution of the M gene by analyzing its sequence in different species. Results Phylogenetic tree revealed host-specific lineages and evolution rates were different among species. Selective pressure on M2 was stronger than that on M1. Selective pressure on M1 for human influenza was stronger than that for avian influenza, as well as M2. Site-by-site analyses identified one site (amino acid position 219) in M1 as positively selected in human. Positions 115 and 121 in M1, at which consensus amino acids were different between human and avian, were under negative selection in both hosts. As to M2, 10 sites were under positive selection in human. Seven sites locate in extracellular domain. That might be due to host's immune pressure. One site (position 27) positively selected in transmembrane domain is known to be associated with drug resistance. And, two sites (positions 57 and 89) locate in cytoplasmic domain. The sites are involved in several functions. Conclusion The M gene of influenza A virus has evolved independently, under different selective pressure on M1 and M2 among different hosts. We found potentially important sites that may be related to host tropism and immune responses. These sites may be important for evolutional process in different hosts and host adaptation.


Background
The influenza virus is a common cause of respiratory infection all over the world. The influenza A virus can infect not only humans but also avian, swine, and equine species. The virus has a negative single-stranded RNA with eight gene segments, namely PB2, PB1, PA, HA, NP, NA, M, and NS. The subtype of influenza A virus is determined by the antigenicity of two surface glycoproteins, hemaglu-tinin (HA) and neuraminidase (NA). The subtypes currently circulating in the human population are H1N1 and H3N2. Influenza A viruses cause epidemics and pandemics by antigenic drift and antigenic shift, respectively [1]. Antigenic drift is an accumulation of point mutations leading minor and gradual antigenic changes. Antigenic shift involves major antigenic changes by introduction of new HA and/or NA subtype into human population.
All known HA and NA subtypes are maintained in avian species, and all mammalian influenza A viruses are thought to be derived from the avian influenza A virus pool [1]. In avian species, influenza A viruses are in an evolutionary stasis [1]. In contrast, all gene segments of mammalian viruses continue to accumulate amino acid substitutions [1]. Today, the emergence of an influenza pandemic is of great global concern. If a novel influenza A subtype acquires the ability to spread between humans efficiently, it could cause the next pandemic [1]. This ability is acquired by reassortment between human and nonhuman influenza A viruses or by the accumulation of mutations in the non-human influenza virus. It is necessary to understand the evolutionary processes of influenza A viruses in various hosts so that we have better knowledge about the emergence of this pandemic virus. We conducted the present study to investigate the evolution of the M gene among different species. Although there are numerous studies on the evolution of the HA gene [2][3][4][5][6][7], only a few studies on the evolution of the M gene have been conducted [8].
The M gene is intriguing because it encodes both matrix and membrane proteins, and has multiple functions. The M gene (1027 bps) encodes two proteins, namely M1 (at nucleotide position 26 to 784) and M2 (at nucleotide position 26 to 51 and 740 to 1007) [9]. M1 is a matrix protein that lies just beneath the viral envelope in the form of dimers and interacts with viral ribonucleoprotein (vRNP) complex, forming a bridge between the inner core components and the membrane proteins [10][11][12][13]. vRNPs harbor the determinants for host range [1,14,15]. M1 contacts with both viral RNA and NP, promoting the formation of RNP complexes and causing the dissociation of RNP from the nuclear matrix [16][17][18][19][20][21]. M1 plays a vital role in assembly by recruiting the viral components to the site of assembly and essential role in the budding process including formation of viral particles [22,23]. M2 is a membrane protein which is inserted into the viral envelope and projects from the surface of the virus as tetramers [24,25]. The M2 protein comprises 97 amino acids -24 in the extracellular domain, 19 in the transmembrane domain, and 54 in the cytoplasmic domain. Extracellular domain of M2 is recognized by hosts' immune system [26][27][28]. Transmembrane domain of M2 has ion channel activity, which involved in uncoating process of the virus in cell [29]. Amantadine inhibits virus replication by blocking the acid-activated ion channel. The cytoplasmic domain of M2 interacts with M1 and is required for genome packaging and formation of virus particles [30][31][32][33][34][35][36].
The molecular mechanism of how the host range of influenza A viruses is determined is still not fully understood. The M gene may be involved in determining host tropism. Besides, novel vaccines targeting M1 or M2 proteins to confer cross-subtype protection have been shown to be promising [37][38][39][40][41][42][43]. Therefore, understanding of evolution of the M gene is of great importance and practical relevance.

Phylogenetic Tree
The phylogenetic trees for the M gene of all the sequence data we analyzed are shown in Figure 1. We defined "lineage" as an aggregate of large branches. The phylogenetic analysis revealed seven host-specific lineages: 1) human lineage (Hu1) consisting of H1N1 between 1918 and 1954 (Spanish Flu and its progeny viruses), H2N2 between 1957 and 1967 (Asian Flu and its progeny viruses), and H3N2 (Hong Kong Flu and its progeny viruses) after 1968; 2) another human lineage (Hu2) consisting of H1N1 (Russian Flu) after 1977; 3) avian lineage (Av1) including viruses mainly from Asia but also from other regions; 4) another avian lineage (Av 2) including viruses mostly from North America; 5) swine lineage (Sw1), located between human and avian lineages, mainly from North America; 6) another swine lineage (Sw2) diverging from Av1 and consisting of swine viruses after 1980, mainly from Europe; and 7) canine/equine lineage (CE) diverging from the root of Av2.
The M gene of all known human influenza A viruses, i.e., H1N1 between 1918 and 1957, H2N2 between 1957 and 1968, H3N2 after 1968, and H1N1 after 1977 was derived from that of the 1918 Spanish Flu. One lineage (Hu1) included three different subtypes (H1N1 between 1918 and 1957, H2N2 between 1957 and 1968, and H3N2 after 1968), which means that the same M gene was maintained in human influenza even after two antigenic shifts in 1957 and 1968. Another lineage (Hu2) included H1N1 after 1977. This M gene was also derived from Spanish Flu, but underwent different evolutionary processes and formed another lineage. Since H1N1 re-emerged in 1977 as Russian Flu, the two subtypes (H1N1 and H3N2) have been co-circulating in human populations and have formed two distinct lineages (Hu1 and Hu2). However, Hu2 exclusively includes H1N1 viruses and all human H3N2 are included in Hu1 ( Figure 1B). On the other hand, both avian influenza lineages (Av1 and Av2) did not show any subtype specificity, and included many dif-ferent subtypes ( Figure 1A and 1B). In avian lineages, even small branches of the phylogenetic tree are shared by different subtypes.
Although strains with the M gene in both avian lineages (Av1 and Av2) have been seen sporadically in humans, they have not been maintained in the population (blue characters in Av1 and Av2, Figure 1A and 1F). Strains with the M gene in swine lineages also infect humans, but these swine viruses have not been established in human populations (blue characters in Sw1 and Sw2, Figure 1A and 1F). All H5N1 viruses that infected humans as well as the H5N1 virus that infected swine possessed share the M gene of the avian influenza lineage (Av1, ( Figure 1E).

Evolutionary Rate
For evolutionary rate analysis, we included the sequences of only host-specific lineages and excluded other sequences such as those of the H5N1 influenza in humans ( Figure 1F. See "Materials and methods"). The profile of the sequences analyzed is shown in Table 1. Evolutionary rates were estimated for each lineage ( Figure 2).
Av2 of avian influenza A viruses showed the slowest evolutionary rate (1.63 × 10 -4 substitutions per site per year). All human and swine Influenza A viruses had a signifi-cantly faster evolutionary rate than avian viruses ( Table  2). In addition, evolutionary rates were significantly different even between lineages of same host. Hu2 has evolved more rapidly than Hu1, and Sw2 has evolved more rapidly than Sw1 ( Figure 2 and Table 2).

Selective Pressures
The selective pressures for the entire sequence (we defined the magnitude of the pressure as "ω") were 0.13 for the entire coding region of the M gene, 0.06 for M1, and 0.45 for M2 ( Figure 3). A higher selective pressure indicates that the gene (or the site) is under stronger selection (positive selection) for amino acid substitution. Lower selective pressure indicates that the gene (or the site) is under stronger negative selection to retain the same amino acid(s) because changes may lead to incompetence or abortion [44,45]. Selective pressure was statistically stronger in M2 than that in M1 for all hosts.
ω of the entire coding region of the M gene for human and swine influenza was significantly higher (no overlap of 95% confidence intervals) than that for the avian influenza ( Figure 3). ω for both M1 and M2 of human influenza are also significantly larger than that for avian influenza ( Figure 3).

Site-by-site Analyses
Site-by-site (by each codon) analyses for human influenza were conducted by SLAC (the entire tree [eSLAC], internal branches [iSLAC], and terminal branches [tSLAC]), and FEL (the entire tree [eFEL] and internal branches [iFEL]) methods [45]. We conducted the analyses by testing hypotheses for the entire tree, internal branches, and terminal branches (See "Materials and methods").
"dN/dS" indicates the magnitude of selective pressure on each codon. When dN/dS on a certain codon is significantly greater than 1, the site is considered to be under significant positive selection. When dN/dS on a certain codon is significantly smaller than 1, the site is considered to be under significant negative selection. Figure   position, i.e., the codon). Figure 5 shows that this site is located at the edge of the structure and is a part of a T-cell and MHC cell epitope. Of ten sites positively selected in M2, seven sites are in the extracellular domain (positions 11, 12, 13, 14, 16, 21, and 23), one site is in the transmembrane domain (position 27), and two sites in the cytoplasmic domain (positions 57 and 89, Table 3).
To define the evolutionary difference for each codon in human and avian influenza, we also calculated site-by-site selective pressures for avian influenza by eFEL. Consensus sequences of human and avian viruses were compared to identify major differences between these two hosts. We identified the sites at which consensus amino acids were different between the human and avian viruses and showed selective pressures ( Figure 6 and Table 4  Selective pressure among hosts Figure 3 Selective pressure among hosts. Selective pressures for the entire sequence (ω) are calculated for the entire coding region of the M gene, and separately for M1 and M2. Error bar shows 95% confidence interval. Selection profile by eFEL and eSLAC  Table 4).

Discussion
The phylogenetic tree showed that the M gene of influenza A viruses has evolved independently in each host. It revealed host-specific lineages, which were compatible with other reports. In previous reports, Av1, Av2, Sw1, Sw2, and CE were named as Eurasian (Old World) avian, North American (New World) avian, classic (old) swine, European (avian-like) swine, and recent (avian-like) canine lineages, respectively [1,8,46,47]. Since the emer-gence of the Russian Flu, both H1N1 and H3N2 have been co-circulating in human populations and undergoing different evolutionary processes, which have resulted in two distinct human influenza lineages, Hu1 and Hu2 ( Figure 1A, B, and 1F). Although reassortment of human influenza A viruses between the same subtype (intratypic recombination) has occurred frequently [48][49][50][51], we found only a few strains that seemed to be generated by reassortment between H1N1 and H3N2 human influenza, including H1N2 strains. These strains were not maintained in human populations. After Spanish Flu, the same M gene has been maintained in human influenza, even after two pandemics (Asian Flu and Hong Kong Flu) that were thought to have been generated by reassortment between avian and human influenza A viruses [1] ( Figure 1A and 1C). In the phylogenetic tree ( Figure 1A), Spanish Flu is located at the root of a human lineage and close to a swine lineage; there is a greater distance between Spanish Flu and the avian influenza A viruses identified around 1918. This result supports the hypothesis that an ancestral virus of Spanish Flu had entered the mammalian population before 1918 [53,54]. It remains to be seen whether this M gene will be retained after further pandemics. It was shown that the M gene of recent human influenza cannot incorporate the HA segment of avian influenza in vitro [55].
3D crystral structure of M1 Figure 5 3D crystral structure of M1. The figure was generated using BioHealthBase. M1 is identified as dimers. Site at position 219 (yellow circles), which is under positive selection for human influenza, is located at the edge of the structure. Figure 6 Consensus sequence. Consensus amino acid sequences of human and avian influenza A virus are shown. The major variable is defined as amino acid variants which are found in 10% or more strains. Different sites are shaded in red. There have been several sporadic infections with viruses from non-human lineages to humans, including the recent H5N1 infections in humans. However, these viruses were not maintained, and therefore, they disappeared from the human population without efficient transmission from human to human. In addition, it is implied that swine can be a "mixing vessel" in which human and avian viruses are reassorted to generate a human pandemic strain [1,56]. However, infections of strains with avian or human M genes in swine were also rare, and most of these viruses were not maintained in the swine population, except for the Sw2 lineage, in which viruses with the avian lineage M gene became established in the swine population.

Consensus sequence
Our phylogenetic analysis showed that viruses were clustered in host-specific lineages. This suggests that the M gene may be host specific and viruses with an M gene from other hosts are difficult to replicate. It is possible that the M gene determines the host range through the interaction between M1 and vRNPs [13,14,57]. An M gene that can match with host-specific vRNPs may be needed to replicate and transmit in a certain host. In addition, many studies have shown the interaction between M1 protein and host proteins, such as RACK1, MAPK, and core histone [13,[58][59][60]. The M gene may be directly and/or indirectly linked to host tropism of the virus.
The evolutionary rate of the M gene was low in avian viruses compared to human and swine viruses ( Figure 2 and Table 2). This result is rational because birds are considered to be a natural host for the influenza A virus [1]. The avian influenza A virus may have already been adapted to the host and not subject to pressure to induce further amino acid changes. This is also supported by the result showing that ω of the M gene was the lowest in avian influenza (Figure 3). Additional amino acid changes might be required in mammalian hosts to allow the viruses to adapt to these relatively new hosts. This stronger selective pressure on human and swine influenza may make human and swine influenza evolve more rapidly than avian influenza (Figures 2 and 3).
Interestingly, evolutionary rates were significantly different between lineages of the same host ( Table 2). The evo-  lutionary rates of Hu2 and Sw2 were faster than Hu1 and Sw1, respectively. The evolution of the M gene might not only be controlled by host species. One possible explanation is that strains in a lineage that appeared more recently such as Hu2 or Sw2, have to evolve more rapidly in order to be adapted better to the host than strains in other preexisting lineages (Hu1 or Sw1), which have already adapted to some extent. Social factors at the time when new lineages appeared such as the growth of the population and globalization may also facilitate a faster evolution. This may be the reason why the evolutionary rates of Hu2 and Sw2 are higher than those of Hu1 and Sw1, respectively ( Figure 2). However, reason of difference between evolutionary rates of Av1 and Av2 is unclear.

Summary of site-by-site analyses
The selective pressure is stronger in M2 than in M1 ( Figure  3) and more sites under positive selection were identified in M2 than in M1 (Table 3 and Figure 7). Among them, most of the sites (7 out of 10) under positive selection in M2 are located in the extracellular domain (Table 3 and Figure 7). Infection of influenza A virus induces the host's immune response to M2, especially to the extracellular domain [26][27][28]. It has been shown that antibodies recognizing the extracellular domain including the sites under positive selection confer protective immunity [37][38][39]. The host's immune response may make stronger selective pressure on M2 than that on M1. However, of course, selective pressure is much higher in the HA segment, the major antigenic component, than in the M2 gene [61], and this M2 gene is thus more conserved than the HA gene [42].
M1 is thought to play a vital role in the assembly and budding process [12,22,23]. Even minor mutations in M1 may cause a critical deficiency in virus replication. This could also explain why M1 is under strong negative pressure and why the selective pressure on M1 is smaller than that on M2 ( Figure 3). Nevertheless, the selective pressure on M1 of the human influenza was stronger than that of the avian influenza ( Figure 3). M1 of human influenza should be under stronger selective pressure than that of avian influenza to be better adapted.
Position 219 in M1 is under positive selection in human influenza. It was also reported that this site was positively selected using a different method of calculation [62]. However, this site is under negative selection in avian influenza (Figure 7). M1 is recognized by cytotoxic T cells [40,63,64] and the C-terminal of M1 determines antigenicity [65,66]. The site, located at the edge of structure ( Figure 5), is part of the T-cell and MHC epitope. M1 may also be under selective pressure from the host's immune response, although this is weaker than M2. Besides, the Cterminal of M1 is important for binding to vRNPs [16]. This site might play an important role in the interaction with vRNPs, being associated with host range. Therefore, it is under positive selection only in the human and not in avian influenza virus.
Positions 115 and 121 in M1, which are under significant negative selection in both human and avian influenza, had different consensus amino acids between these two hosts ( Figure 7). These results indicate that these sites may be important for host tropism and are therefore under negative selection. In addition, position 137 also has different consensus amino acids between the hosts, though this site is not under significant negative selection in human influenza (the site is under negative selection in avian influenza). The two domains in M1 have been reported to affect the disposition of viral RNA. One domain resides in a palindromic stretch of basic amino acids (position 101 to 105) [17,18] and the other domain is located at position 148 to 162 containing a zinc finger motif [19,20]. The three sites (positions 115, 121, and 137) are located between these two domains. These sites might affect the disposition of viral RNA and be involved in the determination of host range.
Position 27, which is a site in the transmembrane domain, is positively selected in M2. This site is associated with amantadine resistance [67]. The selective pressure on the site may be due to drug pressure. However, we could not show any positive pressure on position 31, which is associated with the recent spread of amantadine resistance [68]. Details on drug pressure and possible mechanism for recent surge of amantadine-resistant strains will be described in another manuscript (in preparation).
The cytoplasmic domain of M2 is important for interaction with M1, genome packaging, and formation of virus particles [33][34][35][36]. Two sites are under positive selection in the cytoplasmic domain of M2 (positions 57 and 89, Table 3). In particular, position 57 showed different consensus amino acids between human and avian influenza (Figure 7). These results indicate that the amino acids in these sites have frequently changed, and these sites are likely to be involved in several functions of M2. The M2 cytoplasmic tail (position 45 to 69) has been shown to be a binding domain for M1 [35]. Position 82 to 89 is important for infectious virus production [35]. Another study showed that vRNP packaging is mediated by amino acids at position 70 to 89 of the M2 gene [69]. The M2 gene must, therefore, have evolved with several functions.
In conclusion, the M gene of the influenza A virus has evolved with different selective pressures on M1 and M2 among different hosts. We found potentially important sites that may be related to host tropism and immune responses. These sites may be important for evolutionary processes in different hosts and host adaptation. How-ever, Dunham et al. concluded that it is difficult to predict what specific genetic changes are needed for mammalian adaptation by comparing evolution of avian and swine influenza A viruses [47]. Further studies to clarify the specific role of each site identified in the present study are needed.

Sequence Data
All data were obtained from the influenza sequence database (Influenza Virus Resource on: http:// www.ncbi.nlm.nih.gov/genomes/FLU/FLU.html, accessed on July 21, 2008) [70]. All sequencing data for the strains with a full-length M gene of any subtypes of influenza A from different host species including avian, canine, equine, human, and swine were included. Sequences derived from laboratory strains and duplicate strains verified by the strain name were excluded. A total of 5489 sequences were obtained [accession numbers are listed in additional file 1]. After excluding sequences containing ambiguous nucleotides, minor insertions, minor deletions (data for full length of coding region were used) or premature termination codons, a total of 5060 sequences were used in the analysis. Sequencing data were obtained together with information of the host, subtype, isolation year, and isolation place. The sequencing numbers for the influenza of each host are listed in Table 1. A multiple alignment of the nucleotide sequences, which did not contain any gaps, was constructed using ClustalW.

Phylogenetic Tree Analysis
A phylogenetic tree was inferred by RAxML [71]. The sequences data only for the coding region were used; i.e., at nucleotide position 26 to 1007. The basic sequential algorithm of RAxML is described elsewhere [72]. RAxML is one of the fastest and most accurate sequential phylogeny programs [73]. In this method, a rapid bootstrap search is combined with a rapid maximum likelihood search on the original alignment. The tree was constructed using Web-servers, RAxML BlackBox: "http://phy lobench.vital-it.ch/raxml-bb/" [71]. The tree is colorcoded according to hosts, subtypes, geographical information, or temporal information using FigTree (ver.1.1.2).

Dataset of Influenza for Each Host
Datasets for each host (avian, canine/equine, human, and swine hosts) were constructed. Sequences only from the host-specific lineage in the phylogenetic tree were used. For example, the H5N1 influenza A viruses that had infected humans were excluded from the analyses because humans were accidental hosts infected with the viruses of an avian lineage. Identical nucleotide sequences in the same dataset were removed before further analyses.
The number of base substitutions per site from an average of all sequence pairs was calculated to define the diversity of sequences in each dataset (Table 1) using the maximum composite likelihood method in MEGA (ver. 4) [74].

Evolutionary Rate
The evolutionary rate of each lineage was calculated. To calculate the rate, at least one sequence of each subtype in each year was selected from each dataset. Evolutionary rate was analyzed for the selected sequences as the number of substitutions per site per year compared to the oldest strain in each lineage with a linear regression model. The significance of the correlations was estimated using the Pearson correlation. Differential between slopes of the tangent of simple regression lines were tested by analysis of covariance. The analyses were conducted using SPSS (ver.17).

Consensus Sequence
Consensus amino acid sequences were determined as the sequence of amino acids that were identified most frequently at each position in a dataset, for human and avian influenza. Amino acid substitutions that were identified in more than 10% of the strains were regarded as major variants.

Evaluation of Pressure (ω)
Phylogenetic trees for each dataset by hosts were constructed using the maximum-likelihood method implemented in PhyML-aLRT [75] with the GTR model (four rate categories, all parameters estimated from the data).
Selective pressure for each host population was calculated using the trees. Selective pressure was analyzed by HyPhy [76]. All analyses in HyPhy were conducted after identifying the best fit model from every possible time-reversible model (e.g., F81 and HKY85) according to Akaike's information criterion [45,77].
Global estimates (ω) of relative rates of non-synonymous (dN) and synonymous (dS) substitutions averaged over the entire alignment were compared to calculate the overall strength of selection [45].

Site-by-site Selective Pressure (dN/dS)
Positive selection sites for human influenza were detected using two methods: single likelihood ancestor counting (SLAC) and fixed-effects likelihood (FEL). FEL was also conducted for avian influenza. The relative rates of nonsynonymous and synonymous substitutions were compared. Sites where dN/dS > 1 and dN/dS < 1 were inferred as positively and negatively selected, respectively. The details of the two methods is described elsewhere [45,78,79]. It was shown that many recent non-synonymous substitutions, i.e., those in the terminal branches of the tree, were not represented on internal branches [80]. At codons where internal substitutions are seen, the strength of selection along the terminal branches is high. Analyses were conducted exclusively by testing hypotheses for the entire tree, internal branches, and terminal branches.
Briefly, in the SLAC method, the nucleotide and codon model parameter estimates are used to reconstruct the ancestral codon sequences at the internal nodes of the tree. The single most likely ancestral sequences are then fixed as known variables, and applied to infer the expected number of non-synonymous or synonymous substitutions that have occurred along each branch, for each codon position. SLAC is a substantially modified and improved derivative of the Suzuki-Gojobori method [44]. The FEL method is based on maximum-likelihood estimates. The FEL method estimates the ratio of non-synonymous to synonymous substitutions on a site-by-site basis for the entire tree (eFEL) or only the interior branches (iFEL). iFEL is essentially the same as eFEL, except that selection is only tested along the internal branches of the phylogeny [80].
In the present study, we used a newly developed softwares RAxML [71] and HyPhy [76] for phylogenetic analyses. Markov Chain Monte Carlo (MCMC) and PhyML are widely used and are considered useful for manipulating data sets (hundreds) [2,[81][82][83][84]. However, they cannot process huge data sets in order of thousands on an ordinary desktop computer. Therefore, we used a PC with Windows operating system and 3 GB RAM). PAML, which is a common software package for phylogenetic analyses such as the calculation of selective pressure using maximum likelihood [85], also failed occasionaly in analyzing large data sets. We therefore used the recently developed software, which overcomes these problems. The accuracy of this software has been confirmed [73,79,86]. Through site-by-site analyses, we identified more sites negatively or positively selected than those in a study by Suzuki [61]. This was ascribed to a difference in the number of data sets and/or algorithms, as we analyzed ten times more than the number analyzed in his study.