Accelerated Evolution of H7N9 Subtype Influenza Virus under Vaccination Pressure.

No avian H7N9 outbreaks have occurred since the introduction of H7N9 inactivated vaccine in the fall of 2017. However, H7N9 is still prevalent in poultry. To surveil the prevalence, genetic characteristics, and antigenic changes of H7N9, over 7000 oropharyngeal and cloaca swab specimens were collected from live poultry markets and farms in 15 provinces of China from 2017 to 2019. A total of 85 influenza virus subtype H7N9 strains were isolated and 20 representative strains were selected for genetic analysis and antigenicity evaluation. Results indicated the decreased prevalence of low-pathogenic H7N9 strains while highly-pathogenic H7N9 strains became dominated since the introduction of vaccine. Phylogenetic analysis showed that strains from 2019 formed an independent small branch and were genetically distant to strains isolated in 2013-2018. Analysis of key amino acid sites showed that the virus strains may adapt to the host environment evolutionally through mutation. Our analysis predicted additional potential glycosylation sites for HA and NA genes in the 2019 strains. Sequence analysis of HA gene in strains isolated from 2018 to 2019 showed that there were an increased nucleotide substitution rate and an increased mutation rate in the first and second nucleotides of coding codons within the open reading frame. The hemagglutination inhibition (HI) assay showed that H7-Re1 and H7-Re2 exhibited a lower HI titer for isolates from 2019, while H7-Re3 and rLN79 showed a high HI titer. The protective effect of the vaccine decreased after 15 months of use. Overall, under vaccination pressure, the evolution of influenza virus subtype H7N9 has accelerated.

Introduction

Influenza virus subtype H7N9 is a zoonotic influenza virus with public health significance. Firstly detected in patients with pneumonia in China in March 2013, influenza virus subtype H7N9 initially appeared in live poultry markets, causing infections in humans in several provinces (Chen et al. 2013). There have been five H7N9 outbreaks since 2013 (Lam et al. 2013, 2015; Liu et al. 2013; Xie et al. 2015; Lu et al. 2018). They usually start in the autumn and winter each year and end in the summer of the following year (Gao 2018). According to the World Health Organization (WHO), there have been a total of 1,568 human cases of H7N9 infection, with 616 deaths. On February 19, 2017, the Chinese Center for Disease Control and Prevention has announced that H7N9 exhibits an insertion mutation at the hemagglutinin (HA) cleavage site, which is the characteristic of a highly pathogenic avian influenza virus (HPAIV) compared to subtypes isolated from previous epidemics (Quan et al. 2018). In order to effectively control the epidemic situation of H7N9 subtype influenza, China has carried out a comprehensive policy of H7N9 subtype influenza vaccine immunization for poultry since the autumn of 2017 (http://www.moa.gov.cn/govpublic/SYJ/201707/t20170711_5744436.htm).

H7N9-inactivated vaccine (H7-Re1) has been introduced in the poultry industry in the fall of 2017 following the emergence of H7N9 mutant strains with high pathogenicity in the fifth outbreak. In December 2018, H7-Re2 replaced H7-Re1 and vaccine rGD76 was put in use. In July 2020, H7-Re3 replaced H7-Re2, and rGD76 was replaced by rLN79. Within three years, the vaccine strain was replaced twice. Although vaccination programs have successfully controlled the spread of H7N9 between poultry and humans (Wu et al. 2019), the virus remains prevalent in poultry (Jia et al. 2019). Continuous monitoring is required to check whether vaccination has induced changes in the viral gene sequence and antigenicity (Lee et al. 2004). In this study, virus samples from live poultry markets and poultry farms in 15 provinces of China from 2017 to 2019 were isolated and identified, and 20 strains of influenza virus subtype H7N9 were selected for genomic sequence analysis and antigenic evaluation.

Materials and Methods

Virus Isolation

Oropharyngeal and cloaca swab specimens were collected from poultry in live poultry markets as well as sent by farmers from 2017 to 2019 (Zhang et al. 2020). Each sample was placed in 1 mL of phosphate-buffered saline (PBS) containing penicillin (2000 U/mL) and streptomycin (2000 U/mL) before being centrifuged at 10,000 ×g for 5 min. The resulting supernatant was inoculated into the allantoic cavity of 10-day-old specific pathogen-free embryonated chicken eggs and incubated for 48 to 72 h at 35 °C. All the experiments were conducted in a biosafety level 3 laboratory.

RNA Extraction, RT-PCR, and DNA Sequencing

RNA was extracted from the suspension of virus isolates using RNeasy Mini Kit (Qiagen) by following the protocol recommended by the manufacturer. Two-step RT-PCR was conducted using universal primers as previously reported by Hoffmann et al. (2001), and every gene segment was amplified under standard conditions. PCR products were purified using QIAamp Gel Extraction Kit (Qiagen). DNA sequences were analyzed using ABI 3730 DNA Analyzer (Applied Biosystems).

Phylogenetic Analysis

Seqman and Editseq (DNA Star Lasergene) were used for sequence splicing and open reading frame analyses. EasyCodeML, BioEdit and MEGA were used for sequence format conversion and key amino acid site analyses. DAMBE was used for sequence saturation detection. PhyloSuite and ModelGUI were used for sequence alignment and evolution model selection. MrBayes was used for phylogenetic tree drawing. FIG Tree was used for phylogenetic tree viewing and editing. BEAST and BEAUti were used to calculate the mutation rate of nucleotide sequences. Tracer was used for viewing mutation rate results. Online freeware NetNGlyc 1.0 Server was used to predict potential glycosylation sites (http://www.cbs.dtu.dk/services/NetNGlyc/). The reference sequences of influenza virus subtype H7N9 (human and avian) were downloaded from GenBank and Global Initiative of Sharing All Influenza Data (GISAID) (Supplementary Table S1).

Estimation of Nucleotide Substitution Rates

The HA gene sequences of all strains of H7N9 subtype (avian origin) from 2013 to 2018 were downloaded from NCBI and GISAID. After sequence comparison, unique sequences with complete length were chosen and analyzed with the sequences of 20 isolated strains in this paper (Supplementary Table S2). The sequences were divided into six outbreak groups (first wave, second wave, third wave, fourth wave, fifth wave and 2018–2019). Bayesian-Markov Chain Monte Carlo (MCMC) algorithm was used to estimate the nucleotide substitution rate of HA and the mutation rate of amino acids encoded in the open reading frame.

Hemagglutination Inhibition Assay

To understand antigenicity changes in the isolated strains from 2017 to 2019, the standard positive sera of H7N9 vaccine strains (Re-1, Re-2, Re3, rLN79) and 20 strains of viral antigens were tested for hemagglutination inhibition (HI). Both antigen and antiserum were purchased from the Harbin Weike Biotechnology Development Company. Hemagglutination inhibition test was performed with standard positive serum and 4 HA units of H7N9 viruses were prepared according to the hemagglutination test results of 20 isolates. The highest serum dilution that produced complete inhibition of HA activity was regarded as HI endpoint titers.

Results

Epidemiological Investigation

Over 7000 samples were collected from 15 provinces of China, including Guangdong, Guangxi, Hebei, Shandong, Liaoning, Shaanxi, Hunan, Hubei, Sichuan, Jiangxi, Yunnan, Fujian, Henan, Chongqing and Jilin provinces. A total of 85 influenza virus subtype H7N9 strains were identified. Among them, 70 strains were isolated in 2017, 7 were isolated in 2018, and 8 were isolated in 2019. The isolation rate of goose-, duck- and chicken-derived strains was 28.2%, 18.8%, and 53%. The sequences of those isolates were deposited in GISAID.

Nucleotide and Amino Acid Sequence Analysis

A total of 20 influenza virus subtype H7N9 strains which have wide coverage and strong representativeness, were selected for further analysis based on isolation time, location, and source host (Supplementary Table S2). There were 16 strains with 4 inserted amino acids near the cleavage site of HA protein (PEVPKGKRTAR↓GLF and PEVPKRKRTAR↓GLF), which represent the molecular characteristics of HPAIV. Among the 20 strains, whilst the isolates from 2017 were of LPAIV, those from 2018 to 2019 all showed HPAIV characteristics (Ke et al. 2017; Zhu et al. 2017). The occurrence of T160A, G186V and Q226L amino acid mutations of the HA protein enhance the virus binding ability to sialic acid α-2, 6 galactose receptors (Dortmans et al. 2013). Based on our analysis, T160A, G186V, G186A and Q226L mutations occurred in 13, 19, 1 and 3 strain(s), respectively. In the E128 strain, three mutations (T160A, G186A and Q226L) occurred concurrently. Analysis on NA protein receptor binding site showed that R292K mutation occurred in 2 HPAIV strains isolated in 2019, suggesting that the drug-resistant mutation of this strain may pose a greater potential threat than that of other strains (Table 1). Similar to most early H7N9 strains, the absence of amino acids in the neck of NA protein may induce the virus adaptation to different host environments (Desmet et al. 2013). The absence of amino acids in the neck of NA protein at positions 69–73 in all the 20 isolates may indicate that strains exhibit an enhanced ability to infect and replicate in mammals (Castrucci and Kawaoka 1993).

Table 1

Statistics of key sites found in the 20 H7N9 strains

Gene	Site	Effect	Number
HA (H3 number)	PEIPKG↓GLF	Multiple consecutive basic amino acids are a typical feature of highly pathogenic avian influenza	3
	PEVPKG↓GLF		1
	PEVPKRKRTAR↓GLF		14
	PEVPKGKRTAR↓GLF		2
	T160A	Enhance the ability of the virus to bind to SAα2-6 Gal receptor	13
	G186V		19
	G186A		1
	Q226L		3
NA (N2 number)	H276Y	Resistance to neuraminic acid inhibitors	20
	R292K		2
	Stalk deletion (69–73)	Increase replication and virulence in mammals	20
M2	S31N	Resistance to M2 ion channel blockers	20
PB2	A588V	Increased adaptability and virulence to mice	6
	E627K	Enhance virus replication	20
	D701N	Enhance adaptability in mammals	20

Prediction analysis of potential glycosylation sites in HA protein showed that, except for E656 strain, there were four potential glycosylation sites (30NGT32, 46NAT48, 249NDT251, 493NNT495/497NNT499) in all isolates from 2017 to 2018, which were similar to those from the previous four outbreaks. However, there were five potential glycosylation sites (30NGT32, 46NAT48, 141NGT143, 249NDT251 and 497NNT499), with an additional potential glycosylation site of 141NGT143 in the 2019 isolates (Supplementary Table S3). The loss of glycosylation site at position 150 of HA protein may change the virus’s ability in recognizing mammalian cell receptors (Gao et al. 2013). Prediction analysis of NA protein potential glycosylation sites revealed that, except for E656 and F692 strains, there were seven potential glycosylation sites (42NCS44, 52NTS54 (or 52NPS54, 52NIS54), 63NET65, 66NIT68, 82NLT84, 142NGT144 /142NET144 and 197NAS199) in all isolates from 2017 to 2018. Notably, with an additional potential glycosylation site of 354NNT356, there were eight potential glycosylation sites in 2019 isolates (Supplementary Table S4). NA proteins can hydrolyze receptor-specific glycoproteins on the surface of host cells to facilitate the release of virus particles. Changes in potential glycosylation sites suggest that the 2019 strains may have adapted to the host environment.

Phylogenetic Analyses of Surface Genes

Taking A/Anhui/1/2013 as the root, the H7N9 strains, based on genetic variations in their HA and NA genes, are divided into the Yangtze River Delta branch and Pearl River Delta branch. According to the phylogenetic tree, it can be inferred that the strains isolated in 2018–2019 have evolved from the previously prevalent 2017 strains. Notably, both HA and NA gene of the 2019 strains form an independent small branch (Fig. 1, Fig. 2), indicating that the 2019 strains have evolved and become genetically distant from strains isolated in 2013–2018. In addition, all the isolates except E65 located in the Yangtze River Delta branch form an independent branch with some Guangdong avian influenza virus strains, same results revealed in HA and NA phylogenetic trees. Since the isolates were collected from different provinces (Fig. 1), a cross-regional genetic exchange may have occurred.

Fig. 1

Phylogenic tree of H7N9 influenza viruses using HA gene sequences. The tree was rooted by using A/Shanghai/1/2013(H7N9), which was collected in February 2013.

Fig. 2

Phylogenic tree of H7N9 strains using NA gene sequences. The tree was rooted by using A/Shanghai/1/2013(H7N9), which was collected in February 2013.

Nucleotide Substitution and Amino Acid Codon Mutation Rates of HA Gene

The HA nucleotide substitution rate of isolates from the first and second waves was unchanged. Whilst those from the third and fourth wave showed a decline; and those from the fifth wave showed a slight increase. Strikingly, the HA nucleotide substitution rate of 2018–2019 strains increased 10 times than that of isolates from the fifth wave, indicating an acceleration of HA nucleotide substitution in influenza virus subtype H7N9 (Table 2).

Table 2

HA nucleotide substitution rate

	The nucleotide substitution rate of the HA genes (site/year)
	First wave	Second wave	Third wave	Forth wave	Fifth wave
Time	-2013	2013 ~ 2014	2014 ~ 2015	2015–2016	2016–2017	2018–2019
Rate	5.755E−3	5.734E−3	3.807E−3	1.821E−3	1.917E−3	1.963E−2
95%	3.85–7.77 (E−3)	3.98–7.71 (E−3)	1.72–6.28 (E−3)	5.27E−5–3.67E−3	2.92E−4–3.60E−3	8.65E−3–3E−2

The mutation rate of the first nucleotide within a codon in the entire HA open reading frame of the first wave isolates was found to be increased by nearly 2 times compared to that of the second wave isolates. The mutation rate then remained roughly unchanged, with a slight increase or decrease in isolates from the third, fourth and fifth waves. Meanwhile, the mutation rate of the second nucleotide within a codon in the entire HA open reading frame was observed to be increased by 2 times in the second wave isolates compared to that in the first wave isolates. The mutation rate remained the same in the third wave isolates and subsequently decreased in isolates from the fourth and fifth waves. For the third nucleotide within a codon in the entire HA open reading frame, the mutation rate was found to be continuously decreased in isolates from the first, second and third waves; but slightly increased in isolates from the fourth wave, and then remained the same in isolates from the fifth wave (Table 3).

Table 3

Codon nucleotide mutation rate

	Codon nucleotide mutation rate
	First wave	Second wave	Third wave	Forth wave	Fifth wave
	-2013	2013–2014	2014–2015	2015–2016	20,162,017	2018 ~ 2019
First nucleotide	0.354	0.614	0.71	0.603	0.695	0.939
Second nucleotide	0.248	0.575	0.574	0.413	0.379	0.518
Third nucleotide	2.398	1.811	1.716	1.984	1.926	1.544

In the reaction of H7-Re1 standard positive sera with the virus antigens of isolates, the HI titer of 2017 isolates was observed to be 4log2 ~ 8log2. Whilst isolate F11 showed the lowest HI titer of 4log2; the remaining 2017 isolates had HI titers of 5log2 ~ 8log2. The HI titers of 2018 isolates were found to be 4log2 ~ 5log2; while those of 2019 isolates were observed to be 3log2 ~ 5log2. Meanwhile, H7-Re1 standard antigen and standard positive serum showed a HI titer of 9log2. The HI titer of 2018 ~ 2019 isolates was observed to be 6log2 ~ 9log2 with H7-Re2, and 2019 isolates showed a significant antigenic difference. H7-Re3 and rLN79 had good antigenic matching with the epidemic strains in 2019 (Table 4).

Table 4

Hemagglutination inhibition endpoint titers of 20 influenza virus subtype H7N9 strains

Antigen	Antiserum (log2)
	H7-Re1	H7-Re2	H7-Re3	rLN79
H7-Re1	9	12	9	10
H7-Re2	10	11	8	8
H7-Re3	3	6	10	10
rLN79	4	6	10	10
E65(2017)	7	8	9	9
E112(2017)	8	8	9	9
E122(2017)	5	6	8	9
E128(2017)	8	8	9	10
E194(2017)	5	9	8	10
E282(2017)	6	6	9	10
E594(2017)	5	6	8	10
E656(2017)	5	6	10	9
E664(2017)	6	8	9	9
E743(2017)	6	8	9	10
F11(2017)	4	9	9	10
F450(2018)	4	8	9	9
F690(2018)	5	9	10	10
F691(2018)	5	9	10	9
F692(2018)	5	9	10	10
G285(2019)	4	9	10	10
G323(2019)	5	6	9	10
G375(2019)	4	6	9	9
G502(2019)	3	6	10	10
G503(2019)	3	6	9	10

Discussion

Influenza virus subtype H7N9 infection in humans was first reported in China in March 2013. The virus had caused five outbreaks subsequently. The spread of H7N9 from the early epidemic area of Yangtze River Delta to Pearl River Delta has been reported to be accompanied with the genetic evolution of the virus (Yang et al. 2017). In the early stage during the fifth outbreak, highly pathogenic H7N9 mutant strains that originally derived from low pathogenic strains in poultry have been found in samples collected from the Guangdong live poultry market, causing the epidemic to spread further (Qi et al. 2018a). Since the implementation of comprehensive H7N9 immunization prevention and control policy for poultry by the Chinese government in September of 2017, the spread of H7N9 infection between poultries and humans has been prevented successfully, as indicated by the absence of new outbreaks (Shi et al. 2018). According to published data, there have been a total of 23 H7N9 strains being isolated from live poultry markets and poultry farms nationwide in 2018–2020, 21 of which have been found to exhibit highly pathogenic molecular characteristics (Shi et al. 2018; Jia et al. 2019; Jiang et al. 2020). This indicates that under the influence of vaccine, the prevalence of low pathogenic H7N9 strains has been gradually decreased, and subsequently substituted by high pathogenic H7N9 strains.

Strains of the Yangtze River Delta branch are derived from the first wave of popular strains in Jiangsu, Shanghai, Zhejiang and other regions; while strains of the Pearl River Delta branch are mainly evolved from the early strains of the Guangdong Province (Qi et al. 2018b). One of our 20 isolates from 2017 to 2019 has HA and NA genes being attributed to the Yangtze River Delta branch and the Pearl River Delta branch, respectively. At the same time, isolates that located in the Yangtze River Delta branch form an independent branch with Guangdong poultry strains, suggesting that the virus strains may have cross-regional gene exchange. Although the 2018–2019 strains have a closer genetic relationship with the 2017 strains, the 2019 strains have an independent small branch, suggesting that the 2019 strains are more genetically distant from strains isolated from previous years.

The insertion of four consecutive basic amino acids near the cleavage site of HA protein is a molecular feature of the highly pathogenic (HP) influenza virus subtype H7N9. Whilst 16 strains exhibiting HP molecular characteristics with HA protein cleavage sites of PEVPKGKRTAR↓GLF and PEVPKRKRTAR↓GLF were 2018–2019 isolates; 4 strains with low pathogenic molecular characteristics were 2017 isolates. This indicates that the gradually decreased LPAIV prevalence has been dominated by HPAIV.

Double or triple mutation of HA protein receptor binding sites (T160A, G186V and Q226L) has been shown to affect virus-host specificity by enhancing the ability of the virus to bind sialic acid α-2,6 galactose receptor (Rogers and Paulson 1983; Gog et al. 2007; Dortmans et al. 2013). The key factor for the decreased binding affinity of the new influenza virus subtype H7N9 for sialic acid α-2,3 in 2013 is attributed to the deletion of the 150th glycosylation site of HA protein (Gao et al. 2013). H276Y and R292K mutations in NA protein can induce resistance to neuraminidase inhibitory drugs (such as oseltamivir phosphate and peramivir). Adaptation to different host environments, as well as enhanced infection and replication ability of the virus in mammalian hosts, can be achieved through alteration of glycosylation sites and deletion of amino acids in the neck of NA protein. Three isolates have T160A, G186V(A) and Q226L mutations simultaneously, while 13 strains have T160A and G186V mutations at the same time. Meanwhile, two HPAIV strains have R292K mutation, suggesting that the strains may be resistant with a greater potential threat compared to other strains. All 20 isolates have an amino acid deletion in the neck of NA protein at positions 69–73. Isolates from 2019 have an additional potential glycosylation site respectively in both HA and NA genes, suggesting that the virus has evolved in the host environment after two years’ vaccination. It is conducive for influenza virus subtype H7N9 to adapt to the host environment during evolution.

The nucleotide substitution rate of 2018–2019 strains was found to be 10 times higher than that of isolates from the fifth wave, indicating that the highly pathogenic influenza virus has been mutating rapidly. The generation of antigenic variants under vaccination-mediated immunity has been well investigated in human influenza viruses by a mouse model and H5 HPAIVs by a chicken model (Hensley et al. 2009; Nguyen et al. 2017). Therefore, the mutation of the virus under vaccine pressure is a reasonable speculation. Notably, in 2018–2019 strains, the mutation rates of the first and second nucleotides in codons were found to be greatly increased compared to those in the fifth wave strains. Since the incorporation of different amino acids is mainly determined by the first or second nucleotide within a codon, our results are consistent with the observation that, potential glycosylation sites are found in HA and NA genes of the 2019 isolates as a result of amino acid mutation. Various research findings have indicated the acceleration of H7N9 evolution and mutation since 2017, so prevention and control of the virus should be strengthened (Zhang et al. 2020).

The antigenicity of the H7N9 subtype influenza virus has changed, and the HI titer of Re-1 against the 2017–2019 strains is decreasing year by year, suggesting that the protective effect of the Re-1 vaccine is declining; the HI titer of rLN79 and Re-3 against 2019 isolates maintain at a high level. However, it is unknown whether the new vaccine strain will induce further changes in the virus antigenicity with a greater threat to poultry and humans. It is also unclear whether the protective effect of the H7-Re3 and rLN79 will be decreased over time, so continuous monitoring is required.

In summary, the highly pathogenic H7N9 strains are mainly prevalent at present. H7N9 has been undergoing intense evolution since the application of the vaccine. Phylogenetic analysis showed that the 2019 strains form an independent small branch, suggesting that the HA and NA genes differed greatly compared to those of isolates from previous years. An additional potential glycosylation site is found in the HA gene, as well as in the NA gene of the 2019 isolates. From 2018 to 2019, the nucleotide substitution rate of HA gene was found to be increased, and the mutation rate of the first and second nucleotides of the coding codon was observed to be increased. Lee’s research showed that commercial vaccine was not able to prevent virus shedding when chickens were challenged with antigenically different isolates, the virus evolution has accelerated after the use of the vaccine (Lee et al. 2004), therefore stringent measures should be implemented to prevent, control and monitor the spread of the virus. The antigenicity of the isolates in this study has changed, as indicated by the decreased HI titer of Re-1, and an average HI titer of Re-2 exceeding 6log2. The protective effect of the Re-1 vaccine was observed to be decreased after 15 months of use, suggesting that close monitoring of the Re-2 vaccine efficacy is required. Hence, urgent attention and diligent surveillance of avian influenza virus subtype H7N9 are becoming increasingly essential.

Below is the link to the electronic supplementary material.

Supplementary file1 (PDF 63 KB)

Supplementary file2 (XLSX 17 KB)