Abstract
Fast-spreading variants of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) energize the COVID-19 pandemic. B.1.1.7 (VOC-202012/01) has become the predominant variant in many countries and a new lineage (VOC-202102/02) harboring the E484K escape mutation in the B.1.1.7 background emerged in February 20211. This variant is of concern due to reduced neutralizing activity by vaccine-elicited antibodies2,3. However, it is not known whether this single amino acid change leads to an altered immune response. Here, we investigate differences in the immune transcriptome in hospitalized patients infected with either B.1.1.7 (n=28) or B.1.1.7+E484K (n=12). RNA-seq conducted on PBMCs isolated within five days after the onset of COVID symptoms demonstrated elevated activation of specific immune pathways, including JAK-STAT signaling, in B.1.1.7+E484K patients as compared to B.1.1.7. Longitudinal transcriptome studies demonstrated a delayed dampening of interferon-activated pathways in B.1.1.7+E484K patients. Prior vaccination with BNT162b vaccine (n=8 one dose; n=1 two doses) reduced the transcriptome inflammatory response to B.1.1.7+E484K infection relative to unvaccinated patients. Lastly, the immune transcriptome of patients infected with additional variants (B.1.258, B.1.1.163 and B.1.7.7) displayed a reduced activation compared to patients infected with B.1.1.7. Acquisition of the E484K substitution in the B.1.1.7 background elicits an altered immune response, which could impact disease progression.
Introduction
Over the past few months, several variants of concern (VOC) carrying specific mutations thought to enhance viral fitness have emerged. Specifically, B.1.3514 and P.1 were of particular concern because they carry the mutation E484K within the receptor binding domain (RBD), which has been demonstrated to enhance escape from neutralizing antibody inhibition in vitro2,3 and may be linked with reduced vaccination efficacy5. In February 2021 Public Health England (PHE) published a report of B.1.1.7 genomes with acquisition of the E484K spike mutation1. Since the PHE announcement a total of 1332 B.1.1.7+E484K genomes have been uploaded to GISAID from England and 14 other countries, including 203 from the United States where it was first detected in Oregon (as of 05/24/2021).
The new B.1.1.7+E484K variant (VOC 202102/02) is fast spreading, prompting interest in the response of both unvaccinated and vaccinated individuals to disease from this variant. Recent in vitro data demonstrated that introduction of the E484K mutation into the B.1.1.7 background led to a more-substantial loss of neutralizing activity by vaccine-elicited antibodies compared with the mutations in B.1.1.7 alone2,3, suggesting that this variant represents a threat to the efficacy of the BNT162b2 vaccine. However, it remains unclear whether the B1.1.7+E484K variant elicits a modified immune response in patients compared with the parent B.1.1.7 variant.
To obtain insight into the biological significance of the E484K mutation in the B.1.1.7 backbone, we investigated the transcriptome of a total of 40 hospitalized patients infected with either the B.1.1.7+E484K variant (n=12) or the B.1.1.7 parent strain (n=28). We also investigated the impact of the BNT162b vaccine on both variants.
Results
Immune transcriptome progression after B.1.1.7 infection
Transcriptomes from 28 patients infected with the B.1.1.7 variant and 12 patients infected with the B.1.1.7+E484K variant were compared (Table 1). First, we investigated the temporal progression of the immune transcriptomes in 28 patients hospitalized with infection of the B.1.1.7 variant (Table 1, Fig. 1a, Supplementary Table 1). RNA-seq experiments were conducted on PBMCs isolated between days 1-5 (group A), days 10-14 (group B) and days 15-30 (group C) after reporting symptomology (Supplementary Data 1-4). The presence of SARS-CoV-2 was confirmed by PCR, followed by whole viral genome sequencing. PCA plots demonstrate a separation of the three cohorts, with group C shifting towards the non-COVID controls (Fig. 1b). Expression of a total of approximately 3,071 genes was significantly elevated within five days of reporting symptoms and these clustered in immune-relevant pathways (Fig. 1c). Expression of immune-related gene classes, including the highly activated JAK/STAT pathway, declined up to 95% after 15 days upon onset of COVID-19 (Fig. 1c and d, Supplementary Fig. 1).
Comparison of early response immune transcriptomes of patients infected with B.1.1.7+E484K to those of patients infected with B.1.1.7
Next, we investigated whether the E484K escape mutant in the backbone of B.1.1.7 (B.1.1.7+E484K) leads to an altered immune response. For this, we compared the immune transcriptomes of patients infected with either the B.1.1.7+E484K or the B.1.1.7 parent variant within the first five days after the onset of COVID-19 symptomology. All patients had had been hospitalized and there were no significant differences in age distribution, gender, BMI, vaccination status, treatment, ICU and deaths between the two groups (Table 1). RNA-seq was conducted on seven B.1.1.7+E484K patients and 14 B.1.1.7 patients (Fig. 2; Supplementary Data 2, 5, 6). PCA plots showed separation between RNA-seq samples from non-COVID controls6 and the two variants on first and second principal components (PC1 and PC2) (Fig. 2a). A higher number of differentially expressed genes (DEGs) were identified between patients infected with the B.1.1.7+E484K variant compared to controls (4399) than found in the comparison between patients infected with the B.1.1.7 variant and controls (3071) (Fig. 2b). Hallmark Gene Set Enrichment Analysis (GSEA) established that the significantly up-regulated genes in both cohorts were enriched in immune response pathways, including IL-JAK/STAT signaling and IFNα/γ responses (Fig. 2c). Pathway analyses identified genes that are induced in lung tissue from COVID-19 patients and in cells infected by SARS-CoV-2 (Fig. 2d). Genes preferentially activated in the B.1.1.7+E484K variant are part of the interferon and JAK/STAT pathways (Fig. 2d-e; Supplementary Fig. 2-4). These experiments demonstrate that the acquisition of the E484K escape mutation in the B.1.1.7 background elicits a significantly different immune response than the parent variant itself.
Differences in temporal progression of the immune response in patients following B.1.1.7+E484K as compared to B.1.1.7 infection
Next, we utilized temporally matched samples from patients infected with B.1.1.7 (n=5) and B.1.1.7+E484K (n=4) to compare temporal progression of the immune transcriptome at days 1-5 versus 10-14 following symptomology (Supplementary Data 7-9). PCA plots demonstrated that while the transcriptomes from B.1.1.7 patients diverged on both PC1 and PC2 (Fig. 3a, left panel), the transcriptomes from B.1.1.7+E484K patients diverged only on PC2 (Fig. 3a, right panel). At days 1-5, 294 genes were preferentially expressed in B.1.1.7+E484K patients compared to B.1.1.7 with genes mapping to interferon and JAK/STAT pathways (Fig. 2e; Supplementary Fig. 3-4). The pronounced interferon and inflammatory responses seen at days 1-5 sharply declined between days 10-14 (Fig. 3b, c). Notably, the decrease was more pronounced in the B.1.1.7 variant with a more sustained immune response in the B.1.1.7+E484K variant. Since B.1.1.7+E484K induced a more expansive immune transcriptome than the B.1.1.7 parent variant, we directly compared their transcriptomes between days 10-14 after symptomology (Fig. 3d-e). The transcriptomes separate on the PCA plot (Fig. 3d) and the significantly higher expressed genes in B.1.1.7+E484K are enriched in cytokine responsive signaling and viral, including SARS-CoV-2, infection pathway (Fig. 3e-f). Of the 72 genes enriched in COVID-19 pathway (Fig. 3f), the top 50 induced gene in the B.1.1.7+E484K variant compared to B.1.1.7 at the second time point are related to interferon and JAK/STAT signaling (Supplementary Fig. 5).
Immune transcriptomes in vaccinated patients infected with B.1.1.7 and B.1.1.7+E484K
Next, we investigated the impact of the BNT162b vaccine on the immune transcriptome in the seven patients that had received at least one dose between 6 and 24 days prior to demonstrating COVID-19 symptomology and one fully vaccinated patient (Supplementary Data 10-12). Three vaccinated patients infected with B.1.1.7 died in the hospital (vaccination days 8, 10, 14 prior to symptoms). The transcriptomes from the three B.1.1.7+E484K patients who received their first dose between 8 and 22 days prior to developing symptoms were compared with six unvaccinated B.1.1.7+E484K patients (Fig. 4a). Expression of 19 genes was significantly reduced in the vaccinated cohort as compared to the unvaccinated one (Fig. 4b). These genes map to interferon, JAK/STAT and TNFα pathways (Fig. 4c-d). Differences were reflected in a significant difference on the PC2 but not PC1 component of the PCA plot (Supplementary Fig. 6). A direct comparison of the transcriptomes from three B1.1.7 and three B.1.1.7+E484K patients that had received at least one dose highlighted the different responses (Fig. 4e-g).
Immune transcriptomes in patients infected with other SARS-CoV-2 variants
We also examined the immune transcriptome of patients infected with the variants B.1.177 (n=1), B.1.1.163 (n=2) and B.1.258 (=2) within the first five days of symptomology (Supplementary Table 1; Fig. 5a; Supplementary Data 13-16). COVID-19 patients are distributed by the type of variant on the PCA plot (Fig. 5b). The B.1.1.7 variant elicited stronger interferon responses as indicated by elevated JAK/STAT and inflammatory pathways compared to the other variants (Fig. 5c). Key components of ISGs are expressed highest in patients infected with B.1.1.7+E484K and B.1.1.7 compared to the other variants (Fig. 5d).
Discussion
As the COVID-19 epidemic has evolved, the impact of variants on disease progression and response to vaccine has become a key area of investigation2,3. The B.1.1.7 variant, first detected in Kent (UK) in September 2020, has spread rapidly and has become the predominant form in many countries7,8. The RBD mutation, E484K, first detected in the B.1.351 variant, has evolved independently in infected individuals in different genetic SARS-CoV-2 backgrounds. Here we were able to explore the impact of this escape mutation on the background of the B.1.1.7 variant1,9, utilizing samples from cohorts of hospitalized patients with the same demographics that resided in specific geographic regions and developed infection within the same timeframe, enabling a relatively controlled focus on the impact of this single amino acid substitution.
The clinical course of SARS-CoV-2 infection has been well-characterized, ranging from asymptomatic to fatal disease10. Here we concentrated on the immune transcriptome response of cohorts of more severely ill hospitalized patients and followed response over three weeks from the time of symptom development. In these patients we found that the host transcriptome response to SARS-CoV-2 infection evolves rapidly with the most profound changes within five days of symptomology, dropping significantly over the next ten days, with further drop, but not normalization, within four weeks. In a previous study we documented that asymptomatic and minimally symptomatic SARS-CoV-2 infected individuals from the same general geographic area showed no significant differences to uninfected individuals by six weeks following infection6. This means that studies of transcriptome responses following SARS-CoV-2 infection need to be controlled for time of infection and symptom development for valid comparisons. There was no significant difference in deaths between B.1.1.7 and B.1.1.7+E484K infection. Five and three individuals died in B.1.1.7 and B.1.1.7+E484K groups, respectively. Our study provided the opportunity to examine whether or not transcriptomes from those that survived were significantly different than those that expired from infection. No differences were found. All individuals who died were over the age of 70 and had underlying medical conditions. Drug treatment including corticosteroid and antibody treatment could theoretically alter immune response. In this study the majority of this hospitalized population received corticosteroid treatment, which may have modified immune response, but as there were no significant differences in the percentage treated between B.1.1.7 and B.1.1.7+E484K patients, this was not the reason for the significant differences in transcriptome response found. Only one individual within the cohort received therapy with convalescent plasma before the first blood sampling. Interestingly, the PCA demonstrated that this individual’s transcriptome mapped to a location different than the others. Due to this being a single sample, we cannot make any generalizations from this study about how antibody treatment may modify transcriptome response; however, this is an area of interest for future investigation if sufficient numbers of temporally and variant matched samples are available.
Prior vaccination, even a single dose, significantly modified transcriptome response, even as it did not prevent hospitalization. However, the significant differences in transcriptome response between B.1.1.7 and B.1.1.7+E484K infected individuals persisted even if these vaccinated individuals, documenting the reproducibility of the finding across clinical settings. Infection following first vaccination is well documented with some individuals also becoming severely ill after full vaccination11. Our study also included one individual who was hospitalized with infection following full vaccination but there were no significant differences in transcriptome response with the other vaccinated patients. Antibody titers were not available for any of the vaccinated individuals, but they did not fall into any of the currently recognized groups with inadequate immune response following vaccine.
A small number of samples from individuals infected with (B.1.351) and without (B.1.177, B.1.1.163, B.1.258) the E484K mutation (B.1.351) were available for comparison. Cohort sizes were not sufficient for statistical analyses of significance but temporally-matched transcriptomes did show variation from B.1.1.7 and B.1.1.7+E484K. A recent bulk RNA-seq immune transcriptome study also identified transcript categories, including interferon-activated genes, induced in COVID-19 patients but the variant was not defined12. As samples become available from patients with other variants, the degree of transcriptome variation will be able to be defined.
In summary, our study documents that the E484K mutation is sufficient to alter transcriptome response. This RBD escape mutation has evolved independently in several SARS-CoV-2 backgrounds, including B.1.1.79 with a single immunocompromised patient study demonstrating that it can occur within 10 days of infection13. Here, our significant differences in transcriptome response did not correlate with significant differences in clinical course but this could be due to the selection bias of only hospitalized patients being included in this study. The fact that variants can induce different host immune response has implications for the understanding of both innate and acquired immune response mechanisms. Further investigation of the reasons why these differences occur may help in design of vaccines and antibody treatments that may have more universal coverage of variant.
Limitations
The findings in this report are subject to several limitations. Our study focused on hospitalized patients in a specific geographic area. Our study focused on elderly patients with limited comparison to younger individuals. We have investigated the effect of the BNT162b vaccine but not of any other vaccine.
Methods
SARS-CoV-2 virus sequencing
RNA was extracted from patient’s blood using a Maxwell RSC simply RNA Blood purification kit according to the manufacturer’s instructions (Promega, USA). Library preparation and sequencing was performed as described14. In short, cDNA was obtained by using reverse transcriptase with random priming. Following cDNA synthesis, primers based on sequences from the ARTICnetwork were used to generate 400 bp amplicons in two different PCR pools. After merging of pools and amplification, libraries were constructed using QIASeq FX DNA Library UDI Kit following the manufacturer’s instructions (Qiagen GmbH, North Rhine-Westphalia, Germany).
Sequencing was done with Illumina NextSeq® 500/550 using 149-bp paired-end reads with 10-bp indices (Illumina, California, USA). Obtained viral sequences were assembled using CLC Genomics Workbench v20.0.3 (Qiagen GmbH, North Rhine-Westphalia, Germany). SARS-CoV-2 isolate Wuhan-Hu-1 served as the reference genome (Accession NC_045512.2). SARS-CoV-2 variants were identified by uploading FASTA-files on freely accessible databases (http://cov-lineages.org/).
Extraction of the buffy coat and purification of RNA
Whole blood was collected, and total RNA was extracted from the buffy coat and purified using the Maxwell RSC simply RNA Blood Kit (Promega) 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 μg of total RNA and cDNA was synthesized using SuperScript 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).
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 Trimmomatic15 (version 0.36) and STAR RNA-seq16 (version STAR 2.5.4a) using 150 bp paired-end mode was used to align the reads (hg19). HTSeq17 (version 0.9.1) was to retrieve the raw counts and subsequently, R (https://www.R-project.org/), Bioconductor18 and DESeq219 were used to normalize the counts across samples. Additionally, the RUVSeq20 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 ggplot221. Genes were categorized as significant differentially expressed with an adjusted p-value (pAdj) below 0.05 and a fold change > 2 for up-regulated genes and a fold change of < −2 for down-regulated ones with no zero values in all data set and then conducted gene enrichment analysis (GSEA, https://www.gsea-msigdb.org/gsea/msigdb).
Data Availability
The RNA-seq data of Non-COVID individual were obtained from GSE162562. The RNA-seq data of COVID-19 patients infected by B.1.1.7 or B.1.1.7/E484K variant will be uploaded in GEO before publishing the manuscript.
Statistical analysis
For significance of each GSEA categoty, significantly regulated gene sets were evaluated with the Kolmogorov-Smirnov statistic. Demographic data were analyzed by Chi-square on GraphPad Prism software (version 9.0.0). All graphs were generated using GraphPad.
Ethics approval
This study was approved by the Institutional Review Board (IRB) of the Office of Research Oversight / Regulatory Affairs, Medical University of Innsbruck, Austria (EC numbers: 1064/2021).
Data availability
The RNA-seq data of Non-COVID individual were obtained from GSE162562. The RNA-seq data of COVID-19 patients infected by B.1.1.7 or B.1.1.7/E484K variant will be uploaded in GEO before publishing the manuscript.
Author contribution
LK: recruited patients, collected material, analyzed data; HKL: analyzed data, wrote manuscript; MW: collected and prepared material; AM, AZ, LK Sr, JS and NK: recruited and diagnosed patients; PAF: analyzed data, wrote manuscript; LH: analyzed data, wrote manuscript. All authors read and approved the manuscript
Competing interests
The authors declare no competing financial interests.
Acknowledgments
Our gratitude goes to the patients who contributed to this study to advance our understanding of COVID-19. This work was supported by the Intramural Research Programs (IRPs) of National Institute of Diabetes and Digestive and Kidney Diseases (NIDDK) and utilized the computational resources of the NIH HPC Biowulf cluster (http://hpc.nih.gov). RNA-sequencing was conducted in the NIH Intramural Sequencing Center, NISC (https://www.nisc.nih.gov/contact.htm).