For this study, protein reference sequences from HCoV-NL63, MERS-CoV, SARS-CoV-1, and SARS-CoV-2 and from five variants of concern (VOCs) from SARS-CoV-2 were acquired from the NCBI public database ( Supplementary Tables 1 – 3 ). Figure 1 illustrates the universe of coronaviruses studied here.

Acquired sequences were applied for the prediction of 9-mer peptides in NetCTL 1.2 software (https://services.healthtech.dtu.dk/service.php?NetCTL-1.2/), using peptide-MHC-I affinity, peptide cleavage and transport by the transporter associated with antigen protein complex (TAP) to determine an overall combined score (CS) (17). Peptides with a CS ≥0.75 are presumed to be immunogenic for CD8⁺ T cells and kept for the evaluation of prediction parameters as well as for in silico prediction of cytokine induction and biochemical parameters.

Peptide prediction was performed for six of the most prevalent MHC-I haplotypes in the Brazilian population, namely supertypes HLA-A*01 (HLA-A*01:01, HLA-A*01:02, HLA-A*25:01, HLA-A*30:04, HLA-A*32:01, HLA-A*36:01, HLA-A*43:01, HLA-A*80:01), HLA-A*02 (HLA-A*02:01, HLA-A*02:02, HLA-A*02:03, HLA-A*02:04, HLA-A*02:05, HLA-A*02:06, HLA-A*02:07, HLA-A*02:09, HLA-A*02:14, HLA-A*02:17, HLA-A*68:02, HLA-A*69:01), HLA-A*03 (HLA-A*03:01, HLA-A*11:01, HLA-A*31:01, HLA-A*33:01, HLA-A*33:03, HLA-A*34:02, HLA-A*66:01, HLA-A*68:01), HLA-A*24 (HLA-A*24:02, HLA-A*24:03, HLA-A*24:04, HLA-A*23:01, HLA-A*30:01, HLA-A*30:02, HLA-A*30:03), HLA-A*26 (HLA-A*26:01, HLA-A*26:02, HLA-A*26:03, HLA-A*26:04), and HLA-B*07 (HLA-B*07:02, HLA-B*07:03, HLA-B*07:04, HLA-B*07:05, HLA-B*15:08, HLA-B*35:01, HLA-B*35:02, HLA-B*35:03, HLA-B*51:01, HLA-B*51:02, HLA-B*53:01, HLA-B*54:01, HLA-B*55:01, HLA-B*55:02, HLA-B*56:01, HLA-B*56:02, HLA-B*67:01, HLA-B*78:01) (18, 19).

To assess the HLA-A*02-restricted peptide signature, a customized microarray was performed by PEPperPRINT^© (Heidelberg, Germany) as described previously (20) with modifications. Briefly, 2,555 peptides from HCoV-NL63, MERS-CoV, SARS-CoV-1, and SARS-CoV-2 were printed on a PEPperCHIP^© microarray slide and incubated with a HLA-A*02:01 dimeric protein (Ig : DimerX - BD Biosciences, California, USA) at 1µg/mL or 10µg/mL diluted in staining buffer (PBS + 10% bovine serum albumin). After 16 hours, the slide was washed and incubated with anti-murine IgG1/Cy3 secondary antibody (BD Biosciences, San Jose, CA, USA) diluted in staining buffer (1:5,000). After 45 minutes of incubation, the secondary antibody was removed. The slide was air-dried and digitized using the Affymetrix 428 Array Scanner device (Thermo Fisher, California, USA). Peptide reactivity to HLA-A*02 was analyzed using PepSlide^® Analyzer (PEPperPRINT ^©, Heidelberg, Germany). Quantification was performed using 16-bit grey-scale images and results were expressed in fluorescence intensity.

HLA-A2-restricted peptides derived from HCoV-NL63, MERS-CoV, SARS-CoV-1, and SARS-CoV-2 with Combined Score (CS) ≥0.75 were examined for IFN-γ and IL-4 induction in silico by MHC-II axis with IFN epitope (https://webs.iiitd.edu.in/raghava/ifnepitope/index.php) and IL4pred (https://webs.iiitd.edu.in/raghava/il4pred/index.php) software, applying Support Vector Machine method (21, 22). Peptides were also assessed for physicochemical parameters in IL4pred (https://webs.iiitd.edu.in/raghava/il4pred/index.php) (21).

Peptides with HLA-A2 reactivity in vitro that were shared between two or more species of coronaviruses included in this study were identified using Venn diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/).

Logo sequence analysis was constructed using WebLogo 3 (https://weblogo.threeplusone.com/) to assess conserved residues in HLA-A*02:01 reactive peptides that are either lost or gained in Alpha, Beta, Delta, Gamma, and Omicron SARS-CoV-2 variants (VOCs) when compared to the reference Wuhan sequence (23).

Correlations between the peptide potential of IFN-γ or IL-4 induction and physicochemical properties were assessed by Spearman’s rank correlation test on GraphPad Prism v.8.0 software (GraphPad Software, California, USA). Correlation networks were assembled using Cytoscape v3.9.0 software (24).

Orange3-3.34.0 software was used for the generation of tSNE and heatmap analysis of peptide potential of IFN-γ and IL-4 induction and physicochemical properties using normalized values (25).

Additional graphical and statistical analyses were created using GraphPad Prism v.8.0 software (GraphPad Software, California, USA).

Non-parametric distribution of the dataset was confirmed by the Shapiro-Wilk test, whereas differences between groups were assessed by Mann-Whitney, Wilcoxon, or Kruskal-Wallis followed by Dunn’s post-test. Statistical significance was considered at p<0.05 and is indicated by asterisks in the graphs.

To assess the ability of human coronaviruses (HCoV) and their viral peptides to being immunogenic, a prediction of peptide and human MHC binding was performed using protein sequences from HCoV-NL63, MERS-CoV, SARS-CoV-1 and six strains of SARS-CoV-2 (Wuhan reference sequence, and VOCs Alpha, Beta, Delta, Gamma, and Omicron). The results in Figure 2 show the ability of these peptides to bind with restriction to the six most prevalent MHC-I worldwide (HLA-A*01, HLA-A*02, HLA-A*03, HLA-A*24, HLA-A*26, and HLA-B*07). The results demonstrate a predominance of HLA-A*02-restricted peptides in HCoV ( Figures 2A, B upper panel), which indicates that this haplotype responds majorly to coronaviruses. However, it was also possible to observe that HLA-A*02-restricted peptides presented the lowest values of combined prediction score, with a median of 0.976 ( Figure 2B lower panel), in comparison with the results obtained for other MHC-I.

Furthermore, Wuhan original strain showed the highest binding scores when compared to other variants ( Figure 2B ). Considering these results and the predominance of HLA-A*02 in the population, peptides restricted to this haplotype were selected for further scrutiny (18, 26). These results clearly indicate that, while mutations occur along time, variants are probably losing immunodominant regions within the virus genome as compared to the original strain, suggesting loss of viral immunogenicity along time.

In order to confirm the in silico results for the prediction of loss in viral immunogenicity, a analysis of peptide signatures was carried out by the analysis of peptide-HLA-A*02 reactivity on a microarray chip containing selected SARS-CoV-2 peptide sequences. The results showed a dose-dependence binding for the majority of the peptides tested, with higher means of reactivity at 10µg/mL of HLA-A*02:Ig than at 1µg/mL ( Figure 3A upper panels), which corroborates the in silico results. Moreover, a comprehensive analysis of the HLA-A*02:Ig-reactive peptides in all four coronavirus species revealed abundant reactive peptides, with lower percentages for SARS-CoV-1 and SARS-CoV-2 in comparison to HCoV-NL63 and MERS-CoV ( Figure 3 left panel). However, despite the low percentage of reactive peptides, SARS-CoV-1 presented the highest mean of in vitro HLA-A*02:Ig reactivity, whereas SARS-CoV-2 presented the lowest values of overall reactivity ( Figure 3 right panel), indicating possible HLA-A*02 reactive peptide impairment in inducing CD8⁺ T cell responses as compared to other coronaviruses.

Considering the phylogenetic proximity of coronaviruses studied, we searched for possible overlapping reactive peptides amongst two or more species of coronaviruses. To obtain a panoramic view of overlapping HLA-A*02-reactive peptides, a Venn diagram was built using in silico tools. Colored regions specify the occurrence of shared reactive peptides between HCoV, whereas white regions mark viruses for which there was no overlapping reactive peptides. Interestingly, the results showed that a minority of reactive peptides presented full overlapping sequences for two or more coronaviruses, as shown in Figure 3E . The highest rate of overlap was observed between SARS-CoV-1 and SARS-CoV-2 (11.01%) ( Figures 3D, E ).

To determine which SARS-CoV-2 viral proteins display immunodominant peptides, epitope mapping of HLA-A*02:Ig reactive peptides was assessed for every viral protein. Figure 4A brings a representative scheme of all SARS-CoV-2 proteins and pie charts display the total number of peptides as well as the percentage of reactive peptides. Considering all assessed SARS-CoV-2 variants, a higher percentage of reactive peptides was observed for ORF7b (91.7%), ORF7a (72.2%), ORF3a (70%), ORF 8 (66.7%), N (66.7%) and ORF6 (62.5%), followed by S protein (58.4%), as shown in Figure 4A . Dark green regions highlight proteins with more than 60% of reactive peptides.

However, whereas ORF7b, ORF3a, ORF7a, ORF8, and N proteins were the top 5 main sites of HLA-A*02:Ig reactive peptides for the Wuhan, Delta and Gamma variants, the Omicron variant showed a slight pattern shift, with N protein representing main hotspot as shown in lollipop graphs ( Figure 4B ).

Panoramic assessment of HLA-A*02:Ig-reactive peptide signatures for SARS-CoV-2 VOCs (Alpha, Beta, Delta, Gamma and Omicron) revealed the occurrence of peptide gain and loss, in comparison to the Wuhan original strain. Peptide gain was defined as reactive peptides that are encoded in VOCs strains but were not identified in the Wuhan original strain. Conversely, peptide loss was identified when a specific reactive peptide was derived from the Wuhan sequence but was not found amongst the VOCs. Overall, peptide loss was greater than peptide gain for the 5 variants taken together ( Figure 5A ). The major sites in which such losses occurred were identified as ORF8 and S protein for all assessed VOCs ( Figure 5B right panel). However, peptide loss was observed mainly for S protein in Delta, Gamma, and Omicron variants ( Figure 5C ). Peptide gain was less common and occurred mostly in ORF1ab, but also in ORF3 and M protein ( Figure 5B left panel and 5C).

Next, we aimed at assessing distinct amino acid patterns in peptides that were either lost or gained in SARS-CoV-2 VOCs using logo sequence analysis. For this analysis, the size of the amino acid symbols is directly correlated to the frequency of occurrence of each residue in specific positions within the peptide sequence. Hydrophilic residues are indicated in blue, hydrophobic residues in green, neutral residues in black, and leucine (L) and valine (V) are highlighted in red. In this case, the totality of lost (TL) and gained (TG) peptides showed a similar pattern on conserved amino acid residues, but it was possible to distinguish a conspicuous loss of peptides with leucine (L) in position 2, especially for Beta and Gamma ( Figure 5D right panel). For all SARS-CoV-2 variants evaluated, gained peptides (TG) display aspartic acid (D) as a predominant residue in the peptide sequence ( Figure 5D left panel). As expected, the Omicron variant presented the highest shuffling of peptide residues in all 9 positions, both for gained and lost peptides ( Figure 5D ). These results indicate that mutations in VOCs are causing non silent amino acid sequence changes on CD8 T cell epitopes, which elucidate and reflect on the VOC immunogenicity loss.

Aiming at understanding the implications of HLA-A*02-restricted peptide presentation, we performed an in silico assessment of cytokine induction by these peptides, focusing on the induction of Th1 (IFN-γ) and Th2 (IL-4) mediators. Our results highlighted a dominant tendency of reactive peptides to induce IL-4 (pink) in comparison to IFN-γ (blue) with robust statistical significance for HLA-A*02-restricted peptides derived from all evaluated coronavirus included in the study ( Figure 6A ) (p<0.0001 = ****). The assessment of cytokine induction by HLA-A*02-restricted peptides derived from different proteins of SARS-CoV-2 was displayed in radar charts in which the axis indicates all structural and non-structural viral proteins and the extent of the spokes indicate the score of cytokine production (SVM score). The observed results confirmed the predominance of IL-4 induction by structural and non-structural proteins of Wuhan original strain as well as for all five SARS-CoV-2 variants ( Figure 6B ). The only exception was detected for ORF8 from the Alpha variant, with a strikingly higher potential for IFN-γ induction ( Figure 6B ).

Considering that parameters such as hydrophobicity, hydrophilicity, isoelectric point, and molecular weight can affect the potential of an epitope to elicit different shades of the immune response, the next step was to apply tSNE analysis for the assessment of physicochemical properties of HLA-A*02-restricted peptides derive from SARS-CoV-2 (Wuhan reference sequence and assessed VOCs, Alpha, Beta, Delta, Gamma and Omicron) regarding their potential to induce IL-4 and/or IFN-γ. In this strategy, every dot represents a peptide, and the color pattern is associated with the in silico score obtained for each assessed parameter in each of the tSNE graphs. In order to grant a simultaneous overview of all physicochemical parameters in all VOCs and for the Wuhan reference virus, heat map graphs were also assembled. This strategy allowed the identification of a cluster of hydrophobic peptides with increased potential of IFN-γ induction in comparison to other peptides ( Figures 7A, B ). These peptides also presented medium/low potential for IL-4 induction ( Figures 6A, B ) and were derived both from structural and non-structural proteins from all SARS-CoV-2 variants included in this study and the Wuhan original strain ( Figure 6A ).

Furthermore, correlation networks based on r Spearman’s test showed a positive correlation (full lines) between peptide induction of IFN-γ and hydrophobicity and a negative correlation (dotted lines) between IFN-γ induction and peptide hydrophilicity, whereas the opposite was observed for IL-4 induction ( Figure 7C ). These results may bring some light on the design of vaccine strategies for coronaviruses, considering the inclusion of viral peptides able of inducing IFN-γ, which may be lost by new VOCs and are essential for promoting a protective and long-lasting anti-SARS-CoV-2 immune response.

Omicron or B1.1.529 variant has become dominant as SARS-CoV-2 continues to spread worldwide. This could be attributed to the emergence of subvariants of concern with superior viral fitness and transmissibility as well as evasion from previous immune mechanisms induced by either disease or vaccination. Therefore, it is essential to address this by performing a detailed analysis of the HLA-A*02-restricted peptide signature of the new Omicron subvariants. For that, we carried out peptide prediction for six new Omicron subvariants (BA.1, BA.2, BA.3, BA.4, BA.5 and BQ.1) and assessed the rate of peptide gain and loss in comparison to the original Wuhan strain ( Figure 8A ). Overall, all the subvariants presented significant changes as compared to the original strain, apart from BA.3, which remained most similar to the Wuhan strain ( Figure 8A ). A more detailed viral protein mapping of gained and lost peptides revealed ORF1ab as a major site of peptide changes on BA.1, BA.2, BA.4, BA.5, and BQ.1 subvariants in comparison to the Wuhan strain ( Figure 8B ). Meanwhile, all Omicron subvariants had immunodominant peptide loss rather than gain on S protein, but especially for BA.1 (n=7), BA.4 (n=6) BA.5 (n=6), and BQ.1 (n=6) ( Figure 8B ).

Next, we assessed the percentage of peptides that presented in vitro reactivity to the HLA-A*02:Ig molecule and were derived from the Wuhan reference sequence but were lost on Omicron subvariants. The results are shown as pie charts assembled with the total number of lost peptides derived from each SARS-CoV-2 protein from Omicron subvariants in comparison to the Wuhan reference sequence. Data analysis showed that most peptides lost in the ORF1ab region were reactive in the in vitro assay. Interestingly, 20% of lost peptides were derived from the M protein ( Figure 8C ). As for the S protein, the first subvariants, BA.1 and BA.2, have lost more HLA-A*02:Ig-reactive peptides in this protein sequence, while BA.3 lost only one HLA-A*02:Ig-reactive peptide in this region ( Figure 8C ). These results indicate that the loss of immunodominant epitopes may follow a specific dynamics in which changes in S protein immunodominance are observed early on during viral divergence.