Spread of Gamma (P.1) Variants with Spike Protein Changes Near Key Activation Site and Deletions in the Front Region Drives Ongoing COVID-19 Transmission in Amazonas, Brazil

To better understand the recent evolution of VOC gamma spreading in the Amazonas state, we sequenced the virus genome from 1,137 patients representing a random sample of SARS-CoV-2 positive cases diagnosed from 1st January to 6th July 2021. The Amazonas state health surveillance foundation sent SARS-CoV-2 positive samples from different municipalities for sequencing at FIOCRUZ Amazônia, part of the local health genomics network (REGESAM) and the consortium FIOCRUZ COVID-19 Genomics Surveillance Network of the Brazilian Ministry of Health (http://www.genomahcov.fiocruz.br/↗). Furthermore, we also sequenced the virus from 51 individuals living in Manaus who were part of a cohort of patients receiving the CoronaVac vaccine (CovacManaus https://www.ipccb.org/covacmanaus↗). This resulted in a total of 1,188 SARS-CoV-2 high-quality (<1% of undetermined "N" bases) whole-genome (>98% coverage) sequences, representing 0.6% of all laboratory-confirmed SARS-CoV-2 cases in the Amazonas state from 1st January to 6th July 2021 (n = 202,773) (Fig. 1A).

Our genomic survey confirms that VOC Gamma was the most prevalent lineage representing 99.7% (1,185 of 1,188 genomes) of all samples sequenced in the Amazonas state across the study period. We identified over 127 distinct amino acid (AA) substitutions and deletions in the S protein in addition to those that define lineage P.1, most of them present in less than 10 samples (Table S1). AA deletions covering in the NTD region (Δ144, Δ141-144 and Δ138-143) and substitutions at the S1/S2 junction (N679K, P681H and P681R), however, were particularly prevalent and sharply increase from January to May 2021 (Fig. 1B). The P.1+NTDdel variants increased from 1.6% in January to 12.4% in May 2021, while the P.1+N679K and P.1+P681H variants increased from 0% to 25.4% and 17.6%, respectively, in the same period. During June-July 2021, the proportion of P.1+N679K genomes continued to increase up to 76.9%, the relative frequency of P.1+P681H genomes increased up to 37.1% in June and then decrease to 13.5% in July, while the frequency of variants P.1+NTDdel dropped to 2.1%. Inspection of P.1 sequences available at EpiCoV database in the GISAID (https://www.gisaid.org/↗ 4) on July 22nd, 2021, also detected an increased frequency of variants P.1+NTDdel, P.1+N679K, P.1+P681H/R in other Brazilian states, but in a much lower prevalence than in the Amazonas (Fig. 1C).

To determine if the most frequent P.1 mutations detected in Brazil resulted from independent convergent mutations events, we combined the Amazonian P.1 sequences generated in this and previous studies with P.1+NTDdel (Δ144, Δ143-144, and Δ141-144), P.1+N679K and P.1+P681H/R sequences detected in other Brazilian states that were available at the EpiCoV database in GISAID (https://www.gisaid.org/↗). The Maximum Likelihood (ML) phylogenetic analysis supports that NTD deletions around position Y144 arose multiple (n > 20) times during the evolution of lineage P.1 in Brazil, in agreement with our previous observations (7), as well as mutations S:N679K (n > 5), S:P681H (n > 4) and S:P681R (n > 3) (Fig. 2). This analysis further revealed four well-supported (aLRT = 76–99%) monophyletic P.1+ sub-clades in Amazonas; some of which received the Pango Lineage designation P.1.3 to P.1.6 (Fig. 2). Lineage P.1.3 comprises most P.1+Δ141-144 sequences from the Amazonas state (n = 29/33, 88%). Lineage P.1.4 comprises most P.1+N679K sequences from the Amazonas state (n = 187/197, 95%) and three P.1+N679K sequences from Rio de Janeiro. Lineage P.1.5 comprises the remaining P.1+N679K sequences from the Amazonas state (n = 10/197, 5%) and two P.1+N679K sequences from Roraima and São Paulo states. Lineage P.1.6 comprises most P.1+P681H sequences from the Amazonas state (n = 208/209, 99%) and one P.1+P681H sequence from Rio de Janeiro. The fourth lineage, designated P.1+Δ144_AM, comprises about half of P.1+Δ144 sequences from Amazonas (n = 12/29, 41%). Our analyses also revealed two major well-supported (aLRT = 76–99%) P.1+ lineages designated P.1.7 and P.1.8 that comprises most P.1+P681H (n = 227/234, 97%) and P.1+P681R (n = 13/20, 65%) sequences detected outside the Amazonas state, respectively. Analysis of the mutational profile reveals a variable number of lineage-defining mutations ranging from one in P.1.4 to nine in P.1.8 (Table 1 and Fig. S1). Interestingly, lineages P.1.4 and P.1.5 displayed the same lineage-defining AA substitution (N679K) but different nucleotide mutations (T23599G and T23599A). Most P.1+ sub-clades displayed only one AA lineage-defining mutation in the S protein, except P.1.8 that displayed three mutations (T470N, P681R, and C1235F).

The analysis of temporal structure revealed that the overall divergence of the new P.1+ lineages is comparable with the divergence of other contemporaneous P.1 genomes circulating in the Amazonas, thus supporting a homogenous evolutionary rate among lineages (Fig. 3A). Bayesian time-scaled reconstruction estimated the emergence of the lineages P.1.3 at 7th Mar 2021 (95% HPD: 21st Feb – 12th Mar 2021), P.1.4 at 13th Feb 2021 (95% HPD: 9th Jan – 10th Mar 2021), P.1.5 at 15th Mar 2021 (95% HPD: 21st Feb – 30th Mar 2021), P.1.6 at 3rd March 2021 (95% HPD: 3rd Feb – 24th Mar 2021), P.1.7 on 16th Feb 2021 (95% HPD: 19th Jan – 5th Mar 2021), P.1.8 at 11th May (95% HPD: 1st May – 19th May 2021), and P.1+Δ144_AM at 21st Apr 2021 (95% HPD: 6th Apr – 28th Apr 2021) (Fig. 3B to E). Of note, the estimated median tMRCA of most P.1+ lineages circulating in Amazonas (13th Feb - 15th Mar) coincided with the declining phase of the second COVID-19 epidemic wave in the state.

As expected, phylogeographic reconstructions traced the Amazonas state as the most probable source location (PSP = 1) of lineages P.1.3, P.1.4, P.1.5, P.1.6, and P.1 +Δ144_AM. These P.1 sublineages remained mostly restricted to the Amazonas state, except for a few sporadic disseminations to the Rio de Janeiro (lineages P.1.4 and P.1.6), Roraima (lineage P.1.5) and São Paulo (lineage P.1.5) states (Fig. 3B, C, and E). The state of São Paulo was pointed as the most probable (PSP = 0.68) epicenter of lineage P.1.7 dissemination to several states from Southern (Paraná, Rio Grande do Sul and Santa Catarina), Southeastern (Rio de Janeiro), Central-Western (Goiás and Mato Grosso do Sul), Northeastern (Alagoas, Maranhão, Paraíba, and Pernambuco), and Northern (Pará and Tocantins) Brazilian regions (Fig. 3D). We found evidence of large local transmission clusters of lineage P.1.7 in Goiás and Maranhão and of secondary disseminations from Goiás to Maranhão, Tocantins, Pará, Paraná and Santa Catarina, and from Maranhão to Amazonas, Pará, Paraná, Goiás and Rio de Janeiro. The state of Rio de Janeiro was pointed as the most probable epicenter (PSP = 1) of lineage P.1.8. Despite this variant's overall low estimated prevalence, we found evidence of its dissemination to the Amazonas and Santa Catarina states and of subsequent local transmission in Santa Catarina (Fig. 3E).

The increasing relative frequency of lineages P.1.4 and P.1.6 suggest that those variants displayed an effective reproduction number (Re) >1 in the Amazonas state between March and July 2021. To confirm this hypothesis, we used the birth-death skyline (BDSKY) model to estimate the Re of lineages P.1, P.1.4, and P.1.6. Although the high posterior density (HPD) intervals were relatively large, the trajectory of the median Re values was consistent with each lineage's relative prevalence changes (Fig. 4). Lineage P.1.4 displayed a roughly constant median Re > 1 (1.08–1.12) overall the study period; while lineage P.1.6 displayed a median Re > 1 (1.09–1.25) from February to early June, but decreased to 0.91 in June-July, coinciding with the drop of its relative frequency from 37.1% in June to 13.5% in early July. Lineage P.1 displayed a median Re > 1 (1.96) during November-December 2020, a Re ∼1 during January-March 2021, and a Re < 1 (0.77 to 0.95) from April to July 2021. In June-July 2021, the ratio Re(P.1.4)/Re(P.1) was 1.41–1.42 and the ratio Re(P.1.6)/Re(P.1) was 1.17–1.38.

The distribution of real-time RT–PCR cycle threshold (Ct) scores from single or successive cross-sectional samples is a consistent proxy for virus load in the underlying population and may be shaped by changes in both the infecting virus variant and the epidemic trajectory (8). To test if individuals infected with P.1+ variants might have a higher virus load in the upper respiratory tract than those infected with P.1, we compared the Ct values among SARS-CoV-2 cases diagnosed in Amazonas during the endemic phase of transmission (March to July 2021) when the number of SARS-CoV-2 cases remained roughly stable. Our analyses revealed that as the relative prevalence of P.1+ lineages increase, the mean Ct value of SARS-CoV-2 positive cases progressively reduced from 26.8 (95% CI:26.5–27.0) in March to 24.4 (95% CI:24.2–24.7) in July, which correspond to an ∼5-fold increase in the mean viral load of infected subjects (Fig. 5A). Median Ct values in June-July 2021 were significantly (P < 0.0001) lower than in March-May 2021 (Fig. 5B).

Furin is an endoproteinase whose substrate is a polybasic structural motif of the type R-X-K/R-R↓ (where ↓ represents the peptide bond to be cleaved) (9, 10). In the S protein of SARS-CoV-2, furin recognizes the motif 681-PRRAR↓S-686 which is found in a long loop region between residues 675–692. The S1/S2 furin cleavage site (681PRRAR|S686) is of great importance for SARS-CoV-2 Spike protein structural changes during the initial steps of host cell infection and influence the invasion success and transmission of the virus (11, 12). The specificity of furin for a polybasic substrate is partly due to the high negative charge distributed close to the active site (Fig. 6A) and mutations S:P681H/R, that exchange neutral/nonpolar residues with basic amino acids close to the cleavage site, enhance the S1/S2 cleavability (13–17). We hypothesize that mutation S:N679K may also benefit the enzyme-substrate coupling.

To test this hypothesis, metadynamics simulations that modeled the dissociation of the substrate-enzyme complex were used to estimate the relative binding affinity between the wild-type (WT)/mutant peptides mimicking the S protein loop and the furin enzyme. The dissociation was freely conducted along the CVs-defined pathway, providing the potential of mean force (PMF) barrier between the WT complex and the pre-complex states. For metadynamics simulations in which the employed CV is protein – peptide distance, the PMF is equivalent to the free energy of binding (ΔG). The free energy landscape (FEL) was described through the PMF, depicting the phase space sampled from the CVs. It is worth noting that these values may differ quantitatively if the full S protein head was used. Additionally, the calculated values refer to energies obtained within the determined Gaussian parameters to the set of compounds under study. Thus, these energies should be used as a measure to compare the binding free energies among them, and not as absolute free energies. As it can be seen, the three PMFs display a well-defined energy minimum and metastable states during the dissociation process (Fig. 6B). From the PMF it is inferred that all mutants present more favorable binding energy, and therefore, higher affinity, compared to the WT peptide. Our analyses predict that replacement of a nonpolar residue (P) by a polar one (H) provides a modest increase in affinity compared to the WT peptide, whereas having a positively charged residue at either position 679 (K) or 681 (R) leads to a dramatic enhancement in binding affinity to the furin binding site (Fig. 6B).

The expansion of P.1+ lineages in the Amazonas state may be partly associated with their higher ability to infect individuals who acquired immunity through vaccinations. To test this hypothesis, we compared the frequency of distinct P.1 variants in 38 fully vaccinated breakthrough cases (documented infection occurring 14 days after the second dose) and 104 unvaccinated individuals diagnosed in the Amazonas between April and July 2021. Our analyses revealed that P.1+ lineages were not overrepresented among breakthrough SARS-CoV-2 infections (71%) compared with unvaccinated controls (93%) infected in the same period (Fig. 7). Significant difference (P < 0.0001) in the overall distribution of P.1+ lineages among groups was mainly explained by the much higher frequency of lineage P.1.6 in the unvaccinated (63%) than in the fully vaccinated (32%) group. Although this may indicate a higher susceptibility of unvaccinated individuals to infection with lineage P.1.6 with respect to the fully vaccinated ones, this is a limited sample size and observations should be interpreted with caution.

We analyzed nasopharyngeal and pharyngeal swabs (NPS) collected between 01st January and 06th July 2021 from residents in the Amazonas state positively tested by real-time RT-PCR as a routine diagnostic for SARS-CoV-2. In total, 1,188 samples were submitted to nucleotide sequencing at the Fiocruz/ILMD under the auspices of the Fiocruz COVID-19 Genomic Surveillance Network, the Amazonas State Health Foundation – Dra. Rosemary Costa Pinto (FVS-RCP), and the Brazilian Ministry of Health. This study was approved by the Ethics Committee of the Amazonas State University, which waived signed informed consent (CAAE:25430719.6.0000.5016).

The SARS-CoV-2 genomes were recovered using a previously described sequencing protocol (33) or a commercial kit Illumina COVIDSeq Test (Illumina), with the following modifications in the original manufacturer instructions. The reverse transcription step was optimized to use half of the volume described in the Illumina COVIDSEQ Test protocol. We also included a longer incubation time (65°C for 5 min) and increasing sequential temperature steps 25°C 15'/37°C 15'/45°C 15'/50°C 15'/70°C 15' and 4°C/∞. With these modifications we were able to sequence samples with Ct higher than 25, around 28–30. We also included custom primers at a final concentration of 20 nM to minimize the dropout observed in the primers target regions when only standard Illumina COVIDSeq Test primer set was used. (COVIDSEQ_3732_FNF 5′ - GTTGTTAATGCAGCCAATGTTTACCTTAAA, COVIDSEQ_4186_FNR CAACTTGCTTTTCACTCTTCATTTCCAAA, COVIDSEQ_20883_FNF 5′-TGCTAATTCCATTGTTTGTAGATTTGACACTA, COVIDSEQ_21285_FNR 5′ - CTGAAGTCTTGTAAAAGTGTTCCAGAGG, COVIDSEQ_24259_FNF 5′ - AACATCACTAGGTTTCAAACTTTACTTGCT, COVIDSEQ_24604_FNR 5′ - ATGCAAATCTGGTGGCGTTAAAAAC).

The new generation sequencing libraries were clustered with MiSeq reagent kit v3 (600-cycles) on 2 × 150 cycles paired-end runs (Illumina) or with NextSeq 1000 on 2 × 50 cycles. FASTQ reads were generated by the Illumina pipeline at BaseSpace (https://basespace.illumina.com↗). Consensus sequences were generated using DRAGEN COVID LINEAGE 3.5.1 to 3.5.3, according to the most up-to-date version of this app on each sequencing run. Subsequently, we evaluated the consensus files for quality using Nextclade tool v1.5.2 (https://clades.nextstrain.org↗), those with more than 1% of ambiguities "Ns" had the FASTQ files imported into Geneious Prime 2021 for trimming and assembling using a customized workflow employing BBDuk and BBMap tools (v38.84) and the NC_045512.2 RefSeq as a template with carefully visually inspection. Using both approaches, we generated consensus sequences with mean depth coverage higher than 800X, excluding duplicate reads. Whole-genome SARS-CoV-2 consensus sequences were initially assigned to viral lineages according to the nomenclature proposed by Rambaut et al. (34), using the Pango Lineage web application (https://pangolin.cog-uk.io↗) and later confirmed using phylogenetic analyses as explained below.

SARS-CoV-2 P.1 sequences obtained here were aligned with high quality (<1% of N) and complete (>29 kb) lineage P.1 Amazonian sequences that were available in the EpiCoV database in the GISAID (https://www.gisaid.org/↗) on May 31st, 2021. Sequences were aligned using MAFFT v7.475 (35) and then subjected to maximum-likelihood (ML) phylogenetic analysis using IQ-TREE v2.1.2 (36) under the general time-reversible (GTR) model of nucleotide substitution with a gamma-distributed rate variation among sites, four rate categories (G4), a proportion of invariable sites (I) and empirical base frequencies (F) nucleotide substitution model, as selected by the ModelFinder application (37). The branch support was assessed by the approximate likelihood-ratio test based on the Shimodaira–Hasegawa-like procedure (SH-aLRT) with 1,000 replicates. The temporal signal of the P.1 and P.1+ sequences was assessed from the ML tree by performing a regression analysis of the root-to-tip divergence against sampling time using TempEst (38).

We performed a time-scaled Bayesian phylogeographic analysis for SARS-CoV-2 genomes of lineages P.1.3 to P.1.8 plus P.1+Δ144 sampled in Brazil using the Bayesian Markov Chain Monte Carlo (MCMC) approach implemented in BEAST 1.10.4 (39) with BEAGLE library v3 (40) to improve computational time. The Bayesian tree was reconstructed using the GTR+F+I+G4 nucleotide substitution model, the nonparametric Bayesian skyline (BSKL) model as the coalescent tree prior (41), a strict molecular clock model with a uniform substitution rate prior (8 × 10⁻⁴ substitutions/site/year) and a reversible discrete phylogeographic model (42) with a continuous-time Markov chain (CTMC) rate reference prior (43). MCMC was run sufficiently long to ensure convergence (effective sample size [ESS] > 200) in all parameters estimates as assessed in TRACER v1.7 (44). The maximum clade credibility (MCC) tree was summarized with TreeAnnotator v1.10 and visualized using FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/↗).

To estimate the Re trajectories of lineages P.1, P.1.4, and P.1.6 through time in the Amazonas, we used the birth-death skyline (BDSKY) model (45) implemented within BEAST 2 v2.6.5 (46). The sampling rate (d) was set to zero for the period before the oldest sample and estimated from the data afterward. The BDSKY prior settings were as follows: Become Uninfectious Rate (exponential, mean = 36); Reproductive Number (log-normal, mean = 0.8, sd = 0.5); Sampling Proportion (beta, alpha = 1, beta = 100). Origin parameter was conditioned to the root height, and the Re was estimated piecewise over five-time intervals. The molecular clock was as in the time-scaled trees analysis, and the HKY+G4+F substitution model was used. MCMC chains were run until all relevant parameters reached ESS > 200, as explained above.

The furin cleavage site in the Spike protein of SARS-CoV-2 is found in a disordered region, exposed to the solvent. Due to the high flexibility of this loop, the experimentally resolved structures lack atomic coordinates in this region, even when mutations are introduced to express the protein in its pre-fusion conformation. Starting from the crystallographic structure of an inhibitor bound to human furin (PDB ID: 6HLB) (47), the 679-NSPRRARS-686 loop region was built analogously in the enzyme-interacting conformation (Fig. S2A). To mimic the complete loop linked to the Spike protein, the modeled motif was used as a folding seed (rotamers kept fixed) and the remaining loop residues were modeled using the Rosetta Remodel protocol (48) (Fig. S2B). The final model comprising the residues 675-QTQTNSPRRARSVASQSI-692 (Fig. S2C) was minimized using the FastRelax protocol of Rosetta (49) using the REF2015 energy function and maintaining the catalytic residues of furin fixed. Constraints to Ca2+ ion and its coordinating residues were applied with SetupMetalsMover (50). The variants N679K, P681H and P681R loops were built by using Point mutant ("pmut") scan application of Rosetta onto a modeled native loop (51). The resultant mutant structures were then geometry-optimized following the same protocol with the FastRelax application. The relative binding free energies between furin and the peptides mimicking the Spike furin-loop variants were estimated by metadynamics simulations in an explicitly aqueous environment (52) (as detailed in the Supplemental Material).

Descriptive statistics, tests for normal distribution (D'Agostino & Pearson and Anderson-Darling), and the nonparametric Mann-Whitney test were used to compare the Ct of SARS-CoV-2 RT-PCR positive samples from the upper respiratory tract of patients over time. Only Ct values from samples analyzed with the same RT-PCR diagnostic assay were compared. The Chi-Square test was used for testing the association between the vaccination status and frequency of P.1 variants. The threshold for statistical significance was set to P < 0.05. Graphics and statistical analyses were performed using GraphPad v9.02 (Prism Software, United States).

All the SARS-CoV-2 genomes generated and analyzed in this study are available at the EpiCov database in GISAID (https://www.gisaid.org/↗). The list of accession IDs may be found in the attached supplemental material file.