Co-circulation and persistence of multiple A/H3N2 influenza variants in China

ABSTRACT The spread of influenza A/H3N2 variants possessing the hemagglutinin 121 K mutation and the unexpectedly high incidence of influenza in the 2017–2018 northern hemisphere influenza season have raised serious concerns about the next pandemic. We summarized the national surveillance data of seasonal influenza in China and identified marked differences in influenza epidemics between northern and southern China, particularly the predominating subtype and the presence of an additional summer peak in southern China. Notably, a minor spring peak of influenza caused by a different virus subtype was also observed. We also revealed that the 3C.2a lineage was dominant from the summer of 2015 to the end of the 2015–2016 peak season in China, after which the 3C.2a2 lineage predominated despite the importation and co-circulation of the 121 K variants of 3C.2a1 and 3C.2a3 lineages at the global level. Finally, an analysis based on genetic distances revealed a delay in A/H3N2 vaccine strain update. Overall, our results highlight the complicated circulation pattern of seasonal influenza in China and the necessity for a timely vaccine strain update worldwide.


Introduction
Influenza activity differs between southern and northern China, with two peaks (winter and summer) in southern China and only a single peak (winter) in northern China [1]. However, between May and August of 2017, the abrupt increase in hospitalization fatality rates caused by seasonal influenza A/H3N2 in Hong Kong Special Administrative region of China received broad media coverage worldwide, leading to serious public health concerns [2]. Most striking was that the reported number of influenza cases in the 2017-2018 influenza season ranked second in China after the 2009 pandemic H1N1 [3].
Seasonal human A/H3N2 influenza virus is characterized by periodic and rapid antigenic variation [4][5][6] that may reduce vaccine effectiveness. Since the 2016-2017 influenza season in the northern hemisphere, a number of novel subclades have been detected during routine influenza surveillance worldwide [7,8]. In particular, a novel genetic variant with a N121 K mutation in the hemagglutinin (HA) protein emerged in Hong Kong and comprised more than 35% of the A/H3N2 viruses tested in May 2017 [2]. This mutation was located in epitope D, and variants possessing this mutation were reported to have a less vaccine effectiveness in Denmark in the 2016-2017 influenza season [9]. More importantly, a recent study suggested that most of the global population would be susceptible to the 121 K variants [10]. To date, this variant has been widely spread across Europe [7,9,11], Asia [2,8,12] and North America [13].
To provide a more comprehensive study of the genetic diversity and antigenic variation of A/H3N2, and to trace the origin and spread of the novel 121 K variant in China, we present a description of the national surveillance data of influenza-like illnesscombined with the sequence data of 1471 full-length genomes of A/ H3N2 influenza A virus isolated across China since 2015.

Sample collection
Throat or nasal swabs were collected from outpatients of the sentinel hospitals, with clinical evidence of influenza-like illness (ILI), defined as a person with sudden onset of fever ≥ 38°C and cough or sore throat. Samples collected at sentinel sites were sent to the network laboratories within 48 h. Real-time PT-PCR was performed to determine the influenza positive samples as well as the type and subtype of the influenza virus. The viruses were isolated with Madine-Darby canine kidney (MDCK) cells or embryonated chicken eggs by the network laboratories in a biosafety level 2 facility. Influenza positive isolates were sent to provincial CDC and the Chinese National Influenza Centre for additional laboratory tests.

RNA Extraction, RT-PCR and sequencing
Viral RNA was extracted from 200 μl of the infected allantoic fluid of embryonated chicken egg or cell culture supernatant with QIAamp viral RNA mini kit (Qiagen, Inc.). For each virus strain, the RNA was eluted in 50 µl nuclease-free water. Complementary DNA was synthesized by reverse transcription reaction, and gene amplification by PCR was performed, whole genome sequencing of influenza A virus [14] was performed on the MiSeq high-throughput sequencing platform (Illumina, Inc., San Diego, CA, USA).This resulted in the generation of full-length genome sequences of 1417 influenza A/H3N2 viruses isolated across China since 2015.

Phylogenetic analysis
Full-length genome sequences of human A/H3N2 influenza viruses isolated since 2015 were downloaded from the Influenza Virus Resource at the National Center for Biotechnology Information (NCBI) (http:// www.ncbi.nlm.nih.gov/genomes/FLU/FLU.html) and the GISAID database on 30 November 2017. Repetitive sequences in the two databases were removed by matching strain names using in-house scripts. Sequences of low quality were also removed. The remaining sequences were combined with those generated in the present study. Eight datasets corresponding to the eight gene segments of influenza A virus were first aligned using Muscle [15] and then adjusted manually in Bioedit [16]. Phylogenetic analysis of the aligned HA and neuraminidase (NA) datasets were performed using RAxML [17] under the GTRGAMMA nucleotide substitution model [18] and with random starting trees. To assess nodal support 1000 bootstrap replicates were performed and all other parameters were set to default. Trees were visualized using FigTree (http://tree.bio.ed.ac.uk/software/figtree/).
To investigate the potential origin of the Chinese 121 K variants, the HA genes of global 121 K variants were downloaded from the GISAID database. These dataset comprised 10,088 HA gene sequences, including 165 from China, and was aligned using Muscle and manually adjusted using Bioedit. We used multidimensional scaling (MDS)was to analyze the relationships among the sequences by computing a distance matrix among the sequence dataset. We employed MDS rather than phylogentic trees as a visualization tool because of the very large size of the data set. Indeed, MDS has been previously used to visualize genome evolution [19] as well as antigenic variation of influenza viruses [4]. To reduce the computational burden, we estimated the pairwise sequence for the complete dataset distance matrix using Phylip [20]: this was used the input in the MDS analysis to identify the major groupings among these sequences. Based on the MDS analysis, the complete dataset was then sub-divided into two smaller sub-datasets. We then repeated the phylogenetic analysis described above on these sub-datasets using the same parameters, but with 100 bootstrap replicates. The trees generated in this case were also visualized using FigTree.
To help assess vaccine effectiveness and the timing of vaccine strain update, pairwise HA1 protein distances between vaccine strains and the circulating strains were also calculated using Phylipand visualized using in-house scripts.

Antigenic variation analysis of H3 HAs
The crystal structure of A/H3N2 HA (A/Victoria/361/ 2011; PDB code:4O5N) was used to show the position of the A-E epitopes in H3 HA surface as antigenic sites. Compared with the HA of A/Victoria/361/2011, distinct residues which led to antigenic variation from six representative A/H3N2 virus strains were highlighted. The six HA sequences included were:
Based on the ILI surveillance data obtained from Chinese National Influenza Centre, we found there was a clear seasonality of A/H3N2 activity both in northern ( Figure 1A) and southern China ( Figure  1B). However, a single influenza peak typically occurring during winter months was observed in northern China, and the positive rate of A/H3N2 was high in the 2014-2015 and 2016-2017 influenza seasons ( Figure 1A). In the 2015-2016 influenza season, the B/Victoria lineage dominated over other subtypes in northern China. In particular, a second minor peak was observed in February following the earlier major A/H3N2 peak, as well as in the 2014-2015 and 2016-  In contrast to the situation in northern China, influenza epidemics in southern China tended to have one summer peak and one winter peak ( Figure  1B). The summer peak was clearly observed in 2015 and 2017 with A/H3N2 being the dominant subtype, but not observed in 2016, probably caused by the delayed end of the B/Victoria epidemics originating in early 2016. A/H3N2 also dominated over other subtypes in the 2016-2017 influenza season, which began in August, approximately eight weeks earlier than usual ( Figure 1B). In particular, like northern China, a minor spring peak was observed in southern China in 2016 and 2017 ( Figure 1B).
Southern China differs from northern China in a number of other aspects ( Figure 1A Figure 1C). In addition, these strains were isolated from 15 of 16 provinces (autonomous regions and municipalities) in northern China (n = 172) and 14 of 15 provinces (autonomous regions and municipalities) in southern China (n = 1239), as well as in the Hong Kong and Macao Special Administration Region (n = 6) ( Figure S1).
Consistent with the genetic analysis, homology modelling showed that the Chinese A/H3N2 viruses have evolved rapidly in the known epitopes during this time period (Figure 2). For example, compared to the vaccine strain A/HongKong/4801/2014( Figure  2D), A/Singapore/INFIMH160019/2016(representative strain of clade 3C.2a1) possessed the distinct N121 K and N171 K amino acid substitutions in epitope D, whereas A/Yunnan/Jinghong_1653/2017 (representative strain of clade 3C.2a3) obtained a T135Ksubstitution in epitope A and N121 K in epitope D ( Figure 2F), In addition, A/Guangdong/Dong-guan_F20161100/2016 (representative strain of clade 3C.2a2, Figure 2G) possessed a T131 K and aR142 K in epitope A. Further study on the effects of specific amino acid substitutions at these antigenic sites may have important implications for improved vaccine development.
Similarly, based on our phylogenetic analysis (Supplementary file 2) and two recent reports [21,22], we also determined the NA genotypes of these viruses (  (Table S1).
The HA1 sequence divergence between vaccine and circulating strains The major epitopes of seasonal influenza viruses are located in the HA1 protein [23]. To help assess vaccine effectiveness and the optimal timing for vaccine strain updates, we estimated the distances between the vaccine strains and the Chinese strains since 2015 using the full-length HA1 protein sequences (Figures 3 and  4). Genetic distances between strains were visualized by a colour gradient from blue to red, with blue representing distance of 0 and red representing distance of 0.1    (Figure 4). A delay of vaccine strain update was also observed in other regions in the northern hemisphere, such as North America ( Figure S3) and Europe ( Figure S4), regions in the southern hemisphere, such as Australia ( Figure S5), and even in the tropics, such as Southeast Asia ( Figure S6).

Persistence of A/H3N2 in southern China
Although a small number of viruses with the 121 K mutation have been previously isolated in China, they did not cluster with the lineage 3C.2a (data not shown). In our analysis, a total of 153 Chinese strains possessed the 121 K mutation: 6 in genotype 3C.2a, 102 in 3C.2a1, 45 in 3C.2a2 and one in 3C.2a3. The earliest strain was A/Guangdong-Zhanjiang/SY486/2015, which was isolated on June 1, 2015. Since then, only 15 121 K strains were identified in China before the 2017 off-season in our routine surveillance, whereas 124 121 K variants were reported in the 2017 off-season.
To investigate the origin of the Chinese 121 K variants, we performed an MDS analysis of the worldwide 121 K strains (n = 10,088). This indicated that they could be divided into two major groups ( Figure S7) that generally corresponded to the two named 121 K subclades, 3C.2a1 and 3C.2a3, respectively. We then extracted the sequences of the two groups (3C.2a1, n = 7456; 3C.2a3, n = 2624) and performed separate phylogenetic analyses using RAxML ( Figure 5 and supplementary files 3-4). This revealed that many Chinese 121 K variants were not located at the main trunk of the trees, but rather appeared as terminal branches (i.e. tips) ( Figure 5). This suggests that the majority of the Chinese 121 K variants were imported from external sources multiple times, such as Europe, Australia and other Asian countries ( Figure 5).

Discussion
The influenza peak in the summer of 2017 (May-August) in southern China, especially in Hong Kong,  In addition, we further reported that even in winter months the dominant influenza subtype between southern and northern China was not always the same. In particular, a minor spring influenza peak caused by a virus different from the dominant subtype in winter was found. Together, these results suggest a heterogeneous and complex influenza circulation pattern in China which might be shaped by such factors as herd immunity and climatic variation.
Compared with influenza A/H1N1 and B, numerous studies have revealed that A/H3N2 experienced more rapid antigenic variation and genetic variants of A/H3N2 were reseeded from the East Asia and Southeast Asia (E-SE) regional circulation network [24]. In recent years, several novel A/H3N2 variants have emerged globally, such as 3C.2a and its decedents 3C.2a1 and 3C.2a3 [7]. In the 2014-2015 influenza season, 3C.3a was found to be dominant in China. However, 3C.2a became the dominant genotype in the subsequent off-season of 2015 and the 2015-2016 influenza season. Although 3C.2a1 and 3C.2a3 have been reported to have increasing incidence and have become the dominant subtypes in some regions in 2017, such as Hong Kong [10], 3C.2a2predominated in mainland China, suggesting a diversifying circulation pattern in China. Furthermore, viruses of 3C.2a2were more closely related to A/HongKong/ 4801/2014 than the proposed vaccine strain by WHO in February 2018, A/Singapore/INFIMH160019/2016, which mainly targeted clade 3C.2a1 viruses. Therefore, the efficacy of this vaccine strain update in China warrants further investigation because the prevalence of the 121 K variants in China was not high during our surveillance period.
The amino acid substitutions in the epitopes of the HA protein of human influenza viruses would likely result in changes in antigenicity, helping the viruses evade human immune responses. Our analysis revealed that the clades 3C.2a1 (N121 K and N171 K), 3C.2a3 (N121 K and T135 K), and 3C.2a2 (T131 K and R142 K) possessed amino acid substitutions in the known epitopes A or D of A/H3N2, and experimental data from recent reports have revealed antigen variation among these clades [25]. In particular, the two clades possessing the 121 K substitution -3C.2a1 and 3C.2a3 -had reduced vaccine effectiveness [9]. Therefore, the proposed vaccine strain A/Singapore/ INFIMH160019/2016 from clade 3C.2a1 might have higher vaccine efficacy to the 121 K variants than A/ HongKong/4801/2014.
Attempting to predict the genetic and antigenic variation of human influenza has been one of the most important and controversial topics in influenza research [26] and several sequence-based machine learning methods have been proposed [6,27,28]. Such efforts could be potentially helpful for selection of appropriate candidate vaccine strain from circulating strains, which plays a pivotal role in the vaccine effectiveness. Unfortunately, our sequence-based analysis revealed a distinct delay (approximately one year later) of the A/H3N2 vaccine strain update in China, which has also been noted previously [6]. Similar delays have been observed in other regions, such as North America, Europe and Australia. Although our analysis was sequence-based and this might not reflect the real antigenic relationships between vaccine and circulating strains, considering the fact that the influenza vaccine effectiveness was not high in both China [3] and the USA [29], we believe that a novel and more efficient vaccine selection strategy is urgently required to account for the different influenza circulation patterns across continents and even major populous countries.
Apart from the delayed vaccine strain update, the egg-adaptive mutations occurring during the egg passage of the A/H3N2 vaccine seed strains would alter glycosylation and impair the neutralizing antibody response, resulting in reduced vaccine effectiveness [30]. For example, the L194P mutation on influenza HA was one of the most commonly observed egg-adaptive mutations and would significantly alter HA antigenicity [31]. A recent study reported that introduction of theG186 V mutation would only minimally alter HA antigenicity, but prevent the occurrence of L194P by means of mutationalincompatibility [31].
Influenza infection rates show clear seasonal variation and the peak of influenza season is concentrated in winter in temperate regions, usually between October and March in the northern hemisphere, the socalled peak season, whereas the remaining months are called often termed the off-season [32]. However, in the tropics, influenza infection remains at low levels throughout the year, with a peak often during the rainy season [33,34]. Despite the climatic variation, a recent analysis of the influenza surveillance data of Australia during 2007-2016 revealed distinct continental synchronicity of human influenza virus epidemics [35]. To date, a number of models have been proposed to account for the seasonality of seasonal influenza, such as global migration [5,36,37] and both global [38] and local persistence [39]. We found that the 121 K variant was seldom isolated in 2016, while most of the Chinese 121 K variants were observed in the summer of 2017 in southern China. Combined with the phylogenetic analysis, we suggest that these variants were imported from external sources to China multiple times. A previous study based on the ILI surveillance data during 15 years in Shenzhen, southern China also revealed multiple viral introductions [40]. Both studies suggest that China, even southern China, could also be a sink region in the global persistence model of human A/H3N2 influenza virus.
In sum, our national surveillance of ILI reveals the complex patters of circulation of different influenza subtypes and the co-circulation of multiple A/H3N2 3C.2a derived variants with minor antigenic variations in China. Not only do these data reveal the difficulties in vaccine strain selection, which clearly merit further investigation, but our analyses highlight the multiple introductions and circulation of the 121 K variants in this country as well as the worrying delay in vaccine strain updates both in China and globally.