Genomic monitoring unveil the early detection of the SARS-CoV-2 B.1.351 (beta) variant (20H/501Y.V2) in Brazil.

Sao Paulo State, currently experiences a second COVID-19 wave overwhelming the healthcare system. Due to the paucity of SARS-CoV-2 complete genome sequencing, we established a Network for Pandemic Alert of Emerging SARS-CoV-2 Variants to rapidly understand and monitor the spread of SARS-CoV-2 variants into the state. Through analysis of 210 SARS-CoV-2 complete genomes obtained from the largest regional health departments we identified cocirculation of multiple SARS-CoV-2 lineages such as B.1.1 (0.5%), B.1.1.28 (23.2%), B.1.1.7 (alpha variant, 6.2%), B.1.566 (1.4%), B.1.544 (0.5%), C.37 (0.5%) P.1 (gamma variant, 66.2%), and P.2 (zeta variant, 1.0%). Our analysis allowed also the detection, for the first time in Brazil, the South African B.1.351 (beta) variant of concern, B.1.351 (501Y.V2) (0.5%), characterized by the following mutations: ORF1ab: T265I, R724K, S1612L, K1655N, K3353R, SGF 3675_F3677del, P4715L, E5585D; spike: D80A, D215G, L242_L244del, A262D, K417N, E484K, N501Y, D614G, A701V, C1247F; ORF3a: Q57H, S171L, E: P71L; ORF7b: Y10F, N: T205I; ORF14: L52F. The most recent common ancestor of the identified strain was inferred to be mid-October to late December 2020. Our analysis demonstrated the P.1 lineage predominance and allowed the early detection of the South African strain for the first time in Brazil. We highlight the importance of SARS-CoV-2 active monitoring to ensure the rapid detection of potential variants for pandemic control and vaccination strategies. Highlights Identification of B.1.351 (beta) variant of concern in the Sao Paulo State. Dissemination of SARS-CoV-2 variants of concern and interest in the Sao Paulo State. Mutational Profile of the circulating variants of concern and interest.

INTRODUCTION

The emergence of a new severe acute respiratory syndrome‐related coronavirus 2 (SARS‐CoV‐2) lineage, namely P.1 (gamma variant), in the Brazilian city of Manaus and its rapid dissemination despite the relatively high seroprevalence in the city, demonstrated the importance of the SARS‐CoV‐2 genomic surveillance with the identification of mutation constellations which may impact viral infectivity, immunological evasion, and phenotypic characteristics. Moreover, the importance of these lineages is related to the significant pressure which they exert on the healthcare system of the affected region due to the high morbidity and mortality leading to an overall collapse of the intensive care units, a situation which occurred when P.1 emerged in the Amazonian city of Manaus.

Apart from data originating from the Brazilian Amazon, the circulating lineages in the Sao Paulo State remain largely unknown. Sao Paulo State incorporates the most populous area of Brazil and is currently experiencing a second wave of SARS‐CoV‐2 infections with an estimated total number of cases 3 265 930 up to May 31, 2021, and 111 304 deaths with an increase of approximately 50% of the confirmed SARS‐CoV‐2 cases. Such a situation gives significant pressure on the healthcare system of the most industrialized Brazilian area. Due to this, the Network for Pandemic Alert of Emerging SARS‐CoV‐2 Variants was established to implement large‐scale genomic surveillance to monitor and characterize the circulating SARS‐CoV‐2 lineages in the Sao Paulo State. A total of 210 samples were obtained from large cities which constitute the Regional Health Departments of the State up to March 2021. In this study, we also demonstrated for the first time the presence of a B.1.351 (beta) lineage in Brazil.

MATERIALS AND METHODS

Study design

We randomly selected and sequenced 210 SARS‐CoV‐2 samples with quantitative reverse transcription‐polymerase chain reaction positive results from laboratories that were responsible for SARS‐CoV‐2 testing in the São Paulo State and were comprising the SARS‐CoV‐2 Diagnostic Network established by the Butantan Institute. For viral genotyping samples with cycle threshold (C t) values of up to 30 were only used. Moreover, RNA extraction was performed with the Extracta kit AN viral (Loccus) in an automated extractor (Extracta 32; Loccus) following the manufacturer's guidelines. SARS‐CoV‐2 molecular diagnosis was carried out using the kit Gene Finder™ COVID19 Plus RealAmp kit (OSang Healthcare, Co., Ltd.) in all laboratories comprising the network which reduces variations related to C t values.

The SARS‐CoV‐2 full genome sequencing was performed by Oxford Nanopore Technologies platform using the multiplex PCR amplicon sequencing approach developed by the ARTIC Network (https://www.protocols.io/view/ncov-2019-sequencing-protocol-v3-locost-bh42j8ye). Samples were first reverse‐transcribed using LunaScript RT SuperMix (New England Biolabs). Amplification was carried out using Q5 Hot Start High‐Fidelity 2X Master Mix (New England Biolabs) and an Artik v3 primer set, targeting the whole SARS‐CoV‐2 genome (https://github.com/artic-network/artic-ncov2019/blob/master/primer_schemes/nCoV-2019/V3/nCoV-2019.tsv). Sequencing libraries were prepared using the Oxford Nanopore Ligation Sequencing Kit (SQK‐LSK109) and Native Barcoding Expansion kits (NBD104 and EXP‐NBD114) following previously published protocol. The libraries were loaded on an R9.4 flow cell (FLO‐MIN106) using a MinION MK1B device (ONT).

Data analysis

Regarding the SARS‐CoV‐2 genome obtained from the city of Sorocaba, raw files were basecalled using Guppy v4.4.2 (https://nanoporetech.com/nanopore-sequencing-data-analysis) and barcode demultiplexing was performed using qcat v1.1.0 (https://github.com/nanoporetech/qcat). We used Genome Detective and Coronavirus Typing Tool to obtain consensus sequences by de novo assembling. For the remaining genomes, assembly against the SARS‐CoV‐2 reference (Genbank refseq NC_045512.2) was performed by a pipeline of BWA for mapping, samtools for read indexing and to compile per base variation, bcftools for variant calling, and seqtk for creation of consensus genome.

Phylogenetic and phylodynamic analysis

A representative subset of 3852 genomes obtained from GISAID was obtained following the Nextstrain guidelines (https://nextstrain.github.io/ncov/customizing-analysis.html). These sequences were separated by downsampling as implied by Nextstrain using the most suitable parameters like space and sequence distribution. The 210 full‐length new genomes were appended to this subset for further analyses. Sequence alignment was performed using MAFFT v7.475 and manually curated to remove artifacts using Aliview. Approximate maximum likelihood phylogenetic trees were estimated using FastTree v2.1.10 with local branch support by resampling of site likelihoods 1000 times in conjunction with the Shimodaira–Hasegawa test.

Beast v1.10.4 was used to infer the date of origin of the Sorocaba haplotype within the B.1.351 clade. A Skyride tree prior was employed. Rate prior was incorporated within an uncorrelated lognormal model of rate variation across branches. The parameter ucld.mean (the rate variation prior mean) was set as a uniform [8.99E−4; 1.66E−3] substitutions/site/year following Ghafari et al. Two Markov chain Monte Carlo chains were run in parallel for 100 M generations, with convergence assessment checked individually in Tracer v1.7.1, and until all parameters had effective sample sizes more than 200.

Mutational patterns analysis

Mutational profile was investigated by using the Nextclade tool (https://clades.nextstrain.org/) to describe substitutions. Subsequently, we compared the set of nonsynonymous mutations with the profiles given in the PANGO lineages resource (https://cov-lineages.org/lineages.html) to attribute genomes to lineages.

Ethical statement

This study was approved by the Institutional Ethics Committee of the Faculty of Medicine of Ribeirão Preto (Process CAAE: 38975620.1.1001.5440).

RESULTS

The 210 newly sequenced genomes from the Sao Paulo State classified by Pango lineage demonstrated that the majority of sequences belonged to the P.1 (gamma) lineage (66.2%) followed by B.1.1.28 (23.2%). Underrepresented lineages were B.1.566 (1.4%), B.1.544 (0.5%), and C.37(0.5%). The P.2 (zeta) Brazilian variant was also detected in 1.0% of the cases. The UK B.1.1.7 (alpha) variant of concern (VOC), B.1.1.7 (6.2%) was also detected with a relatively high percentage.

During this genomic surveillance, we could identify and characterize for the first time in Brazil a unique isolate classified as B.1.351 (0.5%) which belongs to the South African lineage (beta variant). This isolate was detected in the southeastern part of Sao Paulo State, and the infected individual reported no recent travel history within and/or outside Brazil. The B.1.351 VOC was characterized by K417N, E484K, N501Y in the spike region and 19 mutations and 2 deletions as follows: ORF1ab: T265I, R724K, S1612L, K1655N, K3353R, SGF 3675_F3677del, P4715L, E5585D; spike: D80A, D215G, L242_L244del, A262D, D614G, C1247F; ORF3a: Q57H, S171L, E: P71L, ORF7b: Y10F, N: T205I, ORF14: L52F (Figure 1D). The dissemination of SARS‐CoV‐2 lineages in the State according to each region is shown in Figure 1A.

Figure 1

Genomic characterization of SARS‐CoV‐2 variants circulating in Sao Paulo State. (A) Map of Brazil and Sao Paulo State (on the right side) showing the frequency of SARS‐CoV‐2 lineages obtained in this study. (B) Approximate maximum likelihood phylogenetic tree including 4069 genomes from GISAID, including the 217 new isolates obtained in this study plus 3852 SARS‐CoV‐2 reference strains collected up to March 21, 2021. (C) Local subtree showing the closest B.1.351 genomes to the African‐like genome from the State of São Paulo (Brazil) after divergence dating

We then explored the genetic relationship of the newly sequenced SARS‐CoV‐2 genomes to those of other isolates by phylogenetic inference. Figure 1B highlights the main SARS‐CoV‐2, VOCs (P.1, B.1.1.7, and B.1.351), and variants of interest (VOIs) (P.2) circulating into the State. To get more insight regarding the early detection of the SARS‐CoV‐2 B.1.351 VOC we performed Bayesian analysis including all B.1.351 VOCs available in GISAID up to March 28, 2021. In Figure 1C, the cluster including the Brazilian B.1.351 (beta variant) is expanded. Although genetic data alone cannot provide evidence of directionality of transmission, our estimates suggest this variant might have reached Brazil through some importation events mediated by returning travelers from European countries, reinforcing how the high connectivity of countries can mediate the introduction of new SARS‐CoV‐2 strains. Furthermore, our estimates suggest that the most recent common ancestor of the Brazilian B.1.351 sequenced genome originated between mid‐October 2020 and end‐December 2020 (95% highest posterior density).

Genomic characterization of the B.1.1.7 (alpha), P.1 (gamma), and P.2 (zeta) lineages

To better characterize the genomes recovered with greater impact in the pandemic, we focused on describing SARS‐CoV‐2 VOC in more detail. For this, we use an analysis where we recover the variations of the genome in relation to the reference strain (NC_045512.2) through the use of bioinformatics pipelines. The UK variant, that is the alpha variant (B.1.1.7, 20I/501Y. V1), was characterized by 14 mutations that define this isolate, 6 of which in spike (ORF1ab: C3267T, C5388A, T6954C; spike: A23063T, C23271A, C23604A, C23709T, T24506G, G24914C; Orf8: C27972T, G28048T, A28111G; nucleocapsid: C28977T). Our strains belonging to this lineage (6.2%) present without exception all these defining mutations. However, a fact that draws attention is the rise of the A23403G (D614G) mutation being present in all samples.

Furthermore, we characterized the P.1 (gamma) (20J/501Y.V3) VOC containing 15 mutations in its genome, 10 of which were located in the spike region. P.1 defining mutations related to each genomic region were the following: ORF1ab: S1188L, K1795Q, E5665D; spike: L18F, T20N, P26S, D138Y, R190S, K417T, E484K, N501Y, H655Y, T1027I; Orf8: E92K; and nucleocapsid: P80. We additionally observed the D614G mutation in the spike, as well as the mutation G25088T (V1176F). Finally, we characterize the P.2 VOI that was distinguished by five mutations: ORF1ab: C100U; Orf8: C28253U; nucleocapsid: G28628U, G28975U, and C29754U; and spike: G23012A (E484K). We also highlight the presence of the two mutations—D614G and V1176F—in the spike gene in all analyzed P.1 and P.2 variants of this study.

DISCUSSION

In this study, we report the dissemination of SARS‐CoV‐2 variants in the Sao Paulo State. After the notification of the first confirmed case of COVID‐19 in Brazil in the city of São Paulo in February 2020, São Paulo State experienced two SARS‐CoV‐2 waves as the second started in December 2020 and is currently ongoing. To date, the average number of daily confirmed COVID‐19 cases increased 1.8‐fold, from 100 601 positive cases confirmed by RT‐PCR during December 2020 to 187 810 cases during March 2021 (unpublished data from The Network of COVID‐19 Diagnosis of Sao Paulo State, coordinated by Butantan Institute). Such uncontrolled growth probably reflects VOC spread, as observed during the P.1 (gamma variant) emergence in the Amazonian city of Manaus, , and requires a robust SARS‐CoV‐2 genomic surveillance.

Here, we demonstrated that two main lineages dominate the current SARS‐CoV‐2 scenario in the Sao Paulo state P.1 (gamma variant) and B.1.1.28 which were highly represented (89.4%). Therefore, a dynamic process in this Brazilian region related to the substitution of the initial SARS‐CoV‐2 lineages (B.1.1.28 and B.1.1) by a mixture of SARS‐CoV‐2 VOC/VOIs including P.1 (gamma), P.2 (zeta), B.1.1.7 (alpha), and B.1.351 (beta) which actually represent 72.9% of circulating SARS‐CoV‐2 lineages is probably occurring. In this context, the increased transmissibility, which has already been associated with VOC suggests a higher number of positive cases, hospitalizations and mortality observed not only in the Sao Paulo State but in Brazil in general.

The mutational and phylogenetic analysis were consistent with the identification of B1.351 (beta) VOC for the first time in Brazil. Based on the molecular, phylogenetic/phylogeographic, and epidemiological data, we suggest that B.1.351 presence in Brazil may be related to an introduction from travelers originating from another country and leading to a local B.1.351 transmission. This is highly possible due to the actual widespread presence of B.1.351 in almost 70 countries (https://cov-lineages.org/global_report_B.1.351.html). Furthermore, a performed deep analysis of the mutational profile of this strain demonstrates that in comparison with the South African reference strain there were present three additional specific mutations in the spike genomic region (A262D, D614G, and C1247F). This suggests that the detected B.1.351 VOC may have experienced molecular evolution locally with the acquisition of characteristic mutations, though at this point the actual transmission chain within the country is still under investigation. Therefore, more genomic data surveillance is necessary to robustly evaluate the dissemination and transmission route of the VOC B.1.351 lineage in Brazil.

Other SARS‐CoV‐2 VOCs were also identified, including the UK variant, which shows that the SARS‐CoV‐2 lineages in the Sao Paulo State are presented as a complex mixture suggesting that more extensive genomic surveillance is urgently needed. In this respect, extensive VOC characterization and collaborative efforts are essential to evaluate the current SARS‐CoV‐2 scenario not only in the Sao Paulo State but also in Brazil due to the general worsening of the pandemic in a nationwide aspect.

In conclusion, the performed study emphasizes the importance of the SARS‐CoV‐2 genomic surveillance to monitor the SARS‐CoV‐2 dissemination in the Brazilian regions and states. Given the vast extension of the country, a collaboration of different sequencing networks and combining SARS‐CoV‐2 genomic data will be of crucial importance to understand in more detail the nation‐wide Brazilian SARS‐CoV‐2 pandemic, especially related to VOC emergence and improve the responsiveness to further possible waves and vaccination strategies to SARS‐CoV‐2.

CONFLICT OF INTERESTS

The authors declare that there are no conflict of interests.

AUTHOR CONTRIBUTIONS

Sample screening: Svetoslav N. Slavov, Marta Giovanetti, Antonio J. Martins, Vincent L. Viala, Evandra S. Rodrigues, Elaine V. Santos, Claudia R. S. Barros, Elaine C. Marqueze, Bibiana Santos, Flavia Aburjaile, Debora B. Moretti, Erika Freitas, David Schlesinger, and Simone Kashima. Sequencing: Svetoslav N. Slavov, Marta Giovanetti, Vincent L. Viala, Evandra S. Rodrigues, Elaine V. Santos, Claudia R. S. Barros, Elaine C. Marqueze, Bibiana Santos, Flavia Aburjaile, João P. Kitajima, Erika Freitas, David Schlesinger, Maria C. Elias, and Simone Kashima. Bioinformatic analysis: José S. L. Patané, Vincent L. Viala, Rafael dos Santos Bezerra, Marta Giovanetti, and Vagner Fonseca. Conceptualization: Svetoslav N. Slavov, José S. L. Patané, Rafael dos Santos Bezerra, Marta Giovanetti, Antonio J. Martins, Vagner Fonseca, Vincent L. Viala, Raul M. Neto, Debora B. Moretti, Rodrigo T. Calado, João P. Kitajima, Luiz C. Junior de Alcantara, Maria C. Elias, Sandra C. Sampaio, Simone Kashima, and Dimas T. Covas. Writing original draft: Svetoslav N. Slavov, José S. L. Patané, Marta Giovanetti, Luiz C. Junior de Alcantara, Maria C. Elias, Sandra C. Sampaio, and Simone Kashima. Funding acquisition: Raul M. Neto, Ricardo Haddad, Rodrigo T. Calado, Luiz C. Junior de Alcantara, Maria C. Elias, Sandra C. Sampaio, Simone Kashima, and Dimas T. Covas.

Supporting information.

Click here for additional data file.