Ongoing evolution of SARS-CoV-2 drives escape from mRNA vaccine-induced humoral immunity

Abstract

Since the COVID-19 pandemic began in 2020, viral sequencing has documented 131 individual mutations in the viral spike protein across 48 named variants. To determine the ability of vaccine-mediated humoral immunity to keep pace with continued SARS-CoV-2 evolution, we assessed the neutralization potency of sera from 76 vaccine recipients collected after 2 to 6 immunizations against a comprehensive panel of mutations observed during the pandemic. Remarkably, while many individual mutations that emerged between 2020 and 2022 exhibit escape from sera following primary vaccination, few escape boosted sera. However, progressive loss of neutralization was observed across newer variants, irrespective of vaccine doses. Importantly, an updated XBB.1.5 booster significantly increased titers against newer variants but not JN.1. These findings demonstrate that seasonal boosters improve titers against contemporaneous strains, but novel variants continue to evade updated mRNA vaccines, demonstrating the need for novel approaches to adequately control SARS-CoV-2 transmission.

Keywords: COVID-19; JN.1; SARS-CoV-2; XBB.1.5 Booster; breadth; infectivity; neutralizing antibodies; spike; vaccination; variants.

Figure 1:. SARS-CoV-2 variants of concern harbor spike protein mutations that are enriched in the RBD.

(A) Scaled stacked area chart of variants emerging during the SARS-CoV-2 pandemic. Arrows denote when each vaccine was made available. Adapted from nextstrain.org. Data as of February 28^th, 2024. (B) Crystal structures of pre-fusion stabilized SARS-CoV-2 spike trimer (PDB ID 6XR8) highlighting mutations found in each strain mapped onto their locations on the wild-type spike trimer. From left to right WT, Delta (B.1.617.2), Omicron (BA.1), BA.5, and JN.1 strains. Mutations highlighted for each strain are from corresponding table S1.

Figure 2:. Neutralization assays with primary vaccine sera revealed regions of vulnerability within the natural mutation landscape.

(A) Schematic of high throughput neutralization assay used in these studies. (B) (Top) Schematic illustrating the time period from which variants were selected and the vaccinee sera to be tested below. (Bottom) Titers that achieve 50% of pseudovirus neutralization (pNT₅₀) are plotted for all individual mutants found in SARS-CoV-2 through Delta (1.617.2) for 24 COVID naive donors who received the primary vaccination regimen (2 doses of Moderna or Pfizer). The solid black line represents geometric mean pNT₅₀ for reference. The following abbreviations are used NTD = N-Terminal Domain, RBD = Receptor Binding Domain, S1 = S1 Subunit, S2 = S2 subunit. pNT₅₀ of each spike mutant was compared to pNT₅₀ of wild type using a One Way ANOVA with Friedman’s test and Dunn’s multiple comparisons test. (* = P<0.0332, ** = P<0.0021, ***= P<0.0002, ****= P<0.0001, red = significantly higher than WT). (C) Two-way hierarchical clustering of pNT₅₀ values of sera obtained from vaccinated donors (rows) across each variant (columns) relative to individual donor WT titer. pNT₅₀ values are plotted in a heat map colored according to neutralizing activity relative to WT. Donors with WT neutralizing titers at the maximum of our assay were excluded. Mutations where a significant majority of donor titers reached the limit of detection were also excluded. Clustering was performed using pheatmap package v1.0.12 in R-Studio.

Figure 3:. Combinations of individual mutations yield neutralization resistance.

(A-B) Primary series vaccinee sera was tested via pseudovirus-based neutralization assays against each individual mutation or combination of mutations found in the (A) (Left) Beta strain or (B) (Left) Delta strain. Each variant was compared to the wild type using a One-way ANOVA with Friedman’s Test and Dunn’s multiple comparisons test. (* = P<0.0332, ** = P<0.0021, ***= P<0.0002, ****= P<0.0001). Primary series vaccinee sera was tested against combinations of (A) (Right) Beta or (B) (Right) Delta RBD mutant pseudoviruses and compared to the complete set of mutations for each strain. Each group was compared to wild type using a One-way ANOVA with Friedman’s Test and Dunn’s multiple comparisons test. (* = P<0.0332, ** = P<0.0021, ***= P<0.0002, ****= P<0.0001).

Figure 4:. Individual Omicron variants harbor novel sites of vulnerability, while recent strains exhibit significant escape from boosted sera.

(A) (Top) Schematic illustrating the time period from which variants were selected and the vaccinee sera to be tested below. (Bottom) Primary or Boosted pNT₅₀ for donor across each variant was normalized to geometric mean titer across all donors for WT and log transformed. Length of each vertical line illustrates pNT₅₀ change. (B) Two-way hierarchical clustering of pNT₅₀ values of sera obtained from boosted donors (rows) across each variant (columns) relative to individual donor WT titer. pNT₅₀ values are plotted in a heat map colored according to neutralizing activity relative to WT. Donors with WT neutralizing titers at the maximum of our assay were excluded. Mutations where a significant majority of donor titers reached the limit of detection were also excluded. Clustering was performed using pheatmap package v1.0.12 in R-Studio. (C) Log-transformed pNT₅₀ geometric mean values for each variant relative to WT for primary series were plotted against matching boosted values. Linear regression (R² = 0.5109; Slope = −0.098; P<0.001) (D) Log-transformed pNT₅₀ geometric mean values for each variant relative to WT colored according to their location in the spike protein comparing primary to boosted vaccinee sera. Statistics represent Wilcoxon two-tailed paired, non-parametric t-test with **** = P<0.0001. Selected RBD mutations from the delta and beta strains connected with lines showing their respective values for primary or boosted donor serum.

Figure 5:. Individual Omicron variants reveal new sites of resistance that are observed in more recent strains to significantly escape from boosted sera.

(A) (Top) Schematic illustrating the time period from which variants were selected and the vaccinee sera to be tested below. (Bottom) pNT₅₀ are plotted for all individual spike mutations occurring after Delta (1.617.2) through XBB.1.16 for 24 COVID naive donors that received primary vaccination (2 doses of Moderna or Pfizer) and a booster shot. The solid black line indicates the geometric mean of the titers against each pseudovirus. The following abbreviations are used NTD = N-Terminal Domain, RBD = Receptor Binding Domain, S1 = S1 Subunit, S2 = S2 subunit. Each pseudovirus was compared to SARS-CoV-2 wild type using a One Way ANOVA with Kruskal Wallis Test and Dunn’s multiple comparisons test (* = P<0.0332, ** = P<0.0021, ***= P<0.0002, ****= P<0.0001, red = significantly higher than WT). (B) Log-transformed geometric mean pNT₅₀ values normalized to WT for mutations that appeared in strains between WT and Delta as compared to those appearing post-Delta to XBB.1.16. Statistics represent Mann-Whitney two-tailed unpaired, non-parametric t-test with ** = P < 0.0021. In this case P = 0.0010. (C) The corresponding spike crystal structures (PDB 6xr8) from (B) with mutations colored according to neutralization resistance to boosted sera. (Left) Highlighting mutations between WT and Delta and (Right) highlighting mutations post-Delta to XBB.1.16. Log transformed values for each mutation in the spike are plotted on a three-color scale with a spectrum of −0.9 to 0.9. (D) Log transformed pNT₅₀ geometric mean values relative to WT for boosted donor sera against all RBD variants in strains that appeared after Delta plotted against the date of first submission to GISAID. (E) pNT₅₀ values for boosted donor sera against variants that appeared post-Delta to XBB.1.16, with Alpha, Beta, and Delta for reference. Strains are plotted in the order that they first appeared. The solid black line indicates the geometric mean of the titers against each pseudovirus. Each pseudovirus was compared to SARS-CoV-2 WT using One Way ANOVA with Kruskal Wallis Test and Dunn’s multiple comparisons test (* = P<0.0332, ** = P<0.0021, ***= P<0.0002, ****= P<0.0001).

Figure 6:. Bivalent and XBB.1.5 boosters improve neutralization against highly infectious variants.

(A) Schematic illustrating the variants and vaccinee sera to be tested below. (B) Geometric mean pNT₅₀ values of longitudinal serum samples over the course of CDC-recommended vaccination schedule against each variant. Strains are plotted in the order in time that they first appear. (C) pNT₅₀ of individual serum samples were normalized to geometric mean titer across all donors for WT and log transformed. Longitudinal samples from donors after WT, Bivalent, or XBB.1.5 boosters are plotted. Length of each vertical line illustrates pNT₅₀ change. (D) Log transformed pNT₅₀ geometric mean values relative to WT for donor sera against variants of concern through JN.1 plotted against the date of first submission to GISAID. Solid lines represent linear regressions for WT Booster (R²=0.63, slope=−0.298, P=0.002), Bivalent Booster (R²=0.57, slope=−0.27, P=0.0044) and XBB.1.5 booster titers (R2=0.14, slope =−0.076, P =0.23). (E) One-way hierarchical clustering of pNT₅₀ values of sera obtained from boosted donors (rows) across each variant of concern (columns) relative to individual donor WT titer. Heatmaps maintaining clustering order of rows across Bivalent and XBB.1.5 boosters are plotted to the right. pNT₅₀ values are plotted in heat maps colored according to neutralizing activity relative to WT. Clustering was performed using pheatmap package v1.0.12 in R-Studio. (F) Pseudovirus infectivity (defined as infectious units per genome copy) was normalized to WT infectivity for each variant spike across three technical replicates (n=3). Pseudovirus infectivity relative to WT was measured for major SARS-CoV-2 variants by calculating fold change in slope of linear regressions fit across known dilutions. Bars and error bars depict mean and standard error of the mean. Each pseudovirus was compared to SARS-CoV-2 wild type using a One Way ANOVA with Holm-Sidak correction for multiple comparisons. (* = P<0.033, ** = P<0.0021, ***= P<0.0002, ****=P<0.0001).