Abstract
Background The possibility of association between SARS-CoV-2 genomic variation and immune evasion is not known among persons with Omicron variant SARS-CoV-2 infection.
Methods In a retrospective cohort, using Poisson regression adjusting for sociodemographic variables and month of infection, we examined associations between individual non-lineage defining mutations and SARS-CoV-2 immunity status, defined as a) no prior recorded infection, b) not vaccinated but with at least one prior recorded infection, c) complete primary series vaccination, and/or d) primary series vaccination and ≥ 1 booster. We identified all non-synonymous single nucleotide polymorphisms (SNPs), insertions and deletions in SARS-CoV-2 genomes with ≥5% allelic frequency and population frequency of ≥5% and ≤95%. We also examined correlations between the presence of SNPs with each other, with subvariants, and over time.
Results Seventy-nine mutations met inclusion criteria. Among 15,566 persons infected with Omicron SARS-CoV-2, 1,825 (12%) were unvaccinated with no prior recorded infection, 360 (2%) were unvaccinated with a recorded prior infection, 13,381 (86%) had a complete primary series vaccination, and 9,172 (58%) had at least one booster. After examining correlation between SNPs, 79 individual non-lineage defining mutations were organized into 38 groups. After correction for multiple testing, no individual SNPs or SNP groups were significantly associated with immunity status levels.
Conclusions Genomic variation identified within SARS-CoV-2 Omicron specimens was not significantly associated with immunity status, suggesting that contribution of non-lineage defining SNPs to immune evasion is minimal. Larger-scale surveillance of SARS-CoV-2 genomes linked with clinical data can help provide information to inform future vaccine development.
Background
Since its emergence, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has been evolving, with evidence for increased transmissibility and ability to evade existing immune responses elicited by prior infection or COVID-19 vaccination.1,2 Multiple mutations that increase transmission and immune evasion were identified when the Delta and subsequent Omicron variants (and sub-variants) emerged. However, few in-depth analyses characterizing the genomic variation within a circulating variant and identifying potential mutations associated with immune evasion have been conducted. In addition, during the Omicron wave, individual immunity could be due to prior infection or multiple levels of vaccination, due to the availability of vaccine boosters in this timeframe, and any immune evasion characteristics found among mutations could be tied to a specific type of immunity, derived from, for example, prior infection, immunization dose(s), or both.
To address these data gaps, we conducted a retrospective cohort study of persons with nucleic-acid amplification test (NAAT)-confirmed SARS-CoV-2 infection with the Omicron variant confirmed by whole genome sequencing. Patients were members of a large, integrated health system that had detailed epidemiologic and clinical data on all participants. The study assessed overall SARS-CoV-2 genetic variation over time as well as the relationship between mutations and the individual’s source(s) of potential COVID-19 immunity. The key assumption in this study is that after adjusting for sociodemographic characteristics and timeframe, an association between mutation presence and vaccine- or infection-derived immunity could imply greater immune evasion for that mutation.
Methods
Setting
Kaiser Permanente Northern California (KPNC) is an integrated health system that serves over 4.5 million members in Northern and Central California and provides comprehensive preventive and curative care in inpatient and outpatient settings across 266 medical offices and 21 hospitals. Members receive most clinical services, including laboratory testing, outpatient, and inpatient care in KPNC facilities. Members have similar sociodemographic characteristics to the population of Northern and Central California.3 SARS-CoV-2 NAAT for outpatients and hospitalized patients is conducted at a regional laboratory or local hospital laboratories using various platforms and specimens are sent to the California Department of Public Health (CDPH) for whole genome sequencing.
Study design
We conducted a retrospective cohort study of persons with NAAT-confirmed SARS-CoV-2 infection and whole genome sequencing-confirmed Omicron variant infection and assessed patterns among non-lineage-defining mutations as well as the association between SARS-CoV-2 genetic variation and immunity source (if any). This immunity was defined by prior reported SARS-CoV-2 infection or vaccination. Epidemiologic and clinical data were obtained from the KPNC Virtual Data Warehouse, a common data model into which standardized data are extracted from clinical and administrative databases including an integrated electronic health record database (Epic, Verona, WI, USA), and from other sources.4,5 The study was approved by the KPNC and University of California, Berkeley Institutional Review Boards with waivers of the requirement for informed consent. This activity was reviewed by the Centers for Disease Control and Prevention (CDC) and conducted consistent with applicable federal law and CDC policy (see e.g., 45 C.F.R. part 46.102(l)(2), 21 C.F.R. part 56; 42 U.S.C. §241(d); 5 U.S.C. §552a; 44 U.S.C. §3501 et seq.).
Inclusion criteria
We included all persons with an incident NAAT-confirmed SARS-CoV-2 infection and whole genome sequencing-confirmed Omicron variant infection from January 1, 2022 to October 31, 2022 who met the following criteria: 1) KPNC membership for at least one year prior to SARS-CoV-2 diagnosis; 2) aged ≥5 years and <90 years; 3) no prior receipt of COVID-19 vaccine or the completion of a primary COVID-19 vaccine series of two doses of BNT162b2 [Pfizer/BioNTech] or m-RNA-1973 [Moderna/National Institutes of Health] or one dose of Ad.26.COV2.s [Janssen] more than 14 days prior to index date; 4) incident infection defined as the first NAAT-confirmed SARS-CoV-2 infection during the study period. The index date was defined as the date of NAAT-confirmation of SARS-CoV-2 infection.
SARS-CoV-2 whole genome sequencing and processing
KPNC is a partner in COVIDNet with CDPH and submits a majority of NAAT-confirmed SARS-CoV-2 specimens to COVIDNet for genomic surveillance; a subset of specimens undergo whole genome sequencing at CDPH or affiliated laboratories based on CDPH priorities, including prioritizing hospitalized cases and a convenience sample of the first specimens received per week (number varies by burden).6 COVIDNet consists of a set of (primarily University of California-affiliated) laboratories tapped to assist CDPH because they had surplus sequencing capacity. Illumina sequencing was used by most of the sites, with about 10% of genomes being sequenced by Oxford Nanopore. For this study, we only included Illumina-sequenced genomes because the error rate of Nanopore is too high to accurately identify low-frequency mutations. Varying sets of primers were used by different institutions, but the most common for the Omicron genomes were the ARTICv3 and ARTICv4 primers. The University of California, San Francisco lab used their own custom set of primers with more primers and shorter amplicons. After sequencing, the CDPH Viral and Rickettsial Disease Laboratory links lineage information with patient identifiers and shares the results with KPNC. The raw sequence files are shared with the University of California, Berkeley (UCB); KPNC and UCB collaborate further to link demographic, clinical, and epidemiologic data from KPNC with subsequent sequence analysis data from UCB.
The raw sequencing data were processed through a SARS-CoV-2 analysis pipeline that was modified for this work as follows. Adapter removal and trimming were performed using bbduk (http://sourceforge.net/projects/bbmap/). The reads were then aligned to the Wuhan reference genome (GISAID ID: EPI_ISL_402125, GenBank accession number MN908947) using minimap27 followed by primer trimming using iVar.8 Samtools was used to create a pileup file from which a consensus sequence was generated. This consensus genome was created with iVar using a minimum depth of 10 reads and majority rule for base calling. We next used iVar to call variants from the pileup file where we set the sensitivity threshold for calling a mutation to 2%. This process calls mutations for any loci where at least 2% of the reads are non-reference. Using this very low threshold allowed us to capture even very low frequency mutations in the sequencing data. This threshold resulted in many spurious mutations from sequencing and alignment errors and therefore we later set a threshold of 5% when processing mutations. This original list of variants was then annotated with the gene and amino acid change (if there was one) or insertion or deletion.
After calling with iVar using a low-frequency threshold, mutations were compiled into a master frequency matrix listing all mutations that appeared in any genome in rows, and all the genomes that each mutation appeared in columns (https://github.com/staciawyman/SC2_GenVar). The cell values were frequency of the alternate reads at the location. A mutation had to have an allele frequency of at least 5% to be included to rule out sequencing errors. This matrix was partitioned into variant-specific parts to facilitate analysis. We constructed an Omicron matrix including all individuals that had the Omicron lineage9 between January 1, 2022 and October 31, 2022 including non-synonymous SNPs, insertions and deletions (indels). Pipeline code is available at https://github.com/staciawyman/SC2_GenVar.
Outcome measure
The primary outcome measure was, among persons with an incident NAAT-confirmed SARS-CoV-2 infection, the presence of individual mutations (nonsynonymous SNPs, insertions and deletions) with ≥5% allelic frequency and ≥5% and ≤95% population frequency, as well as mutation groups created by grouping SNPs with a pairwise Pearson correlation coefficient > 0.7. Additionally, patterns in the prevalence of these mutations were assessed as follows: as noted earlier, pairwise Pearson correlations with other mutations in each infection; presence within Omicron subvariant (BA.1, BA.2, BA.4, BA.5, BCD, BE, BF, Bdot), and prevalence over time.
Primary exposure: SARS-CoV-2 immunity status
The primary exposure of interest in this analysis was immunity level, defined in mutually exclusive categories: i) no recorded vaccination or recorded prior infection; ii) recorded prior infection but no recorded vaccination; and iii) recorded vaccination (primary series complete and/or more than primary series), regardless of prior infection status.
Statistical analysis
Adjusted relative prevalence, shown as an adjusted rate ratio (aRR), and 95% confidence intervals (CIs) for the estimates for our exposure and adjustment variables were obtained from modified Poisson regression models with robust standard errors using PROC GENMOD with a log link function in SAS (Version 9.4, Cary, North Carolina); statistical significance was defined as a p-value <0.05. All p-values were adjusted for multiple comparisons to control the false discovery rate at 5%, using the Benjamini-Hochberg procedure.10 To adjust for confounding, all models included age category, sex, race/ethnicity, Charlson comorbidity index score11, and index month. Adjusting for these demographic and clinical factors accounted for differential vaccination status of these subgroups, increasing vaccination rates over time, and the inherent uncertainty as to where in a population a new mutation might arise.
Role of funding source
The funders had no role in study design, data collection, data analysis, interpretation, or writing of the report.
Results
From January 1, 2022 to October 31, 2022, 312,249 SARS-CoV-2-positive specimens were submitted for sequencing, 18,603 (6.0%) underwent sequencing (Supplemental Figure 1), and 18,476 (5.9%) had reported lineages consistent with the Omicron variant. After applying selection criteria, 15,556 persons with SARS-CoV-2 Omicron variant infection were included in the analysis: 9,172 (58%) vaccinated with a booster; 4,209 (26%) vaccinated with a complete primary series [two doses of Pfizer or Moderna vaccine]; and 2,185 (14%) with no vaccination (Supplemental Figure 1). The prevalence of Omicron subvariants in our specimens, aggregated to the week level, showed increasing diversity over time (Figure 1).
Cohort flowchart
SARS-CoV-2-positive specimens tested by Kaiser Permanente Northern California, proportion by subvariant, by week from January 1, 2022 to October 31, 2022
After applying selection criteria to the sequences, we identified 78 non-lineage defining mutations (all non-synonymous SNPs) in 1,456 SARS-CoV-2 SNPs included in this analysis, shown in Supplemental Figure 2. Identified mutations were in the spike gene (n=44), nucleocapsid gene (n=8) open reading frame 1a (ORF1a) (n=19), ORF1b (n=4), ORF3a (n=1), ORF6c (n=1), ORF10 (n=1). The population frequency within our genomes of individual identified mutations ranged from 5.0% to 95.0% and 46 out of 79 (58.2%) mutations were present in <50% of the overall population. Although we called low-frequency mutations with our pipeline, none of them met the criteria for inclusion in our analysis.
Pearson correlations of SNP presence in specimens, before and after clustering.
Of the 78 identified SNPs that met our inclusion criteria, many co-occurred on the same specimens. Supplemental Figure 2 shows the pairwise correlations of SNP presence in our dataset. Given the strong pairwise correlations observed, we formed groups of SNPs that had a pairwise Pearson correlation > .7, reducing the dimension to 33 SNPs/SNP groups, also shown in Supplemental Figure 2. These mutation groups are largely monophyletic and segregate together across the phylogeny of Omicron genomes (Supplemental Figure 3).
Phylogenies of 15,045 Omicron genomes colored by group. Mutations in group 52 (top) represent a clade defined by the BA.2.12.1 subvariant (in blue). Mutations in group 4 (bottom) fall into a subclade of BA.4 (in yellow).
The demographics of persons with SARS-CoV-2 infection included in this analysis, stratified by vaccination and prior infection status, are shown in Table 1, and the number of total reported infections and those sequenced over time is shown in Supplemental Figure 4. Significant differences were found between groups, thus supporting their inclusion in our covariate adjustment set. Figure 2 details some of the relationships between SNP groups, time, and subvariant among the specimens in our data set.
Persons with SARS-CoV-2 infection tested by Kaiser Permanente Northern California stratified by subvariant by week from January 1, 2022 to October 31, 2022
SNP/group prevalence by subvariant and week in 2022
While certain mutations/groups appear essentially only in a single subvariant (for example, 99.1% of S:H69L mutations appear in specimens classified as BA.1 in Figure 2c), the majority of mutations/groups appear in many subvariants, suggesting that they are not simply proxy measures of subvariant. Further, the presence of mutations/groups do not follow clear temporal patterns (Figure 2b). While some mutations such as Group 12 arise, become plentiful, then wane, others disappear and re-appear over time without a clear pattern (for example, S:D1259H). Smoothed over time, grouped SNP prevalence by month is shown in Supplemental Figure 5.
Proportion of specimens with a given SNP or group, by month of 2022
Using modified Poisson regression models with robust standard errors, we found several of our 33 SNP groupings (Table 2) and/or 78 individual mutations (Supplemental Table 1) that were significantly associated with immunity status; however, after correcting p-values to adjust for multiple comparisons, none of the associations were statistically significant. These findings were not meaningfully changed when people with a booster vaccine were/were not aggregated with the complete primary series group (Supplemental Table 2, Supplemental Table 3). The SNPs/groups with the largest absolute associations with vaccination (Table 2) included some that were strongly associated with a single variant (for example, ORF10:L37F, which occurred almost exclusively in BA.5 subvariants) as well as some that frequently appeared in multiple variants (for example, ORF1b:N498I and group 45).
When looking at the evolutionary relationships with the genomes colored by vaccination status (Supplemental Figure 6), vaccination status did not show a correlation with a clade/subvariant or distance from the root (number of mutations with respect to the reference SARS-CoV-2 genome).
Phylogeny of 15,045 Omicron genomes sampled between 1/1/2022 and 9/25/22. Branch lengths represent divergence from Wuhan reference genome and nodes are colored by vaccination status. Vaccinated category includes primary series and primary plus booster.
Discussion
Characterizing the continuing evolution of SARS-CoV-2, along with active surveillance of infections, is critical to targeting public health measures and vaccine development. Large-scale global whole genome sequencing efforts have contributed substantially to our understanding of the genetic variation of SARS-CoV-2 and its associations with transmissibility and immune evasion. Our current study adds to this evidence base by assessing within-variant (Omicron) variation in the genome beyond the already-identified subvariants such as BA.1 and BA.2 using a unique whole genome surveillance system linking genomic data with detailed data on demographic, clinical and epidemiologic characteristics in a well-defined population.
Understanding the potential for immune evasion based on mutations beyond the spike protein is critical for vaccine development and understanding the mechanisms of the human immune response. In this detailed analysis of genomic variation within the Omicron variant over a 10-month period in California, we identified only 78 mutations that met our inclusion criteria among 15,566 genomes. These mutations occurred across the entire genome, not just in the spike protein. Spike protein mutations previously implicated with immune evasion were excluded based on population frequency. We encourage future research to examine the entire viral genome for potential immune evasion or increased replication or infectivity.
Our results showed that there appear to be some mutations that are tightly grouped with certain subvariants, while a number of other mutations appear across multiple subvariants. While we did not find evidence for immune evasion among our identified mutations/groups after correction for multiple testing, it is possible that we are underpowered to detect small effects over a short time period.12 Many of the SNPs with the highest adjusted risk ratios appeared in a small number of specimens (Supplemental Table 1). While one would expect more variability from the smallest table cells if we were simply observing statistical noise, it is also possible that a larger sample could confirm these stronger adjusted risk ratios.
Our analysis has several limitations. First, not all possible demographic and clinical variables were included in our models, and if those variables are associated with immunity status and SNP presence, unmeasured confounding could bias our analysis. Second, there is a risk of misclassification of immunity status; though the KPNC surveillance system for vaccination is quite robust, prior infections (say, diagnosed by a home test) may have gone unreported. Third, we focused on a set of specimens from Northern and Central California, limiting the generalizability of our findings to this context. Finally, there may be hidden bias as to which specimens were successfully sequenced based on the nucleic acid viral load in the specimen (e.g., specimens with lower viral loads might not be successfully sequenced and thus would have been excluded from this analysis). Fourth, although our surveillance system has a very robust system for assessing individual-level vaccination status, there is a small risk of misclassification of persons as having incorrect immunity status.
As part of the global effort to conduct near real-time genomic surveillance for SARS-CoV-2, there has been tremendous progress in developing national genomic surveillance systems with nearly real-time public reporting of circulating SARS-CoV-2 variants and uploading of consensus sequence data to publicly available genomic databases.13 14 15 16 However, these efforts can be improved. First, enhancing methods to define consensus sequences to also assess for minor variation with <50% allelic frequency would better define the true genomic variation of circulating SARS-CoV-2, including capturing within-host genetic diversity.17,18 Second, current SARS-CoV-2 genomic surveillance systems in the United States have limited metadata on human demographics, clinical factors, epidemiologic factors, and immunity status (e.g. COVID-19 vaccination and documented prior infection), particularly when compared to the United Kingdom.19 This limits the ability to correlate noted patterns in SARS-CoV-2 genetic variation and real-world clinical and epidemiologic outcomes in the United States. We present one potential analysis that uses this data linkage to better understand how SARS-CoV-2 genetic variation might be associated with human disease, in this case evidence for immune evasion associated with specific mutation clusters. Improving the clinical metadata associated with SARS-CoV-2 genomic data should be a priority as we continue to improve and expand genomic surveillance. We have demonstrated this capability through partnership and collaboration between a large healthcare system and a state public health department, namely KPNC and CDPH/COVIDNet.
In summary, we demonstrate that linking SARS-CoV-2 genomic sequence data with detailed clinical data revealed no significant associations between genetic variation among Omicron SARS-CoV-2 mutations and immunity status. Further research using larger sample sets that link SARS-CoV-2 genetic data with clinical data is warranted for the potential to reveal whether SARS-CoV-2 non-lineage defining genetic variation might contribute to the fitness of SARS-CoV-2 by either immune evasion or other mechanism.
Funding Statement
This study was supported by the Rockefeller Foundation (SKW, JS), National Cancer Institute of the National Institutes of Health under Award Number U01CA260584 as part of the Serological Sciences Network (SeroNet) (JS), the Physician Researcher Program of The Permanente Medical Group Delivery Science and Applied Research Program (JS), The Packard Foundation (SKW), and The Sergey Brin Family Foundation (SKW). CDPH/COVIDNet genomic surveillance work was funded by Centers for Disease Control and Prevention, Epidemiology and Laboratory Capacity for Infectious Diseases, Cooperative Agreement Number 5 NU50CK000539.
Declaration of interests
All authors: No conflicts of interest identified. All authors have submitted the ICMJE Form for Disclosure of Potential Conflicts of Interest. Conflicts that the editors consider relevant to the content of the manuscript have been disclosed.
Data Availability
De-identified data that underlie the results reported in this article will be made available upon reasonable request.
Data sharing statement
De-identified data that underlie the results reported in this article will be made available upon reasonable request.
Disclaimer
The content is solely the responsibility of the authors and does not represent the official views or opinions of the National Institutes of Health, Kaiser Permanente, California Department of Public Health or the California Health and Human Services Agency.
*Additional Disclaimer
Use of trade names and commercial sources is for identification only and does not imply endorsement by the California Department of Public Health or the California Health and Human Services Agency.
Acknowledgments
The authors are grateful to all Kaiser Permanente members without whom this study would not have been possible. In addition, we wish to acknowledge our partnership with the California Department of Public Health and participation in the California SARS-CoV-2 Whole Genome Sequencing Initiative (COVIDNet). We thank the VRDL and COVIDNet Teams: Summer Adams, Allison Bailey, Matt Bacinskas, Nikki Baumrind, Elizabeth Baylis, John Bell, Ricardo Berumen III, Yocelyn Cruz, Mojgan Deldari, Alex Espinosa, Sabrina Gilliam, Madeline Glenn, Bianca Gonzaga, Melanie Greengard, Jill Hacker, Kim Hansard, Monica Haw, Thalia Huynh, Chantha Kath, Ruth Lopez, Sharon Messenger, Alexa Quintana, Chris Preas, Clarence Reyes, Maria Salas, Hilary Tamnanchit, Serena Ting, Kathy Jacobson, MD, and Carol Glaser, MD.
Footnotes
Declaration of interests: All authors: No conflicts of interest identified. All authors have submitted the ICMJE Form for Disclosure of Potential Conflicts of Interest. Conflicts that the editors consider relevant to the content of the manuscript have been disclosed.
References
