The HA_Ind protein sequences, reported from various geological locations of India during 2009–2018, are deposited in the NCBI (Table 1), and the same were retrieved for the present study (https://www.ncbi.nlm.nih.gov/genomes/FLU/Database/nph-select.cgi?go=database).

To understand the mutational and evolutionary drift among the HA_Ind protein sequences that circulated during the aforementioned period, MSA was carried out using ClustalW (Thompson et al., 1994). The evolutionarily closest sequence has been identified using the pair-wise distance matrix method. The resulting MSA was used to find the evolutionarily conserved regions in the examined sequences.

Following MSA, phylogenetic analysis was performed on HA_Ind proteins reported during 2009–2018. In order to understand the evolutionary relationship of these HA_Ind proteins of H1N1 strains from India, along with the reported pandemic Californian strain (A/California/04/2009), a phylogenetic tree was constructed using the maximum parsimony method. Parsimony analysis was performed in PAUP (v.4.0) using a heuristic search approach along with the following settings: 1) characters unordered with equal weight, 2) random taxon addition, and 3) branch swapping with the tree bisection–reconnection (TBR) algorithm. Resampling was performed with 1,000 replicates by bootstrapping to check the reliability of the results (Felsenstein, 1985; Tamura et al., 2004; Victoria Martínez et al., 2008). The selected HA_Ind sequences revealed a close relationship with the pandemic HA_Cal protein of the pandemic strain (A/California/04/2009) reported in 2009. Hence, the HA_Ind sequences, evolved as the closest members to the HA_Cal protein, were clustered as one clade and used for further studies (Table 2).

Mutational analysis was carried out using the ClustalW alignment tool of BioEdit (version 7.2.5) (Hall, 1999) with a bootstrap value of 1,000 to generate a global alignment for the selected HA_Ind proteins compared with the HA_Cal protein (Table 2) to investigate whether there is any prevalence of phenotypic variation in the reported HA_Ind protein sequences. The algorithm computed a distance matrix between each pair of sequences based on pairwise sequence alignment scores.

In addition to sequence comparison, the effects of mutations on the structure of HA were also investigated. For the HA structure comparison analysis, the following were selected: 1) one of the Indian HA proteins (Acc No: QCP70896), reported in 2018 (will be referred to as HA_Ind-2018, hereafter) and 2) the reference HA_Cal protein.

The complete crystallographic structure of the reference HA_Cal protein (with 566 amino acids) is not reported in the PDB website, and the reported HA_Cal structure (PDB ID: 3LZG) has only 506 amino acids. Hence, the complete HA_Cal structure was predicted using the ab initio modeling strategy implemented in Robetta (Raman et al., 2009; Song et al., 2013). The sequences of HA_Cal and HA_Ind-2018 (GQ280797 and QCP70896, respectively, in FASTA format) were retrieved from the NCBI databank (http://www.ncbi.nlm.nih.gov/) and utilized for homology modeling. The superimposition of both crystal and modeled HA_Ind-2018 and HA_Cal structures is shown in Supplementary Figure S5.1 and Supplementary Figure S5.2. The modeled HA protein structures were analyzed by WHAT IF and SAVES (http://nihserver.mbi.ucla.edu/SAVES/) servers and were visualized by UCSF Chimera (v.1.15) software (Pettersen et al., 2004). The robustness of the generated models was ensured by the RAMPAGE server (Ramachandran map).

Receptor-binding site analysis was performed on the HA_Ind protein to examine the mutation-mediated variation that emerged at the binding sites when compared to the HA_Cal protein. Wei Hu et al. reported the highly conserved receptor-binding domains of the HA protein of the influenza A (H1N1) virus (Hu, 2010), and their information was used while characterizing both HA_Ind and HA_Cal strains.

The analysis of the EBS gains importance as it provides the hotspot for membrane fusion between the host and pathogens. Analysis of the conserved EBS is essential to understand its dominance over the recognition of the antibody. Apart from the canonical/native epitope sites, the identification of mutation-derived new epitope sites is also essential to explain the exact viral–host interaction during the membrane fusion mechanism. Epitope sites of both HA_Ind and HA_Cal proteins of the influenza A (H1N1) virus were analyzed in the SVMTriP web server (http://sysbio.unl.edu/SVMTriP) using default parameters to explain both the conserved EBS and the newly emerging EBS due to mutation. The potential antigenic sites in the HA sequence were examined using a string kernel-based support vector machine (SVM), SVMTriP (Yao et al., 2012). This SVM model calculates the similarity using the BLOSUM62 matrix for the tripeptides or trimers from the input sequences given in FASTA format. Finally, the predicted epitopes, within the default limit of 20 sites, are ranked according to their scores.

One of the most influential post-translational modifications is N-glycosylation, which affects antigenicity, biological activity, cell–cell interactions, protein solubility, protein folding, localization, and trafficking. Here, the N-glycosylation sites across functional domains of the HA protein are mapped to locate both known and mutation-derived new sites as well. The NetNGlyc 1.0 server (https://services.healthtech.dtu.dk/service.php?NetNGlyc-1.0) was used with default parameters to analyze the N-glycosylation sites that are conserved among the Indian isolates of influenza A (H1N1) viruses. The NetNGlyc 1.0 server predicts all possible sequence patterns, “N-X-S/T” (any amino acids except P at the X position) within HA protein sequences as potential N-glycosylation sites, based on an artificial neural network approach. The most probable N-X-S/T patterns with the highest percentage of occurrence are filled out using the cutoff value of 0.5. The locations of the predicted N-glycosylation sites in the monomer of the HA_Ind protein are numbered according to the full-length HA_Cal protein sequence.

The ProtParam (https://web.expasy.org/protparam/) tool implemented in the ExPASy server is capable of predicting various physicochemical properties from the sequence, such as molecular weight, pI, amino acid composition, atomic composition, extinction coefficient, estimated half-life, instability index, aliphatic index, and grand average of hydropathicity (GRAVY) from the sequence. Here, variation in the amino acid composition of HA_Ind proteins reported during 2009–2018 is analyzed in comparison with the HA_Cal protein using ProtParam with default parameters to understand the genetic susceptibility and evolution of Indian isolates compared with the pandemic strain (A/California/04/2009).

A high degree of conformational plasticity may present a barrier to the development of beneficial antibodies. The GOR IV web server (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?page=npsa_gor4.html) was used with default parameters to 1) understand the degree of conformational plasticity by analyzing the secondary structure (alpha helix, extended strand, and random coil) and also 2) further illustrate the variable and invariable structural changes in HA_Ind proteins reported during 2009–2018 compared with the selected HA_Cal protein.

Electrostatic interactions (EIs) play a vital role in determining biomolecular functions. In particular, the EIs, which govern biomolecular sensing, are highly regulated by the nature of electrostatic potential. Hence, analysis of effective biomolecular sensing requires a thorough characterization of the distribution of ESP over the biomolecular surface boundaries. Here, the electrostatic charge distribution over the surface of both HA_Ind and HA_Cal proteins is calculated at an ionic strength of 0.15 M and visualized/analyzed using the Adaptive Poisson–Boltzmann Solver (APBS) plugin, integrated in VMD software (version 1.9.3). The ESP, represented as an isoelectrostatic potential map, depicts red and blue contours with values of −5.0 and +5.0 kBT/e⁻, respectively.

About 512 HA_Ind protein sequences, reported during 2009–2018 from various geographical locations in India until October 2020, are available in the NCBI repository database under the subsection “Influenza virus.” We collected the available HA_Ind protein sequences in FASTA format (Supplementary File S1) and performed the BLAST search for the identification of reference strains, followed by MSA. It is interesting to observe that the HA_Cal protein of the A/California/04/2009 strain disclosed a very close relation with the 10 Indian isolates (given in Table 1 ) reported from 2009 to 2018. Hence, these strains are selected for further comparative analyses to explore the evolution of HA_Ind proteins. The 3D structures of the HA_Ind protein are unavailable, whereas the crystal structures of the HA_Cal protein, available in PDB (for example, PDB ID: 3AL4, 3LZG, 3UBE, 3UBJ, 3UBN, 3UBQ, 3UYW, 3UYX, 3ZTN, 4JTV, 4JTX, 4JU0, 4M4Y, 5GJS, 5K9O, 5WKO, 6URM, and 6WJ1), lack a complete structure (Xu et al., 2010; Corti et al., 2011; Xuan et al., 2011; Xu et al., 2012; Hong et al., 2013; Zhang et al., 2013; Zhang et al., 2013; Joyce et al., 2016; Wang et al., 2016; Lang et al., 2017; Cheung et al., 2020; Wu et al., 2020). Due to the lack of a complete structure of the HA protein, we modeled the complete structure of both HA_Cal and HA_Ind-2018 using the full-length sequence (566 aa) and compared it with the reported crystal structure (PDB ID: 3LZG). The superimposition of both crystal and modeled structures reveals similar architecture, as shown in Supplementary File S5. Hereafter, the three-dimensionally modeled structures of the entire sequence of HA_Cal and HA_Ind-2018 are used for further comparative structural analyses.

Knowledge on the extent of genetic reassortment, antigenic shifts, and drifts in HA surface proteins of the H1N1 influenza virus isolates reported in India has become an indispensable concept as it discloses the most important factors related to its virulence. By examining the evolution of sequences, we tried to highlight how the selective pressure on the viral protein changes over time, leading to alterations in antigenicity, which further discloses variation in host specificity toward their receptor. A total of 512 HA protein sequences reported from the Indian strains of H1N1 viruses circulated during 2009–2018 (HA_Ind) were retrieved from the NCBI flu database (Supplementary File S1). At first, an exhaustive MSA was performed on these selected HA_Ind and HA_Cal sequences using ClustalW. In ClustalW, the sequences expressing variations due to mutation during the evolution of virus strains are aligned in accordance with the evolutionary distance and are further analyzed for phylogenetic relationships in a year-wise manner. A comparison of all HA_Ind proteins with reference to HA_Cal revealed the presence of new mutations as well. The MSA and phylogenetic tree (constructed using PAUP) revealed that HA_Ind proteins (evolved from 2009 to 2018, as given in Table 2) share a close evolutionary relationship with the HA_Cal protein (Supplementary File S2) and were chosen for further investigation.

Mutational analysis on the selected 10 HA_Ind sequences was performed in comparison with HA_Cal disclosed additional phenotypic variations at 84 positions, out of which 16 mutations gain importance as they share more conservation. Table 3 lists all the observed mutations in the selected HA_Ind proteins. In HA_Ind, positions 114, 180, 273, 300, 516, 202, 468, 214, 220, 100, 338, and 391 (the grey highlighted positions in Table 3) have a greater probability of mutation than HA_Cal. For a better understanding, the frequency (cumulative occurrence) of a mutation, in comparison with the reference HA_Cal, is calculated and depicted in the frequency plot (Figure 1). The residues disclosing more than 50% of mutational occurrences are in red font. Similar colors in the plot depict similar mutations, but the order of mutation from one residue to another is differentiated with the “*” symbol. For instance, the frequency of mutation of proline to serine is observed 10 times in HA_Ind proteins, whereas the mutation of serine to proline in HA_Ind proteins is only observed four times.

A year-wise analysis of mutational occurrence reported among the HA_Ind proteins from 2009 reveals that 1) all selected HA_Ind proteins possess mutations such as P100S, T214A, and I338V; 2) the mutations S220T and E391K are reported from 2010 onward; 3) the residues D114, K180, A273, K300, and E516 of HA_Cal remain conserved among the selected HA_Ind sequences reported until 2013 and the same sites disclosed mutations from 2014 onward such as D114N, K180Q, A273T, K300E, and E516K; 4) the mutations A13T, S101N, S179N, and I233T are reported since 2016; and 5) the mutations S202T and S468N are reported since 2013 and 2012, respectively. All these observations witness the occurrence of additional mutations that evolved over the successive period. It should be noted that the mutations of residues from T to A (T → A, S → N, D → N, K → Q, K → E, P → S, S → T, I → V, A → T, and E → K) play a significant role in the emerging diversity of HA_Ind proteins.

Wei Hu et al. reported seven receptor-binding sites (RBS) that are highly conserved in the HA_Cal protein. The variations analyzed at the RBS of HA_Ind in comparison with the reported RBS of HA_Cal (Hu, 2010) are given in Table 4. The RBS of the HA_Ind strain (2009) is reported with 22 single amino acid mutations at the following positions: 1) RBS3 (S160T, N164I, and S173P), 2) RBS4 (H197P, H198L, S200P, and Q210K), 3) RBS5 (V218G, Y223S, S224R, and A232E), 4) RBS6 (V266G, M274L, A278V, S280P, I282F, and I283F), and 5) RBS7 (T294P, T295P, T298P, A302V, and N304T), in comparison with the HA_Cal protein. Few more mutations reported in HA_Ind strains during 2010–2018 are discussed in the following paragraphs.

RBS1, RBS2, and RBS3 reveal mutations such as V33G (in 2012), S91R (in 2018), and S160G (2012), respectively. In RBS4, the mutation S202T that emerged in 2013 was conserved among the strains reported in successive years. RBS4 also witnessed an additional mutation (S200P) in the 2018 strain. In RBS5, all strains have inherited the A220T mutation along with few more such as 1) I233T in 2016, 2) A232G and I233T in 2017, and 3) I233T in 2018. In RBS6 and RBS7, mutations such as A273T and K300E were identified between 2014 and 2018. Specifically, the 2018 strain that emerged with seven mutations at the RBS (Figure 2) may imply that HA_Ind is more prone to mutation than HA_Cal.

In general, except for the 2009 strain, the RBS analysis reveals that the HA_Ind strains circulated during 2010–2018 were significantly conserved except for the few aforementioned mutations. Despite the observed significant conservation at the sequence level, the emerging single mutation posed a challenge to the inhibitors in sensing the receptor-binding sites and hence prompted the scientific community to design sequence-specific receptor-binding agents for further inhibition.

Here, we present a systematic analysis of the HA proteins of H1N1 to understand the adaptation and divergence among Indian strains due to frequent mutational events. Particularly, the analyses focused on the impact of missense mutations on receptor-binding domains, antigenic site alteration, N-glycosylation site prediction, amino acid compositional variability, and associated variability in secondary structure. In line with this, the present analysis also exclusively emphasizes on Indian strains and compares them with a recognized reference pandemic strain to understand the challenges behind the failure of successful medication in the Indian context. Accordingly, about 512 Indian strains were retrieved along with the A/California/04/2009 strain from the flu database of the NCBI reported during 2009–2018 and were analyzed by MSA to explore the evolutionarily conserved genetic regions. Site-specific variations observed in the aligned sequences, reflecting the rate of mutations or degree of evolution, were further analyzed to figure out the evolutionarily conserved regions with reference to the strain A/California/04/2009. A one-to-one relation between the aligned sequences was visualized using PAUP-generated phylogenetic tree, which displayed the closest evolutionary relationship between Indian and Californian strains of the influenza A (H1N1) virus (Supplementary File S2). Variations in sequences due to missense mutations at various positions play a vital role in altering the structure and function of different domains of the HA protein, such as the receptor-binding domain, epitope-binding domain, and N-glycosylation site.

Alignment of the amino acid sequences of HA_Ind reported during 2009–2018 showed about 84 amino acid substitutions when compared to the reference HA_Cal protein. About 24 substitutions were observed in HA_Ind-2018, in which 16 were highly concerned (Table 3). Analyses reveal that among the 16 mutations, seven mutations were found in the receptor-binding sites (Table 4), four were in antigenic sites (Table 5), and three were involved in the formation of N-glycosylation sites (Table 6). The HA_Ind strains are characterized by the mutations P100S, T214A, S220T, I338V, and E391K, i.e., possible beneficiary mutations that got fixed in the strains reported during 2009–2018 (Table 3). The literature suggests that T214A substitution in HA genes decreases the binding affinity with the host receptor (de Vries et al., 2013). We observed six new amino acid mutations (S91R, S138N, S200P, K319T, I421M, and E523D) in HA_Ind-2018. The mutations S91R and S200P were found to be unique in HA_Ind-2018, and these substitutions were abundant in the complete HA population (in 2018) compared to the pandemic HA_Cal. The substitutions A13T, S101N, D114N, I312V, S468N, and E516K were also observed in HA_Ind and are also reflected in the recent studies (Biswas et al., 2019; Prasad et al., 2020; Siddiqui et al., 2020). In accordance with our results, another research group carried out a mutational examination of H1N1 with random samples and observed that viruses circulated during 2017 have 18 detected substitutions in HA_Ind (Jones et al., 2019). They also reported I233T, S179N, S181T, and I312V as new substitutions, among which S181T and I312V were presented as unique mutations in HA_Ind isolates (Jones et al., 2019). Interestingly, we did not find I312V in 2017. The observed amino acid substitutions (S91, S200, S202, A214, and I233) have been found in receptor-binding sites envisaged to vary during the adaptation process to α2-6-linked sialic acid receptors in humans (Maines et al., 2009). The I223T amino acid substitution is linked with increased binding affinity to human α2-6-linked sialic acid receptors (Al Khatib et al., 2019). Substitutions S200P and S202T are responsible for enhanced receptor-binding avidity by altering the receptor-binding affinity, whereas the A214T substitution is linked to the decreased binding avidity (de Vries et al., 2013). A previous study suggested that S202T is one of the responsible substitutions involved in increased mortality and morbidity (Adam et al., 2019). Studies also support our observations that mutations of HA like P100S, T214A, S220T, I338V, and E391K are conserved mutations specific to the dominant variant(s) of influenza A (H1N1) viruses during post-pandemic circulation in India (Morlighem et al., 2011; Jones et al., 2019; Siddiqui et al., 2020). It is also evident from research that the substitutions S181T and I312V in HA could lead to altered glycan specificity (Jones et al., 2019). The substitution K180Q triggers conformational variation in ligand binding, which might trigger the failure of specific ligand-binding properties as well (Jones et al., 2019). The mutation S179N, associated with glycosylation, is responsible for the increased pathogenicity of the viral particle by preventing the antigenic sites of immune recognition (Al Khatib et al., 2019).

Out of 84 mutational sites, about 12 most probable conserved mutational sites at amino acid positions 100, 114, 180, 202, 214, 220, 273, 300, 338, 391, 468, and 516 have been observed in the last five consecutive years (Table 3). The HA_Cal protein possesses seven characteristic receptor-binding sites (Hu, 2010) and has been compared against all HA_Ind strains. Indian strains expressed mutations mainly like serine-to-threonine, alanine-to-threonine, and lysine-to-glutamine at various binding sites (RBS 4–7) over a period of time, i.e., mutations S220T (at RBS 5), S202T (at RBS 4), A273T (at RBS 6), and K300E (at RBS7) were reported since 2010, 2013, 2014, and 2018, respectively (Table 4). The results suggest that these mutations may trigger the alteration in the RBD and become resistant to available therapeutic options. In comparison, all HA_Ind proteins emerge with more point mutations than the selected HA_Cal, which play a significant role to evade from the known immune defense mechanism and further become life-threatening as well. Epitope mapping gains prime importance in the design of therapeutic monoclonal antibodies or vaccines, and any sequence-level mutations at the antigenic epitope sites hinder or delay the design of effective novel vaccines. In line with this, the studied Indian strains, which revealed significant mutations in the epitope-binding domains of HA proteins (Supplementary File S4), also delayed the successful identification of a vaccine for all Indian strains. Hence, faster evolution of the epitope-binding domain renders more complexity in the eradication process of the influenza A (H1N1) virus. Similar to epitope-binding sites, variations at N-glycosylation sites also increase complexity in the design of inhibitors (Zhang et al., 2004; Wei et al., 2010). Mutations in the glycosylation sites aid the HA protein (Table 6) to co-evolve with the host protein for the successful initiation of further infections.

The glycosylation of the influenza strain can disturb its host specificity, virulence, and contagious nature either directly by changing the biological activity of surface proteins (Schulze, 1997), or indirectly by 1) attenuating receptor-binding sites (Gao et al., 2009), 2) masking antigenic sections of the protein (Abe et al., 2004), 3) obstructing the HA protein precursor via its cleavage into the disulfide-linked subunits HA1 and HA2 (Ohuchi et al., 1989), and 4) regulating the catalytic activity or preventing proteolytic cleavage of the stalk of NA (Matsuoka et al., 2009). A report revealing the destabilization of the coiled coil of the HA protein due to the buried hydrophilic residue, Thr59, also endorses the sequential and structural level of distortions raised by the mutations threonine-to-serine observed in this study (Lin et al., 2018). Hence, the examined mutation-mediated structural diversification of the HA protein gains importance.

The ESPS characterized for the Indian isolate reported in 2018 revealed significant changes in the electrostatic surface, which is also presumed to render strong binding of HA proteins with the host receptors. The specific mutations observed in HA_Ind-2018 (for example, S91R, S181T, S200P, I312V, K319T, I421M, and E523D) may increase the fitness of the virus in a new environment and host, which may render a reduced efficacy toward the available treatment. The mutation-mediated adaptability and efficacy of HA_Ind proteins of the influenza A (H1N1) virus need to be studied critically.