A bioinformatic prediction of antigen presentation from SARS-CoV-2 spike protein revealed a theoretical correlation of HLA-DRB1*01 with COVID-19 fatality in Mexican population: An ecological approach.

SARS-CoV-2 infection is causing a pandemic disease that is reflected in challenging public health problems worldwide. Human leukocyte antigen (HLA)-based epitope prediction and its association with disease outcomes provide an important base for treatment design. A bioinformatic prediction of T cell epitopes and their restricted HLA Class I and II alleles was performed to obtain immunogenic epitopes and HLA alleles from the spike protein of the severe acute respiratory syndrome coronavirus 2 virus. Also, a correlation with the predicted fatality rate of hospitalized patients in 28 states of Mexico was done. Here, we describe a set of 10 highly immunogenic epitopes, together with different HLA alleles that can efficiently present these epitopes to T cells. Most of these epitopes are located within the S1 subunit of the spike protein, suggesting that this area is highly immunogenic. A statistical negative correlation was found between the frequency of HLA-DRB1*01 and the fatality rate in hospitalized patients in Mexico.

INTRODUCTION

The coronavirus disease (COVID‐19) was declared as a pandemic by the World Health Organization (WHO) in March of 2020. It is estimated that by June the 10th of 2020 there were over 6.19 million confirmed cases and 370,000 deaths worldwide.

COVID‐19 is a disease generated by the novel severe acute respiratory syndrome‐coronavirus‐2 (SARS‐CoV‐2), with a wide range of clinical manifestations, like fever (88.7%), cough (67.8%), fatigue (38.1%), and acute respiratory distress syndrome in severe cases. Interestingly, the molecular and clinical manifestations of the disease vary between asymptomatic, mild‐symptomatic, and severe patients, requiring hospitalization in some cases to prevent fatal outcomes.

Currently, the SARS‐CoV‐2 genome has been characterized as a new betacoronavius, which shares around 87% of genomic identity with bat‐SL‐CoVZC45 and bat‐SL‐CoVZXC21 viruses. A recent analysis by Zhou et al. reported that there is a 96.2% identity with BatCoVRaTG13 and a 79.5% identity with SARS‐CoV. The genomic characterization of the virus not only provides information about its taxonomy and probable origin but also offers opportunities to perform deeper analysis using bioinformatics tools.

The angiotensin‐converting enzyme‐2 (ACE‐2) receptor and the transmembrane serine protease 2 are essential components of the human host for the virus entry into the upper respiratory epithelial cells. The virus recognizes ACE‐2 through the viral spike glycoprotein (S), and this event leads to the virus–cell membrane fusion. The S glycoprotein is found as a homotrimer of three identical monomers, each one of which is divided into two subunits: S1 and S2. The first subunit folds in four domains: A, B, C, and D. The B domain possesses a receptor‐binding domain that recognizes ACE2, hence it is important for viral entry. The S2 subunit sequence has two tandem domains, namely, HR1 and HR2, that play an essential role in the viral fusion to the membrane. Furthermore, analysis of the spike protein showed that it is conserved among SARS‐CoV and SARS‐CoV‐2 with 76.3% of identity and 87.3% of similarity.

Several studies focused on viral diseases have shown that clinical severity is closely associated with some individual factors, such as genetic background and immune response. The human leukocyte antigen (HLA) is responsible for the antigen presentation to T cells and, therefore, a key component for adaptive immune response initiation. The HLA genes are the most polymorphic genes in the human genome, and these polymorphisms influence the ability to present different sets of epitopes to T cells. Some HLA molecules are more efficient than others presenting certain antigens, which may lead to a better induction of immune responses. This fact has already been proven for some viral diseases like A H1N1 influenza  and HIV.

It has been previously reported an association between SARS‐CoV infection and HLA‐B*07:03, HLA‐Cw*08:01, HLA‐B*46:01, and HLA‐B*54:01. Specifically, it has been reported that the individuals who are HLA‐B*46:01 positive have a higher risk of severe infection, whereas the frequency of HLA‐DRB2*03:01 is lower among patients with COVID‐19.

Mexico is one of the top 10 countries with higher mortality, and its number of cases and deaths keeps increasing significantly. Some of the most common haplotypes reported in Mexico′s less affected states are HLA A*02‐B*35‐DRB1*08‐DBQ1*04, A*68‐B*39‐DRB1*04‐DBQ1*03:02, and A*02‐B*15‐DRB1*08‐DBQ1*04, according to the Allele Frequency Net Database website (www.allelefrequencies.netwww.allelefrequencies.net).

On the other hand, up to now, Mexico City is the region with the highest number of reported cases. The studies regarding allele frequency in this city have reported that its haplotype is largely composed of Native American haplotypes, specifically 63.85 ± 1.55% American, 28.53 ± 3.13% European, and a less apparent 7.61 ± 1.96% African. Individually, some studies have reported that the most frequent alleles in Mexican population are HLA‐A*02, ‐A*24, ‐A*68, ‐B*35, ‐B*39, ‐B*51, ‐DRB1*04, ‐DRB1*08, ‐DRB1*07, ‐DQB1*0302, ‐DQB1*0301, and ‐DQB1*0201. Nonetheless, there are no studies related to the HLA association with the susceptibility or the resistance against COVID‐19 in the Mexican population. The understanding of the relationship between viral infection, HLA, and disease susceptibility is important to drive towards vaccine development and molecular epidemiology research that can contribute to novel therapies.

So far, the control of the COVID‐19 pandemic remains a challenge, resulting in thousands of new cases and deaths reported daily. It is necessary to find prophylaxis and specific treatments to contain this uncontrolled infection and to reduce the global morbidity and mortality. The generation of a vaccine that targets this virus remains as the primary solution, however, the lack of knowledge regarding the immune response, such as the HLA‐virus interactions, makes it a challenging task.

In addition, the genetic variations among different populations and their possible link with SARS‐CoV‐2 viral responses remain unknown. In this report, we analyze which epitopes of the SARS‐CoV‐2 spike protein are highly immunogenic and able to be presented by HLA Class I and II in different populations using bioinformatic tools. We also demonstrate an ecological correlation between HLA allele frequency and the predicted fatality rate in hospitalized patients of 28 Mexican states.

METHODS

Study design

A bioinformatic epitope prediction of the S glycoprotein was performed. This gave information about the most immunogenic peptide‐HLA matches and the HLA alleles that are more likely to present these epitopes efficiently. Also, an ecological study was made to look for correlations between the HLA allele frequencies and the predicted fatality rate of hospitalized patients with COVID‐19 to May 29th, 2020.

Bioinformatic epitope prediction

Bioinformatic analyses were performed to predict HLA Class I and II epitopes using the sequence of the SARS‐CoV‐2 spike protein. The sequence for the SARS‐CoV‐2 S glycoprotein was obtained from the GenBank with the accession number QHR63290.2 in FASTA format. This sequence was then submitted to the TepiTool server from the IEDB Analysis Resource database (http://tools.iedb.org/tepitool/). The epitope prediction was performed for the 27 most frequent HLA‐A and ‐B alleles that cover for most populations (Table S1). Once the total epitope list was obtained, it was submitted to the T cell Class I pMHC immunogenicity predictor server (http://tools.iedb.org/immunogenicity/) to get the immunogenicity score, which is predicted according to the amino acid residues of the peptide.

The peptide‐HLA pairs with a positive immunogenicity score and a predicted IC50 level lower than the established cut‐off (Table S2) from the complete list were chosen, considering that the lower the IC50 value, the higher the binding affinity. The 10 more immunogenic peptide‐MHC combinations from this list were selected.

The epitopes for HLA Class II molecules were also predicted using the same sequence as before and submitting it to the IEDB MHC Class II epitope prediction tool (http://tools.iedb.org/mhcii/) using the IEDB recommended 2.2 algorithm and the most common HLA‐DP, DQ, and DR alleles (Table S1). The predicted epitopes with an SMM‐predicted IC50 value higher than 50 were excluded and the sequences were ordered by the percentile rank. The MHC‐II prediction tools use a core of nine amino acids to predict the best peptide binding affinity, even when the Class II molecules bind peptides with 15 amino acid length, so the 10 SMM cores with the minor percentile rank—what means the highest affinity binding—were selected.

Structural modeling

To provide a graphical representation of the location of the epitopes, we used the structural model the full‐length SARS‐CoV‐2 spike glycoprotein (ID:6VSB_1_1_1). The full‐length SARS‐CoV‐2 structural model is available at CHARMM‐GUI13 COVID19 Archive.

The three‐dimensional (3D) structure was obtained and analyzed using PyMOL® software (Schrödinger LLC. Molecular Graphics System [PyMOL] Version 1.80 LLC, New York, NY, 2015). The basic local alignment search tool online (https://blast.ncbi.nlm.nih.gov/Blast.cgi) was used to assess the position of the predicted peptides in the glycoprotein and the protein sequence was adjusted manually using the PyMOL tools.

Analysis of HLA alleles frequency and fatality rates

We selected 28 states of Mexico considering the homogeneity in epidemiological reports and registered the allele frequency of the main capital city of each state. All the states were included except for Mexico State, Baja California Sur, and Tamaulipas because no information was found.

We used the Allele Frequency Net Database (http://www.allelefrequencies.net/default.asp) and searched for populations in North America′s geographical region and used Mexico′s (132) database. The total state samples reported on the databases was of 5840.

For HLA Class I, the subgroup alleles was not reported for 26 of the states. However, Mexico City Mestizo and Veracruz Xalapa did contain subgroup data.

In the selection of the Class II molecules, the HLA‐DPA1*03:01, DPB1*04:02, HLA‐DPA1*01:03, DPB1*02:01, HLA‐DPA1*02:01, DPB1*01:01, and HLA‐DQA1*05:0 alleles were not found in the database of any population. All the frequency data is summarized in Table S3 organized per city.

Fatality rate

We used national public data reporting all individuals with a result for SARS‐CoV‐2 in Mexico to July 8th, 2020 (SARS‐CoV‐2 Mexico database). This database is compiled by the Ministry of Health (available at https://www.gob.mx/salud/documentos/datos-abiertos-152127). We considered the following information: age, sex, state of birth, date of birth (if applicable), and type of healthcare facility where the patients were assisted—IMSS, ISSSTE, SSA, private hospital, and others. There is also information about comorbidities—diabetes, hypertension, obesity, asthma, immunosuppression, chronic kidney disease, and cardiovascular disease, as well as smoking status and hospitalization status. The registration options were yes, no, not known, or not specified. Finally, it is specified whether the patients were attended in sentinel units. The primary care sentinel institutions test for SARS‐CoV‐2 to one of every 10 patients with an acute respiratory infection, while the nonsentinel institutions perform tests according to physician criteria. The 100% of patients with severe acute respiratory infection who require hospitalization are tested in both institutions, sentinel, and nonsentinel.

The total database contained 684,804 records. We only included records of Phase 3 (614,370). We excluded 60,520 patients who were admitted for hospitalization after July 1st to allow the presentation of the outcome “death,” since the median from hospitalization to death was 7 days. Of the 553,850 remaining records, 307,421 had a negative or pending result, 173,724 were not hospitalized, 123 did not have information of the state of birth, and 1026 were indigenous people. We excluded 448 pregnant women because the immune response is expected to be different. Finally, nine records were eliminated because the date of death was before the admission date. Hence, our final sample was 71,099 records.

Statistical analysis

To create a predictive model of the hospitalized fatality rate—number of deaths caused by COVID‐19, we performed a stepwise approach with all the variables reported in the SARS‐CoV‐2 Mexico database in a Poisson model. All the variables that were significantly associated with death were kept in the model: age, sex, diabetes, hypertension, obesity, chronic kidney disease, type of healthcare, being a sentinel unit or not, and admission date. We explored if a multilevel model, using state of birth as a second level, would be a better fit for the data, but the LR test was not significant (p = 1). Hence, the state of birth variable was included in the Poisson model. Afterward, the predictive risk of death in each state was calculated. Then, a factorial analysis was performed with the 21 HLA types to determine groups that explained the variance between them and selected the representative HLA allele of each factor as the one with the maximum correlation within the factor. We selected seven factors that explained 85.2% of the variance and selected the HLA with the highest correlation within each factor as follows: Factor 1 HLA‐A*68:01, Factor 2 HLA‐A*11:01, Factor 3 HLA‐DRB1*07:01, Factor 4 HLA‐A*01:01, Factor 5 HLA‐B*57:01, Factor 6 HLA‐DRB1*01:01, and Factor 7 HLA‐B*58:01 (Table S4). Afterward, a Spearman's rank correlation was performed between the seven HLA allele frequencies and the risk of death at state level. A p < .05 was considered statistically significant. The analyses were performed in Stata v14 and figure were created using Graphpad Prism version 6.0®.

RESULTS

To assess the best Spike protein epitope‐HLA Class I matches, its sequence was analyzed looking for epitope predictions in the most frequent HLA‐A and HLA‐B alleles. The 10 most immunogenic peptides with a higher affinity binding to its restricted HLA are shown in Table 1.

Table 1

HLA Class I epitope prediction

	Peptide/protein residues (predicted immunogenicity score)	HLA Restriction
1	GTHWFVTQR/1096‐1104	HLA‐A*31:01	HLA‐A*68:01	HLA‐A*11:01	HLA‐A*03:01
	(0.3513)
	Predicted IC50	9.5	14.5	29.6	379.3
2	RSFIEDLLF/813‐821	HLA‐B*58:01	HLA‐B*57:01	HLA‐A*32:01
	(0.2744)
	Predicted IC50	7.5	24.6	62.7
3	FIAGLIAIV/1218‐1224	HLA‐A*02:03	HLA‐A*02:06	HLA‐A*02:01	HLA‐A*68:02
	(0.272)
	Predicted IC50	3.2	6.3	8.5	13.7
4	EVFNATRFA/338‐346	HLA‐A*68:02
	(0.2182)
	Predicted IC50	12
5	QYIKWPWYI/1205‐1213	HLA‐A*23:01	HLA‐A*24:02
	(0.2162)
	Predicted IC50	4.3	6.9
6	NTQEVFAQV/775‐783	HLA‐A*68:02
	(0.1788)
	Predicted IC50	5.2
7	WTAGAAAYY/256‐264	HLA‐A*26:01	HLA‐A*68:01	HLA‐A*01:01	HLA‐A*30:02
	(0.1525)
	Predicted IC50	9.9	27.4	31.1	36.4
8	SANNCTFEY/160‐168	HLA‐B*35:01
	(0.1327)
	Predicted IC50	14.1
9	YLQPRTFLL/267‐275	HLA‐A*02:01	HLA‐A*02:03	HLA‐A*02:06	HLA‐B*08:01	HLA‐A*23:01	HLA‐A*24:02	HLA‐A*32:01
	(0.1305)
	Predicted IC50	4.1	7.8	9.1	23.9	125.3	201.3	202.7
10	VVFLHVTYV/1057‐1065	HLA‐A*02:03	HLA‐A*02:06	HLA‐A*02:01	HLA‐A*68:02
	(0.1278)
	Predicted IC50	9.3	11.9	21.2	24.5

Abbreviation: HLA, human leukocyte antigen.

Although the most immunogenic peptide from this list is GTHWFVTQR, the match with the highest affinity was between the peptide FIAGLIAIV and HLA‐A*02:03. Of note, here we analyzed the most frequent Class I A and B alleles, so this analysis reveals epitopes that can be used for vaccine development and the HLA alleles that best present epitopes of this particular protein.

The best epitopes and HLA Class II alleles were also predicted, as shown in Table 2. The prediction tool for HLA Class II uses a core of nine amino acids to predict the binding efficiency of peptides to the pocket of the molecules, even if this core is in the middle of different peptides of 15 amino acids. Interestingly, among this whole set of peptides, only seven HLA molecules resulted with a high binding affinity: HLA‐DPA1*01:03/DPB1*02:01, HLA‐DPA1*02:01/DPB1*01:01, HLA‐DPA1*03:01/DPB1*04:02, HLA‐DQA1*05:01/DQB1*03:01, HLA‐DRB1*01:01, HLA‐DRB1*07:01, and HLA‐DRB1*09:01.

Table 2

HLA Class II epitope prediction

	SMM core	Peptides	HLA‐restriction	Percentile rank
a	FELLHAPAT	LSFELLHAPATVCGP	HLA‐DRB1*01:01	0.03
		VLSFELLHAPATVCG	HLA‐DRB1*01:01	0.03
		VVLSFELLHAPATVC	HLA‐DRB1*01:01	0.03
		SFELLHAPATVCGPK	HLA‐DRB1*01:01	0.09
		VVVLSFELLHAPATV	HLA‐DRB1*01:01	0.09
		FELLHAPATVCGPKK	HLA‐DRB1*01:01	0.71
b	VVVLSFELL	QPYRVVVLSFELLHA	HLA‐DPA1*03:01/DPB1*04:02	0.24
		PYRVVVLSFELLHAP	HLA‐DPA1*03:01/DPB1*04:02	0.25
		YRVVVLSFELLHAPA	HLA‐DPA1*03:01/DPB1*04:02	0.25
		PYRVVVLSFELLHAP	HLA‐DPA1*02:01/DPB1*01:01	0.3
		QPYRVVVLSFELLHA	HLA‐DPA1*02:01/DPB1*01:01	0.3
		YQPYRVVVLSFELLH	HLA‐DPA1*02:01/DPB1*01:01	0.3
		YRVVVLSFELLHAPA	HLA‐DPA1*02:01/DPB1*01:01	0.3
		PYRVVVLSFELLHAP	HLA‐DPA1*01:03/DPB1*02:01	0.36
		QPYRVVVLSFELLHA	HLA‐DPA1*01:03/DPB1*02:01	0.36
		YQPYRVVVLSFELLH	HLA‐DPA1*01:03/DPB1*02:01	0.36
		YRVVVLSFELLHAPA	HLA‐DPA1*01:03/DPB1*02:01	0.36
		RVVVLSFELLHAPAT	HLA‐DPA1*02:01/DPB1*01:01	0.63
		VVVLSFELLHAPATV	HLA‐DPA1*02:01/DPB1*01:01	0.68
		RVVVLSFELLHAPAT	HLA‐DPA1*03:01/DPB1*04:02	0.85
		YQPYRVVVLSFELLH	HLA‐DPA1*03:01/DPB1*04:02	2.2
c	FLVLLPLVS	FVFLVLLPLVSSQCV	HLA‐DRB1*01:01	0.24
		MFVFLVLLPLVSSQC	HLA‐DRB1*01:01	0.24
		VFLVLLPLVSSQCVN	HLA‐DRB1*01:01	1.3
		FLVLLPLVSSQCVNL	HLA‐DRB1*01:01	1.8
d	VLSFELLHA	RVVVLSFELLHAPAT	HLA‐DRB1*01:01	0.24
e	GYQPYRVVV	GYQPYRVVVLSFELL	HLA‐DPA1*02:01/DPB1*01:01	0.3
f	FGAGAALQI	SGWTFGAGAALQIPF	HLA‐DRB1*09:01	0.33
		TSGWTFGAGAALQIP	HLA‐DRB1*09:01	0.34
		GWTFGAGAALQIPFA	HLA‐DRB1*09:01	0.35
		WTFGAGAALQIPFAM	HLA‐DRB1*09:01	0.67
		WTFGAGAALQIPFAM	HLA‐DQA1*05:01/DQB1*03:01	1.6
g	FVFLVLLPL	MFVFLVLLPLVSSQC	HLA‐DPA1*03:01/DPB1*04:02	0.34
		FVFLVLLPLVSSQCV	HLA‐DPA1*03:01/DPB1*04:02	0.36
		MFVFLVLLPLVSSQC	HLA‐DPA1*01:03/DPB1*02:01	5.2
h	RVVVLSFEL	GYQPYRVVVLSFELL	HLA‐DPA1*01:03/DPB1*02:01	0.36
		GYQPYRVVVLSFELL	HLA‐DPA1*03:01/DPB1*04:02	6.2
i	FTISVTTEI	AIPTNFTISVTTEIL	HLA‐DRB1*07:01	0.4
		PTNFTISVTTEILPV	HLA‐DRB1*07:01	0.51
		IPTNFTISVTTEILP	HLA‐DRB1*07:01	0.52
		TNFTISVTTEILPVS	HLA‐DRB1*07:01	0.52
		NFTISVTTEILPVSM	HLA‐DRB1*07:01	2.5
		FTISVTTEILPVSMT	HLA‐DRB1*07:01	2.6
j	TNFTISVTT	IAIPTNFTISVTTEI	HLA‐DRB1*07:01	0.47

Abbreviation: HLA, human leukocyte antigen.

To track down and illustrate the specific location of the peptides in the SARS‐CoV‐2 spike glycoprotein, the corresponding 3D model was obtained. In this model, the different predicted epitopes (Tables 1 and 2) were searched in the protein structure considering its subunits and domains (Figure 1). Notably, HLA Class I peptides WTAGAAAYY, SANNCTFEY, and YLQPRTFLL—7, 8, and 9—are located in the A domain, which is highly conserved among other coronavirus species, suggesting that these could also be epitopes for other coronaviruses. On the other hand, it was found that the Class II epitopes FELLHAPAT, VVVLSFELL, FLVLLPLVS, VLSFELLHA, and FTISVTTEI—a, b, c, d, and h—and the HLA Class I EVFNATRFA—4 are preferentially found in the B domain.

Figure 1

Localization analysis of immunogenic peptides of SARS‐CoV‐2 spike glycoprotein by three‐dimensional modeling. (A) Structure of the SARS‐CoV‐2 spike glycoprotein with S1–S2 subunits. The S1 domains consist of A, B, C, and D. The S2 subunit consists of the fusion peptides and domains HR1 and HR2. (B) The predicted epitopes for HLA Class I are shown in red (C) and the suggested peptides for HLA Class II in blue. The peptides are marked individually, listed from 1 to 10 for Class I and a–j for Class II, corresponding to the immunogenicity Tables 1 and 2. HLA, human leukocyte antigen; SARS‐CoV‐2, severe acute respiratory syndrome‐coronavirus‐2

HLA allele analysis and correlation with a predicted fatality rate in hospitalized patients

After factorial analysis, we found a significant negative correlation between the frequency of the HLA‐DRB1*01:01 allele and the predicted fatality rate in hospitalized patients (R = −0.44, p = .02; Figure 2). No other significant correlations were observed (Table 3).

Figure 2

Spearman's correlation of HLA‐DRB1*01:01 frequency and fatality rate. The correlation is shown as a dot plot graph with the regression tendency line. The frequency of this HLA allele in Mexico was obtained and a correlation was performed with the predicted risk of death associated with SARS‐CoV‐2 infection in hospitalized patients. According to the bioinformatic prediction, the HLA‐DRB1*01:01 molecule can efficiently present eleven of the S protein predicted epitopes (LSFELLHAPATVCGP, VLSFELLHAPATVCG, VVLSFELLHAPATVC, SFELLHAPATVCGPK, VVVLSFELLHAPATV, FELLHAPATVCGPKK, FVFLVLLPLVSSQCV, MFVFLVLLPLVSSQC, VFLVLLPLVSSQCVN, FLVLLPLVSSQCVNL, and RVVVLSFELLHAPAT; Table 2). HLA, human leukocyte antigen; SARS‐CoV‐2, severe acute respiratory syndrome‐coronavirus‐2

Table 3

Correlation between the representative HLA alleles (7) resulted from factorial analysis and fatality rate in Mexico states (n = 26)

	R	p
F1: HLA‐A*68	.15	.45
F2: HLA‐A*11:01	−.3	.12
F3: HLA‐DRB1*07:01	.11	.6
F4: HLA‐A*02:01	.05	.79
F5: HLA‐B*57:01	.35	.07
F6: HLA‐DRB1*01:01	−.44	.02
F7: HLA‐B*58:01	−.14	.49

Note: Spearman's rank correlation

DISCUSSION

Determining HLA interactions with epitopes for optimal presentation is crucial for understanding the immunological response to SARS‐CoV‐2. Here, we present a group of epitopes of the spike protein that can be efficiently presented to CD8 and CD4 T cells and are probably related to the virus′s immune‐mediated elimination. These peptides can be either used for the peptide‐based design of vaccines or in further analysis of the immunogenicity and structure of this relevant protein.

COVID‐19 vaccine development includes five clinical‐Phase I vaccine candidates, 11 preclinical‐vaccine candidates, and 26 research‐stage vaccine candidates. , Recently, the full proteome of the SARS‐CoV‐2 has been characterized through in silico analysis to show the prediction of the most immunogenic epitopes from each viral protein for 438 MHC alleles—either Class I or Class II. , This knowledge has been considered in the design of two of the Phase I‐vaccine candidates, which are LV‐SMENP‐DC and pathogen‐specific aAPC. Nevertheless, most of the other vaccine candidates have been designed based on the spike protein of the SARS‐CoV‐2 due to better immunogenic and protective potential. The S protein is the main target for COVID‐19 vaccine development. Even though the S gene sequences of SARS‐CoV‐2 have a 93.2% nucleotide sequence identity to the bat coronavirus RaTG13 and less than a 75% nucleotide sequence identity with the SARS‐CoV, three out of the five Phase I‐vaccine candidates—which are mRNA‐1273, Ad5‐nCoV, and INO‐4800s, have been designed using this protein as the main target. ,

Remarkably, our structural analysis of the protein shows a higher abundance of epitopes in the A and B domains of the S1 subunit of the virus, indicating that, in the case of this part of the protein being processed by the host cells, it could represent a highly immunogenic region. In this analysis, we did not look for B cell epitopes in the structure of the protein. We cannot confirm that the specific target of the presented epitopes could interfere with its viral function, as would be the case of neutralizing antibodies.

HLA peptide groove sequence determines which epitopes from an antigen are presented to the immune system to elicit an effective response. The high rate of polymorphisms in the HLA locus can indicate a different ability to respond to certain antigens by different individuals. Furthermore, some HLA alleles can be more efficient in presenting certain antigens, thus also in protecting from certain infections. Our analysis from the most representative HLA alleles revealed those that present more effectively the spike protein antigens of SARS‐CoV‐2, hence, one can hypothesize that their presence in an individual might confer an enhanced ability to defend against the virus.

To assess this, we analyzed the frequency of these alleles and their relation to the disease dynamics in different states of Mexico. Although it would be interesting to extrapolate these results to several countries with different epidemiological behaviors of the disease, epidemiological reports would be highly heterogeneous and data at an individual level associated with risk of death would be needed to adjust the fatality rate.

While there is a myriad of factors related to the lethality of the disease, little is known about the involvement of the immune system in this regard. It has been proposed that many patients develop an exaggerated immune response against the infection, accompanied by a cytokine releasing syndrome or autoinflammatory syndromes. Also, Grifoni et al. showed that T helper cell responses (initiated by HLA Class II molecules) seem to be protective against the infection through a strong correlation with the production of virus‐specific antibodies, and also that they are highly represented by S‐protein specific clones.

A significant negative correlation was found between the frequency of the Class II HLA‐DRB1*01 allele and the fatality rate in hospitalized patients from the states that were included. Remarkably, this correlation was weak, suggesting that other important factors apart from HLA could be involved in the protection. Therefore, it is plausible that the correlation we found based on bioinformatic predictions, would mean that these alleles could show some degree of protection against lethal outcomes of the disease. Although, the frequency of this specific allele is low in the different states, so the overall effect in fatality rates might be small. Thus, further experimental studies are needed to reinforce these outcomes.

HLA‐DRB1*01 alleles have been previously associated with multiple sclerosis resistance. Nevertheless, its role in the susceptibility to viral diseases remains poorly understood. A recent report demonstrated, using molecular docking, that this molecule can interact with the VYQLRARSV epitope from the ORF‐7a protein of the SARS‐CoV‐2 virus. Our results revealed an ecological negative correlation of this allele and that it can present a set of epitopes. Previous reports have identified that this allele can present at least nine epitopes of the M protein and 11 of the N protein (Table S4), revealing that this molecule can be highly relevant for SARS‐CoV‐2 immunity.

A remarkable characteristic of this study is that we narrowed it to the S protein, which has been the most used target for vaccine development. Considering that we did not include other viral proteins, we made an exhaustive bibliographic review that allowed us to compile a total of 77T cell epitopes for the M protein and 87 for the N protein that were already evaluated experimentally and included an analysis of the HLA alleles used for its prediction. As shown in Table S4, the HLA‐A*26:01, HLA‐A*03:01, HLA‐A*11:01, HLA‐A*31:01, HLA‐A*32:01, HLA‐A*68:01, HLA‐B*57:01, HLA‐B*58:01, HLA‐A*01:01, HLA‐A*02:01, HLA‐A*02:03, HLA‐A*02:06, HLA‐A*68:02, HLA‐A*23:01, HLA‐A*24:02, HLA‐B*35:01, HLA‐A*30:02, and HLA‐DRB1*01:01 alleles—which resulted in our epitope prediction—can also be effective presenting peptides of other proteins like M and N.

Other studies have reported an association between HLA I alleles and several SARS‐CoV outcomes within specific populations: HLA‐B*07:03 with infection rate in China ; or HLA‐B*46:01  , with severity and HLA‐Cw*08 with infection in Taiwan. Besides, HLA‐DR*03*01 has been associated with a lower frequency of SARS‐CoV infection.

Several limitations need to be acknowledged. First, the association of the frequency of the HLA allele and fatality rate is ecological and cannot be applied at an individual level. Other studies need to be conducted to explore if the association persists at an individual level in hospitalized patients. Second, the predictive model of the fatality case was conducted using only data from hospitalized patients. Given that different comorbidities can lead to hospitalization, we cannot exclude the possibility of collider bias. That is, the conditioning of analysis on hospitalization can produce biased associations between the risk factors and the outcome “fatality rate” in this case. Third, we do not rule out the possibility of misclassification since the information on comorbidities is self‐reported. However, our aim was not to make an inference of the fatality rate at an individual‐level factor, but rather to create a predictive model that was as less biased as possible. Fourth, there may be other state characteristics that are associated with death, such as the health infrastructure or the number of available specialized medical staff that are not considered in the model. Finally, the HLA allele frequencies do not include minorities like the indigenous population, who might have different HLA alleles frequencies.

CONFLICT OF INTERESTS

All authors declare not to have any conflict of interests.

AUTHOR CONTRIBUTIONS

José P. Romero‐López conceived the idea and directed the project. Martha Carnalla‐Cortés and Nancy López‐Olmedo performed all the epidemiological analysis and correlations. Diana L. Pacheco‐Olvera performed the molecular modeling and epitope localization. Julia Moreno‐Manjón, Jacqueline Oliva‐Ramírez, and Brian Bernal‐Alferes participated in manuscript writing and revision. Ethel García‐Latorre and María L. Domínguez‐López revised the manuscript. Arturo Reyes‐Sandoval and LuisJiménez‐Zamudio participated in the discussion of the results. The authors did not receive any funding for this study.

Supporting information.

Click here for additional data file.

Supporting information.

Click here for additional data file.

Supporting information.

Click here for additional data file.

Supporting information.

Click here for additional data file.