Genetic diversity of influenza A viruses circulating in Bulgaria during the 2018–2019 winter season

Introduction Influenza viruses evolve rapidly and change their antigenic characteristics, necessitating biannual updates of flu vaccines. Aim The aim of this study was to characterize influenza viruses circulating in Bulgaria during the 2018/2019 season and to identify amino acid substitutions in them that might impact vaccine effectiveness. Methodology Typing/subtyping of influenza viruses were performed using real-time Reverse Transcription-PCR (RT-PCR) and results of phylogenetic and amino acid sequence analyses of influenza strains are presented. Results A(H1N1)pdm09 (66 %) predominated over A(H3N2) (34 %) viruses, with undetected circulation of B viruses in the 2018/2019 season. All A(H1N1)pdm09 viruses studied fell into the recently designated 6B.1A subclade with over 50 % falling in four subgroups: 6B.1A2, 6B.1A5, 6B.1A6 and 6B.1A7. Analysed A(H3N2) viruses belonged to subclades 3C.2a1b and 3C.2a2. Amino acid sequence analysis of 36 A(H1N1)pdm09 isolates revealed the presence of six–ten substitutions in haemagglutinin (HA), compared to the A/Michigan/45/2015 vaccine virus, three of which occurred in antigenic sites Sa and Cb, together with four–nine changes at positions in neuraminidase (NA), and a number of substitutions in internal proteins. HA1 D222N substitution, associated with increased virulence, was identified in two A(H1N1)pdm09 viruses. Despite the presence of several amino acid substitutions, A(H1N1)pdm09 viruses remained antigenically similar to the vaccine virus. The 28 A(H3N2) viruses characterized carried substitutions in HA, including some in antigenic sites A, B, C and E, in NA and internal protein sequences. Conclusion The results of this study showed the genetic diversity of circulating influenza viruses and the need for continuous antigenic and molecular surveillance.


InTRoduCTIon
Influenza viruses cause highly contagious infection, that occurs as annual epidemics affecting approximately 5-15 % of the human population, resulting in 3-5 million cases of severe illness and 290 000-650 000 deaths globally [1]. Since 2009, two subtypes of influenza A viruses, A(H1N1) pdm09 and A(H3N2), together with two genetic lineages of B viruses, B/Victoria and B/Yamagata, have been responsible for seasonal epidemics. Influenza A viruses cause moderate to severe illnesses, affect all age groups, have a large animal reservoir and undergo faster evolution than type B viruses. Type B viruses usually cause milder illnesses, most commonly affecting children and, in contrast to type A viruses, have neither an animal reservoir nor pandemic potential.
The genomes of influenza A and B viruses consist of eight segments of single-stranded RNA with negative polarity, each of them encoding at least one protein with specific function(s). The surface glycoproteins -haemagglutinin (HA) and neuraminidase (NA) -elicit protective antibodies and are major components of influenza vaccines. Under host immune pressure, amino acid substitutions accumulate in these glycoproteins, leading to escape from natural or vaccine-induced immune responses and to a decrease in vaccine effectiveness. This phenomenon, known as antigenic drift, necessitates twice-yearly updating of influenza vaccine composition in order to correspond to circulating strains. Amino acid substitutions underlying antigenic drift are mainly located in the so-called antigenic (antibody-recognized) sites of HA. Five distinct antigenic sites have been identified in the HAs of both A(H1N1) (Ca1, Ca2, Cb, Sa, Sb) and A(H3N2) (A-E) viruses [2][3][4]. A relatively small number of amino acid substitutions in sites surrounding the receptor binding site (RBS) may result in altered antigenic phenotype [5][6][7], but other substitutions also have significant impact on antigenicity [8]. Notably, substitutions that create potential glycosylation sequons, NX S/T where X can be any amino acid but P, in HA and NA is another mechanism that enables the virus to evade the immune system. It is suggested that the added hostderived glycans alter antigenicity of HA by masking antigenic epitopes and preventing antibody recognition [9,10].
Laboratories of the World Health Organization (WHO) Global Influenza Surveillance and Response System (GISRS) perform continuous global surveillance of circulating influenza viruses with tracking of their evolutionary dynamics in order to optimize the selection of vaccine viruses. Analysis of the viruses identified allows early detection of novel genetic variants with altered antigenicity, epidemic potential, increased virulence or reduced susceptibility to antivirals. Based on this global surveillance, WHO makes recommendations twice a year, for countries in the Northern and Southern Hemispheres, respectively, for the composition of influenza vaccines in the forthcoming influenza season. The main aims of this study were to investigate the circulation pattern of influenza viruses causing the 2018/2019 season in Bulgaria and to determine their genetic and phenotypic characteristics for comparison with those of the vaccine and other epidemic viruses.

Study population and specimen collection
From October 2018 to May 2019, patients presenting with acute respiratory infection (ARI) from different regions of the country were enrolled in the national influenza surveillance program. ARI was defined in compliance with the European Centre for Disease Prevention and Control (ECDC) with sudden onset of symptoms, at least one of the following four respiratory symptoms: cough, sore throat, shortness of breath, coryza and a clinician's judgement that illness was due to an infection [11]. The diagnosis of each patient was determined by their attending physician based on standard clinical criteria. Nasopharyngeal specimens from the enrolled patients were collected using polyester collection swabs that were preserved in virus transport medium (Deltalab, Spain) during the visit to the doctor or within the first 24 h of admission. Swabs were stored at 4 °C for up to 72 h before shipment to the National Laboratory 'Influenza and ARI' for influenza virus detection and characterization. The laboratory is recognized as a WHO National Influenza Centre. It is the only laboratory in the country that conducts research on influenza viruses and performs year-round diagnostic testing of samples from outpatients and hospitalized patients to monitor the circulation of seasonal influenza viruses. Specimens were processed immediately or stored at −80 °C before testing.

Extraction of viral nucleic acids and real-time RT-PCR
Virus RNAs were extracted automatically from the respiratory specimens using a commercial ExiPrep Dx Viral DNA/ RNA kit and ExiPrep 16DX equipment (BioNeer, Korea) in accordance with the manufacturer's instructions. Detection and typing/subtyping of influenza viruses were carried out by a real-time RT-PCR method using a kit -SuperScript III Platinum One-Step Quantitative RT-PCR (qRT-PCR) System (Invitrogen, USA). All samples were first tested for the presence of influenza A and B viruses. Those that were positive for influenza A were subsequently screened for A(H1N1) pdm09 and A(H3N2). Primers, probes and positive controls were provided by the International Reagent Resource, USA. Amplification was performed with a Chromo 4 thermal cycler (Bio-Rad) in accordance with the protocol recommended by CDC, Atlanta, USA (reverse transcription at 50 °C for 30 min, Taq inhibitor inactivation at 95 °C for 2 min, followed by 45 cycles of denaturation at 95 °C for 15 s and annealing/amplification at 55 °C for 30 s). A Ct value <38 was regarded as positive.

Virus isolation and antigenic characterization
Real Time RT-PCR positive clinical specimens with Ct values <28 were inoculated onto cultures of Madin-Darby Canine Kidney (MDCK) and MDCK-SIAT1 (a cell line expressing increased levels of α2,6-sialyltransferase cells [12]). Cultures were incubated at 35 °C in a 5 % CO 2 atmosphere and observed daily for 7 days for evidence of cytopathology. The presence of virus in culture was confirmed by haemagglutination assay following standard protocols using a 1 % suspension of guinea pig red blood cells. Cell-culture isolates of viruses detected at different stages of the epidemics, in different regions of the country and in patients of different age groups were selected for antigenic and genetic characterization. Antigenic testing was performed by the haemagglutination inhibition (HI) assay, in accordance with the WHO Manual, using vaccine viruses/antigens and their corresponding antisera provided by the WHO Collaborating Centres (WHO-CCs) in London and Atlanta [13]. More detailed HI assays of representative Bulgarian influenza isolates with panels of reference viruses and antisera were performed at the WHO-CCs in London and Atlanta. Virus isolates were identified as antigenically related to the vaccine virus if they showed no more than fourfold reduced HI titre with antiserum raised against the vaccine virus, as compared to the homologous titre. A reduction of at least eightfold in the HI titres was considered a signal of antigenic difference.

Genetic characterization
HA and NA gene sequences of influenza viruses detected in Bulgaria during the 2018/2019 season were determined at the WHO-CC London. Whole-genome sequencing was carried out at WHO-CC, Atlanta. Sequences have been deposited in the EpiFlu database of the Global Initiative on Sharing All Influenza Data (GISAID) [14]. For phylogenetic analyses, sequences of study viruses, reference viruses with known genetic group identities and viruses representing different countries of Europe during the 2018/2019 season, were retrieved from the EpiFlu database of GISAID (Table S1, available in the online version of this article). Multiple alignments for HA and NA sequences were carried out using the muscle algorithm and maximum-likelihood phylogenetic trees were constructed using RaxML v8.2X [15], followed by annotation with amino acid substitutions defining nodes and individual virus gene products (indicated in parenthesis after virus names) using treesub [16]. Trees were visualized using FigTree [17] and highlighted using Adobe Illustrator CC 2015.3 [18]. Phylogenies relate to nucleotide alignments of coding sequences for HA gene products, with signal-peptide encoding and stop codons removed to give mature HA glycoprotein amino acid numbering.

deduced amino acid sequence analysis and prediction of N-glycosylation motifs
The amino acid sequences were generated by translating nucleotide sequences with the standard genetic code using mega (version 6.06) software [19]. The deduced amino acid sequences of the study strains were compared to those of vaccine strains and other reference strains to identify amino acid substitutions. The amino acid identity was calculated using FluSurver [20]. Putative N-glycosylation motifs in the HA and NA were predicted using the NetNGlyc 1.0 web Server to identify sequence motifs N-X-S/T (sequon), where X can be any amino acid except proline. Only threshold values of >0.5 were considered as potential glycosylation sites. The occupancy of glycosylation sites in the HAs of some recent H1 and H3 viruses has been determined [21,22].

Antiviral susceptibility surveillance
Screening of A(H1N1)pdm09 viruses for the presence of point mutations encoding NA H275Y amino acid substitution, known to confer oseltamivir resistance, was carried out using a real-time RT-PCR assay that allowed discrimination of a single nucleotide difference between oseltamivir sensitive and resistant viruses. Two TaqMan probes differing in position 823 of the NA gene were used simultaneously: the first probe contained a cytosine at position 823 and was labelled with VIC (H275), while the second probe contained thymine in the same position and was labelled with FAM (275Y). Primer/probe sequences and protocol were kindly provided by Public Health England, London. Reference influenza viruses A/Denmark/524/2009 (sensitive, H275) and A/ Denmark/528/2009 (resistant, 275Y), provided by WHO-CC, London, were used as positive controls. A phenotypic analysis (MUNANA test) of influenza virus susceptibility to neuraminidase inhibitors (oseltamivir and zanamivir) was performed at WHO-CC, London.

Statistics
Age and gender of patients, the clinical features of their illness and the incidence of each virus were compared using the Chi square or Fisher's exact tests for categorical variables. P-values of <0.05 were considered statistically significant.

RESuLTS
The 2018/2019 influenza epidemic lasted for 7 weeks (weeks 2-9/2019) and peaked during week 4/2019. The ARI incidence rate during the epidemic peak (247.92 per 10 000 population) was higher than the incidence rates during the epidemic peaks

Influenza virus detection
The study population consisted of 1696 patients presenting with ARI symptoms: 392 (23.1 %) of these were persons attending outpatient healthcare centrers; 1304 (76.9 %) were hospitalized patients, of which 71 were treated in intensive care units (ICU). The patients' ages ranged from 14 days to 89 years (average age 28.  (Fig. 4). Antigenic characterization of clade 3C.2a viruses using HI assay was difficult because the majority of these viruses did not agglutinate red blood cells. Ten Bulgarian subgroup 3C.2a1b and two subclade 3C.2a2 viruses analysed by HI with post-infection ferret antisera showed low reactivity with antisera raised against the egg-propagated vaccine virus A/Singapore/INFIMH-16-019/2016 (subclade 3C.2a1).

Amino acid sequence comparisons
Sequences of full-length HA and NA glycoproteins of Bulgarian influenza A(H1N1)pdm09 and A(H3N2) viruses (36 and 28, respectively) were compared to those of eggpropagated vaccine viruses to identify substitutions that might impact vaccine effectiveness.   All studied viruses that fell in subgroups of subclade 6B.1A carried HA1 S183P substitution. Bulgarian representatives of the individual subgroups and 'outlier' group possessed different amino acid substitutions defining these subgroups (Table 1) [23]. Of the amino acid substitutions, three were located in antigenic sites: S74R in site Cb; S164T and N129D in site Sa [24]. Nine conserved potential N-glycosylation motifs were encoded by the HA genes of the 36 Bulgarian A(H1N1)pdm09 viruses: at HA1 positions 10/11, 23, 87, 162, 276, 287 and HA2 positions 154 and 213. The substitution S164T modified the 162-164 sequon, NQS→NQT, within the Sa antigenic site.  Table 2. The substitutions L85I, S105N, I188T, N449D and T452I occur at amino acid positions associated with NA antigenic sites [25,26]. NA catalytic site residues (R118,  (Table 2). Bulgarian isolates contained amino acid substitutions in all internal proteins, most being in PB2. Some of these substitutions were detected in other Balkan countries but the possible functional importance of these substitutions is not known. As with the vast majority of A(H1N1)pdm09 viruses, M2 proteins carried S31N substitution associated with resistance to M2-ion channel blockers (amantadine and rimantadine) [28].   Internal proteins sequences of four Bulgarian A(H3N2) viruses were compared to those of A/Singapore/INFIMH-16-0019/2016 vaccine virus and A(H3N2) viruses detected in Albania (5) and Macedonia (9) ( Table 4). Amino acid substitutions were detected at six PB2, three PA, two NP, one each PB1, M2 and NS1 amino acid positions.  NAs of Bulgarian A(H3N2) viruses all carried I212V and N329S substitutions compared to A/Singapore/INFIMH-16-019/2016. The 17 HA subclade 3C.2a1b viruses and the nine HA subclade 3C.2a2 viruses possessed different amino acid substitutions defining these subgroups ( Table 1). The catalytic site and supporting framework residues within NA were highly conserved, as were seven of eight potential N-glycosylation sequons (at positions 61, 70, 86, 146, 200,  234, 245 and 367), two of which (146 and 367) locate close to the enzyme active site [29]. Substitutions T72N, N234D and N329S identified in some or all Bulgarian viruses resulted in the loss of N-glycosylation sequons.

A(H3N2
Internal proteins sequences of four Bulgarian A(H3N2) viruses were compared to those of A/Singapore/INFIMH-16-0019/2016 vaccine virus and A(H3N2) viruses detected in Albania (5) and Macedonia (9). Amino acid substitutions were detected at six PB2, three PA, two NP, one each PB1, M2 and NS1 amino acid positions.

Antiviral susceptibility surveillance
Of the Bulgarian viruses assessed for susceptibility to NA (oseltamivir, zanamivir) and PA (baloxavir marboxil) inhibitors, none showed evidence of reduced susceptibility. Testing was performed on 280 detected A(H1N1)pdm09 viruses by real-time RT-PCR with respect to the NA H275Y oseltamivir resistance substitution at the NRL in Bulgaria; NA and PA sequences of all study A(H1N1)pdm09 and A(H3N2) viruses were screened for known markers of reduced susceptibility to NA inhibitors and baloxavir marboxil [30,31] and phenotypic testing (by the MUNANA method) was performed on 13 influenza A(H1N1)pdm09 and 5 A(H3N2) influenza viruses.

dISCuSSIon
The 2018/2019 influenza epidemic in Bulgaria was characterized by moderate intensity but more specimens from severe and fatal influenza cases were examined compared to the previous three seasons. According to media reports, greater numbers of influenza deaths were also reported in other Balkan Peninsula countries (Romania, Croatia, Serbia and others) during the same period. The total number and percentage of influenza positive cases were slightly higher (34 vs 33 %) compared to the previous 2017/2018 season, however, among the hospitalized patients, the difference was greater (34.7 vs 31.9 %). Only influenza A(H1N1)pdm09 and A(H3N2) viruses, at a ratio of 2 : 1, were identified in contrast to the previous season when type B viruses (mainly B/Yamagata-lineage) predominated accounting for 80 % of all influenza cases [32]. A similar proportion amongst circulating influenza viruses was observed in most European countries. Cumulative data for the WHO European region showed strong prevalence of influenza A viruses in the 2018/2019 influenza season accounting for 99 % of all positive influenza cases with A(H1N1)pdm09 prevailing over A(H3N2) [33].
Generally, influenza A viruses account for the majority of influenza cases, prevail in most seasons, cause a higher percentage of serious illness and more frequently undergo genetic and antigenic changes in comparison to influenza B viruses [7]. In Bulgaria, influenza A was dominant in all but two ( [34]. During the 2018/2019 epidemic, influenza A(H1N1) pdm09 viruses predominated, were most prevalent among children aged 5-14 years and were responsible for a higher proportion of serious illnesses than A(H3N2) viruses, which most severely affected 15-29 year (16%) and 5-14 year (15%) age groups. Analysis revealed increased proportions of influenza A(H1N1)pdm09 virus infection among hospitalized patients, among those treated in ICUs and among fatal influenza cases compared to proportions related to A(H3N2) infection.
Since the emergence of influenza A(H1N1)pdm09 virus during the 2009/2010 pandemic, its HA gene has undergone significant genetic changes and evolved in eight genetic groups and several clades/subclades. The recent clade 6B.1 viruses emerged during the 2015/2016 season and caused many severe and fatal influenza cases in a number of countries, including Bulgaria [35][36][37]. The most recently circulating viruses formed a subclade designated 6B.1A and seven subgroups have been defined within it: 6B.1A1 -6B.1A7. A(H1N1)pdm09 viruses circulating in Bulgaria in 2018/2019, in contrast to previous seasons, were genetically diverse falling in four of these subgroups and a 6B.1A outlier group carrying several HA and NA amino acid substitutions. A total of 24 different substitutions (three located in antigenic sites Sa and Cb and one associated with egg-adaptation of the vaccine virus) were identified in HA and 18 in NA compared to vaccine virus, A/Michigan/45/2015. One to two amino acid changes in HA antigenic site Sa dramatically affect antibody recognition [38]. Gao et al. showed that N449D is one of three NA substitutions resulting in loss of reactivity with some human monoclonal antibodies [39]. Subgroup 6B.1A5 viruses prevailed globally and carried more substitutions in antigenic sites than other subgroups. Despite these changes in functionally important sites of HA and NA, subclade 6B.1A A(H1N1) pdm09 viruses remained antigenically similar to vaccine virus as assessed by use of post-infection ferret antisera raised against the A/Michigan/45/2016 vaccine virus [23]. However, recent subclade 6B.1A A(H1N1)pdm09 viruses with HA S183P substitution [58 % of the sequenced Bulgarian A(H1N1)pdm09 viruses] were antigenically distinguishable from the vaccine virus by panels of post-vaccination human sera, and this resulted in the choice of A/Brisbane/02/2018 as a vaccine component for the Northern Hemisphere 2019/2020 influenza season [40].
A(H1N1)pdm09 HA1 D222G/N/S or Q293H substitutions have been associated with cases of severe disease and fatalities earlier in the pandemic period [41]. D222G/N polymorphism has been reported in recently circulating viruses [42,43] and two of the recent Bulgarian viruses (A/Bulgaria/013/2019 and A/Bulgaria/043/2019) carried D222N substitution. The D222G substitution was shown to cause a shift from α2,6sialic acid receptor specificity to mixed α2,3/α2,6-sialic acid receptor specificity, adduced thereby to facilitate lung infection [44,45].
Among the seasonal influenza viruses, A(H3N2) undergo the fastest evolution and diversification [7,46]. Since the 2013/2014 Northern Hemisphere influenza season two genetic clades, 3C.2a and 3C.3a, have been in circulation. The majority of clade 3C.2a viruses showed poor or total loss of ability to agglutinate red blood cells, so it was difficult or impossible to assess the correlation between amino acid/ glycosylation changes and the antigenic properties of these viruses using HI assay resulting in the implementation of virus neutralization assays [23]. During the recent influenza epidemics, several subclades (3C.2a1, 3C.2a2, 3C.2a3 and 3C.2a4) and subgroups (3C.2a1b and 3C.2a1a) have been identified within clade 3C.2a with subgroup 3C.2a1b viruses being predominant in the 2018/2019 season [47]. The majority (67.9 %) of Bulgarian A(H3N2) isolates belonged to subgroup 3C.2a1b, the remainder (32.1%) to subclade 3C.2a2. Bulgarian subgroup 3C.2a1b viruses showed amino acid substitutions at 17 positions in HA (inclusive of three associated with egg-adaptation) and 15 in NA, compared to the 2018/2019 vaccine virus, while subclade 3C.2a2 viruses were genetically closer to the vaccine virus with substitutions at 13 HA (inclusive of three associated with egg-adaptation) and 7 NA positions. Of the A(H3N2) HA amino acid substitutions seven at positions in HA1 were located within antigenic sites A, B (T160K substitution associated with egg-adaptation), C and E, and these changes could potentially alter the antigenicity of the HA. Antigenic sites B and A, localized on the top of HA around the RBS are main targets of human neutralizing antibodies and single substitutions, which occurred at seven positions in these sites (145 in site A and 155, 156, 158, 159, 189, 193 in site B) have been responsible for major antigenic cluster transitions from 1968 to 2003 [5]. In our study no substitutions were found in any of these seven key positions. Attachment and loss of oligosaccharides in HA and NA can also alter the biological properties of the viruses and the glycosylation in the antigenic regions of HA has been reported to associate with increased virulence [48]. For example, the N-glycan at HA1 position 158 of recent A(H3N2) viruses is in the vicinity of the residues at positions 156, 158 and 159 that are responsible for major antigenic changes. Non-HA substitutions, e.g. V263I in NA and K196E in NS1, associated with increased clinical severity in some reports [49], were not identified in Bulgarian isolates. Overall, our studies show increased amino acid substitution in antigenically important sites and higher levels of glycosylation in H3-HA compared to H1-HA (Table 3) supporting other authors' observations that H3 viruses evolve more rapidly and display more diversity compared to H1 viruses, which necessitates more frequent updates for this component of influenza vaccines [7,50,51].
Analysis of internal proteins in A(H1N1)pdm09 and A(H3N2) viruses revealed the presence of the number of amino acid substitutions, mainly in the PB2 protein, but their biological significance and clinical relevance is unknown. Antiviral susceptibility testing of influenza viruses circulating in Bulgaria, carried out by means of genotypic and phenotypic methods, showed that all studied type A viruses were susceptible to the within country licensed NA inhibitors oseltamivir and zanamivir. These investigations are necessary to allow rapid detection of viruses with reduced susceptibility or resistance to available antiviral drugs. Global surveillance of influenza virus susceptibility to NA inhibitors has shown very low levels (<1 %) of antiviral resistance among recently circulating influenza viruses [52].

Conclusion
The 2018/2019 influenza season in Bulgaria was caused by genetically diverse A(H1N1)pdm09 and A(H3N2) viruses that co-circulated. The increased numbers of serious and fatal cases of influenza during the season could be explained by predominance of recently emerged subclade 6B.1A influenza A(H1N1)pdm09 viruses and low vaccination coverage (2-3 % of the population). Continuous monitoring of circulating influenza viruses is necessary in order to evaluate the burden of influenza infection and to adopt an effective national strategy for the control and prevention of influenza.