This retrospective cohort study was conducted with the approval of the institutional Ethics Committee, Hamadan University of Medical Sciences (IR.UMSHA.REC.1399.005). Blood samples were collected from 142 COVID-19 patients who were admitted to Sina University Hospital of Hamadan in the northwestern part of Iran between July and August 2020. As controls, 143 ethnically matched healthy volunteers were recruited among blood donors during the same period and from the same geographic area. Among them, 135 were negative for IgG antibody (Quanti SARS-CoV-2 Anti-spike IgG antibody, Pishtazteb Co. Tehran, Iran) and without any symptoms related to COVID-19, and the 8 others not being tested did not have any symptoms since the beginning of the pandemic.

Diagnosis and confirmation of SARS-CoV-2 infection were carried out based on the presence of viral RNA in the nasopharyngeal swab samples (laboratory confirmed disease) and/or observation of the radiological changes in CT scan as well as known clinical presentations of COVID-19 for the suspected cases. In terms of disease severity, the patients were classified into three subgroups, namely, 46 hospitalized cases with the moderate form of the disease but not requiring admission at the intensive care unit (ICU) and not requiring supplemental oxygen (“mild/moderate” subgroup), 52 patients with severe COVID-19 who require supplemental oxygen at the ICU (“severe” subgroup), and 44 patients with critical COVID-19 at the ICU that required invasive mechanical ventilation (“critical” subgroup). Classification of the patients was based on the local guidelines and radiological findings as follows: The radiologist evaluated all five lobes of both lungs for the presence of inflammatory abnormalities including ground-glass opacities, mixed ground-glass opacities, and consolidation, according to the method presented by Li et al. (26) and Li et al. (27). In terms of the percentage of involvement, a score of 0.0 to 4.0 was considered for each lobe: 0 (0%), 1 (1%–25%), 2 (26%–50%), 3 (51%–75%), or 4 (76%–100%). Then, the total severity score (TSS) was calculated by summing the points of the five lobes, which ranges from 0 to 20. A TSS equal to or less than 3 was considered as mild involvement; 4 to 7, moderate; and ≥ 8, severe (28). Moreover, in the presence of imaging criteria of acute respiratory distress syndrome (ARDS) including ground-glass attenuation associated with traction bronchiolectasis or bronchiectasis, airspace consolidation associated with traction bronchiolectasis or bronchiectasis, crazy-paving pattern, and honeycombing, or in the presence of complications such as pneumothorax, a TSS equal to or more than 8 was considered as critical disease (27). Due to their similarities, the severe and critical subgroups were also combined as Severe/Critical for some analyses.

Age and sex information were recorded for both controls and patients. The main clinical characteristics and information of comorbidities (e.g., diabetes, renal disease, liver disease, hypertension, cardiovascular disease, malignancies, and other infectious diseases) were documented using medical records for all patients ( Supplementary Table 1 ).

Primarily, genomic DNA was extracted from EDTA containing peripheral blood samples by implementing an improved salting-out method. In the next step, genotypes of HLA-A, -B, -C, and -DRB1 loci for all COVID-19 patients and HLA-A, -B, and -DRB1 loci for healthy controls were determined by polymerase chain reaction with a sequence-specific primer (PCR-SSP) method using low-resolution HLA-A-B-C and HLA-A-B-DR SSP kits (Olerup SSP^®A-B-C and Olerup SSP^®A-B-DR SSP Combi Trays, Stockholm, Sweden) according to the manufacturer’s protocols. Unfortunately, HLA-C locus was not typed for controls. Specific HLA-A, -B, -C, and -DRB1 allele families (1st-field, which will be referred to as “alleles” for reason of simplicity) were determined by SCORE software, v5.00.80.02 T/07 provided by the company (29). The HLA-DRB1 data of patients have already been reported recently (12) in comparison with a different set of controls.

To study the difference in binding repertoires of viral peptides between HLA molecules, we extracted the whole proteome of a B.4 SARS-CoV-2 variant (30) submitted to GenBank (MT994849.1). This variant was identified to be the major viral cluster during March–July 2020 in the area (31). A total of 9,660 different 9-mer and 9,594 different 15-mer peptides were obtained, the two lengths representing the most common binders of HLA Class I and Class II molecules, respectively. According to a list of 2nd-field HLA alleles we prepared, the binding affinity of each corresponding Class I molecule (HLA-A, -B, and -C) to each 9-mer peptide and that of each Class II molecule (HLA-DR) to each 15-mer peptide were predicted by applying the state-of-the-art prediction tools netMHCpan 4.1 (32) and netMHCIIpan 4.0 (32, 33), respectively.

The list of our 2nd-field HLA alleles includes all those considered as “commonly” distributed in worldwide population (observed in at least five populations) that we defined in a previous study (34). In order to adapt the prediction results for 2nd-field alleles to our 1st-field data of patients and controls, we assigned each observed 1st-field allele to the most probable 2nd-field one, using as a reference a set of 2nd-field HLA-A, -B, and -DRB1 allele frequency data recently reported for 90 Iranians from Yazd province in the center of Iran (35), who share similar ethno-linguistic background with Hamadan people ( Supplementary Table 2 ). To avoid sampling bias, we also performed a random assignment procedure, during which each observed 1st-field allele was assigned to a random common 2nd-field allele sharing the same 1st-field number.

Data of two indices, i.e., IC₅₀ and %Rank, measuring binding affinity were retrieved from the raw output of netMHCpan and netMHCIIpan, respectively. Both indices may be used to determine if a specific peptide can be considered as a binder to an HLA molecule. Based on the thresholds suggested for IC₅₀ and %Rank to define weak (IC₅₀: 500 nm; %Rank: 2 for Class I molecules, 10 for Class II molecules) and strong binders (IC₅₀: 50 nm; %Rank: 0.5 for Class I molecules, 2 for Class II molecules), we computed for each HLA molecule the numbers of predicted weak and strong bound peptides derived from SARS-CoV-2 proteome. We kept the results from both indices and both thresholds to control possible bias introduced by the choice of these thresholds, which may actually vary among HLA molecules (36, 37). In reality, IC₅₀ has been used in most of the aforementioned prediction studies focusing on HLA–COVID-19 associations, whereas %Rank seems to be more realistic according to recent arguments (38, 39). Number of distinct weak and strong binders predicted for each HLA molecule was then computed according to each of the two indices. To measure the overall HLA capacity to bind SARS-CoV-2-derived peptides, we further computed, for each patient and control, the cumulative numbers of weak and strong binders, namely, the overall viral peptide repertoire sizes of her/his HLA molecules encoded by (1) both HLA-A alleles (n_A), (2) both HLA-B alleles (n_B), (3) all four HLA-A and -B alleles (n_AB), and (4) both HLA-DR (n_DR) alleles. For each patient, since the HLA-C genotype was available, those numbers for (5) both HLA-C alleles (n_C) and (6) all six HLA-A, -B, and -C (n_ABC) alleles were equally included.

By the use of the GENE[RATE] tools available on HLA-net (40), we first tested Hardy–Weinberg equilibrium in the control group for each locus, then estimated allele frequencies as well as two-locus haplotype frequencies. We further computed, for each locus, two slightly different measures of genetic diversity (41), i.e., heterozygosity index (h) and frequency of homozygotes.

A general comparison of HLA allele distribution between the patient and control groups was performed by computing pairwise F statistics (F_ST) using Arlequin (42) v3.5.2.2, the significance of which is accessed by a procedure of 10,000 permutations. Considered as two samples, patients and controls were also compared to the Yazd Iranian sample (43). Comparisons of the frequencies of specific HLA genotypes, alleles, and haplotypes between patient and control groups and among patient subgroups were performed using Fisher’s exact test (44), and their corresponding odds ratio (OR) was estimated with 95% confidence interval. p-value less than 0.05 was considered as statistically significant, and the Benjamini–Hochberg method for multiple comparisons was used to control the false discovery rate (45).

Distributions of the numbers of predicted bound peptides (n_A, n_B, n_AB, n_DR, n_C, and n_ABC) were summarized by using kernel density estimation (46), and their difference was compared by using two-tailed Wilcoxon test (47). The relationship between severity and the overall repertoires of HLA molecules was further studied through generalized linear models (48).

To study the influence of age, the widely observed risk factor of COVID-19, comparisons between patients and controls and between patient subgroups were also performed for individuals under 60 years old.

All analyses were performed with R (49) v4.1.2 implemented in RStudio (50) unless otherwise specified. Data visualization was accomplished by using the ggplot2 package (51) v3.3.5.

Significant difference was found for age (p < 0.0001) between patients and controls, which was not the case for sex (p > 0.05), though more male patients were found among patients (77/142 vs. 68/143). Since the healthy controls had no comorbidities, comparison of the main clinical features and demographic factors was only performed between the three patient subgroups ( Supplementary Table 1 ). Significant differences were also observed for age (p < 0.001) and most clinical indices and comorbidities, but not for sex, either. It is interesting to note that there is a much higher proportion of patients with olfactory dysfunction (OD) as well as a decreased proportion of patients with negative PCR results in the critical subgroup, the former being in contrast with previous reports that show that OD appeared to be more prevalent in patients with mild-to-moderate symptoms (52).

At the within-population level, no deviation from Hardy–Weinberg equilibrium was observed in the control group for any locus we analyzed ( Supplementary Table 3 ). At a first glance, the patient group had lower heterozygosity (h) for all three HLA loci ( Table 1 ), and the differences became more visible when displaying the distribution of HLA frequencies in the two groups ( Figure 1 ). With the observed alleles ranked by frequency for each locus, patients had more alleles with high frequencies and less alleles with intermediate frequencies than controls, implying an excess of homozygosity. Accordingly, a higher proportion of homozygotes was found among patients than controls for HLA-A and -B, which was not the case, however, for HLA-DRB1 ( Table 1 ), and none of these differences in numbers was significant according to Fisher’s exact test. Interestingly, when looking at genotypes of specific alleles, we noticed that the homozygotes of some specific alleles were unevenly distributed in patients and controls, with an obvious concentration of A*03, A*24, and B*35 homozygotes in the patient group ( Supplementary Figure 1 ).

Among patient subgroups, heterozygosity did not show any consistent difference. In contrast, proportions of HLA-A, -B, -C, and -DRB1 homozygotes all increased with severity, especially when severe and critical subgroups were combined ( Table 1 ). Homozygotes of specific alleles were not compared among subgroups because of low counting numbers.

At the among-population level, when we looked at the pairwise F_ST between patients, controls, and Yazd Iranians ( Supplementary Table 4 ), the highest values were found between patients and Yazd Iranians for the three loci, and all with significant p-values. The controls were genetically closer to Yazd Iranians, much clearer for HLA-A, where the F_ST was not significant, but less obvious for HLA-DRB1.

Figure 2 visualizes odds ratios with 95% confidence interval estimated for HLA-A, -B, and -DRB1 alleles and results of Fisher’s exact test between patients and controls (detailed results are listed in Supplementary Table 5 ). Among HLA alleles, A*03 (OR = 2.06, p = 0.0025), B*35 (OR = 1.49, p = 0.0494), and DRB1*16 (OR = 3.13, p = 0.0237) were significantly more frequent in the patient group (susceptible), whereas A*32 (OR = 0.38, p = 0.0388), B*55 (OR = 0.24, p = 0.0033), B*58 (OR = 0.12, p = 0.0376), and DRB1*14 (OR = 0.42, p = 0.0300) were significantly more prevalent in the control group (protective). After multiple testing corrections, only the result for A*03 (p _C = 0.0403) remained significant. Concerning severity, A*03, A*32, B*27, B*39, B*55, and DRB1*16 showed differences among the control group and patient subgroups, the frequencies of which did not, however, seem to be associated with severity ( Supplementary Figure 2 ).

In the same way, Figure 2 also shows the odds ratios with 95% confidence interval for HLA-A~B, B~DRB1, and A~DRB1 haplotypes more frequent than 1% in at least one group (detailed results are listed in Supplementary Table 6 , and a full list of estimated two-locus haplotype frequencies is available in Supplementary Table 7 ). The majority of the haplotypes with significant results were composed of one or two of the aforementioned alleles with different frequencies in the two groups, such as A*02~B*35, A*03~B*35, and B*35~DRB1*11. Only the result for B*35~DRB1*11 (p _C = 0.0010) remained significant after multiple testing corrections. Concerning patient subgroups, no haplotype showed an interesting association between frequency and severity (results not shown).

When HLA alleles were ranked according to the numbers of predicted bound peptides derived from SARS-CoV-2 proteome, none of the corresponding alleles detected by frequency comparisons showed any particularity, without extremely low or high values, respectively ( Supplementary Figure 3 ; values for each HLA molecule are available in Supplementary Table 8 ). For example, the only allele with significant result after corrections, namely, HLA-A*03, was medium ranked in IC₅₀-based lists and low ranked in %Rank-based lists. Actually, the ranking lists depended strongly on indices and thresholds chosen and differed considerably between the results from %Rank and IC₅₀.

In contrast, by comparing cumulative numbers of different SARS-CoV-2-derived bound peptides at the individual level, patients and controls did show remarkable differences, which were more consistent among indices and thresholds chosen. Figure 3 visualizes density distributions of these numbers for each of the two groups, taking %Rank index and threshold for weak binders as example, and density charts with different indices and/or thresholds are found in Supplementary Figures 4 – 6 . Indeed, a higher proportion of patients carried HLA-A molecules predicted to bind only 500 to 750 different viral peptides, whereas more controls carried HLA-A molecules predicted to bind more than 750 different viral peptides. When considered as a risk factor, n_A less than 750 had an odds ratio of 2.04 (95% CI: 1.20 to 3.50; Fisher’s exact test p = 0.0084). In contrast, the HLA-B molecules showed less difference between the two groups, and the controls only showed a very slightly higher proportion of values around 700. Considering HLA-A and HLA-B together (n_AB), the patients also show less bound peptides (peak at approximately 1,400) compared to controls (peak at approximately 1,600). As for HLA-DR molecules, the distribution for patients seemed more concentrated compared to controls. Wilcoxon test confirmed significant difference for n_A (p = 0.0002) but not for n_B, n_AB, and n_DR ( Table 2 ). More interestingly, when only younger (under 60 years old) individuals were included in the comparisons did numbers of predicted bound peptides differ more apparently between patients and controls for n_A and n_AB ( Figure 3 and Supplementary Table 4 ). In this case, the odds ratio for n_A less than 750 increased to 3.01 (95% CI: 1.54 to 5.93; Fisher’s exact test p = 0.0014).

Among the control group and patient subgroups, a smaller HLA-A overall repertoire also seemed to differ, especially when severe and critical subgroups were combined ( Figure 4 for %Rank-based weak binders). Generalized linear models revealed a significant association between n_A and severity (p = 0.0337, Table 3 ). This is also compatible with IC₅₀-based weak and strong binders (p = 0.264 and p = 0.311, respectively), and with %Rank-based strong binders for which the result was marginally significant (p = 0.053; see Supplementary Table 10 ). Visible but not always significant differences were also observed for n_AB, n_C, and n_ABC, and for individuals under 60 years old ( Table 3 and Supplementary Table 10 ).

Despite minor differences, the repetition of these analyses after the random assignment procedure we designed did not change the results of comparisons (results not shown).