Genetic study of COVID-19 spread in Esteio, Rio Grande do Sul, Brazil

In December 2019, the causative agent of Coronavirus disease-2019 (COVID-19) pandemic named Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), emerged in Wuhan, Hubei, China [1]. As of 28 April, 2021, there are 148.963.836 confirmed cases and 3.140.213 million deaths in 192 countries [2]. Unprecedented international efforts of viral sequencing have allowed the submission of ~1.3 million genomes in the Global Initiative on Sharing All Influenza Data (GISAID) up to date [3], which are now available for studies of genomic epidemiology to follow the evolutionary history and dynamics of SARS-CoV-2 through space and time. In this sense, some important studies were already conducted in highly-affected countries, including USA [47], Italy [8], Netherlands [9], Australia [10, 11], and Brazil [1214].

By using a nomenclature developed to capture local and global patterns of genetic diversity of the virus, two main lineages (A and B) were identified, both originated in Wuhan and with simultaneous spreading around the world [15]. The dynamics of viral transmission in the Brazilian territory was investigated through the sequencing of ~500 genomes until the end of April, 2020. It was determined that: (i) B.1 and derived lineages were predominant at the beginning of the pandemic; (ii) >100 independent international introductions occurred in the country; (iii) a significant movement of the virus among the Brazilian regions was observed after international travel restrictions; and (iv) non-pharmacological measures were able to reduce the reproduction number (R₀) from >4 to 1 [12].

The genetic diversity of SARS-CoV-2 has been extensively studied, evidencing the presence of recurrent mutations such as S: D614G, S:E484K, S:N501Y across the world [1618], related to increased pathogenicity and transmissibility (higher viral loads, increased replication on lung epithelial cells, and enhanced binding affinity) [1924]. Furthermore, the E484K mutation was associated with immune evasion from neutralizing antibodies produced in response to currently available vaccines [2527]. In addition to SARS-CoV-2 mutations, co-infection with other pathogens (e. g. Staphylococcus aureus, Haemophilus influenzae, rhinovirus/enterovirus, respiratory syncytial virus, and seasonal coronaviruses) may be associated with poor clinical outcomes, although these rates appear to be limited [2831].

The viral mutations are not the only factor affecting the COVID-19 pathology and SARS-CoV-2 infectious capacity. Human host factors such as: (i) rare genetic variants governing interferon immunity [32], (ii) DNA polymorphisms in key host factors (e. g. Angiotensin-converting enzyme 2 [ACE2] and transmembrane protease serine 2 [TMPRSS2]) [33, 34], (iii) heritage and ethnicity [35], (iv) the presence of comorbidities (hypertension, diabetes, obesity, and immunological diseases) [36, 37] were already associated to increased disease severity, although more integrative studies are still needed to identify the relative contribution of each of these factors. By analyzing 27 candidate genes and Human leukocyte antigen (HLA) alleles in 954 admixed Brazilian exomes, 395 nonsynonymous variants were found. Of these, six were previously associated with the rate of infection or clinical prognosis of COVID-19. Seventy were identified exclusively in the Brazilian sample, and seven (10%) of these were predicted to affect protein function using in silico analysis [38].

As of March 10, 2020, a 60-year-old man who had been in Italy, became the first confirmed case in the southernmost state of Brazil (Rio Grande do Sul - RS) [39], which is the most populous state in the South Region of Brazil and the fifth in the whole country (~11.5 million inhabitants) [40]. As of April 28, 2021, Brazil has ~9.7% of worldwide cases (~14.4 million) and is the third worst-hit country [2]. The RS State reported ~956,030 cases and 24,458 deaths, with ~8% of cases requiring hospitalization [39]. The municipality of Esteio, located in the metropolitan region of RS capital, reported 9272 cases (total population: 83,202) and 302 deaths [41]. As Esteio is a commuter town, many workers move to and return from the state capital every day. Importantly, the case fatality rate (CFR) in Esteio (3.26%) was slightly higher compared to the RS state (2.56%), and both are greater than previous CFR estimates (~1%) [42, 43].

Thus, we aimed to characterize the main circulating lineages in Esteio (RS, Brazil) and their relationship with global, national and regional lineages using phylogenetics and phylodynamics inference from 21 SARS-CoV-2 genome sequences, including the investigation of putative viral mutations related to poor outcomes. Additionally, due to our typical subtropical climate and therefore high occurrence of respiratory infections, we investigated the occurrence of co-infections with other viral pathogens in these samples. The choice of a small municipality as the target of this study was important since we could more easily and precisely follow infected individuals, allowing a more detailed surveillance on the spread of the virus and detection of variability.

In the present work, we accessed SARS-CoV-2 mutations, circulating lineages and phylogenetic patterns of SARS-CoV-2 from a time- and age-representative set of patients admitted in a municipal healthcare system from the Southern region of Brazil, the third country most affected by the Covid-19 pandemic. As the study was conducted in a small municipality, we were able to track two clusters of viral mutations in epidemiologically-linked patients, highlighting the importance of viral dissemination in small areas of the community.

The SARS-CoV-2 spike (S) glycoprotein mediates the interaction with the ACE2 receptor in the host cells and it is the primary target of neutralizing antibodies [45]. There are structural unique spike features that contribute to its pandemic capacity: (i) a flat sialic acid-binding domain enables faster viral surfing over the epithelial surface before receptor interaction; (ii) tight and nearly perfect binding to the ACE2 entry receptor; (iii) the capacity to use furin and other proteases for cell entry [46]. A mutation in S (D614G) was recently associated with higher viral loads [19], increased replication on human lung epithelial cells [23], and younger age of patients [22], but not with the disease severity [19, 22]. This mutation is associated with an abolition in the hydrogen bond between the aspartate originally located at 614 position and a threonine residue in the 859 neighbouring protomer of the spike trimer, thus increasing the probability of RBD being found in the open state [47]. This promotes the binding with the ACE2 receptor leading to enhanced infectivity [23]. Importantly, this replacement was detected in 90.5% of our samples, compatible with the dominance of this variant in global context until December 14, 2020 (221,700 sequences [48];). Also concerning the S protein, we identified the V1176F variant in 38.1% of our samples and in all of them there was a co-occurrence with D614G. Importantly, 349/535 (65.2%) of the sequences isolated in the world including this replacement were from Brazil [48]. More recently (April 28, 2021), 1167 (53.7%) of worldwide sequences were from Brazil, representing mostly the spread of B.1.1.28 and P.1 lineages that harbor this mutation. This substitution is located on the C-terminal portion (S2 subunit), more specifically on the heptad repeat 2 (HP2), which is a major target of adaptive evolution in MERS-CoV-related viruses and carry sites associated with expanded host range in other coronaviruses [49].

Another replacement in the RDB of S protein (E484K) was also assigned. Until mid-December 2020, 157 genomes have this mutation globally, 114 (72.6%) isolated from South Africa, where a new lineage (B.1.351 or 501Y.V2) characterized by three RBD mutations (K417N, E484K and N501Y) recently emerged [50]. E484K also emerged independently in multiple lineages, including P.1, P.2 and B.1.1.33 firstly identified in Brazil between late 2020 and early 2021 [51]. As of 28 April, 2021, 38,436 genome sequences harbor this mutation, which are found in ~10% of the sequences generated to date on average (https://outbreak.info/situation-reports?muts=S%3AE484K↗) In Brazil, >90% of the sequenced samples since February 2021 carry this mutation, which is present in P.1 and P.2 emergent lineages (https://outbreak.info/situation-reports?muts=S%3AE484K&selected=BRA&loc=BRA↗). Recent evidence showed that E484K replacement enables viral escape from neutralizing monoclonal antibodies or polyclonal sera [2527] facilitating reinfection by emerging lineages harboring this mutation as reported in Brazil [52]. Importantly, two sequences (9.5%) from this study had four of five shared mutations with a new lineage (subsequently called P.2) reported in the Rio de Janeiro state (synonymous C28253T (F120F; ORF8), missense G28628T (A119S) and G28975T (M234I) in N protein and C29754U (3 UTR)), in addition to E484K RBD replacement [53]. Our sequences from Esteio were sampled in mid-October, as the first cases of the RJ lineage (Additional File 2). Phylodynamics inferences pointed that this lineage emerged in early July, approximately four months before the detection of its first genomes [53]. Moreover, its identification in a small municipality from the RS state (located 1.5 thousand kilometers from RJ) demonstrates that it emerged months before October and is already widely distributed in the Brazilian territory, but went unnoticed so far by the lack of appropriate genomic surveillance in Brazil.

The substitution of a negatively charged amino acid (Glutamic Acid) for the positively charged Lysine has a profound impact upon a highly flexible loop at the RBD. More specifically, it creates a strong ion interaction between lysine and amino acid 75 of hACE-2 (the main SARS-CoV-2 receptor). This link is not present in the wild type E484. Shifting key elements of RBD responsible for interactions may have a major impact in the risk of immune evasion. Also E484K could be related to enhanced infectivity, which may be associated with the rapid dissemination of these escape mutants [54]. The publication of a recent reinfection case of SARS-CoV-2 harboring E484K and the presence of this mutation in COVID-19 patients during the current second wave in Northern Brazil [55] are highly suggestive that this mutation is critical for viral evolution and thus must be investigated thoroughly [52].

The ORF6 accessory protein plays a critical role in antagonizing host antiviral responses and viral replication. Therefore, it potentially inhibits both type I interferon (IFN) production and downstream signaling [56, 57]. We have found the ORF6:I33T mutation in 42.9% of our samples, raising its potential association with immune suppression. The N protein packages the genomic RNA, playing a fundamental role during viral self-assembly [58]. It is also associated with replication-transcription complexes [59], and used as a target for diagnostic and immunogenic applications. Interestingly, a tri-nucleotide mutation in the Nucleocapsid gene (RG203KR) resulting in a double amino acid change was observed in high frequency in this study (~100), and are characteristic of B.1.1 and derived lineages that spread rapidly around the world [48]. The missense mutation in ORF1ab (L3930F) was reported in other 429 sequences (125 from Brazil) and also in two sequences of lineage B.1.1.248 from the Philippines [60]. Interestingly, the I33F and I292T mutations that were found in ORF6 and N, respectively, have been considered dominant mutations in the SARS-CoV-2 sequences from Brazil [61, 62]. The co-occurrence of these mutations represents the signature of Clade 2 (subsequently named B.1.1.33 lineage), which were one of the three most prevalent Brazilian viral lineage groups in the beginning of the pandemic, highly widespread in 16 states from Brazil [12]. It is important to emphasize that functional studies are necessary to characterize the effect of each viral mutation on transmissibility and pathogenicity. It is expected that most mutations will not have a great impact on viral evolution and its relationship with the human host. However, some of them may increase the viral fitness or represent some viral advantage in the pathogen-host interaction and, consequently, may become fixed in the population. Further studies to elucidate the interaction between human gene variants (especially in the ACE2 receptor) and SARS-CoV-2 mutations are necessary to establish the possible impacts of spike amino acid replacements with regard to ACE2 binding and function [63].

The knowledge of ACE2 physiological functions and specific features could explain how comorbidities like hypertension, diabetes, obesity, and immunological diseases can enhance the severity of symptoms. Thus, modulation of ACE2s function might promote pulmonary inflammation, thrombosis, obesity-induced hypertension, and cardiac failure, which are especially unfavorable to COVID-19 patients [36]. The expression of other cell receptors potentially involved with COVID-19 infection depending on age, gender, and characteristics such as obesity, smoking, and polymorphisms can contribute to patterns of severe symptoms [64]. Moreover, the hypertension-associated elevated immunological activity is a noteworthy factor that promotes the increased risk of hypertensive patients for critical COVID-19 outcomes. The increased immune activation can be observed in these patients and might elucidate the hyperinflammatory response (cytokine storm) [37].

The Brazilian population was mainly formed through an admixture process that comprises mostly European, sub-Saharan African, and Native-American ancestry [6567]. The broad spectrum of symptoms, days to symptoms onset, and the unpredictability of outcome in COVID-19 patients can also be linked with the vastly admixed Brazilian inhabitants. Moreover, recent reports identified HLA alleles previously associated with SARS-CoV-2 counteraction [38], and showed positive selection in genes associated with obesity, type II diabetes, lipid levels, and waist circumference [44]. Regarding COVID-19, epigenetics, specific variants, ACE2 and TMPRSS2 polymorphisms, ethnicity, as well as inborn immunity errors, have been reported worldwide. This suggests that different host genetic backgrounds might contribute to discrepancies in SARS-CoV-2 aggressiveness [3235].

We observed a higher CFR in the municipality of Esteio when compared to the RS state and the majority of Brazilian states, which might be linked to the emergence of new viral variants. However, COVID-19-related mortality is determined by both intrinsic factors of the infected individuals (age, comorbidities, and genetic characteristics) [68] and extrinsic aspects such as the access to healthcare assistance (hospital beds, mechanical ventilators, medicines). Additionally, the Southern Brazilian states (Rio Grande do Sul, Santa Catarina, and Paran) have important determinants of mortality: older population than other regions and highest historical incidence of SARS in the country [69]. Their proximity to the states of So Paulo and Rio de Janeiro (that account for ~30% of the national population) also facilitates travel between the two regions and the rapid dissemination of emerging lineages.

Importantly, CFR is highly influenced by the underreporting of confirmed cases and deaths. States with low testing capacity tend to generate higher CFRs, and recently many deaths with an undetermined cause have been reported in Brazil, which also affects the quality of the records [70]. Therefore, the analysis of lethality should take this combination of factors into account [71]. In the case of Esteio, it is the municipality that tested more proportionally to the number of cases (38,416 tests per 100,000 inhabitants), thus these underreporting biases should be less significant, in contrast to the observed in other Brazilian states.

Compared to three other studies conducted in Brazil [12], in the states of Minas Gerais [14] and Pernambuco [13] in the early phase of pandemic which showed the introduction of viral lineages from other continents (mainly Europe) by international returning travelers , this study suggests a minor role of international lineages in the ongoing viral transmission in Esteio. We speculate a trend towards the perpetuation and diversification of the lineages found in this study (B.1.1.248 and B.1.1.33) inside Brazil. The dissemination of these lineages were also reported in the Uruguayan-Brazilian border, driving viral introductions mainly from Southeast and Southern Brazil (especially RS state) to Uruguay [72]. In this study, we found consistent results, mainly regarding B.1.1.33 diffusion from Southern and Southeast to South-American countries (e. g. Argentina, Chile and Uruguay). These lineages have already formed new sublineages (https://cov-lineages.org/lineages.html↗). B.1.1.33 has evolved in 10 new sublineages (N.1 to N.10). Furthermore, B.1.1.248 has evolved in P.1 (lineage first identified in Manaus [55] associated with a constellation of spike mutations like B.1.1.7 [73] and B.1.351 [50]), P.2 (lineage firstly identified in RJ state and also found in this study), and P.3. Importantly, all these three sublineages harbor the E484K mutation, which arose independently in both of them and appear to be evolving under diversifying positive selection [50, 51].

We built a time-resolved phylogeny prioritizing sequences that are genetically and spatially closer, but maintaining a global representativity of viral spread. This allowed us to confidently identify that our sequences fell into two main clades, with a broad presence of Brazilian and local sequences. We also inferred the spatiotemporal diffusion of these main lineages in regional and global context, finding the key role of Southeast in disseminating these lineages across Brazilian states and other continents. We also found a broader clade represented by Southern Brazil sequences and its important contribution in disseminating B.1.1.33 to South American countries. Moreover, we found four broader clades in the B.1.1.248 MCC phylogeny, suggesting multiple introduction events from Southeast followed by community transmission.

Our evolutionary rate estimates for both the broader ML tree and lineage-specific MCC trees (6 to 710⁴ subst/site/year) were slightly smaller than previous findings (8 to 910⁴ subst/site/year) [74, 75]. These differences have contributed to the estimates of the T_MRCA for both B.1.1.248 and B.1.1.33, which were dated to late 2019, contrasting with the first description of these lineages [12]. Furthermore, other probable sources of these inconsistencies are: closely related samples having the same age (phylo-temporal clustering), among-lineage rate variation and non-random sampling [76]. Although it is possible that different branches of the phylogeny have different rates, when we used a model that allows different rates across the tree (uncorrelated lognormal relaxed clock), the T_MRCA estimates have remained unchanged.

An important caveat for the phylodynamic analysis is that samples are not equally distributed geographically or temporally. This is a consequence of episodic sampling efforts prompted by research resource availability, and does not necessarily resemble a representative uniform sample. Unequal temporal distribution implies that some of the conclusions are disproportionately influenced by events in heavily sampled periods (FebruaryMay). Additionally, a large proportion of the samples come from the Southeastern region of Brazil. While this is in fact a heavily hit region and an economic and travel hub for the country, other regions such as the North are underrepresented. Thus, the strong support for the Southeast as prime center for viral dispersion and location of the root of both clades might be somewhat inflated, and epidemiological links between other regions could be downplayed due to undersampling. However, a study from the beginning of pandemic (FebruaryMarch 2020) estimated that the main destinations of the international passengers arriving to Brazil were So Paulo (46.1%), Rio de Janeiro (21%) and Belo Horizonte (4.1%), three capitals from the Southeast and therefore routes for COVID-19 importation [77]. Moreover, during mid-February and mid-March, SARS-CoV-2 spread mostly locally and within-state borders. In contrast, during mid-March and mid-April there was an ignition of the epidemic from the Southeast region to other states [12], which is consistent with our findings.

Some limitations should be considered. Firstly, it was not possible to analyze a larger sample size. Moreover, the low quantity and spatial representativity of sequences from the RS state to contextualize our sequences limited the inference of events of introduction and movement of the virus with municipal and state resolution. Still in this respect, we have observed a dramatic drop in the sequencing efforts from Brazil after April 2020 [12], which made it difficult to measure the main circulating lineages in the country during our investigation period (MayOctober, 2020) and may introduce confounding factors.

Since the E484K mutation identified in this study has been associated with loss of neutralizing activity from convalescent plasma (immune evasion) and enhanced interaction with hACE-2, lineages containing this substitution must be the subject of intense surveillance. More specifically, it is critical that immune strategies such as convalescent plasma and vaccines be tested against these new variants. Attempts to demonstrate activity against S mutants should be a priority effort for all vaccine and monoclonal antibody makers. Second generation immune therapies might have to be directed at more conservative neutralizing binding sites (such as S2 fusion domain) or elicit strong cellular response in order to keep on long term protection. Finally, human genetic factors, patient heritage and health conditions should also be studied in an integrated way for a broader understanding of vaccine effectiveness in different populations.

Our results provide a comprehensive view of SARS-CoV-2 mutations from a time- and age-representative sample from May to October 2020, highlighting two frequent mutations in spike glycoprotein (D614G and V1176F), an emergent mutation in spike RBD (E484K) characteristic of B.1.351 and P.1 lineages, and the adjacent replacement of 2 amino acids in Nucleocapsid phosphoprotein (R203K and G204R). In particular, to our best knowledge, we described the earliest SARS-CoV-2 sequences harboring E484K in Southern Brazil. A significant viral diversity was evidenced by the absence of identical isolates in our samples. Furthermore, we identified patterns of SARS-CoV-2 viral diversity inside Southern Brazil, demonstrating the major role of community transmission in viral spreading and the establishment of Brazilian lineages ignited from the Southeast to other Brazilian regions. Our data show the value of correlating clinical, epidemiological and genomic information for the understanding of viral evolution and its spatial distribution over time. This is of paramount importance to better inform policy making strategies to fight COVID-19.