Predicted CTL responses from pressured epitopes in SARS-CoV-2 correlate with COVID-19 severity

CTL epitope prediction

Putative CTL epitopes in the SARS-CoV-2 proteome were predicted using the well-established NetMHCPanBA method made available by IEDB. Epitope predictions were carried out for the top 10 abundant proteins in the viral proteome, whose complete sequences were documented on the NCBI Virus database for (i) the 4 Variants of Concerns (VOCs) namely B.1.1.7 (Alpha), B.1.351 (Beta), P.1 (Gamma) and B.1.617.2 (Delta), (ii) 4 Variants of Interest (VOIs) namely B.1.525 (Eta), B.1.526 (Iota), B.1.617.1 (Kappa) and C.37 (Lambda) (as documented in the official WHO website as of August 2021) along with the Wuhan strain (A) and (iii) 2 non-virulent variants – B.1.526.1 and P.2. Predicted epitopes binding with their cognate HLA-alleles with IC50 < 50nM were classified as strong binding epitopes. Strong binding epitopes were predicted individually for all HLA Class I alleles documented across various ethnicities around the world in The Allele Frequency Net Database²⁵. A total of 2022 strong binding epitopes exhibiting high affinity binding with 1299 HLA-alleles were identified (Supplementary File-1). This epitope pool was pruned to identify (a) epitopes conserved across all variants and (b) pressured epitopes – that are under negative selection pressure. A cohort from Bangalore consisting of 36 individuals with COVID-19 severity ranging from asymptomatic and mild to moderate and severe was analysed to obtain a theoretical measure of the CTL responses in each case which were then tested for any correlations with disease severity. CTL responses were predicted by identifying the HLA Class-I genotypes of each individual in the cohort, and identifying the pool of high-binding pressured epitopes, which was then consolidated in to a CTL response score. A schematic overview of the workflow is depicted in Figure 1.

• Download figure
• Open in new tab

Figure 1:

CTL Response Analysis Overview for SARS-CoV-2 Reference strain (Wuhan, Pangolin Lineage A) in comparison with SARS-CoV-2 Variants of Concern (VOCs), Variants of Interest (VOCs) and 2 non-virulent variants - B.1.526.1 and P.2. Using bioinformatic epitope prediction tools, CTL epitopes were predicted for the 10 SARS-CoV-2 proteins using the NetMHCPanBA method²⁸. A total of 1827 HLA Class I alleles were used to predict binding affinities with SARS-CoV-2 epitopes. The initial epitope pool thus obtained was filtered based on peptide-HLA binding affinity (IC50 < 50 nM). From this shortlisted pool, epitopes under selection pressure across variants were identified and used for further analysis. COVID-19 blood samples from patients with a wide range of disease severity were collected and their HLA-genotypes were identified via HLA-typing. Finally, we performed correlation analysis between recognition of pressured epitopes and COVID-19 severity in patients.

Immune pressure on CTL epitopes across SARS-CoV-2 variants

After shortlisting strong binding epitopes, we compared these epitope sets across variants of SARS-CoV-2. A pairwise estimate of difference in epitope sets as well as the protein sequence distance for each protein (Figure 2A, B) revealed higher variation in epitope sets compared to whole-protein sequence across the variants. This indicated that most proteins of SARS-CoV-2 are under immune pressure (Figure 2A, B, Supplementary Figure 1), except for Membrane protein and ORF6 (Supplementary Figure 1) which showed relatively higher conservation across variants. Immune pressure at the sequence level was evidenced by significant overlap between regions of high epitope density and higher variation, calculated as the Shannon entropy (Figure 2A, B).

• Download figure
• Open in new tab

Figure 2:

Analysis of CTL epitopes under immune pressure. Multiple sequence alignment of a representative epitope under immune pressure across the 11 variants of SARS-CoV-2 is depicted (top) for (A) Spike glycoprotein and (B) ORF3a protein. Symmetric difference between epitope sets for each pair of variants is plotted against the whole-protein sequence distance (bottom). The points in the scatter plot lie above the diagonal (y=x) for both proteins indicating higher selection pressure on epitopic regions. Sequence variation in the protein is computed as the Shannon entropy at each residue of the protein (middle - red). Distribution of epitopes within the protein is represented as the epitope density at each residue (middle - green). (C) Number of pressured epitopes in each protein of SARS-CoV-2, recognized by the HLA alleles (top 20 based on the number of epitopes recognized) are represented as a heat map (top). Contribution of each protein to the pressured epitope set is depicted by an inverted bar plot. ORF1ab is excluded from this plot since it is a polyprotein. Spike glycoprotein is seen to be the major contributor to the pressured epitope set.

Based on the observed variation across SARS-CoV-2 variants, from the pool of all strong binding epitopes, two subsets were selected: (a) epitopes which are conserved across all variants and (b) epitopes which are under selection pressure presumably due to host CTL responses and hence show sequence variation across variants. These pressured epitopes majorly belonged to the Spike glycoprotein (Figure 2C) while ORF6 lacked pressured epitopes as suggested by its high conservation across SARS-CoV-2 variants. We observed that HLA-B*15:62 followed by HLA-A*02:11 recognize the highest number of pressured epitopes among the HLA-Class I alleles considered. A ranked list of top-responding HLA-alleles is provided in Supplementary Table 1.

Correlating predicted CTL responses with COVID-19 disease severity

Next, we tested if the pressured epitopes suggest the possibility of conferring protection against severe COVID-19 disease. Whole-blood samples from COVID-19 patients (n=36) with varying severity of the disease (Asymptomatic, Mild, Moderate and Severe) were collected from patients in Victoria Hospital, Bangalore, India and their HLA Class I genotypes were identified through their whole-blood RNA-seq data (Figure 3A). The predicted HLA genotypes for each of the samples in this Bangalore cohort are provided in Supplementary Table 2. The CTL response of an individual was calculated as the number of distinct epitopes recognized by the HLA alleles comprising the HLA genotype.

• Download figure
• Open in new tab

Figure 3:

Correlation of predicted CTL responses with COVID-19 severity. (A) HLA-genotypes of COVID-19 patients from the Bangalore cohort, with 4 levels of disease severity - Asymptomatic, Mild, Moderate and Severe were identified using HLA-typing. CTL response measured in terms of the number of pressured epitopes (B) and (C) total number of epitopes recognized by individuals is plotted against the corresponding COVID-19 severity of the individual. Asymptomatic individuals recognize a significantly larger number of pressured epitopes as compared to those with severe COVID-19. (D) Schematic representation of computation of epitope protectivity score: Let [HG_Aj]_1≤j≤7, [HG_Mij]_1≤j≤9, [HG_Moj]_1≤j≤10 and [HG_Sj]_1≤j≤10 be the sets of individuals with asymptomatic, mild, moderate and severe COVID-19 respectively, represented by their HLA-genotypes. Each HLA genotype is associated with a corresponding set of distinct epitopes where n is the number of epitopes recognized by the HLA-alleles comprising the genotype and n_s is the number of individuals with disease severity s. Since epitope sets recognized by individuals with different disease severities overlap in many cases, the frequency of association of an epitope E with a disease severity DS, f_E,DS where DS ∈ {A, Mi, Mo, S} is calculated for each epitope. A protectivity value is then assigned for each disease state ranging from 3 for asymptomatic cases to 0 for severe cases. The protectivity score for an epitope is computed as the sum of the ‘protectivity value’ for each disease severity weighted by the frequency of association of the epitope with the particular disease severity (see Methods). (E) Protective CTL response calculated as the number of protective epitopes recognized by an individual is plotted against the corresponding COVID-19 severity of the individual. The number of protective epitopes recognized by asymptomatic individuals is significantly higher than those with mild, moderate or severe COVID-19. (F) The number of protective epitopes under immune pressure is plotted against disease severity of individuals recognizing these epitopes. Here, the number of pressured and protective epitopes recognized by asymptomatic individuals is significantly higher than those with moderate or severe COVID-19. (G) Protectivity scores of HLA alleles (top) based on the disease severity they associate with, along with the distribution of disease severity of patients expressing the allele (bottom) is shown. HLA alleles with higher protectivity are largely seen in Asymptomatic individuals while those with lower protectivity are observed in individuals with severe disease. Statistical significance for the bar plots was computed using the Wilcoxon paired signed-rank test (NS - Not Significant, * - p-value ≤ 0.05, ** - p-value ≤ 0.01, *** - p-value ≤ 0.001).

The resulting epitope pool of the entire cohort, referred to as the ‘all epitope’ pool, was further classified in 2 ways. Firstly, a set of locally pressured epitopes was shortlisted as a subset of the initially characterized globally pressured epitopes recognized by individuals within the Bangalore cohort. Secondly, epitopes which offered higher protectivity against severe COVID-19 were selected. Epitopes recognized by each HLA genotype within the cohort were assigned a protectivity score based on the severity of disease exhibited by the individual (Figure 3D). A pool of ‘protective’ epitopes was then constructed by shortlisting epitopes whose protectivity score crosses a set threshold corresponding to Asymptomatic and Mild disease states (see Methods).

We observed significantly higher recognition of pressured CTL epitopes by asymptomatic individuals as compared to severe cases (Figure 3B). In contrast, no significant difference was observed across disease states with the same analysis but using the ‘all epitope’ pool (Figure 3C). Our analysis shows that the overall CTL response decreases in the order of disease severity, with highest in asymptomatic and least in severe cases, suggesting that the HLA genotypes in asymptomatic individuals are the most protective in our cohort.

Next, we assigned a protectivity value for each disease severity, ranging from 0 for severe to 3 for asymptomatic disease states. We computed a protectivity score for each epitope as the sum of protectivity values of the 4 disease states, weighted by the frequency of occurrence of the epitope in individuals with the particular disease state. We observed a stronger correlation between the size of the protective epitope pool recognized and COVID-19 severity. The number of protective epitopes recognized by asymptomatic individuals was significantly higher than those recognized by mild, moderate and severe cases (Figure 3E). A pressured-protective epitope pool was then constructed by shortlisting pressured epitopes with high protectivity scores. A similar trend was observed where the size of the pressured-protective epitope pool was significantly higher in asymptomatic individuals as compared to those with moderate and severe disease (Figure 3F). A list of these pressured-protective epitopes is provided in Supplementary Table 3.

Next, we studied the association of individual HLA-alleles with disease severity in COVID-19 patients in the Bangalore cohort. We observed that many HLA-alleles mapped to multiple disease states in the cohort (Figure 3G). Similar to epitopes, a protectivity score was assigned to each HLA allele based on disease severity in the individual expressing the allele. HLA alleles A*02:03, A*02:06, A*03:02, A*29:01, B*07:05, B*38:02 and C*15:05 were seen to confer the highest protection while alleles B*07:06, B*49:01 and C*07:01 offered the least protection against severe disease. A list of HLA genotypes which offer highest protection against COVID-19 is provided in Supplementary Table 4. A gradient of disease states from completely asymptomatic in case of alleles with higher protectivity, to purely severe in case of alleles with least protectivity can be observed (Figure 3G) suggesting a major role of HLA Class-I mediated presentation of a selected epitope pool to CD8+ T cells in determining disease outcome.

SARS-CoV-2 protein sequences and HLA Class-I allele data

A total of 11 SARS-CoV-2 variants including 4 Variants of Concern (B.1.1.7 - Alpha, B.1.351 - Beta, P.1 - Gamma, B.1.617.2 - Delta), 4 Variants of Interest (B.1.525 - Eta, B.1.526 - Iota, C.37 - Lambda, B.1.617.1 - Kappa) as listed in the official WHO webpage as of August, 2021 along with 2 other Variants namely, B.1.526.1 and P.2 and the Wuhan strain (A) were considered for our study. Amino acid sequences of 10 proteins, namely - ORF1ab polyprotein, Spike glycoprotein, ORF3a, Envelope protein, Membrane protein, ORF6, ORF7a, ORF8, Nucleocapsid protein and ORF10 were obtained from NCBI Virus²³ (https://www.ncbi.nlm.nih.gov/labs/virus/vssi/#/).

A total of 1827 HLA Class-I alleles were obtained from 240 ethnic groups with HLA-allele frequency data from the Allele Frequency Net Database²⁵ (AFND) (http://www.allelefrequencies.net/). The ethnic groups were shortlisted based on i) documentation of frequencies of all 3 HLA Class I alleles - A, B and C, ii) High resolution HLA allele data with polymorphisms denoted for all the alleles documented.

Prediction of CTL epitopes

HLA Class I Epitope predictions for the 10 proteins were performed using the NetMHCPan BA 4.1 tool²⁸ made available by IEDB (Immune Epitope Database)²⁴ (http://tools.iedb.org/mhci/). A comprehensive list of 1827 HLA Class I alleles covered by the 240 ethnic groups were considered for pMHC binding predictions. Length of the predicted epitopes was restricted to 9 amino acids. Predicted 9-mers which bind to a given HLA allele with IC50 < 50 nM were considered as strong binding epitopes. Epitope predictions were performed for all 11 variants of SARS-CoV-2.

Immune pressure parameters

For the immune pressure plot, 2 parameters were considered for each protein: the protein sequence pairwise distance and the epitope set pairwise distance between each pair of variants. These parameters are defined as follows:

Let V₁ and V₂ be the two variants being compared and let P₁ = [a_1i]_1≤i≤n and P₂ = [a_2j]_1≤j≤n be the aligned amino acid sequences of a particular protein belonging to the 2 variants with an alignment length n. Let E₁ = [e_1a]_1≤a≤p and E₂ = [e_2b]_1≤b≤q be the sets of distinct epitopes of the protein corresponding to the variants V₁ and V₂ respectively, where p and q are the respective number of epitopes of the protein. These 2 parameters are calculated for each pair of variants and are represented as points in a 2D scatter plot (Figure 2A, B)

Shannon Entropy and Epitope Density

For each residue of a protein, 2 further parameters namely, Shannon Entropy and Epitope Density are calculated as follows:

Let P_V = [p_i]_1≤i≤11 be the set of amino acid sequences for a particular protein across the 11 SARS-CoV-2 variants along with the Wuhan strain considered, where p_i = [a_ij]_1≤j≤n for all 1≤ i ≤ 11 is the amino acid sequence of protein p of variant i with an alignment length of n. Let R_j= [r_ij]_1≤i≤11 be the set of unique amino acids across variants for protein p at residue j and where m_kj is the number of variants possessing amino acid a_kj at position j.

Support from Literature for the CTL epitopes

Out of 577 experimentally validated epitopes from over 25 publications²⁹, 482 high binding epitopes appeared in our study. The total epitope pool of 2022 epitopes in case of the SARS-CoV-2 Wuhan strain could not be covered completely due to binding of some epitopes to indigenous HLA alleles, which have not been validated experimentally in previous studies.

COVID-19 sample collection

Ethical approval for this study was obtained from the Institutional Ethics Committee at Bangalore Medical College and Research Institute, Bangalore, India (BMCRI/PS/02/2021-21), and IISc (1-26062020), Bangalore, India. Written informed consent was obtained from all study participants before sample collection.

COVID-19 patients were recruited for this study at the hospital upon obtaining informed consent from the participants, during September-December, 2020. The presence of the disease was confirmed by qPCR for viral RNA in nasal and throat swabs of the individuals. Patients who had prior episodes of COVID-19 infection, or had chronic infections such as tuberculosis, HIV or HCV were excluded. A total of 36 Covid-19 patients (median age = 52 + 16, both male and female participants), volunteered their samples for this study. The disease severity classification was performed based on the clinical symptoms presented by the patient at the hospital, which are as follows: Asymptomatic (n=7), Mild (n=9), Moderate (n=10) and Severe (n=10). Whole blood samples were collected in PAXgene Blood RNA tubes and stored at -20°C until further use.

RNA isolation and HLA Typing

RNA isolation was performed using Paxgene Blood RNA kit (Cat# 762164) following manufacturer protocol. It was then subjected to quantification and quality control, followed by RNA sequencing. 1μg of total RNA from samples with RIN > 5 was used for poly-A mRNA enrichment with NEBNext Poly(A) mRNA Magnetic Isolation Module (NEB #E7490), followed by library preparation with NEBNext Ultra II Directional RNA Library Prep Kit (NEB #E7760) using manufacturer directions. The paired-end RNA sequencing was performed in an Illumina Hiseq 4000 instrument. The library preparation and sequencing was performed by Biokart India Private Limited, Bangalore. The raw sequencing data was subjected to adapter trimming with the tool TrimGalore (version 0.6.6) and Class-I HLAs were predicted from the RNA sequencing reads using the software arcasHLA v0.2.0³⁰ with standard settings.

Epitope and HLA allele protectivity

Protectivity of epitopes and HLA alleles were calculated based on the HLA genotypes corresponding to the clinical samples with varying disease severities.

Let [HG_Aj]_1≤j≤7, [HG_Mij]_1≤j≤9, [HG_Moj]_1≤j≤10 and [HG_Sj]_1≤j≤10 be the sets of individuals with asymptomatic, mild, moderate and severe COVID-19 respectively, represented by their HLA-genotypes. Each HLA genotype is comprised of a set of 6 HLA alleles. Since some HLA-alleles are associated with more than one disease severity, the frequency of association of a HLA allele H with a disease severity DS, f_H,DS where DS ∈ {A, Mi, Mo, S} is calculated as follows: Each HLA genotype is also associated with a corresponding set of distinct epitopes where n is the number of epitopes recognized by the HLA-alleles comprising the genotype and n_s is the number of individuals with disease severity s recognizing the epitope. Since epitope sets recognized by individuals with different disease severities overlap in many cases, the frequency of association of an epitope E with a disease severity DS, f_E,DS where DS ∈ {A, Mi, Mo, S} is calculated for each epitope as follows: A protectivity value is then assigned for each disease state as follows: The protectivity score for a HLA-allele H is then calculated as follows: The protectivity score for an epitope E is calculated as follows: