SARS-CoV-2 Variants Associated with Vaccine Breakthrough in the Delaware Valley through Summer 2021

Sampling was carried out from March 2020 to September 2021. Over the course of the study, several waves of infection are evident as increased test positivity rates (Fig. 1, light gray curve) and COVID-19-attributed hospitalizations in the city of Philadelphia (Fig. 1, dark gray curve). Widespread vaccination was introduced in winter 2020–2021, reaching over 70% (in adults 18 years old and older) in Philadelphia by September 2021 (Fig. 1, black curve). As is described below, during this period, SARS-CoV-2 variants alpha and delta, designated VBM/VOC, became the overwhelming majority of genomes identified by our sequence surveys (Fig. 1, green and red curves).

Patient samples (nasal or nasopharyngeal swab and saliva) used for SARS-CoV-2 whole-genome sequence analysis are listed in Table S1. To preserve patient confidentiality, samples were deidentified except for collection date and rationale for collection. Surveillance samples (n = 2621) were defined as those acquired from clinical diagnostic laboratories across the Delaware Valley (n = 2485), from hospitalized subjects (n = 116), and asymptomatic subjects testing positive in a university screening program (n = 20) (33). Vaccine breakthrough cases (n = 159) were included only if detected at least 2 weeks after the final vaccine dose (second dose for Moderna and Pfizer-BioNTech mRNA vaccines or single dose for the Johnson & Johnson adenovirus vector platform) and the subject tested positive by clinical laboratory assay. Data were not available on subject immune responses following vaccination.

As a positive control, S gene target failure samples (n = 172) were also compared. The TaqPath COVID-19 Combo kit by Thermo Fisher Scientific targets three regions of the SARS-CoV-2 genome for viral detection; the ORF1ab, nucleocapsid (N), and spike (S) genes. The region of the S gene interrogated by the assay overlies a characteristic deletion in the alpha variant (del 69–70), so samples containing alpha lineage virus are selectively negative for the spike amplicon, while the other two targets are detected. Samples with these characteristics were targeted for sequencing early during the wave of alpha infections to track the variant; here these spike target gene failures serve as positive controls for the statistical model.

Several sequencing strategies were used to acquire whole-genome sequences. The POLAR protocol with Artic primers and Illumina sequencing was used for most samples (34). Smaller numbers of samples were acquired using the Paragon, Illumina RPIP, and Illumina CovidSeq methods. Viral genome sequences were judged to be high quality and included in the study if 95% of the viral genome was covered by at least 5 reads. In all, 2952 high-quality sequences were generated and analyzed. Viral variants were assigned using Pangolin lineage software (35, 36).

Figure 2A shows the proportions of variants detected in surveillance samples over the course of the study from March 2020 to September 2021. Variants detected are summarized by the color code on the bar graph (variant designations used are summarized in Table S2). The numbers of genomes analyzed per week are indicated above each column. Numbers varied both as surveillance sequencing efforts accelerated and as the availability of samples varied.

Members of the B.1 lineage predominated March 2020 until fall 2020, at which time B.1.2 and B.1.243 became predominant. Later the B.1.526 and B.1.1.7 (alpha) variants emerged, with B.1.1.7 becoming predominant by spring 2021. Delta (B.1.617.2 and AY lineages) became detectable in late spring and predominant in early summer. By late summer, delta lineages were the only variants detected.

Wide-spread vaccination was introduced in midwinter 2020–2021, but nevertheless breakthrough infections were detected in some vaccinated individuals, raising the question of whether specific viral features identifiable in sequence data might be associated with vaccine evasion. We sequenced 159 vaccine breakthrough samples collected between February 22 and September 3, 2021, and compared the distributions of lineages or genomic variations to those observed in surveillance samples from the same community (n = 2621).

Challenges in the analysis include the facts that: (i) the distribution of variants in the surveillance samples is changing over time; (ii) sampling is uneven over time; and 3 sampling is subject to stochastic fluctuations. To address these issues, we developed a model based on Bayesian analysis that combines an autoregressive moving average model with multinomial logistic regression. The underlying variant proportions at each week of the study were estimated from the counts of SARS-CoV-2 variants assuming relatively smooth changes over time and accounting for stochastic fluctuations during sampling (Fig. 2B). These estimated surveillance proportions were compared to the variant counts observed in spike gene target failure samples (Fig. 2C) and vaccine breakthrough samples (Fig. 2D). Since the surveillance population proportions were estimated for each week, the weekly vaccine breakthrough or spike gene target failure variant counts could be compared to the corresponding time-matched surveillance estimates.

Examples of the temporal profiles of the modeled lineage succession are shown in more detail in Fig. 3. One benefit of the model is that uncertainty in the surveillance estimates can be included. For example, note that time points with fewer surveillance samples (lighter gray bars) have larger 95% credible intervals for the proportion estimate (colored shading). From this view, it is evident that several further lineages waxed and waned notably over the sampling period, including B.1.1.434 and B.1.526.

As a control, we assessed the variants associated with spike (S) gene target failure samples (Fig. 2C and 4 and Table S3), where an amplicon overlapping the 69–70 deletion in spike failed selectively. To track the spread and variation of B.1.1.7 during the alpha wave of infection, we sequenced 172 S gene target failure samples in the Delaware Valley from January to April 2021. Upon variant assignment, genomes from 96.1% of spike target gene failure samples corresponded to B.1.1.7. Rarely, other lineages were seen that share the spike deletion with B.1.1.7 (B.1.375 and B.1.525) or were likely stochastic failures of S amplification or genomes with novel combinations of mutations.

As expected, we found that the B.1.1.7 lineage was estimated to be highly enriched in the S gene target failure set over what would be expected from its frequency in the surveillance lineage counts (Fig. 4) (odds ratio of 120; 95% credible interval 38–470). The B.1.525 lineage, which also contains the S69-70 deletion, was estimated to be enriched as well (odds ratio 33; 95% credible interval 1.7–380). Note that this enrichment was detected based on counts of only two B.1.525 genomes in the S gene target failure set and 12 surveillance genomes, both emphasizing the sensitivity of the model and highlighting that the detection of enrichment can be easier in lineages rare in the surveillance population. Other lineages showed little association with spike target gene failures, as expected (note that only a single B.1.375 lineage sample was observed and thus it did not reach the threshold for lineage-specific analysis and was grouped in “Other”). Thus, the analysis of S target failures confirmed that the Bayesian autoregressive moving average categorical regression model was effective at identifying overrepresented lineages relative to the time-adjusted expectation from surveillance.

We then applied the Bayesian model to analyzing vaccine breakthrough samples (Fig. 2D and 5 and Table S4). There was less clear enrichment in vaccine breakthroughs than in the S gene target failure set. The delta lineage B.1.617.2 did show signs of enrichment in breakthrough samples, with mean odds ratio of 3 (95% credible interval 0.89 to 11). The one-sided posterior probability of any enrichment for B.1.617.2 in vaccine breakthrough was estimated at 96% with a 73% probability of more than a 2-fold enrichment in odds of appearing in vaccine breakthrough. This enrichment was not observed for all lineages grouped within delta, with both a group-wise estimate for all delta lineages and the specific estimates for the other delta AY lineages, including AY.3, AY.4, AY.14, and AY.24, showing little enrichment over surveillance. These alternative delta lineages also had fewer samples observed, so a lack of power may be a partial explanation. Note that the power to detect enrichment of delta is limited by its rapid spread–observation of a delta vaccine breakthrough case at a time point when almost all surveillance population infections were also delta provides little information. Similarly, lineages that had already waned by the time of widespread vaccination are difficult to assess since there was little chance for them to appear in vaccine breakthrough cases.

We next investigated whether any individual base substitutions in the viral genome were selectively associated with spike gene target failure and vaccine breakthrough samples. A summary of mutations detected is presented in Table S5. We adapted the Bayesian autoregressive moving average logistic regression method to assess the behavior of single mutations and amino acid substitutions. Fig. 6 and Table S6 summarizes results. Many mutations varied in frequency over the course of the study often paralleling the profile of the most abundant lineages containing them (Fig. 6A).

As expected, the 69–70 deletion was found to be highly enriched in the S gene target failure set (odds ratio of 205, 95% credible interval of 65–581). Enrichment among S gene target failures was also seen for another 28 mutations (Table S6), where all the most strongly enriched substitutions were characteristic of the alpha variant and thus due to “hitchhiker” effects. Underrepresentation in spike target gene failures was also seen in mutations in multiple additional open reading frames, reflecting mutations found in variants lacking the 69–70 deletion. Each mutation was analyzed separately without accounting for genomic linkage so this is as expected.

For vaccine breakthroughs (Fig. 6B and Table S7), the most notably enriched substitution was N501Y, which showed an odds ratio of 2.04 (95% credible interval of 1.25–3.18). The N501Y substitution is found in multiple VBM, including alpha, beta, and gamma, and is reported to increase the affinity of spike protein binding to the ACE2 receptor (23, 37, 38) and to diminish binding of some human antibodies to spike (24). The spike substitution P681H was also slightly enriched; this substitution is near the furin cleavage site and may promote efficient proteolysis. Several additional substitutions in spike were potentially enriched (D614G, P681R, D950N), but the 95% credible interval included one, so the evidence for enrichment was weaker (Table S7). Two other closely studied spike substitutions, E484K and L452R, were not notably enriched among vaccine breakthrough cases, and the D253G substitution was modestly depleted. Among non-spike substitutions, the P199L substitution in the nucleocapsid showed notable depletion (Table S7). Stepping back, these findings together emphasize the potential importance of the N501Y substitution in vaccine breakthrough.

The University of Pennsylvania Institutional Review Board (IRB) reviewed the research protocol and deemed the limited data elements extracted with positive SARS-CoV-2 specimens to be exempt from human subject research per 45 CFR 46.104, category 4 (IRB #848605). For hospitalized subjects at the University of Pennsylvania, following informed consent (IRB protocol #823392), patients were sampled by collection of saliva, oropharyngeal and/or nasopharyngeal swabs, or endotracheal aspirates if intubated, as previously described (32). Clinical data were extracted from the electronic medical record. Further samples were collected from asymptomatic subjects detected in a screening program at the Perelman School of Medicine at the University of Pennsylvania and symptomatic subjects tested throughout the PennMedicine clinical network under IRB protocols #843565 and #848608. Human samples were collected at Children’s Hospital of Philadelphia under protocol # 21-018726 approved by the Children’s Hospital of Philadelphia IRB. Human samples were collected at Thomas Jefferson University under protocol number IRB Control # 21E.441 approved by the Thomas Jefferson University IRB. Vaccine breakthrough cases were identified by either report to HUP Infection and Control, Jefferson Infection and Control, CHOP Department of Infection Control and Prevention, the Philadelphia Department of Public Health, or, where applicable, chart review.

Several sequencing methods were used to acquire viral whole-genome sequences.

The POLAR protocol was used to acquire the majority of viral genome sequences (47). Illumina’s NextSeq instrument was used to gather sequence data. In detail, 5 μl of viral RNA, 0.5 μl of 50 μM Random Hexamers (Thermo Fisher, N8080127), 0.5 μl of 10 mM dNTPs Mix (Thermo Fisher, 18427013), and 1 μl nuclease-free water was incubated at for 5 min at 65°C proceeded by a 1 min incubation at 4°C. To perform reverse transcription, 6.5 μl from the previous reaction, 0.5 μl SuperScript III Reverse Transcriptase (Thermo Fisher, 18080085), 0.5 μl of 0.1M DTT (Thermo Fisher, 18080085), 0.5 μl of RNaseOut (Thermo Fisher, 18080051), and 2 μl of 5X First-Strand Buffer (Thermo Fisher, 18080085) was incubated for 50 min at 42°C, followed by an incubation for 10 min at 70°C, and then held at 4°C. To amplify the cDNA, artic-ncov2019 version 3 primers were used (IDT). To perform PCR amplification of the viral cDNA, the following reagents were added to 2.5 μl of the previous mixture: 0.25 μl Q5 Hot Start DNA polymerase (NEB, M0493S), 5 μl of 5X Q5 Reaction Buffer (NEB, M0493S), 0.5 μl of 10 mM dNTPs Mix (NEB, N0447S), either 4.0 μl of pooled primer set 1 or 3.98 μl of pooled primer set 2, and nuclease-free water to bring to a final volume of 25 μl. This PCR amplification of the viral cDNA used the following conditions: 98°C for 30 s for 1 cycle, 25 cycles at 98°C for 15 s and 65°C for 5 min, and then held at 4°C. Amplicons generated by the two primer sets from the same sample were pooled then diluted to a concentration of 0.25 ng/μl. The Nextera library was prepared using the Nextera XT Library Preparation kit (Illumina, FC-131-1096) and the IDT for Illumina DNA/RNA UD Indexes Set A and B (Illumina, 20027213, 20027214, 20027215, 20027216). The Quant-iT PicoGreen dsDNA quantitation assay kit was used to quantify the DNA of each sample (Invitrogen, P7589). The samples were pooled in equal quantities, and the pooled library was quantified using the Qubit1X dsDNA HS assay kit (Invitrogen, Q33230). The library was sequenced on an Illumina NextSeq.

Another sequencing method, used at CHOP, was Paragon. Paragon sequencing was carried out as in (22). Briefly, RNA was extracted from nasopharyngeal swab samples using QIAamp Viral RNA Mini (Qiagen). Whole-genome sequencing was carried out by the Genomics Core Facility at Drexel University. Amplification was performed using Paragon Genomics CleanPlex SARS-CoV-2 Research and Surveillance NGS Panel 1 and 2. Libraries were quantified using the Qubit dsDNA HS (High Sensitivity) assay kit (Invitrogen) with the Qubit Fluorometer (Invitrogen). Library quality was assessed using Agilent High Sensitivity DNA kit and the 2100 Bioanalyzer instrument (Agilent). Libraries were normalized to 5 nM and pooled in equimolar concentrations. The resulting pool was quantified again using the Qubit dsDNA HS (High Sensitivity) assay kit (Invitrogen) and diluted to a final concentration of 4 nM; libraries were denatured and diluted according to Illumina protocols and loaded on the MiSeq at 10pM. Paired-end and dual-indexed 2x150bp sequencing was carried out using MiSeq Reagent Kits v3 (300 cycles).

The Thomas Jefferson University site carried out sequencing using the Illumina RPIP and CovidSeq methods essentially as per the manufacturer’s instructions. A nasopharyngeal swab specimen that was tested positive by PCR for SARS-CoV-2 was used for this study. Vaccine breakthrough specimens were identified by Jefferson Occupational Health Network for Employees and Students (JOHN) and were sequenced with a designation of VBT. Randomly selected residual positive specimens from the Molecular & Genomic Pathology Laboratory at Jefferson during the same period were sequenced for the purpose of epidemiology surveillance and labeled as SURV. RNA was extracted from 200 μl of the specimen and eluted in 110 μl using the bioMérieux EasyMag Extraction System, following EasyMag’s Generic protocol. Whole-genome sequencing for SARS-CoV-2 was subsequently performed at the Molecular & Genomic Pathology Laboratory at Thomas Jefferson University Hospital. With the exception of a few samples, samples were sequenced using Illumina COVIDSeq protocol (Illumina) following manufacturer’s guideline. Remaining specimens were sequenced using the Swift Normalase Amplicon Panels (SNAP) Core kit along with the SARS-CoV-2 Additional Genome Coverage Primer Panel following the manufacturer’s protocol. For data analysis, Illumina’s Local Run Manager’s GenerateFASTQ module was used to generate the fastq files for all specimens. The fastq files were transferred to UPenn’s secure data server for further processing to obtain QC metrics and lineage data.

Samples with VSP numbers lower than VSP00256 (Table S1) were previously described in Everett et al. (32).

To process sequence data, sequence reads are trimmed to remove low quality base calls (< Q20) and aligned to the original Wuhan reference sequence (NC_045512.2) with the BWA aligner tool (v0.7.17) (48), after which alignments are filtered with the Samtools package (v1.10) (49). Sequencing depth is determined for each position in the viral genome and genomes are accepted for analysis when ≥ 95% of genome positions are covered with a read depth of ≥ 5 reads. Variant positions are called with the Bcftools package (v1.10.2-34) (50) requiring PHRED scores ≥ 20 and variant read frequencies ≥ 50% of the total reads. The nature of substitutions is determined by retrieving reading frames from the reference GENBANK record, translating, and determining the native and mutant residues.

Variants were assigned using the Pangolin lineage software (Pangolin version 3.1.11 with the PangoLEARN 2021-08-24 release). Note that these lineages are updated regularly and so are expected to change over time. Point mutations were assigned using a previously described bioinformatics pipeline (31).

Statistical analysis of variant enrichment and mutation enrichment in subsets of the data were separately assessed using a Bayesian model and Markov chain Monte Carlo sampling implemented in Stan (51). The model takes the vector of counts of the variants or mutations seen in the surveillance sampling for a given week, counts_i, and assumes they are multinomially distributed where:with probabilities modeled as multinomial logistic:where p_i,w is the true proportion of variant or mutation i on week w. p_*,w indicates the vector of all lineage probabilities for week w and n is the total number of variants or mutations observed. The underlying proportions for a given week are assumed to be centered around the proportions observed in the prior week (autoregressive) plus a change term that is itself correlated with the change observed in the prior week (moving average):where change_i,w is the change in log odds of variant or mutation i on week w. The initial starting proportions are given flat priors:

The standard deviation for changes, σ, was given a Gamma(1, 2) prior distribution.

To assess if there's enrichment of a particular variant or mutation in the vaccine breakthrough, counts of vaccine breakthrough samples for a given week, vaccineCounts_w, were modeled as:where:andwhere β_i measures the log odds enrichment/depletion of variant or mutation i in the vaccine breakthrough population over the surveillance population, group_i indicates the WHO classification for Pango lineage i and δ_groupi indicates the enrichment or depletion for a CDC VBM/VOC classification containing more than one Pango lineage, e.g., Delta containing B.1.617.2 and other variants with an AY prefix. In these data, Delta, Gamma and Non-VBM/VOC groupings contained more than one lineage above the abundance threshold. For WHO classifications containing only a single lineage, δ_groupi was set to 0. The β and remaining δ were given DoubleExponential(0,1) priors.

This was repeated equivalently in the samples collected as S gene target failures, with the counts, SFailCounts_w, modeled as:where:and

The α were given DoubleExponential(0,1) priors.

For the assessment of lineage enrichment, genomes from lineages that had more than 10 genomes assigned to them were included as their individual lineages, genomes from lineages with 10 or fewer assignments that were listed within VBM/VOC were included as miscellaneous classification for their particular WHO category, e.g., “Other delta” and all other genomes were grouped into an Other category. For the assessment of mutations, each mutation that was found in greater than 5% of genomes was independently modeled as above but replacing the Multinomial distributions with Binomial and removing the inapplicable lineage group δ terms.

All viral genome sequences acquired in this study have been deposited in GISAID and at NCBI under accession numbers listed in Table S1. Analysis code is archived on Zenodo at https://doi.org/10.5281/zenodo.5888338. Sequence processing software and intermediate files used in this study are available at https://doi.org/10.5281/zenodo.5559699. A list of key reagents used in this study is in Table S8.