{'fraction_non_alphanumeric': 0.03569495073447143, 'fraction_numerical': 0.03399862850543184, 'mean_word_length': 5.069660460021906, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Pan-India novel coronavirus SARS-CoV-2 genomics and global diversity analysis in spike protein

The mortality rates due to COVID-19 have been found disproportionate globally and are currently being researched. India mortality rate with a population of 1.3 billion people is relatively lowest to other countries with high infection rates. Genetic composition of circulating isolates continues to be a key determinant of virulence and pathogenesis. This study aimed to analyse the extent of divergence between genomes of Indian isolates (n = 2525 as compared to reference Wuhan-1 strain and isolates from countries showing higher fatality rates including France, Italy, Belgium, and the USA. The study also analyses the impact of key mutations on interactions with angiotensin converting enzyme 2 (ACE2) and panel of neutralizing monoclonal antibodies. Using 1,44,605 spike protein sequences, global prevalence of mutations in spike protein was observed. The study suggests that SARS-CoV-2 genomes from India share consensus with global trends with respect to D614G as most prevalent mutational event (81.66% among 2525 Indian isolates). Indian isolates did not reported prevalence of N439K mutation in receptor binding motif (RBM) as compared to global isolates (0.54%). Computational docking and molecular dynamics simulation analysis of N439K mutation with respect to ACE 2 binding and reactivity with RBM targeted antibodies viz., B38, BD23, CB6, P2B–F26 and EY6A suggests that variant have relatively higher affinity with ACE 2 receptor which may support higher infectivity. The study warrants large scale monitoring of Indian isolates as SARS-CoV-2 virus is expected to evolve and mutations may appear in unpredictable way.


Introduction
The current 2019 coronavirus pandemic  is caused by a positive RNA virus, referred to as severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [1]. Despite all the global efforts, the virus continues to spread and infect a large population and has affected 218 countries, with more than 101,053,721 confirmed cases, 2,182,867 deaths globally (WHO, 29 Jan 2021) [2]. One of the most striking features of COVID-19 is the variable mortality rate across countries [3]. Global case fatality rates (CFR) varied from 0.06 to 18.94%. France recorded the highest CFR of 18.65%. Belgium and Italy had double-digit CFR rates, while the United States, which experienced the largest outbreak, had a CFR of 5.76%. India, the second most populated country, reported 7,761,312 cases in all with 117,306 deaths (WHO, October 23) [4]. Trends in mortality rates in India are intriguing, despite the second largest reported cases of infection ("https://covid19.who.int/table"), second after the USA (more than 76 lakh) in terms of total confirmed cases. Although numerous hypotheses are proposed, one hypothesis includes the variability of circulating viral strains. Monitoring this variability in key viral genomes and immune-dominant antigens may explain the infectivity, pathogenesis, and transmission of SARS-CoV-2 [10].
Genome structure of SARS-CoV-2 is characterized by a size of around 30 kB with a large gene region (ORF) encoding non-structural proteins (NSP) and the genes encoding spike (S) glycoprotein, envelope protein, membrane (M) glycoprotein and nucleocapsid (N) proteins [7]. SARS-CoV-2 transient glycoprotein contains a region, the receptor binding domain (RBD), which specifically identifies the human angiotensin conversion enzyme (ACE2) as its receptor [8]. Antibodies targeting spike proteins, especially against the RBD are neutralizing in function [9]. Mutations that lead to variations in these hot spots may affect infectivity, pathogenesis, and transmission of SARS-CoV-2 [10,72,79].
We report here comparative genomics of Pan India isolates (n ¼ 2525) vis a vis reference genome and isolates from countries showing high fatality rates including France, Italy, Belgium, and the USA. The study also reports deep scanning of mutational events and their global frequency at each amino acid residue in spike protein using 1,44,605 global sequences submitted at GISAID. In-silco studies on impact of mutations on structure and binding affinity to critical human ACE2 receptor [73] and RBM targeted antibodies B38, BD23, CB6, P2B-F26 and EY6A were also carried out using molecular docking and simulation analysis [11].

Multiple sequence alignment and phylogenetic analysis of Indian isolates
Sequences from GISAID were downloaded and consensus sequences were aligned using MAFFT [18]. A maximum-likelihood phylogenetic tree was constructed using IQ-TREE and visualized with iTOL [19,20].

Structure prediction and docking analysis
Impact of mutations on binding affinity towards human ACE2 receptor was performed using docking and MD simulations analysis. The structures of wild type SARS-CoV-2 Spike RBD domain complexed with host ACE2 receptor was retrieved from Protein Data Bank (PDB ID: 6LZG) [21]. Spike protein was mutated using UCSF Chimera 1.10 software [22]. Mutant structures were retrieved from Chimera and energy was minimized for further analysis. All the wild and mutant structures were docked using ZDOCK docking server [23]. Docked structures were subjected to energy minimization and the binding energy was calculated using PDBePISA [24]. The structures of wild type SARS-CoV-2 Spike RBD domain complexed with antibodies B38,BD23,CB6, P2B-2F6, EY6A was retrieved from Protein Data Bank (PDB ID: 7BZ5, 7BYR, 7CO1, 7BWJ and 6ZER) respectively [12,13,14,25].

Modelling of the wild type and N439K variant with respect to ACE-2 binding
The structures of wild type SARS-CoV-2 Spike RBD domain complexed with host ACE2 receptor was retrieved from the Protein Data Bank (PDB ID: 6LZG) [21]. Spike protein was mutated using UCSF Chimera 1.10 software [22]. These structural models were used for molecular dynamics simulations. All molecular dynamics (MD) simulations were performed with GROMACS 2019 [74] software using CHARMM36 forcefield [75]. Explicit TIP3P water model was used to represent the water molecules. To neutralize the systems Na þ and Cl À ions were added. Energy minimization was performed using the Steepest Descent algorithm for 50000 steps. The system was equilibrated under the NVT conditions for 100 ps followed by NPT equilibration of 1000 ps. Temperature coupling was applied to maintain the system temperature at 300 K using the velocity rescaling algorithm. Semi-isotropic pressure was maintained using Parrinello Rahman pressure coupling with a pressure of 1 bar. Atomistic simulations were run for 100 ns for both the systems. Trajectories were viewed and analysed using Visual Molecular Dynamics (VMD) [75] tool. Standard GROMACS tools were used to plot the Root Mean Square Deviation (RMSD) and Root Mean Square Fluctuation (RMSF) of the system [77]. Hydrogen bond analysis was performed in VMD using Hydrogen bonds plugin. Binding free energies were calculated using MM-PBSA program in GROMACS [76]. To understand the diversity in the SARS-CoV-2 isolates from India as compared to reference strain (hCoV-19/Wuhan/WIV04/2019/EPI_-ISL_402124), we analysed a total of 2525 SARS-CoV-2 genomes retrieved from global dataset GISAID (23/10/2020) [26]. The samples were isolated between February-October 2020. A total of 2525 complete or near complete genome sequences with high coverage deposited from India were used in this study. The complete characteristics of 2525 SARS-CoV-2 genomes sequences are provided in additional file 1.
Higher frequency mutations were mapped for key SARS-CoV-2 protein structures, including nucleocapsid, nsp3, nsp12, nsp14, nsp2 and spike protein ( Figure 2). Out of 954 events, 585 events were found in open reading frame1ab (ORF1ab), which is the longest ORF occupying two thirds of the entire genome. ORF1ab is transcribed into a multiprotein and subsequently cleaved into 16 non-structural proteins (NSPs). Of these proteins, NSP12 has the largest number of variants including, P323L as dominant (n ¼ 1998) followed by A97V (n ¼ 328). D614G mutation in spike protein, which is considered as a prevalent global mutation, was present in 2062 of the 2525 (81.66%) sequenced genomes [27,28].
Regardless of fourth highest in the number of reported cases, Gujrat has highest death rate in the country as 6.2% which is double than the global mortality rate evaluated by WHO which 3.2% [2] was. Gujrat followed by Madhya Pradesh and West Bengal both of which have an estimated death rate of 4.2%. Gujrat and Madhya Pradesh showed similar trends in prevalent mutations in spike (D614G) and NSP12 (P323L), where West Bengal showed NSP12 (A97V), N (S202N), NSP2(G339S) mutations predominately ( Figure 3). Statewise prevalence of marker variant differed significantly ( Figure 3). Statewise mutation analysis suggests predominance of mutation in spike D614G with increased transmission and infectivity.

Phylogeny and lineage analysis
To study, SARS-CoV-2 classification for their multitude of distribution across geographical locations of India, vis-a-vis states of India, we explored the lineages spread across country. We further observed the prevalence of different clades within high fatality rate populations.

Viral clade analysis
The global isolates were assigned to major GISAID clades according to the marker variant and existing GISAID clades designations [32,33,34]. . The most represented clades were clade GR observed in 11, 298 complete genomes (GISAID, 18/06/2020) [35]. The clade was more prevalent in genomes sequenced from South America and Europe. The second most frequent clade represented was GH which was observed frequently from the genome sequences in North America and Africa continents [36]. The most represented clade from the Asia region was categorized as ancestral type O, which is similar with originating strain from Wuhan, China [6]. The clades are colour coded in squares and shown above chart in the diagram. As transmission across country is increasing trends in clade predominance was observed during these six months of the pandemic in the country. Ever since the emergence of outbreaks clade "O" is still prevalent, but the clades G, GH and GR were rapidly increasing since ease in social restrictions such as lockdowns and migration of workers across country. The wild type of clade "O" showed gradual decrease in the number of incidences over a time ( Figure 5).
Further to study association of fatality rate and predominant clade, we compared a total of 10,188 complete and high coverage genome sequences from countries reported higher fatality rates across the globe such as France, Italy, Belgium, and the USA, showed a dominant clade G in France, Italy, Belgium, and clade GH in the USA (Figure 4) (Additional file 3). Clade G and GH identified by the signature marker D614G variant in spike protein which is increasing its frequency across globe and specially developed countries. Recently reported as per Korber et al, D614G mutation increased infectivity of the virus [27]. India although reportedly third highest country in the world in terms of infected cases of COVID-19, fatality rate is still less as compared to the many developed world (https://www.mohfw.gov.in/). We observed a dominant clade ancestral type "O" in India was dominant during early phase of pandemic. We observed clade G, GH and GR as frequent clade in Gujrat, West Bengal, and Madhya Pradesh respectively (Figure 4). Enrichment and diminution of specific clades was observed in certain states of India [37]. Overall, within genome sequences from high fatality rates of COVID-19 regions, we observed a mutation D614G as prevalent with its resembling clades G, GH and GR.

Phylogenetic analysis
Phylogenetic analysis based on whole genome alignment of Indian genomes was performed using maximum likelihood tree rooted against reference hCoV-19/Wuhan/WIV04/2019/EPI_ISL_402124 ( Figure 6). Phylogenetic analysis of 864 genomes was done as per the definitions of the PANGOLIN lineage and GISAID clades [32].  Spike glycoprotein (S protein) of SARS-CoV-2 mediates receptor recognition and membrane fusion with the host cell [38,39]. During viral infection, the trimeric S protein is cleaved into S1 and S2 subunits and S1 subunits are released which contains the receptor binding domain (RBD), which directly binds to the peptidase domain (PD) of angiotensin-converting enzyme 2 (ACE2) [40]. Whereas S2 is responsible for membrane fusion. When S1 binds to the host receptor ACE2, another cleavage site on S2 is exposed and is cleaved by host proteases, a process that is critical for viral infection [39,40,41]. Many preprint studies reported different mutations at protein level and predicted its possible effect on binding affinity with host receptor [42,43,44]. Here, we scanned mutation at each single residue in 1273aa long protein sequence and explored their global prevalence. We also mapped mutation at different sites and regions on spike protein to understand its impact on its role in pathogenicity and disease transmission.

SARS-CoV-2 RBD mutation mapping
Total 29 mutations were reported in receptor binding domain of spike protein of SARS-CoV-2. Complete list of mutation, events and its origin are listed in Table 1. RBD domain showed only 0.54% of total mutational events (including n ¼ 1) where RBM showed 0.25% frequency of mutations. The receptor-binding motif (RBM) is the main functional motif in Figure 6. The phylogenetic tree based on the whole genome alignment of 863 genomes sequenced from India. Tree showed the presence of 7 different clades and demonstrates that SARS-CoV-2 is wildly disseminated across 237 distinct geographical location assessed till 19/06/2020. Tree was rooted with hCoV-19/Wuhan/WIV04/ 2019/EPI_ISL_402124 as a reference. Clades were represented as coloured ranges shown in left panel of the diagram. First outer circle represents lineages corresponding which are represented as strip label. Outer circle represents coloured labels representing different geographical regions from India; Maharashtra-Orange, Gujarat-Blue, Delhi-Green, Karnataka-Yellow, Telangana-Red, West Bengal-Purple, Uttar Pradesh-Green, Haryana-Grey.
RBD and is composed of two regions (region 1 and 2) that form the interface between the S protein and hACE2 [40].

SARS-CoV-2 RBM-hACE2 interacting residue mutation mapping
Computer modelling interactions have predicted the receptor binding motif (438-506) which act as main functional motif within RBD domain Figure 7. Prevalence of specific amino acid mutations among spike protein. Bubble map representing amino acid mutation and its prevalence observed among 1,44,605 genomes. Single amino acid position represents reference hCoV-19/ Wuhan/WIV04/EPI_ISL_402124. Circle size shows sum or frequency of events. The marks are labelled by mutations. The mutation of Aspartic acid to Glycine at position 614 was observed 1,23,415 times among the available 1,44,605 SARS-CoV-2 spike protein sequences used in this study. Overall, the S2 domain and specially heptads repeats were rich in mutational events observed and RBM region was found least variable (0.25%) among 1273 bp long protein. contains essential seventeen key residues that are potentially involved in binding with host cell receptor ACE2 (Table 2) [39,40]. We scanned mutation at these seventeen key residues Table 2 which may impact binding of virus and host receptor [47]. Mutations at contacting residues will be important to consider for the vaccine and therapeutic targeted against S protein (RBD region) [9,48]. Among these seventeen residues, twelve residues reported unique variations but at exceptionally low incidence (<2%) which can be due to sequencing errors or error prone replication of viruses within hosts. We observed key contacting residues were conserved amongst 99% of genomes of SARS-CoV-2 sequences.

Interpreting mutation effect of N439K in RBM region
To estimate the functional changes suggested by the RBM mutations, we used the prototype SARS-CoV-2 RBD domain and ACE2 and compared RBM mutants to assess their binding energy to human ACE2. The complex structure of the SARS-CoV-2 S-protein RBD domain and human ACE2 was obtained protein data bank (PDB ID: 6LZG) [49]. Mutant amino acids of the SARS-CoV-2 RBD mutants were directly replaced in the model, was used as a template modelled using SWISS-MODE [50]. We screened and evaluated interaction analysis of mutant N349K as it was most prevalent mutation observed in RBD-RBM region. RBM mutant types (N439K) exhibited significantly lowered ΔG, suggesting a slight  increased affinity to human ACE2; compared to the prototype. The ΔG of these mutant types were around -11.2 kJ/mol, lower than the prototype strain (-11 kJ/mol), where 0.093 kcal.mol-1.K-1 (increase of molecule flexibility) was observed with decrease in two hydrogen bonds.

Simulations predict that N439K variant binds tighter to human ACE2
As described in the Methods section, the structure of wild type Spike Receptor Binding Domain (RBD) complexed with host ACE2 receptor was retrieved from the PDB (PDB ID: 6LZG) [21]. The crystal structure of the Spike RBD complexed with ACE2 receptor reveal that the interactions are dominated by polar contacts mediated by hydrophilic residues [21]. Even a single mutation introduced in these contacts, for example, K353A could abrogate the interactions between the two proteins, highlighting the importance of the polar contacts. Thus, to understand the effect of the N439K variation, the mutant structure was built using Chimera software. Further, we performed 100 ns atomistic molecular dynamics simulations of the wild type and variant complexes. To quantify the structural variations in the protein complexes over time, the RMSD and RMSF was plotted. Both analyses indicated that the wild type of protein complex was more stable than the variant (Figure 8 A and B). On the contrary, a differential plot of the decomposed residue-wise binding free energies (Binding Free Energy of Residue Variant -Binding Free Energy of Residue Wild Type ) indicated a large change at position 439 of the Spike RBD, indicating that the residue contribution to free energy at this position was much more favorable in the variant than the wild type (Figure 8 C). A closer analysis indicates that in the presence of Lysine at position 439 of Spike protein forms water mediated hydrogen bonds with residues Gln 325 and Glu 329 of the ACE2 protein which are more distant in the presence of Asn at position 439 ( Figure 8 D and E). Thus, simulations coupled with binding free energy calculations suggest that the N439K variant binds ACE2 tighter than the wildtype.

Interpreting N439K in RBM region in the context of antibody epitopes
The RBD is the dominant target of neutralizing antibodies to SARS-CoV-1 also, SARS-CoV-2 [51,52]. RBD-targeting antibodies could neutralize binding SARS-CoV-2 with ACE2 [57]. Several antibodies that target the RBD of SARS-CoV-2 have now been reported in very recent studies [53,54,55]. It is unclear to what extent the RBM will evolve to escape neutralizing antibodies in a manner evocative of vaccine escape mutants [58]. Many antibodies have epitopes that overlap the RBD ACE2 contact interface and are therefore strongly constrained by mutation effects RBM region [ Figure 8], such as B38, BD23, CB6,P2B-F26, EY6A and REGN 10933 [12][13][14][15]51,54,56]. To better define the RBM's mutations N439K for antibody escape, we examined structural and binding constraint in the epitopes of antibodies with prototype structures of SARS-CoV-2 RBD.
Structural binding predictions suggests that prevalent mutation in spike RBD region, N439K did not overlap with currently characterized epitopes of neutralizing antibodies recovered from human convalescent patient [47,48,49]. The importance of tracking such mutations within RBM region is demonstrated by a recent study that identified mutations enabling escape from RBD-directed neutralizing antibodies [60]. Our data indicate that the none of the currently characterized mutations within RBM region escape the five RBM targeted neutralizing antibodies used in this study.

Discussion
COVID-19 pandemic since its emergence in China has been showing varied patterns of transmission and infectivity globally [61,62]. Studies on global diversity and real time tracking of mutations will be important in the disease control strategies [59,63,64,65]. India represents an excellent opportunity for such studies for its size, observing a strictest complete lockdown and varied environmental conditions. Many studies have reported the genomic characterization of Indian isolates till date with limited sample size [66]. We present here a complete comprehensive genomic characterization of pan India SARS-CoV-2 genomes till 23/10/2020 with more than eight months of pandemic, with an objective of understanding its global positioning. The study suggests similarity of Indian genomes with predominant types found across different parts of the globe [67]. It was noted that D614G mutant strains in India increased rapidly during eight months of the pandemic and is present in more than 81% of the genome sequences reported from India. The recent studies by Korber et al [27] and WHO collaborating study in China demonstrated D614G strains are 10-fold more infectious than the original Wuhan-1 strain [69]. There are reports suggesting associations of increased prevalence of G614 variant to the higher case fatality rates in 12 different countries [70]. In India, 81.66 % of genomes exhibited G614; however, its association with increased CFR could not be discovered in India. Distinctively, Indian isolates also did not displayed emergence of co-evolving mutations associated with mutant strain G25563T (ORF3a), C26735T (NSP14) and C18877T (M protein), which are observed globally. The study also attempted to look for such correlations in Indian states which exhibited highest mortality. Yet even in such states, the association with D614G could not be found.
Recently, the emergent UK variant strain named as VUI-202012/01 with defined as a set of 17 mutations, more significant N501Y was also observed and reported in India. As per recent studies, the UK variant reported higher binding affinity towards RBD region of spike protein. Preliminary modelling results communicated by the UK on 19 th December suggested that the variant is significantly more transmissible than previously circulating variants, with an estimated increase in reproductive number (R) by 0.4 or greater with an estimated increased transmissibility of up to 70% [78].The variant has proven to transmit more easily and had spread to 10 more countries over the short time. Similarly, another variant GH/501Y.V2,B. 1.351, famously known as South African Variant have been reported. The variant has spread to 31 countries and is shown to be responsible for higher fatality rates. Recently, four cases with the South Africa variant of SARS-CoV-2 were reported in India. Also, one case of Brazil variant known as 20J/501Y.V3 or P.1 lineage was also reported recently. India is reporting a spurt in the number of cases in some states. However, it needs to be ascertained whether this increase is related to emergence of UK, South African, Brazilian variants, or new Indian variants. To this effect, Indian Govt have launched Indian SARS-CoV-2 Genomic Consortia (INSACOG) which will accelerate virus surveillance, genome sequencing and characterization through a multi-laboratory network.
Many therapeutics are currently being pursued with the goal of developing protective immunity against the SARS CoV-2 virus [55,72]. Notably, immunogens and diagnostics are under development targeting spike protein sequence from the Wuhan reference strain. Several spike protein mutant strains are being reported so far. This study analysed 1, 44,605 spike protein sequences reported globally. Despite of a total of 3895 mutational events observed in spike protein sequences, frequency of events and its occurrence is still low in total available sequences except for D641G. RBD being an important region in host-virus interaction and thus for development of therapeutics, frequency of mutational events observed was very low [68]. The prevalence of mutations in RBD was found low (~0.57%) as compared to genome sequences reported globally. Single variant N439K at RBD region was observed in 0.52% (n ¼ 293) of genomes. Notably, N439K mutant is reported to show increased affinity binding affinity with hACE2. Interestingly, Indian isolates did not show presence of this variant. It is further noted that this variant is reported majorly few countries, including Scotland, England, and Romania. This led us to study the impact of observed N439K mutation within RBM region using computational docking and MD simulations approaches. Using a panel of RBM directed monoclonal antibodies viz. B38, BD23, CB6, P2B-F2F, and EY6A [12,13,14,15,51,54,56], study predicts slightly higher affinity to ACE receptor. This needs further confirmatory studies using in-vivo and in-vitro studies.
One of the reasons for lack of association of genomic signatures of Indian isolates with fatality rates could be low sample size. Large scale monitoring of Indian isolates is required for conclusive correlations. The data from this study along with serosurveillance data from India and clinical meta data will be useful to decipher associations between viral clades and severity. The approach and global positioning of Indian isolates presented in this study will be helpful for monitoring the pandemic in India. Several other factors for low susceptibility of Indians to COVID-19 are being also being proposed. Such factors include host innate immunity status, genetic diversity in immune responses, epigenetic, ABO blood group association and universal BCG immunization need to be investigated with the clinical studies.

Conclusions
The study suggests that SARS-CoV-2 virus diversity in India is consistent with global trends with respect to most prevalent mutations. No major mutation event was reported in RBD region of studied Indian isolates. Global spike protein mutation prevalence analysis suggests impact of some mutations on RBD and ACE receptor interaction which may affect virus infectivity. The study findings are cautiously optimistic that viral diversity will not be hindrance to current vaccine development strategies and will be broadly effective against current isolates.

Author contribution statement
Sunil Gairola: Conceived and designed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.
Shweta Alai: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Wrote the paper.
Nidhi Gujar: Performed the experiments. Manali Joshi: Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.
Manish Gautam: Conceived and designed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement
This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Data availability statement
Data associated with this study has been deposited at GISAID under the accession number provided in additional information.

Declaration of interests statement
The authors declare no conflict of interest.

Additional information
Supplementary content related to this article has been published online at https://doi.org/10.1016/j.heliyon.2021.e06564.