Evolution and Vaccination of Influenza Virus

Abstract

In this study, we present an application paradigm in which an unsupervised machine learning approach is applied to the high-dimensional influenza genetic sequences to investigate whether vaccine is a driving force to the evolution of influenza virus. We first used a visualization approach to visualize the evolutionary paths of vaccine-controlled and non-vaccine-controlled influenza viruses in a low-dimensional space. We then quantified the evolutionary differences between their evolutionary trajectories through the use of within- and between-scatter matrices computation to provide the statistical confidence to support the visualization results. We used the influenza surface Hemagglutinin (HA) gene for this study as the HA gene is the major target of the immune system. The visualization is achieved without using any clustering methods or prior information about the influenza sequences. Our results clearly showed that the evolutionary trajectories between vaccine-controlled and non-vaccine-controlled influenza viruses are different and vaccine as an evolution driving force cannot be completely eliminated.

1. Introduction

The rapid growth of the influenza genome sequence data due to the advanced development of sequencing technology in recent years has provided the opportunity for a more comprehensive sequence analysis of the influenza virus. The difficulty in sieving through and making sense of this mountain of data relying solely on phylogenetic approaches has become increasingly limited, in part, due to the poor scalability of the relevant algorithms (Nicholas, 2007). Therefore, a different methodology needs to be utilized to take advantage of the massive amount of available data but at the same time to be able to expose important information or structure within the data. Here, we present an application paradigm in which an unsupervised machine learning approach is applied to the high-dimensional influenza genetic sequences so that the evolution of the vaccine-controlled and non-vaccine-controlled influenza viruses in the past century can be visualized. The main objectives of this study are twofold: (1) to visualize the evolution trajectories of influenza under vaccine pressure and in the wild without using any prior information about the viruses and (2) to provide statistical confidence to support the visualization results. The influenza virus is believed to have originated from a natural reservoir consisting of wild aquatic birds (Webster et al., 1992; Taubenberger and Kash, 2010).

The influenza A virus is divided into subtypes based on differences in the surface proteins hemagglutinin (HA) and neuraminidase (NA), which are targets of the human immune system. Antigenic variants or immunologically distinct strains of A/H1N1, A/H3N2, and Type B have continued to emerge since its introduction into humans (Schweiger et al., 2002). Vaccination is the main strategy in stopping the infection and transmission of the virus in humans (Hannoun, 2013). There are three components in a seasonal flu vaccine: (1) A/H1N1, (2) A/H3N2, and (3) Type B influenza. Each component is designed to fight the specific strain in each subtype that is predicted to be the dominant circulating strain in the upcoming flu season. Over the years, there have been more than 24 vaccine updates for the A/H3N2 strain, more than 17 updates for the Type B strain, and 10 updates for the A/H1N1 strain. Each vaccine update is designed to provide immunity to the new antigenic variant that has emerged from the previous flu season. However, the long-term effects of vaccination on the evolution of the virus itself are not clear.

To shed light on this seemingly unsuspected problem, we used the nucleotide sequences from seasonal human A/H3N2 influenza virus from 1971 to 2016 as an example to demonstrate the evolutionary progress of this influenza virus against each successive vaccine introduction from 1971 to 2016. Figure 1 shows progression of influenza evolution based on the nonsynonymous substitutions (dN) and synonymous substitutions (dS) ratio analysis by using the HA1 domain of the HA gene from A/H3N2 virus. However, some flu seasons did not provide sufficient information of the vaccine strain accession number; thus, the actual plot did not ascend at 2016. The HA1 domain is a hypervariable domain of the HA gene where constant mutational changes can be observed due to the immune pressure generated from the host. A \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$dN / dS$$ \end{document} ratio greater than 1 indicates that the site is under positive selection pressure and is undergoing molecular adaptation. In Figure 1, a constant shift of positively selected sites (blue color: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$dN / dS$$ \end{document} ratio greater 1) could be observed whenever a new vaccine (green square) was introduced, which indicated that a new antigenic variant had emerged.

FIG. 1.

Seasonal human A/H3N2 influenza \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$dN / dS$$ \end{document} ratio analysis against time of vaccine introduction. A constant shift of positively selected site location when a new vaccine was introduced. Horizontal axis represents the position of HA1 domain of the HA gene. Vertical axis represents time progression from 1971 (bottom) to 2016 (top) when each new (green square) and repeated (black square) vaccine was introduced. Red color bars denote the range of positions, with \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$dN / dS$$ \end{document} ratio from 0.8 to 1. Blue color bars denote the range of positions, with \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$dN / dS$$ \end{document} ratio greater than 1. HA, hemagglutinin.

When a repeated vaccine was introduced, the positively selected sites identified from the previous season tended to remain unchanged; even though there were some years that blue sites changed, red regions indicated the possible regions where the actual positively selected sites might have settled. If multiple vaccines are recommended and at least one of them is recommended in the previous year, one still regards that the current year has updated the vaccine recommendation; however, if the same multiple vaccines are recommended for two successive years, they will be regarded as repeated strains, and different strains of this set will be run in the program for these 2 years. Given the results from the \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$dN / dS$$ \end{document} ratio analysis, we compared the evolution trajectories of vaccine-controlled with non-vaccine-controlled influenza viruses and sought to better understand the effect of vaccination on the evolution of the influenza virus. In the present study, we used the human A/H3N2, A/H1N1, Type B, and avian H5 HA sequences as the vaccine-controlled samples. We used the human H5N1 and avian H5N1 HA sequences as the non-vaccine-controlled samples.

2. Background

Influenza viruses have the ability to infect a very broad range of avian and mammalian hosts. Their genomic diversity is acquired through two biological mechanisms: antigenic drift and antigenic shift (Webster et al., 1992). Antigenic drift consists of the accumulated and continual mutations on surface proteins, resulting in the generation of antigenic variants. Of these surface proteins, we are focused on the HA protein. Antigenic shift occurs when complete gene segments are exchanged among different subtypes of influenza viruses within a host cell, resulting in what effectively amounts to a whole new influenza virus genome. Both antigenic drift and antigenic shift allow for the virus to evade the host's immune response and to rapidly adapt to new hosts (Suzuki, 2006; Caron et al., 2009). The evolution of influenza A virus is driven by the high rate of mutations and the ability to reassort gene segments. Because of its high rate of mutation combined with the lack of error-correcting mechanisms during replication, the influenza virus can easily generate different phenotypes that have the ability to survive within its host and to infect others.

To keep track of the evolution of the virus, an annual update to the influenza vaccine composition is needed to provide a vaccine-induced immunity to the general public (Boni, 2008). The main process in influenza vaccine strain selection is to assess the match between the vaccine strain and the currently circulating strains and the potential new antigenic variant (Russell et al., 2008). If the vaccine strain does not match the currently circulating strains or the new antigenic variant that is likely to be the major variant in the upcoming influenza season, the vaccine composition is updated to contain a representative of the new variant (Russell et al., 2008). Each vaccine update is designed to provide immunity to the new antigenic variant that has emerged from the previous flu season. The seasonal influenza vaccine is used to prevent the infection and transmission of the virus, but its effect on the evolution of the virus itself is not clear.

3. Methods

In this study, utilizing the online NCBI influenza database (Bao et al., 2008), we collected HA sequences from human A/H3N2, A/H1N1, Type B, and avian H5 HA sequences that represent the vaccine-controlled samples. We also collected human H5N1 and avian H5 HA sequences that represent the non-vaccine-controlled samples. Table 1 lists the year range and number of HA nucleotide sequences from each sample. To focus on evolution trends in more recent years, we incorporated recent year data collected from 2007 or 2008 onward and randomly chose the same number of samples for each year.

Table 1.

Vaccine-controlled and Non-vaccine-controlled Human and Avian Sequences

Samples	Year	Sequences
Human A/H1N1	1918–2013	2140
Human A/H1N1^a	2008–2016	1440
Human A/H3N2	1968–2009	175
Human A/H3N2^a	2007–2016	1168
Human B (Vict/Yam)	1970–2013	818
Human B (Vict/Yam)^a	2008–2016	920
Human H5N1^b	1997–2012	127
Avian H5 (Mexico)	1994–2002	32
Avian H5 (China)^b	1997–2002	32

Vaccine.

Non-vaccine.

3.1. Influenza evolution visualization

All genetic sequences were first converted into binary strings according to the method outlined in Lam et al. (2012). Nucleotide sequences are represented by strings of characters out of an alphabet of four letters: A, C, G, T. To obtain the binary string, each letter is replaced by a code of 4 bits: 1000, 0100, 0010, and 0001, respectively. In the meantime, nucleotide sequences have some “wild cards” for cases where a single A, C, G, or T cannot be perfectly determined. We, therefore, regarded that they might appear equally possibly. For instance, if Y stands for C or T, then we use 0101 to represent it. All binary strings were collected into a matrix to which Principal Component Analysis (PCA; Jolliffe, 2002) was applied to extract the dominant variation from the dataset. Here, we briefly outline the sequence of steps involved in the PCA. Consider a data matrix \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${X_{m , n}}$$ \end{document} of dimensions m by n, with m being the number of strains and n being the number of sites or positions (in this case, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$n = 987 \times 4 = 3948$$ \end{document} for nucleotide sequences). Each row of X corresponds to a strain of virus, and each column of X corresponds to a particular position. We first center the columns of the data matrix X with \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\hat X = X - { \textstyle \frac { 1 } { m } } e { e^T } X$$ \end{document} , where e is a column vector of all ones, and then we obtain the sample covariance matrix C from \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\hat X$$ \end{document} by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$C = { \textstyle \frac { 1 } { { ( m - 1 ) } } } { \hat X^T } \hat X$$ \end{document} . C is a square symmetric \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$n \times n$$ \end{document} matrix whose diagonal entries are the variances of the individual sites across strains, and the off-diagonal terms are the covariances between different sites.

The PCA algorithm is then applied to matrix C. The result is then visualized by plotting the top two or three principal components of the projected data. Since each strain is encoded as a binary string and PCA works at the binary data level, the pairwise distance relationship between the strains in a reduced space can be understood as follows: Let \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\parallel {s} - {t}{ \parallel _H}$$ \end{document} denote the pairwise Hamming distance between two strains \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$s , t$$ \end{document} (number of differences in genetic sequences). Let \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\parallel s - t{ \parallel _{bin{ \kern 1pt} \;1}}$$ \end{document} , \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\parallel s - t{ \parallel _{bin{ \kern 1pt} \;2}}$$ \end{document} denote the distance between the binary encodings of the two sequences (1-norm and 2-norm, respectively), and let \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\parallel s - t{ \parallel _{proj}}$$ \end{document} denote the 2-norm distance in lower-dimensional space after projection onto the leading principal components. Every single change in the genetic sequence alphabet corresponds to changes to 2 bits in the binary encoding. Hence, we have the relationship between the distance in the lower-dimensional space shown on the plots, with the Hamming distance among the original sequences: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\parallel {s} - {t} \parallel _{proj}^2 \le \; \parallel {s} - {t} \parallel _{bin \;{ \kern 1pt} 2}^2 = \parallel {s} - {t}{ \parallel _{bin \;{ \kern 1pt} 1}} = 2 \parallel {s} - {t}{ \parallel _H}$$ \end{document} .

3.2. Quantification

To provide statistical support to the graphical results obtained, we performed a statistical analysis based on a method that combined a multi-class scatter matrix computation and class labels randomization. The projected data points served as the viruses' two-dimensional (2D) coordinates, and the year of isolation of each virus served as the class label. The multiclass scatter matrix involves the computation of between-class matrix (B) and within-class matrix (W; Box 1). These computed matrices were not used explicitly as we only sought the trace of B and W. These are just the scalar scatter values: sum of squared distances between points and their respective centers. The class separateness measure \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \lambda _o}$$ \end{document} is the ratio of trace B over trace W. A large \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \lambda _o}$$ \end{document} indicates that the classes or clusters are well separated between each other and that elements within a cluster are strongly related or share the same property. This is basically an estimate on how well a multi-class Fisher's linear discriminant could separate the classes (Alpaydin, 2010). A class label randomization algorithm (Algorithm 1) provided the “distance measure” as a surrogate for the probability of observing the observed \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \lambda _o}$$ \end{document} by chance. This is because the area under the tail of the randomized \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\lambda$$ \end{document} distributions beyond the observed separateness values was below a rounding error of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${10^{ - 16}}$$ \end{document} , which made the computation of p value not possible. The larger the distance, the less likely the observed \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \lambda _o}$$ \end{document} is generated by chance. We have also computed the \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\lambda _o^*$$ \end{document} by using the new data collected from 2007 or 2008 onward.

Box 1:

Virus isolation year as class label

C: Number of classes

N_i: Number of data points in class \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$i = 1 , 2 , \ldots C$$ \end{document}

• \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\lambda = { \frac { tr ( B ) } { { tr } ( { W } ) } } $$ \end{document}

• B: Between-class scatter matrix

– \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\sum \nolimits_i^C ( {u_i} - M ) ( {u_i} - M{ ) ^T}$$ \end{document}

– \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$M = \frac { 1 } { c } \sum \nolimits_i^C { u_i } $$ \end{document} “global mean of dataset”

• W: Within-class scatter matrix

– u_i: Mean of class i

Algorithm 1: Estimate Separateness Measures
Let \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ { \lambda _o } = { \textstyle { \frac { tr ( { B_o } ) } { tr ( { W_o } ) } } } $$ \end{document} be the observed separateness value.
Repeat \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$j = 1:K2$$ \end{document}
Repeat \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$i = 1:K1$$ \end{document}
Generate a randomization of the class labels
Compute the within-cluster scatter W
Compute the ratio \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ { \lambda _i } = { \textstyle { \frac { tr ( B ) } { { tr } ( { W } ) } } } = { \textstyle { \frac { tr ( T ) - tr ( W ) } { { tr } ( { W } ) } } } $$ \end{document}
Compute the mean \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\mu$$ \end{document} and std \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\sigma$$ \end{document} for all \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \lambda _{i = 1 , ..K1}}$$ \end{document}
Compute the distance \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ { d_j } = { \textstyle \frac { { \mu - { \lambda _0 } } } { \sigma } } $$ \end{document}
Compute the mean \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\bar d$$ \end{document} and std \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\hat d$$ \end{document} of all \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${d_{j = 1..K2}}$$ \end{document}
Report the distance of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \lambda _o}$$ \end{document} from the mean in the form of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\bar d \pm \hat d$$ \end{document}

4. Results

The application of the high-throughput unsupervised method to the high-dimensional influenza virus genetic sequence data has made possible the visualization of the evolution of the influenza virus in the span of almost half a century. In this study, we present the graphical results from visualization of vaccine-controlled and non-vaccine-controlled influenza viruses based on their genetic sequences alone. The human influenza A/H3N2 has the highest number of vaccine updates among the three vaccine-controlled influenza viruses circulating in humans. Given the observation about the constant shifting of positively selected sites whenever a new vaccine was introduced, we sought to visualize the evolution trajectories of vaccine-controlled and non-vaccine-controlled influenza samples. We also set out to compute the separateness values of class or clusters for both vaccine-controlled and non-vaccine-controlled samples by using the multi-class scatter matrix computation method for both the before and after class labels randomization process. We performed 1000 runs of Algorithm 1 on these samples and listed the results in Table 2. The observed separateness values \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \lambda _o}$$ \end{document} of vaccine-controlled samples are consistently higher than those of the non-vaccine-controlled samples. This suggested that the vaccinated samples have a very good separability by isolation years.

Table 2.

Class Separateness Results: Vaccine- and Non-vaccine-controlled Human and Avian Samples

Sample (human)	\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \lambda _o}$$ \end{document}	\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\lambda_o^*$$ \end{document}	Distance	Distance
A/H3N2 (1968–2009/2007–2016^b)	30.5	2.72	978.3 ± 0.031	1031 ± 0.031
B: Victoria (1970–2013/2008–2016^b)	26.3	4.09	1310 ± 0.02	499.1 ± 0.0154
B: Yamagata (1970–2013/2008–2016^b)	25.3	0.68	1327.8 ± 0.019	105.2 ± 0.0029
A/H1N1 (1918–2013/2008–2016^b)	24.7	5.74	617.2 ± 0.04	2787 ± 0.82
H5N1 (1997–2002)^a	1.0	—	34.8 ± 0.029	—

Sample (avian)	\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \lambda _o}$$ \end{document}	\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\lambda_o^*$$ \end{document}	Distance	Distance
Avian H5 Mexico (1994–2002)	1.7	—	12.23 ± 0.11	—
Avian H5N1 China (1997–2002)^a	0.268	—	3.16 ± 0.0.6	—

Non-vaccine.

Vaccine.

In Figure 2, we observed that the human A/H3N2 viruses clustered around vaccine seed strains chronologically since their introduction into humans in 1968. The evolution trajectory is directional, going from lower left to lower right in the figure. In Figure 3a, two separate lineages of human Type B influenza are co-circulating and each lineage shows the same observational characteristics as the A/H3N2; type B viruses are also clustered around vaccine seed strains. Next to it, Figure 3b shows the trajectories of two lineages (Yamagata and Victoria) that eventually got completely separated. Even though it is only the situation in the United States, separateness is quite clear. Black circles as vaccine strains are shown only in one of the lineages because only one virus in the United States was chosen for a vaccine strain, which is from Yamagata lineage. Since the classical H1N1 virus was replaced by the pandemic swine H1N1 strain, we focused on its evolution trend from 2008 onward. H1N1 behaved actively as two jumps (or “V” shapes), which are seen in Figure 4. One was in 2013, whereas the other was in 2015. These two discontinuities indicated abrupt changes in gene sequences from their correspondingly previous years, so the evolution trends suddenly jumped in completely different directions. However, Hamming distance plot shows that even the mutation rate between 2008 and 2016 was considerably low, ∼3%, let alone the mutation rate between two consecutive years. There is, hence, not sufficient evidence to speculate about mutation. A vaccinated avian sample was used (avian H5) to further understand the evolution characteristic of vaccine-controlled influenza.

FIG. 2.

Seasonal human A/H3N2 influenza virus evolution trajectory. Each arrow points to a vaccine seed strain (red dot). The directional evolution can be seen as traveling from lower left to the top and then coming down to the lower right.

FIG. 3.

Seasonal human Type B influenza virus evolution trajectory. Two separate lineages (Victoria and Yamagata) are evolving simultaneously: (a) top to lower left and to lower right and diverging further (b). Vaccine introductions are indicated by year labels in (a) and shown as black circles in (b).

FIG. 4.

Seasonal human H1N1 influenza virus evolution trajectory in three dimensions. Vaccine strains are black crosses.

In late 1993, an outbreak of avian H5 influenza in poultry in Mexico was detected and a long-term vaccination program was implemented in the hope of bringing the outbreak under control and of eradicating the virus (Lee et al., 2004; Escorcia et al., 2008). The vaccination program was in effect for more than 13 years, but an increase in the respiratory signs of disease was observed in vaccinated chickens (Escorcia et al., 2008). In other words, the vaccine strain used in the vaccination program no longer matched the circulating strain in the field. The vaccine strain (A/Ck/Mexico/CPA-232/1994) was isolated in 1993 and has been in use for the duration of the program for more than a decade. Using the available genetic HA sequences from these vaccinated chicken, we produced a three-dimensional (3D) PCA plot (Fig. 5) to show the evolution of the field isolates from 1994 to 2002. The first observation from Figure 5 is that a directional evolutionary trend similar to other vaccinated samples can be seen in this figure. Second, a chronological pattern is obvious, indicating that the virus had undergone constant evolution or antigenic drift away from the early strains. A split in the evolutionary path can be seen occurring in the 1990s. This split or divergence has been reported in studies by Lee et al. (2004) and Escorcia et al. (2008) based on phylogenetic analyses conducted on the same sequence sample.

FIG. 5.

Vaccine-controlled avian H5 influenza virus evolution trajectory in three dimensions. The vaccine was introduced in early 1990s, and the virus slowly evolved away from the vaccine strain and established two separate lineages.

Figure 6 illustrates the evolution trajectory of the non-vaccine-controlled human H5N1 influenza from 1997 to 2002. We included the human H5N1 virus as the “control” since this subtype is not currently being vaccinated against in humans but is under active research due to its high mortality rate in infected humans. Figure 6 suggests that this subtype has evolved into a few dominant clusters since 1997. Three major evolutionary trends or clustering patterns can be seen originating from the center cluster that contains viruses from 1997. This also implies that this influenza subtype has undergone HA gene diversification. Although it has diversified since 1997, the specific H5 HA gene identified in 1997 has remained present in these days (Wei et al., 2012).

FIG. 6.

Non-vaccine-controlled human H5N1 influenza virus evolution trajectory in three dimensions. The virus has evolved into a few dominant lineages since 1997. Three major evolutionary lineages can be seen originating from the center cluster, which contains viruses from 1997. However, the specific H5 HA gene identified in 1997 has remained present in these days.

Figure 7 shows the evolution of the non-vaccine-controlled avian H5 influenza virus. The overall observation that arises from this figure is that rather than forming a restricted directional trend, the evolution of the virus is characterized by a collection of clusters scattered on the plot. The collection of clusters suggests a diverse pool of the genetic diversity of the virus. For the avian H5 subtype, a less focused evolutionary trend than vaccine-controlled influenza viruses can be observed. The increased genetic diversity since 2000 has been observed by Garcia et al. (1997) and is captured in this figure, with clusters scattered to the left and extended to the upper and lower corner at almost the same time. This clearly suggests the co-circulation of multiple clades or sublineages of the avian H5 subtype. The diverse genetic diversity of the avian H5 represented by multiple clusters across a long time period indicated that the avian subtype in the wild evolves much slower than seasonal human influenza viruses.

FIG. 7.

Non-vaccine-controlled avian H5 influenza virus evolution trajectory in three dimensions. Multiple clusters scattered throughout, sharing almost the same time periods, thus suggesting the co-circulation of multiple clades or sublineages of the avian H5 subtype.

5. Discussion and Conclusions

Vaccination is the principal measure for preventing influenza and for reducing its impact (Wood et al., 2001; Webby et al., 2004). Almost a century ago after the isolation of the first influenza virus, influenza vaccines have been persistent and have evolved to respond to the evolution of the influenza viruses evolving in humans (Gunn et al., 2010; Hannoun, 2013). Antigenic drift of influenza viruses occurs frequently among circulating strains, which leads to new antigenic variants. However, whether the drift mechanism occurs with the presence of vaccine pressure is an important question that needs to be addressed at a different level as vaccination is the primary method in the prevention and protection of humans against the influenza virus. Two studies (Lee et al., 2004; Hensley et al., 2009) have shown that vaccination forces mutations on the HA protein of the influenza virus. These mutations changed the way in which the virus gradually evolved and adapted to a new vaccine-protected environment. Here, we extended the spectrum of analysis to include vaccine-controlled human and avian samples and non-vaccine-controlled human and avian samples to better compare and contrast and understand the evolutionary dynamic of influenza viruses under vaccine pressure.

Using vaccinated and non-vaccinated samples from both human and avian hosts, we hope to minimize potential data selection bias and at the same time to provide a fair comparison across hosts under vaccination pressure. Our method utilized only the genetic composition of the HA sequences alone without using any specific clustering algorithms. As mentioned earlier and shown in Figure 1, genetic sequences contain important signals that are used to detect evolutionary trends between different influenza subtypes under vaccination pressure. The genetic composition combined with the implicit positional information of the HA gene is enough to provide clues that the vaccine-controlled influenza viruses are under pressure to mutate to escape immune responses. Our method takes advantage of the binary coding of each sequence that preserves the positional information of each HA gene.

In this study, we have demonstrated that the evolutionary trajectories for vaccine-controlled influenza are directional and restricted. The restricted directional evolutionary trends and clusters formation around the vaccine strains along the evolutionary paths exhibited by the vaccine-controlled influenza viruses are in sharp contrast to the non-vaccine-controlled influenza viruses. Apart from this distinction, the naturally emerged chronological ordering of vaccine-controlled influenza viruses in both 2D and 3D visualizations is much more noticeable than the non-vaccine-controlled viruses. This natural chronological ordering reflects the active adaptation of the viruses to their changing environment. The class separateness measure exposes the fact that vaccine-controlled influenza viruses that share the same isolation year have the tendency to cluster tightly together with good separateability. Each separate cluster or group represents a distinct genetic diversity of the virus group. In contrast, non-vaccine-controlled influenza viruses that were isolated within the same time period appeared to be more scattered, and the clusters exhibited a much larger within-cluster distance with no narrow restricted bands being observed. These observations suggested that the mutations on the HA gene were not restricted to certain sites alone and that the majority of these mutations were most likely synonymous nucleotide substitutions on the HA gene.

Also, the number of clusters observed are almost identical to the number of vaccine updates for the seasonal human A/H3N2 and influenza B viruses. The very low value of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \lambda _o}$$ \end{document} computed from non-vaccine-controlled influenza viruses has clearly captured the fact that non-vaccine-controlled viruses are not actively evolving by the year. In contrast, the vaccine-controlled influenza viruses have been actively evolving and adapting to the changing environment constantly as new vaccine composition is being introduced year after year. This is clearly reflected in the very high \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \lambda _o}$$ \end{document} value for vaccine-controlled influenza viruses. Evidently, the same observation can be drawn from \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\lambda _o^*$$ \end{document} values. The absolute \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\lambda _o^*$$ \end{document} values are notwithstanding far less than those from longer year spans since \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\lambda _o^*$$ \end{document} describes a shorter time, which makes far less groups.

Another consequential observation is that Yamagata lineage shows far less separateness values than any other \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\lambda _o^*$$ \end{document} . It may confound one if one looks at Figure 3b because Yamagata lineage, at the left, seems to have a very well-spread line. Nonetheless, if one looks carefully, viruses from each year as a matter of fact bounce back and forth along this line, which explains them not being successfully separated. Although our analysis was based on genetic sequences alone, the results suggested that a clear difference existed among influenza viruses evolving in a vaccine-protected environment than in the wild. This difference is shown through the multi-class scatter computation of their evolutionary paths. This quantitative measurement also serves as a basic statistical support to the observed differences in the evolution dynamics between vaccine-controlled and non-vaccine-controlled influenza viruses.

There are other potential factors besides vaccination that can affect the evolution of influenza viruses, such as host-specific immune response, the large difference in life expectancy between humans and avian species, vaccine efficacy and effectiveness, the transmission channel of the virus in a different environment, and geographical regions. These factors have not been considered in this present study because our overall objective is to present a genetic sequence-only approach as the first step in understanding the evolution of influenza viruses in a protected environment. Our approach works directly at the sequence level with no prior assumption about the evolution of the virus. It is a departure from the traditional one-dimensional phylogenetic approach in that we visualize influenza evolution in 2D and 3D space. All phylogenetic methods make or rely heavily on the assumptions about the underlying evolutionary process (Jenkins et al., 2002).

By using methods that avoid making assumptions about the parentage relations among the strains, we can avoid possible misinterpretation of the results. As has been shown in this article, a data-driven approach with no prior assumptions about the evolution of the influenza virus affords us a different perspective in directly visualizing how the virus evolves in a span of more than half a century. This perspective has given us insight into the way we think about the driving forces behind the emergence of human seasonal influenza antigenic variant strains season after season. Perhaps, vaccination did play a role in forcing the virus to undergo a different evolutionary path to continue to establish itself in its occupied host. A definitively scientific conclusion cannot be drawn without a thorough study of the virus in a controlled experiment for an extended period, which should, nonetheless, include multiple influenza epidemics in humans.

Footnotes

Acknowledgments

This research was supported in part by NSF grants IIS 1319749. Influenza research in Srinand Sreevatsan lab is funded by the National Institute of Allergy and Infectious Diseases, National Institutes of Health, and Department of Health and Human Services, under Contract No. HHSN266200700007C. Its contents are solely the responsibility of the authors and do not necessarily represent the official views of the NIH or NSF. The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the article.

Author Disclosure Statement

No competing financial interests exist.

References

Alpaydin

2010. Introduction to Machine Learning, 2nd ed. MIT Press, Cambridge, MA.

Bao

, Bolotov

, Dernovoy

, et al. 2008. The influenza virus resource at the national center for biotechnology information. J. Virol. 82, 596–601.

Boni

M.F.

2008. Vaccination and antigenic drift in influenza. Vaccine, 26(Suppl 3), C8–C14.

Caron

, Gaidet

, de Garine-Wichatitsky

, et al. 2009. Evolutionary biology, community ecology and avian influenza research. Infect. Genet. Evol. 9, 298–303.

Escorcia

, Vzquez

, Mndez

S.T.

, et al. 2008. Avian influenza: Genetic evolution under vaccination pressure. Virol. J. 5, 15.

Garcia

, Suarez

D.L.

, Crawford

J.M.

et al. 1997. Evolution of h5 subtype avian influenza a viruses in North America. Virus Res. 51, 115–124.

Gunn

J.L.

, Craddock

, and Giles-Vernick

2010. Influenza and Public Health: Learning from Past Pandemics. Earthscan, Florence, KY.

Hannoun

2013. The evolving history of influenza viruses and influenza vaccines. Expert Review of Vaccines, 12, 1085–1094.

Hensley

S.E.

, Das

S.R.

, Bailey

A.L.

, et al. 2009. Hemagglutinin receptor binding avidity drives influenza a virus antigenic drift. Science, 326, 734–736.

10.

Jenkins

G.M.

, Rambaut

, Pybus

O.G.

, et al. 2002. Rates of molecular evolution in rna viruses: A quantitative phylogenetic analysis. J. Mol. Evol. 54, 156–165.

11.

Jolliffe

I.T.

2002. Principal Component Analysis. Springer-Verlag, New York, NY.

12.

Lam

H.C.

, Sreevatsan

, and Boley

2012. Analyzing influenza virus sequences using binary encoding approach. Sci. Program. 20, 3–13.

13.

Lee

C.-W.

, Senne

D.A.

, and Suarez

D.L.

2004. Effect of vaccine use in the evolution of mexican lineage h5n2 avian influenza virus. J. Virol. 78, 8372–8381.

14.

Nicholas

2007. Evolution, 1st edn. Cold Spring Harbor Laboratory Press, Long Island, NY.

15.

Russell

C.A.

, Jones

T.C.

, Barr

I.G.

, et al. 2008. Influenza vaccine strain selection and recent studies on the global migration of seasonal influenza viruses. Vaccine, 26, D31–D34.

16.

Schweiger

, Zadow

, and Heckler

2002. Antigenic drift and variability of influenza viruses. Med. Microbiol. Immunol. 191, 133–138.

17.

Suzuki

2006. Natural selection on the influenza virus genome. Mol. Biol. Evol. 23, 1902–1911.

18.

Taubenberger

J.K.

, and Kash

J.C.

2010. Influenza virus evolution, host adaptation, and pandemic formation. Cell Host Microbe, 7, 440–451.

19.

Webby

R.J.

, Perez

D.R.

, Coleman

J.S.

, et al. 2004. Responsiveness to a pandemic alert: Use of reverse genetics for rapid development of influenza vaccines. Lancet, 363, 1099–1103.

20.

Webster

R.G.

, Bean

W.J.

, Gorman

O.T.

, et al. 1992. Evolution and ecology of influenza a viruses. Microbiol. Rev. 56, 152–179.

21.

Wei

, Chen

, et al. 2012. Evolution and adaptation of hemagglutinin gene of human h5n1 influenza virus. Virus Genes, 44, 450–458.

22.

Wood

J.M.

, Nicholson

K.G.

, Zambon

, et al. 2001. Developing vaccines against potential pandemic influenza viruses, Vol. 1219, 751–759. In International Congress Series. Elsevier, Amsterdam, The Netherlands.