Heterologous ChAdOx1-BNT162b2 vaccination in Korean cohort induces robust immune and antibody responses that includes Omicron.

Heterologous ChAdOx1-BNT162b2 vaccination induces a stronger immune response than BNT162b2-BNT162b2. Here, we investigated the molecular transcriptome, germline allelic variants of immunoglobulin loci, and anti-Omicron antibody levels in 46 office and lab workers from the Republic of Korea following ChAdOx1-BNT162b2 vaccination. Anti-spike-specific IgG antibody levels against the ancestral SARS-CoV-2 strain increased from 70 AU/ml to 14,000 AU/ml to 142,000 AU/ml one, three and seven days following the second vaccination. Titers against VOC, including Omicron, were two-fold to three-fold lower, yet higher than those measured following BNT162b2-BNT162b2 vaccination. RNA-seq of peripheral immune cells demonstrated activation of interferon pathways with increased IGHV clonal transcripts encoding neutralizing antibodies. scRNA-seq revealed enriched B cell and CD4+ T cell responses in both ChAdOx1-BNT162b2 and BNT162b2-BNT162b2 recipients, but a stronger clonal expansion of memory B cells with ChAdOx1-BNT162b2. In summary, heterologous ChAdOx1-BNT162b2 provides an innate and adaptive immune response that exceeds homologous BNT162b2 vaccination.

Introduction

The ChAdOx1 nCoV-19 vector (AZD1222) and BNT162b2 mRNA vaccines (hereafter referred to as ChAd and BNT, respectively) have been widely used and shown to induce robust immune responses against the spike protein of SARS-CoV-2. Both vaccines have shown remarkable efficacy in preventing COVID-19 disease. Based on the approval by government agencies and available supply, ChAd was the main vaccine used in the Republic of Korea (South Korea) in the spring of 2021. The manifestation of rare events, including thrombosis and thrombocytopenia syndrome, associated with adenovirus-based vaccines prompted pausing of the distribution of the ChAd vaccine in many countries, including South Korea (Hwang et al., 2021). This promoted partially vaccinated people to complete their vaccination with an mRNA vaccine, Pfizer/BioNTech (BNT) or Moderna.

The ChAd-BNT heterologous vaccination strategy resulted in significantly greater immunoglobulin G (IgG) immune responses aimed against the SARS-CoV-2 spike protein compared with the ChAd-ChAd strategy (Agrati et al., 2021; Barros-Martins et al., 2021; Deming and Lyke, 2021; Groβ et al., 2021; Pozzetto et al., 2021; Schmidt et al., 2021b). Although the immune response to different vaccine regimens has been investigated, there is limited knowledge about the molecular transcriptome responses, elicited by heterologous vaccinations. Here, we hypothesized that the heterologous vaccination elicits unique molecular responses in immune cells by activating transcription of specific germline allelic variants. To test this hypothesis, we investigated immunogenicity and reactogenicity in a Korean cohort receiving the heterologous ChAd-BNT vaccine regimen in comparison with a homologous BNT-BNT cohort (Lee et al., 2022b).

Although many vaccine studies are based on individuals with European and African ancestry (Arunachalam et al., 2021; Sokal et al., 2021), an inclusion of samples from other genetic backgrounds would enable the equitable advancement of genomic vaccinology. A dearth of genomic vaccinology studies on Asian populations prompted this study. Addressing the distinct usage of germline allelic variants of immunoglobulin loci upon COVID-19 vaccination of different populations might enhance our understanding of antibody production. Although there have been COVID-19 vaccine studies of Asian/Korean populations (June Choe et al., 2021; Kim et al., 2021b; Lim et al., 2021a, 2021b), none of them had investigated the molecular immune response.

Here, we use serology and RNA-seq transcriptome analyses to investigate the immune response of office and laboratory workers from a major health care center in the Republic of Korea. Study participants received one ChAd dose followed either by the BNT vaccine. The sequencing depth in this study, and our previous study on individuals receiving two doses of BNT (Lee et al., 2022b) facilitated the identification of an expanded immune response in individuals receiving the heterologous vaccination.

Results

Study design

As current discussion centers around the immune response elicited by heterologous ChAdOx1-BNT162b2 (referred to ChAd and BNT throughout the manuscript) vaccinations as compared to homologous two ChAd or two BNT doses, we investigated the molecular immune responses elicited by a ChAd-BNT heterologous vaccination and identified differences and similarities to a homologous BNT-BNT vaccination from our previous study (Lee et al., 2022b) and a smaller ChAd-ChAd cohort. The study had originally been designed to investigate the immune response after two doses of ChAd and employees (n = 46 (41 Females, 5 Males), median 35 years old, median BMI 22.4) from the Asan Medical Center, Seoul, Republic of Korea received the ChAd vaccine in June 2021 (Figure 1A; Table S1). Owing to a revised vaccination policy the participants received the BNT vaccine as the second dose in August 2021. Four participants had underlying health conditions and two reported mild postvaccine symptoms after the initial ChAd dose. We measured anti-spike IgG antibody levels, circulating cytokines levels and immune transcriptomes (RNA-seq) on peripheral immune cells before the vaccination (referred to as D0) and on days 2-4 (mean 2.8, median 3.0 days) and 7–10 (mean 7.9, median 7.0 days) postvaccination (referred to D3 and D7 throughout the manuscript). For selected samples, we also conducted scRNA-seq. Four additional study participants received two doses of ChAd and although these results are not presented in this study, their RNA-seq data have been uploaded to GEO.

Figure 1

Experimental design and antibody

(A) Schematic presentation of the experimental workflow. The 46 study subjects received ChAd as a first vaccine dose and BNT as a second dose. Blood was collected on the day of vaccination (I-D0 and II-D0) and at days 1-3 (median, day3; D3) (I-D3) and day 7-10 (median, day7; D7) (I-D7) postvaccination as indicated by the colored circles.

(B) Plasma IgG antibody binding the SARS-CoV-2 RBD (spike) from different strains and the SARS-CoV-2 N protein in the ChAd-BNT vaccination groups. p values are from 2-way ANOVA with Tukey’s multiple comparisons test.

(C) Plasma IgG antibody binding the Omicron spike protein at day 7 after the second dose (II-D7).

(D) Neutralizing antibody response to virus spike protein of SARS-CoV-2 original and variants. p values are from two-way ANOVA with Tukey’s multiple comparisons test.

(E) ACE2 inhibition to Omicron spike protein at day 7 after the second dose (II-D7).

(F and G) Relative levels of anti-spike IgG antibody and ACE2 inhibition to the level from the ancestral strain in the serum at day 7 after the second dose (II-D7). p values are from two-way ANOVA with Sidak’s multiple comparisons test. Median: line inside violin plot.

Antibody response after the first and second vaccination

First, we measured circulating antibody responses in plasma samples from the individuals receiving the heterologous ChAd-BNT vaccine regimen using enzyme-linked immunosorbent assay (ELISA) (Figure 1B). On average, anti-spike (WH04) IgG levels were approximately 70 AU/ml within three days (I-D3) following the initial ChAd dose and 14,000 AU/ml within three days (II-D3) following the BNT vaccine dose (Figure 1B). Another 10-fold increase was observed at day seven (II-D7) following the BNT vaccination. A similar pattern was observed for the spike proteins from the Alpha, Beta, Gamma, Iota, Kappa, Zeta, and Delta variants. At day 7 (D7) post the BNT vaccination, antibody concentrations against Omicron (B.1.1.351) were reduced by approximately 65% in the ChAd-BNT group compared to the level against the ancestral strain (Figure 1C).

At this point in the pandemic, a critical question is whether antibodies induced by vaccination can neutralize current variants. Here, we used the degree of serum antibodies blocking the binding of ACE2 to SARS-CoV-2 variants’ spike as a proxy to in vitro neutralization activity. Inhibition of ACE2 binding for the ancestral strain increased approximately 13-fold between days three after the first (I-D3) and second (II-D3) vaccination and another 8-fold at day 8 (II-D8) (Figure 1D). The increase for more recent variants ranged from 3-fold for Alpha to 16-fold for Delta following the second vaccination. The majority of the 46 study subjects did not display any ACE2 inhibition against Omicron after the second dose (Figure 1E), suggesting little neutralizing activity. Only 14 individuals displayed measurable ACE2 inhibition against Omicron. Increased neutralizing activity against different variants, including Omicron, was also detected in a reference cohort (n = 14, median 58 years old) receiving a homologous BNT-BNT vaccine protocol (Lee et al., 2022b). Next, we directly compared anti-SARS-CoV-2 antibody titers between the heterologous ChAd-BNT cohort and the homologous BNT-BNT reference cohort within seven days after the second dose (II-D7) (Figure 1F). The anti-spike IgG levels, including Omicron, were significantly higher in the heterologous cohort. However, although median levels of antibody to variants were higher in the heterologous cohort, no significant differences in ACE2 binding inhibition were observed (Figure 1G).

Immediate immune response

Early responses to vaccination materialize in elevated levels of interferons and other cytokines. To gauge the early vaccine response, we measured plasma levels of a panel of cytokines in all individuals before and after the initial ChAd vaccination and the BNT dose (Figure 2). Out of the 23 cytokines measured, a transient increase of circulating IL-10 and CXCL10 was observed within three days after the first vaccination (ChAd) and returned to baseline by day 7 (I-D7) (Figures 2A and 2B). Upon delivery of the second (BNT) vaccine dose, CXCL10, IL-10, and IL-1Rα levels were statistically elevated (Figures 2C and 2D). CXCL10 (IP-10) expression is rapidly activated following vaccination and viral infections and it has been described as a biomarker associated with COVID-19 severity (Huang et al., 2005; Laing et al., 2020; Sobolev et al., 2016). Its regulation by IFN-γ is mediated by the JAK-STAT signaling pathway (Lee et al., 2021a).

Figure 2

Serum cytokine and chemokine levels

(A) Volcano plots depict differentially expressed cytokines and chemokines at days 0, 3, and 7 after the first dose vaccination. Black dots indicate significantly upregulated and downregulated cytokines (p value < 0.05). Gray dots indicate cytokines not regulated by the vaccination.

(B) CXCL10, IL1-Rα, and IL-10 serum concentrations before and after vaccination (pg/mL).

(C) Volcano plots depict differentially expressed cytokines and chemokines at days 0, three, and 7 after second dose of vaccination.

(D) Concentrations of CXCL10, IL1-Rα, and IL-10 in serum after second dose vaccination. Boxes show median, 25th and 75th percentiles and whiskers show the range. p value are from two-way ANOVA with Tukey’s multiple comparisons test and listed in ∗ p <0.05, ∗∗p < 0.01, ∗∗∗∗p < 0.0001. Median, middle bar inside the box; IQR, 50% of the data; whiskers, 1.5 times the IQR.

Transcriptome response

To further understand the molecular responses elicited by the heterologous vaccination, we analyzed the transcriptomes induced after the first (ChAd) and second (BNT) vaccine dose (Figures 3, S1, and S2; Table S2). Functional enrichment analyses of common and unique DE-Gs to each vaccination regimen revealed differences between the ChAd-BNT and BNT-BNT regimens (Figure S2). Expression of 308 genes was induced within two days (I-D2) after the prime ChAd vaccination (Figure S1A). GSEA analyses linked them to defense response pathways innate immune responses, including interferon and cytokine signaling (Figure S1B). Specifically, genes regulated by the JAK-STAT pathway, including STAT1 itself, interferon-induced (IFI) genes, and the antiviral OAS family were activated within two days after the vaccination (Figure S1C).

Figure 3

Transcriptional signatures elicited by the second vaccination in the ChAd-BNT and BNT-BNT cohorts

(A) Volcano plot of DE-Gs comparing day 0 (II-D0) versus day 2 (II-D2) after the second vaccination.

(B) Genes expressed at significantly higher levels at II-D2 were significantly enriched in Gene Ontology (GO). X axis denotes statistical significance as measured by minus logarithm of FDR q-values. Y axis ranked the terms.

(C) Venn diagram displays the number of significantly induced genes between II-D0 and II-D2 in the ChAd-BNT and BNT-BNT vaccination cohorts.

(D and E) Fold activity of IFN response genes that are significantly induced in both cohorts as compared to day 0. The names of key IFN response genes in COVID-19 patients or vaccination are shown in red.

Next, we analyzed the activation of genes induced in the heterologous group within two days after the second dose (Figure 3). Expression of 318 genes was significantly induced at day 2 as compared to day 0 (Figure 3A) and levels of 84 genes remained elevated at day 7 (Table S3). The 318 BNT-induced genes are enriched in IFN response, complement and cytokine signaling (Figures 3B and S2). Genes with elevated expression at day 7 following the BNT vaccination are associated with proliferative responses. To further investigate the impact of the heterologous vaccination regimen, we compared its immune response with that of the homologous reference cohort (Lee et al., 2022b) after the second BNT dose (Figure 3C). Although expression of 235 genes was significantly higher in the heterologous group, 90 were selectively elevated in the homologous reference group and 83 genes were significantly induced in both cohorts. Thirty-one genes that were preferentially activated in the heterologous group are part of interferon-activated pathways, including CXC chemokines and interferon induced genes (Figure 3D). Fold change differences ranged from two to 90 (log21 ∼ log26.5). Expression of 21 IFN-response genes was higher in the homologous group compared to the heterologous group (Figure 3E), with fold change differences ranging from two to eleven (log21 ∼ log23.5). Fifty-two interferon-regulated genes are preferentially expressed after the second dose and activation of most of them has been observed in peripheral immune cells of COVID-19 patients (Knabl et al., 2022; Lee et al., 2022a). These include classical interferon-response genes as well as antiviral classes, such as the OAS family.

Germline allelic variants of immunoglobulin loci

The preferential increase in anti-spike antibodies and neutralizing antibodies in the heterologous cohort upon receiving the second (BNT) vaccine, led us to dig deeper and interrogate the expression profiles of specific germline variable gene classes. An average sequencing depth of 230 million reads per sample permitted a detailed analysis of germline gene usage. We determined the range of IGHV, IGKV, and IGLV gene usage in the ChAd-BNT heterologous cohort in comparison to the homologous BNT-BNT reference cohort (Lee et al., 2022b) after the second vaccination (Figures 4 and S3). RNA-seq was conducted before the first vaccination (I-D0) and seven days after the second dose (II-D7). Analysis of immunoglobulin heavy chain variable (IGHV) and light chain variable (IGKV and IGLV) genes with first complementarity determining region (CDR1) and CDR2 revealed the use of a broad range of germlines in both cohorts. Although the number of final clonotypes in the heterologous cohort increased approximately 1.3-fold, an almost two-fold increase was observed in the BNT-BNT homologous cohort (Figure 4A). We observed increased frequency of transcription of several VH genes, including IGHV3-30 and IGHV1-61 (Figures 4B and S3A). These particular heavy chain variants have been identified in SARS-CoV-2 infected patients that developed neutralizing antibodies (Andreano et al., 2021; Andreano and Rappuoli, 2021; Zhang et al., 2021; Zost et al., 2020). Similarly, specific IGKV and IGLV allelic variants were induced in the heterologous vaccine group and some variants were induced in both groups (Figures 4C, 4D, S3B, and S3C). Although IGHV3-30 and IGHV4-61 clones were induced in both cohorts, their levels were higher in the BNT-BNT group. IGHV2-70 that had been detected in Korean COVID-19 patients as a potently therapeutic neutralizing antibody (Kim et al., 2021a) was not induced in the two cohorts. Overall, more diverse immunoglobulin genes were induced in the BNT-BNT group compared to the ChAd-BNT group. A similar pattern was observed for IGKV and IGLV genes (Figures 4C and 4D). In contrast to BCR genes, the activation of TCR genes (TRAV, TRBV) in the ChAd-BNT group exceeded that seen in the BNT-BNT group (Figures 4E and 4F). T cell responses are greater in the heterologous group (ChAd-BNT) compared with the homologous group (BNT-BNT, ChAd-ChAd) (Liu et al., 2021b).

Figure 4

SARS-CoV-2-RBD-specific B cell memory between ChAd-BNT and BNT-BNT cohorts

(A) Fold-induction of the number of final clonotypes that are generated between the first dose (I-D0) and the second dose (II-D7) in both cohorts.

(B–D) Boxplots show the induced frequency (%) of variable (V) gene usage for all paired heavy (B), kappa (C), and lambda (D) chains in ChAd-BNT and BNT-BNT groups between day 0 (II-D0) and day 7 (II-D7) after the second vaccination. Those identified in both cohorts are marked in red.

(E and F) V region distribution in samples. The usage rate of V regions from TCRα chain (E) and TCRβ chain (F). Floating bars (min to max); line at mean.

Identifying B cells and T cells through scRNA-seq

Single cell RNA-seq (scRNA-seq) and CyTOF data of PBMCs or whole blood from COVID-19 infection and SARS-CoV-2 mRNA vaccinations reveal an increase of monocytes, B cells and CD4+ T cells (Arunachalam et al., 2021; Liu et al., 2021a; Sureshchandra et al., 2021; Tian et al., 2022). To specifically assess distinct memory responses induced by the heterologous and homologous vaccinations, we performed scRNA-seq on PBMCs at day 2 (II-D2) and day 7 (II-D7) after the second dose and identified cell populations in the immune compartment (Figures 5A, 5B, and S4). According to the expression of canonical cell-type markers, six major cell types were identified: B cells (CD74, CD79A, and MS4A1), CD4 T cells (IL7R, CD3D, and CDL11B), CD8 T cells (CCL5, CD8A, CD8B, CD3D, and CD3G), dendritic cells (DC) (CST3, LYZ, and FCER1A), monocytes (Mono) (CST3, LYZ, and PSAP), and natural killer (NK) cells (GNLY, NKG7, and GZMB). B cell populations, including naive and memory B cells, were elevated in the ChAd-BNT group, whereas T cell populations (CD4 and CD8 T cells) were different depending on participants in both groups between II-D1 and II-D7 after the second dose (two-way ANOVA followed by Tukey’s multiple comparisons test, ∗∗p < 0.01, ∗∗∗p < 0.001) (Figure 5C).

Figure 5

Dissection of immune dynamics using scRNA-seq

(A) The UMAP projection of PBMCs at day 7 (II-D7) after the second dose in the ChAd-BNT group was colored by each cell population.

(B) Violin plots of cell identity markers. Expression levels are shown on the y axis.

(C) Frequencies of total B, CD4 T, CD8 T, and their naive and memory cells at day 2 and 7 (marked as D2 and D7) following second dose vaccination.

(D) Enrichment of DE-Gs comparing heterologous ChAd-BNT vaccination and homologous BNT-BNT vaccination at day2 after the second dose. X axis denotes statistical significance as measured by minus logarithm of FDR q-values. Y axis ranked the Hallmarks.

To understand the difference between cell response to different vaccine regimen, we performed DEG analysis on the B, CD4 T, and CD8 T cells between two groups at day 2 after the second dose (Figure 5D: Table S4). As expected, IFNs responses showed more activation in the B cells of the ChAd-BNT group compared to the BNT-BNT group. However, in the T cells, IFNγ responsive genes were detected higher in the BNT-BNT group compared to the ChAd-BNT group. Interestingly, cell cycle-related modules are identified in the ChAd-BNT group, whereas TNFα and mTOCR1 signaling is stronger in the BNT-BNT group than the ChAd-BNT group.

Discussion

Our study adds to the molecular understanding of findings that a heterologous ChAdOx1-BNT162b2 vaccination schedule promotes a more robust immune response than two doses of either BNT162b2 or ChAdOx1 (Agrati et al., 2021; Deming and Lyke, 2021; Pozzetto et al., 2021). Although ChAd-BNT heterologous vaccination studies (Barros-Martins et al., 2021; Schmidt et al., 2021b), including from South Korea (Kang et al., 2021; Lim et al., 2021b), have been published, this is the first report investigating the molecular immune response using bulk and single cell RNA-seq as well as the antibody response to different variants, including the fast-spreading and now dominant Omicron.

Surveying vaccine-induced genomic responses in peripheral immune cells through RNA-seq approaches aids the identification of transcriptional signatures associated with effective antibody production (Andreano and Rappuoli, 2021; Arunachalam et al., 2021; Lee et al., 2021b). Our sequencing depth of average 230 million reads per sample permitted the identification of gene signatures linked to the heterologous ChAd-BNT and a comparison with a reference cohort receiving a two dose BNT homologous vaccination (Lee et al., 2022b). The BNT vaccine in the heterologous cohort elicited a stronger interferon signature with an activated JAK-STAT pathway than the second BNT dose in the homologous reference cohort. Our analysis permitted the preferential activation of specific variable germline classes as well as CDR3 classes and an increase of specific IGHV clonal transcripts encoding neutralizing antibodies was preferentially detected in the heterologous ChAd-BNT and the homologous BNT-BNT cohort.

The benefits of adding BNT162b2 as the second dose in heterologous vaccination regimens with ChAdOx1 as the primary dose is well established (Barros-Martins et al., 2021; Schmidt et al., 2021a). In extension, recent studies demonstrated that a three-dose heterologous regimen with two doses of CoronoVac (Sinovac) immunization followed by BNT162b2 was associated with a 1.4-fold increase of neutralizing antibodies against Omicron as compared to a homologous two-dose BNT162b2 regimen (Cheng et al., 2022; Pérez-Then et al., 2022). However, neutralizing activity against Omicron was reduced by 7.1-fold and 3.6-fold compared to the ancestral and Delta variants, respectively. In contrast, compared to the ancestral strain, Omicron targeted antibody titers in our heterologous cohort were less than 3-fold and a 4.5-fold reduction was seen in the BNT homologous reference cohort. Further evolution of COVID-19 variants is an open question, but it is likely that additional variants will emerge, with manufacture of targeted vaccines following their identification. Therefore, it is reasonable to examine vaccine regimens for their ability to induce a range of antibody responses that may provide coverage for emerging variants. Here, the ChAd-BNT regimen appeared to be successful for that purpose. To date, up to 2.5 billion doses of ChAd have been delivered with an additional 500 million doses of alternative adenovirus based vaccines from Johnson and Johnson and Sputnik from Russia (Dolgin, 2022). Data presented here supports that from an immunological perspective, a heterologous adenovirus-based mRNA-based regimen is effective and may be a reasonable cost-effective approach for specific populations where public health officials deem that adenovirus-based vaccine side effects are low risk. A heterologous strategy might also foster enhanced immune responses in immunocompromised individuals.

Limitations of the study

There are several limitations to this study. First, our study population was limited to a specific geographic area (South Korea) and a specific genetic population. Second, most study subjects were healthy females of normal weight. Third, our study used an established ACE2 binding as a proxy to neutralization and not a live or pseudo-SARS-CoV-2 virus assay. Fourth, the reference cohort receiving a homologous two dose mRNA vaccination had a different ethnic background (Austrian Alps), was on average 20 years older and had a higher BMI than the Korean cohort.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact Lothar Hennighausen (lotharh@nih.gov).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

Study population, study design, recruitment and Ethics statement

46 SARS-CoV-2 naïve healthy volunteers (Table S1) were recruited for the study under informed consent. Recruitment and blood sample collection took place between June and September 2021. This study was approved (IRB No 2021-0898) by the Institutional Review Board (IRB) of Asan medical center in Korea. Written informed consent was obtained from all subjects. Participant information was coded and anonymized. This study was determined to impose minimal risk on participants. All methods were carried out in accordance with relevant guidelines and regulations. All research has been have been performed in accordance with the Declaration of Helsinki (https://www.wma.net/policies-post/wma-declaration-of-helsinki-ethical-principles-for-medical-research-involving-human-subjects/). In addition, we followed the ‘Sex and Gender Equity in Research – SAGER – guidelines’ and included sex and gender considerations where relevant.

Method details

Quantification of immunoproteins

Plasma samples from all participants were collected from their blood. After thawing, the samples were centrifuged for 3 min at 2000 g to remove particulates prior to sample preparation and analysis. Human XL Cytokine Luminex Performance Panel Premixed Kit (R&D systems, #FCSTM18) was used to measure proinflammatory proteins (IFN-α, IFN-γ, IL-1α, IL-1β, IL-1Rα, IL-2, IL-3, IL-4, IL-6, IL-7, IL-8, IL-10 and TNF-α), cytokines (IL-12p70, IL-13, IL-15, IL1-7A and VEGF) and chemokine (CCL2, CCL3, CCL4, CCL11 and CXCL10). The analysis was performed according to the manufacturers’ protocols.

COVID-19 serology and ACE2-neutralization assay

A multiplexed solid-phase chemiluminescence assay (Meso Scale Discovery, MD) was evaluated for the detection of IgG binding to various SARS-CoV-2–derived antigens (V-PLEX SARS-CoV-2 Panel 17 (IgG) Kit, K15524U, and Panel 23 (IgG) Kit, K15567U) and the quantification of antibody-induced ACE-2 binding inhibition to various variants’ spike antigens (pseudo-neutralization assay) (V-PLEX SARS-CoV-2 Panel 18 (ACE2) Kit, K15535U, and Panel 23 (ACE2) Kit, K15570U). Plates were coated with the specific antigen on spots in the 96 well plate and the bound antibodies in the samples were then detected by anti-human IgG antibodies or ACE2 conjugated with the MSD SULPHO-TAG which is then read on the MSD instrument which measures the light emitted from the tag.

Extraction of the PBMC and purification of RNA

Whole blood was collected, PBMCs were isolated using density gradient medium, and total RNA was extracted from PBMC and purified using the RNeasy Mini Kit (Qiagen, #74104) according to the manufacturer’s instructions. The concentration and quality of RNA were assessed by an Agilent Bioanalyzer 2100 (Agilent Technologies, CA).

mRNA sequencing (mRNA-seq) and data analysis

The Poly-A containing mRNA was purified by poly-T oligo hybridization from 1 mg of total RNA and cDNA was synthesized using Super-Script III (Invitrogen, MA). Libraries for sequencing were prepared according to the manufacturer’s instructions with TruSeq Stranded mRNA Library Prep Kit (Illumina, CA, RS-20020595) and paired-end sequencing was done with a NovaSeq 6000 instrument (Illumina) yielding 200–350 million reads per sample.

The raw data were subjected to QC analyses using the FastQC tool (version 0.11.9) (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). mRNA-seq read quality control was done using Trimmomatic (Bolger et al., 2014) (version 0.36) and STAR RNA-seq (Dobin et al., 2013) (version STAR 2.5.4a) using 150 bp paired-end mode was used to align the reads (hg19). HTSeq (Anders et al., 2015) (version 0.9.1) was to retrieve the raw counts and subsequently, Bioconductor package DESeq2 (Love et al., 2014) in R (https://www.R-project.org/) was used to normalize the counts across samples and perform differential expression gene analysis. Additionally, the RUVSeq (Risso et al., 2014) package was applied to remove confounding factors. The data were pre-filtered keeping only genes with at least ten reads in total. The visualization was done using dplyr (https://CRAN.R-project.org/package=dplyr) and ggplot2 (Wickham, 2009). The genes with log2 fold change >1 or < -1 and adjusted p value (pAdj) < 0.05 corrected for multiple testing using the Benjamini-Hochberg method were considered significant and then conducted gene enrichment analysis (GSEA, https://www.gsea-msigdb.org/gsea/msigdb).

For T- or B-cell receptor repertoire sequencing analysis, trimmed fastq files from bulk RNA-seq were aligned against human V, D and J gene sequences using the default settings with MiXCR (Bolotin et al., 2015, 2017). CDR3 sequence and the rearranged BCR/TCR genes were identified. The diversity of BCR/TCR genes was investigated by the Chao1 index (Chao, 1984).

Single-cell RNA sequencing (scRNA-seq) and data analysis

The isolated PBMCs were frozen in freezing media (Thermofisher) and stored at −80°C until use. Single-cell suspensions were then immediately loaded on the 10X Genomics Chromium Controller with a loading target of 20,000 cells. Libraries were generated using the Chromium Next GEM Single Cell 3′ Kit v3.1 (Dual Index) according to the manufacturer’s instructions. Libraries were sequenced using the NovaSeq 6000 instrument (Illumina).

The raw reads were aligned and quantified using the Cell Ranger with Feature Barcode addition (version 5.0, 10X Genomics) against the GRCh38 human reference genome. The quality control, normalization, dimension reduction, cell clusters, UMAP projection, and cell type annotation were performed using Seurat (version 4.0) (Hao et al., 2021).

Statistical analysis

Differential expression gene (DEG) identification used Bioconductor package DESeq2 in R. P-values were calculated using a paired, two-side Wilcoxon test and adjusted p value (pAdj) corrected using the Benjamini–Hochberg method. Genes with log2 fold change >1 or < -1, pAdj <0.05 and without 0 value from all sample were considered significant. For significance of each GSEA category, significantly regulated gene sets were evaluated with the Kolmogorov-Smirnov statistic. p-values of cytokines were calculated using two-stage linear step-up procedure of Benjamini, Krieger and Yekutieli on GraphPad Prism software (version 9.0.0). For comparison of RNA expression levels, antibody levels or cytokine levels between two groups, data were presented as SD in each group and were evaluated with a two-way ANOVA followed by Tukey’s multiple comparisons test using GraphPad PRISM (version 9.0). A value of ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001 was considered statistically significant.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Biological samples
Human serum	Seoul, Korea	N/A
Human immune cells	Seoul, Korea	N/A
Human immune cells	Austria, Tyrol	N/A
Critical commercial assays
RNeasy Mini Kit	Qiagen	Cat# 74104
Human XL Cytokine Luminex Performance Panel	R&D systems	Cat# FCSTM18
V-PLEX SARS-CoV-2 Panel 17 (IgG) Kit	Meso Scale Discovery	Cat# K15429U
V-PLEX SARS-CoV-2 Panel 23 (IgG) Kit	Meso Scale Discovery	Cat# K15567U
V-PLEX SARS-CoV-2 Panel 18 (ACE2) Kit	Meso Scale Discovery	Cat# K15535U
V-PLEX SARS-CoV-2 Panel 23 (ACE2) Kit	Meso Scale Discovery	Cat# K15570U
TruSeq Stranded mRNA Library Prep Kit	Illumina,	Cat# RS-20020595
Deposited data
Raw and analyzed data	This paper	GEO: GSE201535
BNT162b2 homologous vaccination data	Lee et al., 2021b	GEO: GSE190747
Human reference genome UCSC, hg19	UCSC Genome Browser	http://hgdownload.soe.ucsc.edu/downloads.html#mouse
Software and algorithms
MSD DISCOVERY WORKBENCH analysis software		https://www.mesoscale.com/en/products_and_services/software
Trimmomatic (version 0.36)	Bolger et al., 2014	http://www.usadellab.org/cms/?page=trimmomatic
STAR (2.5.4a)	Dobin et al., 2013	https://anaconda.org/bioconda/star/files?version=2.5.4a
HTSeq	Anders et al., 2015	https://htseq.readthedocs.io/en/master/
R (3.6.3)		https://www.R-project.org/
Bioconductor		https://www.bioconductor.org/
DESeq2	Love et al., 2014	https://bioconductor.org/packages/release/bioc/html/DESeq2.html
RUVSeq	Risso et al., 2014	https://bioconductor.org/packages/release/bioc/html/RUVSeq.html
dplyr		https://CRAN.R-project.org/package=dplyr
ggplot2	Wickham, 2009	https://ggplot2.tidyverse.org/
GSEA		https://www.gsea-msigdb.org/gsea/msigdb
MiXCR	Bolotin et al., 2015	https://mixcr.readthedocs.io/en/master/
Cell Ranger (5.0)	10x Genomics	https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger
Seurat (4.0)	Hao et al., 2021	https://satijalab.org/seurat/index.html
PRISM GraphPad (9.0.0)		https://www.graphpad.com/scientific-software/prism/