Abstract
We developed a genomic surveillance program for real-time monitoring of SARS-CoV-2 variants of concern in Uruguay. Here, we present the first results, including the proposed qPCR-VOC method, the general workflow and the report of the introduction and community transmission of the VOC P.1 in Uruguay in multiple independent events.
Introduction
Throughout the pandemic, as part of the natural viral evolution, genetic variants of SARS-CoV-2 have emerged. Some of these have shown increased transmissibility, more severe disease with increased hospitalizations, intensive care units admissions and deaths, significant reduction in neutralization by antibodies generated during previous infection or vaccination, reduced effectiveness of treatments or vaccines and/or diagnostic detection failures [1, 2]. These are called variants of concern (VOCs). To date, five VOCs have been identified: B.1.1.7, B.1.351, P.1, B.1.427 and B.1.429 [3]. A robust surveillance workflow for early identification of VOCs is key to accelerate pandemic responsiveness.
Surveillance data from Brazil demonstrated a sharp increase in the number of SARS-CoV-2 cases, hospitalisations and deaths after the emergence of the VOC P.1 in the Amazonas state in November 2020 [4, 5]. The VOC P.1 displays a higher transmissibility than previous local SARS-CoV-2 lineages [4, 5] and rapidly spread through the country, becoming the predominant strain in most country regions by February-March, 2021 [6]. This VOC has also spread worldwide and as of 4 May 2021, the lineage P.1 has been detected in at least 41 countries worldwide [7]. This raises the concern that VOC P.1 might follow in other countries the same epidemic trajectory observed in Brazil.
Uruguay, which shares 600 miles of dry border with Brazil and hosts more than 170,000 people living in twin cities with intense economic and social interaction, has experienced an exponential increase of COVID-19 cases since late February, 2021, placing it at the top of the countries with the higher number of daily confirmed cases and deaths per million people. Despite the Brazilian-Uruguayan border has remained closed to nonresident foreigners since March 2020, there is evidence of high viral flux between both countries [8, 9], implying the risk for imminent introduction of P.1 into Uruguay and the necessity of coming up with an organized strategy to monitorize this and other viral variants circulating in the country.
We developed a comprehensive genomic surveillance program for real-time monitoring VOCs emergence in Uruguay to help sanitary authorities to manage the crisis. We created an inter-institutional working group (IiWG) to run this program aiming to put together: i) a diagnostic network, ii) expertise and resources to handle large-scale sequencing, iii) computational scientists to analyze genomic datasets, and iv) an affordable and decentralizable “in house” qPCR test specifically designed to detect worldwide known VOCs (B.1.1.7, B.1.351 and P.1). Here, we present the results of the first two weeks of the IiWG, including the proposed qPCR-VOC method, the general workflow and the report of the introduction and community transmission of the VOC P.1 in Uruguay in multiple independent events.
Results
Inter-institutional working group (IiWG) for real time genomic surveillance of SARS-CoV-2
This task force is form by four Institutions (Institut Pasteur de Montevideo, Universidad de la República, Zurgen-Sanatorio Americano and Ministerio de Salud Pública), including more than 30 people, among them researchers, physicians and health authorities. Four diagnostic laboratories are included in the network, that are able to process all together more than 3,000 nasopharyngeal samples per day. Of those, between 200 – 300 SARS-CoV2 positive samples are received weekly for qPCR-VOC analysis and 50% of them were further processed for SARS-CoV-2 genome sequencing. Those numbers are projected to increase with time if exponential growth of cases is sustained. Figure 1A shows the IiWG workflow.
Detection of P.1 variant by VOCs qPCR method and sequencing
In the first two weeks of the IiWG work, a total of 251 SARS-CoV-2 positive RNA samples were obtained from hospitalised cases and outpatient cases collected from all over the country (15 out of 19 Uruguayan departments) between 01/11/2021 and 03/26/2021 (Supplementary Table 1) and processed according to Figure 1A. Patients involved in the study presented Cts for initial diagnostic PCR which ranged from 9 to 34.7, both genders were included (95 males, 95 females and 61 unknown) and age ranged from 1 to 85 years old.
Results from qPCR-VOC assay showed that 67/251 samples (27%) were classified as putative P.1/B.1.351 (Supplementary Figure 1). All VOC positives samples, plus 31 additional samples, were sequenced using the ARCTIC Network protocol for the MinION platform (Oxford Nanopore)[10] as described in Supplementary Methods, to validate PCR-VOC classification. For 74 final consensus sequences with high quality that were assigned to SARS-CoV-2 lineages using Pango nomenclature [11], we achieved a 100% agreement between qPCR-VOC and genome sequencing results. As expected, given the proximity to Brazil, samples classified as P.1/B.1.351 by qPCR-VOC, were assigned to lineage P.1 after genome sequencing. qPCR-VOC is a feasible, precise and scalable method to follow known VOCs in real-time.
Given the qPCR-VOC screening and the confirmation of P.1 by genome sequencing, we established the circulation of lineage P.1 in 15 of 19 departments in Uruguay here analyzed (Figure 1D and Supplementary Figure 1 and 2). To estimate the geographic source and the number of independent introductions of P.1 variant into Uruguay, P.1 Uruguayan sequences obtained in this work were combined with other P.1 sequences from Uruguay (n=3) collected in February (Supplementary Table 1), and with 691 P.1 sequences from South America available in EpiCoV/GISAID [12] (Supplementary Table 2). The ML phylogeographic analysis identified at least 12 independent introductions of lineage P.1 in Uruguay from Brazil and at least six local transmission clusters (designated as UY-I to UY-VI) of between three and 24 sequences (Figure 1B). The median TMRCA of these six Uruguayan P.1 clades was estimated with Bayesian analysis to between mid-February and early March 2021 (Figure 1D, Supplementary Figure 3 and Supplementary Table 3), which coincides with a period of increasing mobility in the country and the beginning of the exponential growth phase (Figure 1C).
Conclusion
The rapid emergence of the SARS-CoV-2 lineage P.1 in the South American region requires the need for increased screening of this highly transmissible virus. Here we developed a comprehensive genomic surveillance program and provide a clear example of how multidisciplinary teams are key mechanisms for helping sanitary authorities to better manage the crisis. Thus, we created an inter-institutional working group with diagnostic laboratories, hospitals (public and private), public and academic institutions across Uruguay that allows real-time VOCs surveillance. Our findings revealed that the VOC P.1 was introduced in Uruguay at multiple times over a period of increasing mobility between mid-February and early March, 2021. The introduction of the highly transmissible VOC P.1 coupled to the increasing human mobility probably contributed to the rapid local spread of this variant and to the COVID-19 epidemic worsening in Uruguay in 2021.
Data Availability
Data was uploaded to GISAID
Funding
This work was supported by FOCEM-Fondo para la Convergencia Estructural del Mercosur (COF03/11).
Supplementary Material File
Supplementary Methods
Supplementary References
Supplementary Tables
Supplementary Figures
Supplementary GISAID
Acknowledgment Table
1. Supplementary Methods
Ethics statement
Residual de-identified RNA samples from positive patients were remitted to the Institut Pasteur Montevideo (IPmon). IPmon was validated by the Ministry of Health of Uruguay as an approved center providing diagnostic testing for COVID-19. All samples were de-identified before receipt by the study investigators.
Development of a “drop out” qPCR SARS-CoV-2 approach to detect VOCs
Primers and “drop out” probes designed to target the deletion sites in both Spike (S) and ORF1ab viral genes were used for this assay. Primers and probes targeting ORF1ab deletion (Δ3675-3677SGF) were designed using the nucleotide sequence from Wuhan strain (GenBank acc. NC_045512) and are shown in Supplementary Table 4. Oligonucleotide complementarity and Tm were assessed with OligoAnalyzer (IDT)[1]. Primers and probes targeting S deletion (Δ69/70 HV) were taken from [2]. B.1.1.7 has both deletions whereas P.1 and B.1.351 have only deletion in ORF1ab.
Hence, we expect no amplification with both probes with B.1.1.7 lineage viruses, only amplification of the S gene with P.1 and B.1.351 lineages and amplification of both genomic regions with other lineages. These primer/probe sets were analyzed with BLAST to rule out similarities with other sequences other than SARS-CoV-2. The detection strategy is based on the standard hydrolysis probe system (TaqMan® Technology). The S gene specific probes are labeled with the HEX fluorophore, while the Orf1a probe is custom labeled with CY5 fluorophore. The primers and probes were synthesized by Integrated DNA Technologies (IDT) and purified by high-pressure liquid chromatography. The OneStep Reaction mix was prepared using 5 µl of RNA samples, 5 µl of 4 x TaqMan Fast Virus 1-Step Master Mix (Thermo Cat. No. 4444436), with 0.3 µM of each S primers, 0.1 µM of S probe, 0.3 µM of each Orf1a primer, 0.1 µM of Orf1a probe and molecular biology grade water to a final volume of 20 µl. The RT-qPCR assays were performed using QuantStudio™ 7 (Applied Biosystem).
Genome Sequencing
Sequencing libraries were prepared according to the Eco PCR tiling of COVID-19 virus protocol (Oxford Nanopore Technologies, United Kingdom), based on the method described by Quick J. [3, 4] with some modifications. RNA samples previously screened by a qPCR assay, were reverse transcribed using SuperScript™ II Reverse Transcriptase (Thermo Fisher Scientific Inc., MA, USA) or the LunaScript® RT SuperMix Kit (New England Biolabs, Ipswich, MA, USA). A negative control was included at this point and carried throughout the protocol. V3 Artic Network primers (IDT) were used for SARS-CoV-2 genome amplification using the Q5® High-Fidelity DNA Polymerase (NEB) in two multiplex PCR reactions. At this point a positive control from a previously sequenced sample was added and carried throughout the protocol. Both reactions were pooled and diluted for the end prep reaction using the NEBNext Ultra II End repair / dA-tailing Module (NEB). Blunt/TA Ligase master mix (NEB) was used to ligate unique barcodes (Native Barcoding expansion kit, EXP-NBD196, Oxford Nanopore Technologies, United Kingdom) to the end prepped samples. Barcoded samples were pooled, cleaned up using 0.4X volume of AMPure XP beads (Beckman coulter, USA) and quantified using the Qubit HS dsDNA kit (Thermo Fisher Scientific Inc., MA, USA). AMII sequencing adapters were ligated using the NEBNext Quick Ligation Module (New England Biolabs, Ipswich, MA, USA) and washed with SFB buffer (Oxford Nanopore Technologies, United Kingdom). Final library was eluted on EB buffer (Oxford Nanopore Technologies, United Kingdom) and quantified using the Qubit HS dsDNA kit. Approximately 50 fmol was loaded into a FLO-MIN106D R9.4.1 flow cell and sequenced on the GridION X5 sequencing platform (ONT) until a minimum sequencing depth of 500X was achieved.
Basecalling and demultiplexing was performed with Guppy 4.3.2 [5] using the high accuracy mode. Samples were demultiplexed using both front and rear barcodes using the --requiere_both_barcodes option. Consensus genomes were generated and using the poreCov pipeline [6, 7] implemented in Nextflow [8], which includes the Artic Network workflow [9], Pangolin for PANGO lineage [10] assignment, nextclade for clade assignment [11], Kraken2 and Krona for contamination detection and visualization [12, 13], minimap2 for read mapping [14], medaka for consensus generation, and using singularity as container engine [15]. Amplicons that are not sequenced or whose depth is less than 20x are not included in the consensus sequences and these positions are represented by “N” stretches. Complete sequences with up to 15% of Ns were kept for further analysis. All genomes obtained in this study were uploaded to the EpiCoV database in the GISAID initiative under the accession numbers: EPI_ISL_1992241, EPI_ISL_1992100, EPI_ISL_1991974, and EPI_ISL_2031707 to EPI_ISL_2031762 and additional metadata for Uruguayan genomes is presented in Supplementary Table 1.
Phylogenetic and phylogeographic analysis
Uruguayan P.1 genome sequences (n = 59) were manually curated in specific genome positions of interest, such as P.1 synapomorphies, and were analyzed in the context of additional P.1 sequences from South America, downloaded from EpiCoV database in the GISAID initiative [16]. Downloaded sequences (n = 905, see Supplementary GISAID Acknowledgment Table) were complete, high quality, with full collection date information and sampled before March 31, 2021. Alignment was performed with MAFFT v7.471 [17]. After carefully alignment inspection, several sequences were missing the characteristic indels (del11288-11296 and ins28269-28273) of P.1 and were excluded, retaining 691 sequences from Brazil (n = 590), Chile (n = 72), Colombia (n = 4), French Guiana (n = 15), Peru (n = 2) and Suriname (n = 8) (Supplementary Table 2). Maximum likelihood (ML) phylogenetic analysis of P.1 samples was conducted using IQ-TREE version 1.6.12 under the GTR+F+R3 nucleotide substitution model selected by the built-in ModelFinder option [18]. Branch support was assessed by the approximate likelihood-ratio test based on a Shimodaira-Hasegawa-like procedure with 1,000 replicates. The tree was rooted using the sequence EPI_ISL_833137 with collection date 2020-12-04 and visualized using iTOL [19]. ML phylogeographic analysis and time-scaled Bayesian analysis were next performed to infer the geographical source of Uruguayan P.1 samples and to estimate the time of the most recent common ancestors (TMRCA) of the Uruguayan P.1 clades, respectively. The previously generated ML phylogenetic tree was time-scaled using TreeTime [20] applying a fixed clock rate of 8 x 10-4 substitutions/site/year, based on previous estimates [21, 22], a skyline coalescent model with eight grid points and keeping polytomies. The time-scaled tree was then employed for the ancestral character state reconstruction (ACR) of epidemic locations with PastML [23], using the Marginal Posterior Probabilities Approximation (MPPA) method with an F81-like model. A time-scaled Bayesian phylogenetic tree was constructed using the Bayesian Markov Chain Monte Carlo (MCMC) approach implemented in BEAST 1.10 [24] with BEAGLE library v3 to improve computational time. Bayesian analysis was conducted using: the non-parametric Bayesian skyline (BSKL) model as the coalescent tree prior [25], the GTR+I+G model of nucleotide substitution and a strict molecular clock model with a uniform prior on substitution rate (8-10 x 10-4 substitutions/site/year). Two MCMC chains were run for 50 million generations and then combined to ensure stationarity and good mixing. Convergence (Effective Sample Size > 200) in parameter estimates was assessed using TRACER v1.7 [26]. The maximum clade credibility (MCC) tree was summarized with TreeAnnotator v1.10 [27] and visualized using FigTree v1.4.4 [28].
3. Supplementary Tables
4. Supplementary Figures
5. Supplementary GISAID Acknowledgment Table
We gratefully acknowledge the following Authors from the Originating laboratories responsible for obtaining the specimens, as well as the Submitting laboratories where the genome data were generated and shared via GISAID, on which this research is based.
All Submitters of data may be contacted directly via www.gisaid.org
Authors are sorted alphabetically.
Acknowledgements
Thanks for efforts from different groups that contribute with SARS-CoV-2 genomes to the EpiCoV GISAID initiative (Supplementary GISAID Acknowledgment Table). We thank Marcelo Fiori and María Inés Fariello for sharing the Uruguayan mobility data and Christian Brandt from Institute for Infectious Diseases and Infection Control (Jena University Hospital), for adapting their poreCov Nextflow pipeline to our requirements and to quickly clarify doubts about its implementation.
We also want to thank Dr. Nicolas Nin and Javier Hurtado from Hospital Español.