Eva Prašnikar
Andrej Perdih
Jure Borišek1*- 1Theory Department, Laboratory for Chemical Informatics, National Institute of Chemistry, Ljubljana, Slovenia
- 2Faculty of Medicine, Graduate School of Biomedicine, University of Ljubljana, Ljubljana, Slovenia
- 3Theory Department, Laboratory for Computational Biochemistry and Drug Design, National Institute of Chemistry, Ljubljana, Slovenia
- 4Faculty of Pharmacy, University of Ljubljana, Ljubljana, Slovenia
MHC class I antigen E (HLA-E), a ligand for the inhibitory NKG2A/CD94 receptor of the immune system, is responsible for evading the immune surveillance in several settings, including senescent cell accumulation and tumor persistence. The formation of this ligand-receptor interaction promotes the inhibition of the cytolytic action of immune system natural killer (NK) cells and CD8+ T-cells expressing this receptor. The final outcome of the HLA-E/NKG2A/CD94 interaction on target cells is also highly dependent on the identity of the nonameric peptide incorporated into the HLA-E ligand. To better understand the role played by a nonameric peptide in these immune complexes, we performed a series of multi-microsecond all-atom molecular dynamics simulations. We generated natural and alternative variants of the nonameric peptide bound to the HLA-E ligand alone or in the HLA-E/NKG2A/CD94 complexes. A systematic study of molecular recognition between HLA-E and peptides led to the development of new variants that differ at the strategic 6th position (P6) of the peptide and have favorable in silico properties comparable to those of natural binding peptides. Further examination of a selected subset of peptides in full complexes revealed a new variant that, according to our previously derived atomistic model, can interfere with the signal transduction via HLA-E/NKG2A/CD94 and thus prevent the target cell from evading immune clearance by NK and CD8+ T-cells. These simulations provide an atomistic picture of how a small change in amino acid sequence can lead to a profound effect on binding and molecular recognition. Furthermore, our study also provides new data on the peptide interaction motifs as well as the energetic and conformational properties of the binding interface, laying the structure-based foundation for future development of potential therapeutic peptides, peptidomimetics, or even small molecules that would bind to the HLA-E ligand and abrogate NKG2A/CD94 recognition. Such external intervention would be useful in the emerging field of targeting senescent cells in a variety of age-related diseases, as well as in novel cancer immunotherapies.
1 Introduction
Natural killer (NK) cells are part of the innate immune system and belong to a lymphocyte lineage with cytotoxic and cytokine-producing functions. They are involved in the organism’s machinery to combat various types of tumors and microbial infections. In addition, they also have regulatory functions as they can influence many other cell types, limiting or enhancing the immune response. NK cells must effectively recognize the target cells while ensuring tolerance to “self” cells (Vivier et al., 2008). This is provided by a variety of activating (e.g., KIR2DS, CD94/NKG2C, -E, -H, NKG2D) and inhibitory (e.g., KIR2DL, KIR3DL, ILT2, CD94/NKG2A, -B) receptors on their cell surface, resulting in a dynamic interplay of antagonistic signals that control the cytotoxic action (Vivier et al., 2008; Borrego et al., 2002; Prašnikar et al., 2021a; Prašnikar et al., 2021b). These receptors are specific for the MHC class I molecules and can be generally divided into three families: 1) killer cell Ig-like receptors, 2) immunoglobulin like transcripts, and 3) C-type lectin receptors. These receptors recognize the classical (class Ia) HLA class I molecules and the non-classical (class Ib) HLA-G and HLA-E molecules. Most representatives of the C-type lectin receptors exist as heterodimers of the CD94 protein covalently linked to a member of the NKG2 family (Borrego et al., 2002) (Figure 1).
FIGURE 1. (left) Overview of some of the known NK cell receptors and their corresponding ligands presented on the target cells. (right) Magnified representation of the HLA-E/β2m/peptide/NKG2A/CD94 immune complex comprised of HLA-E heavy chain with three alpha domains (α1–α3), HLA-E light chain β2m, nonameric peptide, CD94 and NKG2A proteins (magenta, salmon, and pink, grey, yellow, green, and blue, respectively).
Among them the NKG2A/CD94 receptor is one of the inhibitory receptors of the NK and CD8+ T-cells that recognizes and acts via the HLA-E ligand expressed on the target cells. The HLA-E/NKG2A axis plays a role in colorectal cancer (Eugène et al., 2020), multiple myeloma (Yang et al., 2021), and other cancers (Abd Hamid et al., 2019; Borst et al., 2020) and is rapidly gaining attention as a potentially druggable immune checkpoint. The NKG2A-blocking antibody monalizumab, which targets this axis, has already shown promising results in clinical trials and (Borst et al., 2020) the HLA-E ligand has been recognized as a potential target for the treatment of multiple myeloma (Yang et al., 2021). Overexpression of HLA-E has also been detected in senescent dermal fibroblasts and senescent cells in human skin sections from the elderly, and blocking the HLA-E/NKG2A interaction enhanced the immune response against the senescent cells in vitro (Pereira et al., 2019). Thus, the HLA-E/NKG2A axis opens up new promises in many therapeutic areas, and detailed understanding of the key events of molecular recognition is essential for the development of potential therapies.
The HLA-E ligand is expressed on the cell surface as a trimeric complex consisting of the heavy chain, light chain [beta-2 microglobulin (β2m)], and the specific nonameric peptide (further also referred as peptide) derived from the signal sequences of other HLA class I molecules (e.g., HLA-A, -B, -C, -G) (Braud et al., 1998b; Monaco et al., 2008). Remarkably, only specific sequences of the peptide allow high-affinity binding to HLA-E (Braud et al., 1998a; Miller et al., 2003; Ruibal et al., 2020) and its stable surface expression. Moreover, only a handful of them permit effective recognition by the NKG2A/CD94 receptor located on the NK and CD8+ T-cells, with the subsequent signal transduction, and resulting target cell protection from elimination (Braud et al., 1998a; Borrego et al., 1998; Valés-Gómez et al., 1999; Michaëlsson et al., 2002; Kaiser et al., 2005; Hoare et al., 2008; Lauterbach et al., 2015a). Furthermore, studies have also shown that not only nonameric peptides, but also longer peptides can bind to HLA-E (Kraemer et al., 2015; Celik et al., 2016). In our previous studies, we proposed this molecular recognition event takes place via a tangled network of hydrogen bonds formed in the HLA-E/β2m/peptide/NKG2A/CD94 immune complexes (Prašnikar et al., 2021a).
Some HLA-E-binding nonameric peptides were reported not to elicit protection from the lysis mediated by the NKG2A/CD94+ immune cells. These are Hsp60sp, R5V mutated leader peptide of HLA-B7 (Michaëlsson et al., 2002), the leader peptide of HLA-Cw7, and HLA-B27 (Braud et al., 1998a; Borrego et al., 1998; Valés-Gómez et al., 1999; Hoare et al., 2008). There are two possibilities proposed for such a result: 1) the peptide does not provide adequate stabilization of the HLA-E ligand on the surface of the target cell due to its weak binding to it, 2) HLA-E with bound peptide is not actually a ligand for the NKG2A/CD94 receptor. The former has been already shown to be responsible for the lack of protection in cytotoxicity assays for the peptides Cw7 and B27 (Brooks et al., 1999).
Of all the peptides studied so far, the HLA-G leader sequence is the most effective mediator of the inhibitory signaling to the effector immune cells (Kaiser et al., 2005; Lauterbach et al., 2015a) and among the different receptors, the inhibitory CD94/NKG2A receptor is the one with the highest affinity for the HLA-E complex (Lauterbach et al., 2015a). Meanwhile, signal sequences of the HLA-A alleles [VMAPRTL(L/V)L] and the HLA-B7 (VMAPRTVLL) seem to have one of the highest affinities for the HLA-E ligand (Braud et al., 1998a; Hoare et al., 2008).
Available structural and biochemical information on the HLA-E ligand and its bound peptides offer a good starting point for carrying out all-atom molecular simulations in order to unravel the molecular background for their diverse binding affinities (Dror et al., 2012). The acquired data could be then used as the basis for the targeted design of potential therapeutic peptides that would bind to the HLA-E expressed on the target cell surface and abrogate the NKG2A/CD94+ immune cell inhibition, over pronounced in some pathological conditions and in senescent cell accumulation. Indeed, therapeutic peptides can represent an attractive alternative to small molecules and proteins having affinities comparable to antibodies and specificities much better than small molecules or antibodies (Shrivastava et al., 2009; Lau and Dunn, 2018). In their development the utilization of computational techniques can play a key role in optimal design of such amino acid-based therapeutics (D’Annessa et al., 2020; Farhadi and Hashemian, 2018).
Here, we report the results of extensive all-atom molecular dynamics (MD) simulations of naturally occurring and alternative designed variants of the nonameric peptides in complex with HLA-E/β2m ligand alone or in the HLA-E/β2m/NKG2A/CD94 complexes to broaden the understanding of the molecular background for their diverse binding affinities which are often associated with only small changes in their structure. Simulations revealed new peptide variants with comparable affinity as natural peptides with highest HLA-E binding affinities. One investigated variant can according to our derived atomistic model (Prašnikar et al., 2021a), interfere with the signal transduction taking place in HLA-E/β2m/peptide/NKG2A/CD94 immune complex thus preventing the target cell from evading the immune clearance via NK and CD8+ T-cells. With this study, we aim to demonstrate and provide an atomistic picture of how small changes in amino acid sequence can significantly affect such molecular recognition outcomes.
2 Materials and Methods
2.1 Structural Models
Altogether we built 12 models consisting of HLA-E heavy and HLA-E light chain [called beta-2-microglobulin (β2m)] with various nonameric peptides incorporated into the HLA-E complex along with the apo form of HLA-E. Next, we constructed four models containing the NKG2A/CD94 heterodimer in addition to the HLA-E heavy and light chains and the corresponding nonameric peptide accounting for a total of 17 simulated models. All HLA-E complexes are built using either the crystal structure of the human non-classical major histocompatibility complex molecule HLA-E with bound class Ia MHC leader peptide derived from HLA-B8, solved at 2.85 Å resolution (PDB entry 1MHE) or the crystal structure of human HLA-E, solved at the resolution of 2.50 Å that contains the leader peptide of the HLA class I histocompatibility antigen, alpha chain G (PDB entry 3BZE).
The first model (i) HLA-EG_wt consists of HLA-E, β2m, and the G peptide (sequence: VMAPRTLFL), while the remaining models differ in nonameric peptide, incorporated into the HLA-E. The second model (ii) HLA-EB27_wt contains a peptide derived from HLA-B27 (sequence: VTAPRTLLL, peptide from PDB entry 1KTL). The rest of the models of the HLA-E complex contain a peptide corresponding to the leader sequence of HLA-B8 and HLA-B7 (sequence: VMAPRTVLL). Namely, the third model (iii) HLA-Ewt contains the wild type of the B7 peptide. Then, single mutations of the B7 peptide were introduced in the following models: R5V for model (iv) HLA-ER5V, T6M for (v) HLA-ET6M, T6A for (vi) HLA-ET6A, T6F for (vii) HLA-ET6F, T6Y for (viii) HLA-ET6Y, T6W for (ix) HLA-ET6W, T6L for (x) HLA-ET6L, and T6I for (xi) HLA-ET6I. We also constructed a model with the wild type B7 peptide with capped peptides’ N-terminal residue, (xii) HLA-Ecap. Capped residues are used to neutralize terminals in peptide chains and for neutralization of peptides’ N-terminal in HLA-Ecap model we used the ACE (-C(=O)-CH3) residue. Finally, we prepared the apo model of HLA-E/β2m, (xiii) HLA-Eapo. In the superscript of all model names the name or the variant of the subjected nonameric peptide is listed.
All simulated HLA-E/β2m/peptide/NKG2A/CD94 complexes are based on the crystal structure of CD94/NKG2A in complex with HLA-E with the incorporated leader sequence of HLA-G, solved at 3.40 Å resolution (PDB entry 3CDG). Peptide G derived from HLA-G, originally present in the crystal structure, was replaced with the peptide B7 (sequence: VMAPRTVLL, peptide from PDB entry 1KPR) to which single mutations were introduced, based on the initial computational assessment of known peptides: T6M in model (xiv) NKG2A/HLA-ET6M, T6Y in (xv) NKG2A/HLA-ET6Y, T6W in (xvi) NKG2A/HLA-ET6W, and T6I in (xvii) NKG2A/HLA-ET6I. All x-ray structures used in our simulation contain arginine residue at the position 107 thus corresponding to the HLA-E*01:01 variant.
2.2 System Preparation and Molecular Dynamics Simulations
For the assignation of protonation states of all ionizable residues PDB2PQR web tool at pH 7 was used (Dolinsky et al., 2004). Carboxylic amino acids were found in their deprotonated states, while histidines were protonated at both imidazole sites or only at Nδ or Nε. The protein complex was embedded in a 10 Å thick layer of TIP3P water molecules (Jorgensen et al., 1983). For the HLA-E complexes, this resulted in a box of 83.144 Å3 × 103.623 Å3 × 99.798 Å3. The assembled systems accommodated approximately 85,000 atoms, including water molecules and 15 Na+ counterions. The dimensions of the HLA-E/β2m/peptide/NKG2A/CD94 complexes are 130.968 Å3 × 123.440 Å3 × 122.079 Å3, including water molecules and 18 Na+ counterions, the systems numbered approximately 193,500 atoms. Tleap module of Ambertools 18 was used to construct disulfide bonds and model topologies (Case et al., 2018).
First, energy minimization was performed, followed by heating of the systems in two successive steps. In the first step, heating was performed from 0 to 100 K over 5 ps and in the second step from 100 to 303 K over the next 120 ps. During the stepwise heating, positional restraints of 200 kcal/mol Å2 and 100 kcal/mol Å2 were applied to all heavy atoms in the first and second steps, respectively. Then an isothermal-isobaric ensemble (NPT)-based MD simulation was performed over the next 10 ns, where all restraints were removed and a Berendsen barostat was used to control the pressure of 1 bar (Berendsen et al., 1984). Next, 1.2 μs of MD simulation was carried out using a canonical ensemble (NVT) with periodic boundary conditions and an integration time step of 2 fs, giving a total simulation time of ∼20 µs During the simulation, temperature control (T = 303 K) was ensured with a collision frequency of 1 ps−1 using the Langevin thermostat (Loncharich et al., 1992). Hydrogen bonds were constrained using the SHAKE algorithm (Ryckaert et al., 1977) and long-range electrostatic interactions were calculated using the particle mesh Ewald method (Harvey and De Fabritiis, 2009), with a cutoff of 10 Å.
The first ∼200 ns of the NVT MD run were considered an equilibration phase, while the remaining 1 μs of the MD trajectories, stripped of water and counterions, was used for subsequent analysis. MD trajectories were visualized using VMD (Humphrey et al., 1996) and PyMol (SchrӧdingerLLC, 2015) software and analyzed with the cpptraj module in Ambertools 18 (Case et al., 2018) and Gromacs 2016 (Van Der Spoel et al., 2005) suite. The cpptraj module in Ambertools 18 (Case et al., 2018) was used for clustering of MD conformations, using a hierarchical agglomerative approach with a distance cutoff of ∼2 Å and a distance metric of mass weighted root-mean-square deviation (RMSD) of backbone atoms (Shenkin and McDonald, 1994).
2.3 Binding Free Energy and Interaction Energy Calculations
The binding free energies (ΔGb) between the peptide and the rest of the protein complexes as well as its per-residue distribution were calculated by the Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) method (Massova and Kollman, 2000) and the Amber18 code (Case et al., 2018). The calculations were performed using the equilibrated part of the trajectories corresponding to the last 1,000 ns of the MD production run. MM-GBSA calculations were performed for 100 equally spaced frames between 600 and 900 ns of the equilibrated part of the MD trajectory. Per-residue decomposition and energy analyses were performed using a salt concentration of 0.1 M and igb flag value set to 5. In these calculations, the conformational entropic free energy contribution was omitted as it was reported not to contribute to the better quality of the results when using the MM-G(P)BSA (Genheden and Ryde, 2015; Borišek et al., 2020). The interaction energies (ΔE) between the peptide and HLA-E/β2m were additionally calculated using the gmx energy module of the Gromacs 2016 (Van Der Spoel et al., 2005).
2.4 Generation of Dynamic Pharmacophores
To gain a deeper understanding of the interaction patterns between different peptide variants and HLA-E ligand, we also derived dynamic pharmacophore models (i.e., dynophores). 1,000 equally spaced frames between 10 ns and 1 μs of the equilibrated parts of the trajectories were exported and analyzed using the DynophoreApp (Bermudez et al., 2015; Sydow, 2015; Bock et al., 2016). The resulting models were then visualized and analyzed in LigandScout (Wolber and Langer, 2005). The LigandScout capabilities were also used to generate apo-site grids of the peptides’ binding pocket on the HLA-E ligand. Buriedness of 0.70, a surface grid of 0.15, and feature probabilities of 1.00 were used for visualization and analysis.
3 Results
All models simulated in this study were firmly based on the available crystal structures of the extracellular domains of the human CD94/NKG2A proteins in complex with HLA-E (PDB entry: 3CDG) or HLA-E ligand alone (PDB entry: 3BZE, and 1MHE). Constructed models include all resolved components consisting of extracellular domains of HLA-E, β2m, and nonameric peptide nested between the α1 and α2 domains of the HLA-E, referred to as ligand part. The full immune complexes consisting of ligand and receptor parts also includes the resolved extracellular domains of the NK receptor counterpart, namely proteins NKG2A and CD94 (Figure 1). Together we studied 17 immune complexes (Table 1) of which 13 models consist only of HLA-E ligand part of the immune complex: HLA-EG_wt, HLA-EB27_wt, HLA-Ewt, HLA-ER5V, HLA-ET6M, HLA-ET6A, HLA-ET6F, HLA-ET6Y, HLA-ET6W, HLA-ET6L, HLA-ET6I, HLA-Ecap, and HLA-Eapo, while the remaining four models also encompass the receptor NKG2A/CD94 counterpart: NKG2A/HLA-ET6M, NKG2A/HLA-ET6Y, NKG2A/HLA-ET6W, and NKG2A/HLA-ET6I, where the superscript indicates the name or the variant of the subjected nonameric peptide.
TABLE 1. List of constructed models of the HLA-E and HLA-E/NKG2A complexes with specific names of the model, nonameric peptide presence in the complex, its sequence, the mutations introduced, and the crystal structure used to build the model on.
3.1 Identification of the Prospective Peptide Mutation Points by the End-Point Binding Free Energy Calculations
In the first part of the study, we evaluated the structural features and dynamic behavior of the HLA-E ligand part of the immune complex by simulating three different peptides (peptide G, B27, and B7) with experimentally reported affinities for HLA-E (Braud et al., 1998a). Per-residue decomposition of the binding free energy of the simulated models HLA-Ewt, HLA-EG_wt, and HLA-EB27_wt using the MM-GBSA method (Genheden and Ryde, 2015; Borišek et al., 2020) (Supplementary Table S1; Figure 2) revealed the binding contribution of each peptide residue and provided the rationale for the selection of the mutation points. In agreement with the literature, B7 and G peptides’ positions 2 and 9 represent the energetically most important anchor points. This has been demonstrated experimentally by mutation studies (Braud et al., 1997; Miller et al., 2003), the conservation of amino acids in naturally occurring HLA-E binding peptides (Ruibal et al., 2020), and solved crystal structures. In the latter, deep and shallow binding pockets were observed to accommodate the side chains of P2, P9, as well as P7, and P3 in addition to P6, respectively (Supplementary Table S2) (O’Callaghan et al., 1998). Moreover, in case of the peptide B27 the anchor role of P2 is weakened upon replacement of Met with Thr which is in line with previous reports of Thr not being a preferred amino acid at this position (Miller et al., 2003). Per-residue decomposition additionally revealed that positions 5 and 7 also make a substantial contribution to binding, while the contribution of amino acids at positions 1, 4, and 6 was significantly smaller (Figure 2).
FIGURE 2. Nonameric peptide residue positions, binding pockets, altering points, binding site, energetic profile and per-residues Root mean square fluctuation (RMSF) values. (A) Wild type peptide B7 with highlighted residue positions P1-P9 and alteration points with binding pockets on the HLA-E/β2m accommodating peptides’ side chains. (B) Per-residue decomposition of the mean MM-GBSA binding free energy calculated over 100 equally spaced frames between 600 and 900 ns of the equilibrated part of the MD trajectories for HLA-E/β2m models including wild type peptides B7 (blue), G (orange), and B27 (grey) and capped peptide cap (red). (C) B7 peptide in its HLA-E binding pocket, colored according to the amino acid hydrophobicity (Sweet and Eisenberg, 1983). (D) Per-residues Root mean square fluctuation (RMSF) values calculated on the equilibrated part of the Molecular Dynamics (MD) trajectories for wild type peptides B7 (light blue), G (orange), and B27 (grey) and capped peptide cap (violet) for residues’ Cα atoms. Marked whiskers denote the obtained standard deviation.
Proline residue at position 4 is conserved in naturally occurring peptides derived from the class Ia HLA signal sequences (O’Callaghan et al., 1998; Chen et al., 2014), which may well indicate an important structural function via an imposed conformational constraint, so its mutation could significantly affect binding and was thus not considered. Hence, to develop a B7 peptide variant with higher affinity for HLA-E ligand, peptide interactions at positions 1 or 6, besides P4 proline, could be optimized as these residues of the peptide have the smallest per-residue binding free energy contributions.
Interestingly, position 1 represents the N-terminal residue of the peptide and its positive charge could be the main reason for its poor contribution to binding observed in the energetic analysis (Figure 2). To test this initial hypothesis that the charge of the N-terminal residue negatively affects the interactions of P1 with the binding pocket in HLA-E/β2m ligand, we decided to cap the N-terminal amino acid with ACE (-C(=O)-CH3) capping residue, and perform another MD simulation. Comparison of the per-residue decompositions of the capped peptide of model HLA-Ecap with the wild type peptide HLA-Ewt (Supplementary Table S1; Figure 2) corroborated the positive charge of the N-terminus as a likely culprit for the antagonistic contribution of P1 to peptide binding.
Furthermore, capping indeed enhanced the interactions of P1 residue in HLA-Ecap with its surroundings to the extent that it can be considered as one of the anchoring residues. In contrast, the RMSF calculation showed that capping of the P1 slightly destabilized the binding of the peptide at the N-terminal site, particularly the residue P1 as it now became more flexible. However, compared with the bound peptides G and B27 in HLA-EG_wt and HLA-EB27_wt, the capped peptide still exhibits an overall low flexibility (Figure 2; Supplementary Table S3).
Subsequent comparison of the most representative clusters extracted from the MD trajectories (the last 1 μs of simulation time) revealed an interesting orientation of the side chain of the Arg62HLA-E residue in HLA-Ecap model. Namely, in HLA-Ewt, HLA-EG_wt, HLA-EB27_wt models containing wild-type peptides, Arg62HLA-E residue is extending toward the N-terminal site of the peptide, while in the HLA-Ecap model its side chain was pointing more upwards extending in the direction of the C-terminal site of the peptide (Figure 3). This observation is in line with a similar experimentally determined orientation that was described in the MHC-I complex with N-terminally extended, N-methylated 10-mer peptide (Li et al., 2019). Such positioning of Arg62HLA-E was proposed to open pocket A of the peptide binding site, allowing elongated peptides to protrude from the binding groove (Li et al., 2019). Thus, based on our results it could be suspected that this distinct conformation of Arg62HLA-E might be the consequence of the capping.
FIGURE 3. (A) Conformations of the Arg62HLA-E residue in most representative clusters of HLA-Ecap, HLA-EG_wt, HLA-EB27_wt, and HLA-Ewt. Compared to other models, the Arg62HLA-E residue points more towards C-terminal in the HLA-Ecap model. (B) B7 peptide in yellow depicted as spheres and HLA-E (magenta, slate and pink)/β2m (gray) ligand complex HLA-Ewt in surface representation. (C) Determined mean binding free energies (ΔGb) between peptides and HLA-E/β2m calculated by MM-GBSA and (D) interaction energies (ΔE) between peptides and HLA-E/β2m calculated with Gromacs. Marked whiskers denote the obtained standard deviation.
The MM-GBSA binding free energy as well as Gromacs interaction energy calculations successfully reproduced the experimentally determined order of binding affinities of all three wild-type peptides B7, G, and B27 for HLA-E/β2m (Supplementary Table S4; Figure 3) (Braud et al., 1998a). This provides the necessary reassurance that these calculations can be used to generally distinguish between the peptides with higher or lower affinity for the HLA-E. Furthermore, binding free energy calculations for the HLA-Ecap model revealed that capping of P1 did not improve the binding affinity of the peptide. Both, binding free energy and interaction energy between the capped peptide and the HLA-E/β2m ligand proved to be less favorable than the energies of the HLA-Ewt with bound B7 peptide and closer to that of the HLA-EG _wt model.
These results suggest that the nonameric peptide’s P6 position holds the most potential for further optimization to enhance its binding to the HLA-E ligand. Finally, we analyzed the peptide binding site on the HLA-E with the grid apo-site analysis performed in Ligandscout (Wolber and Langer, 2005) as another viable approach to search for possible favorable interaction areas in this binding pocket. This analysis, which is unbiased in the respect to the existing bound ligands, revealed positive ionizable, aromatic, and hydrophobic pharmacophore features in close proximity to the P6 site further substantiating the previous analysis (Supplementary Figure S1).
3.2 Mutation of the 6th Position of the Nonameric Peptide
Based on the results obtained, the hydrophobic nature seems to be strongest near the P6 site, where the introduction of additional charged structural elements could lead to a desolvation penalty (Salari and Chong, 2010), which could be problematic for peptide binding. This implies that the introduction of hydrophobic residues I, F, L, W, M, A, Y at this position could be considered good substitutions for more optimal binding.
Therefore, we generated these seven variants of the peptide B7 that differ in the amino acid at P6 position to evaluate the effects of this point mutation on the peptides binding to HLA-E ligand. Additionally, we also simulated another model, where we introduced the R5V mutation to the B7 peptide, which has been shown to abrogate the cytolytic effect of NKG2A/CD94+ NK cells (Michaëlsson et al., 2002) in order to assess R5V peptide binding properties. Determined MM-GBSA binding free energies and Gromacs interaction energies both showed that designed peptide variant T6I in the HLA-ET6I system has a comparable affinity to the initial wild type peptide B7 in HLA-Ewt. Both energy calculations also suggested good affinity for designed peptide T6L in HLA-ET6L (Supplementary Table S4; Figure 4). Subsequent RMSF calculations assessing the peptides flexibility revealed that of all the simulated peptides, the B7 wild type peptide was the most stably bound in the binding pocket (Figure 4; Supplementary Table S3). Further comparison of the RMSFs of the wild-type peptides shows that the peptides are most stably bound at the C-terminal site of the peptides. This can also be observed for all mutant peptides. Moreover, peptide T6I is similarly stable in its binding pocket as the peptide B7, whereas the same is not true for the peptide T6L. Less flexibility of the peptide may imply stronger binding, so these observations further corroborate the observed comparable binding affinities of peptides B7 and T6I. Although the B7 peptide shows higher flexibility at the N-terminal site after capping, the capped peptide still shows low overall flexibility compared to peptides other than the T6I-B7 variant, being in line with the trend observed in the models involving wild type peptides.
FIGURE 4. Binding free energy calculations between peptides and HLA-E/β2m and Per-residues Root mean square fluctuations (RMSF). (A) Mean binding free energies (ΔGb) calculated by MM-GBSA method. (B) mean interaction energies (ΔE) calculated with Gromacs. (C) Per-residues Root mean square fluctuation (RMSF) values calculated on the equilibrated part of Molecular Dynamics (MD) trajectories for peptides: wild type peptide G (G_wt), B27 (B27_wt), B7 (wt), B7 variants R5V, T6I, T6L, T6A, T6M, T6F, T6Y, T6W, and B7 (cap) (left panel). Sequences of subjected peptides in corresponding models colored according to RMSF values from blue to red indicating increased flexibility (right panel). Marked whiskers denote the obtained standard deviation.
Next, we created dynamic 3D pharmacophore (dynophore) models, also known as dynamic pharmacophores, for each of the 12 simulated HLA-E:peptide complexes. In them, dynamic nature of the MD simulations and static structure-based pharmacophore models are combined, resulting in a dynamical characterization of the peptides interaction patterns throughout the simulations (Janežič et al., 2021). In this way, we have added another layer of analysis to the previously described energy-based investigation by now focusing more on the geometric aspect of molecular recognition between the peptides and HLA-E ligand with included dynamical component provided by the MD simulations. According to the dynophore analysis, peptides B7 in HLA-Ewt and B7_T6I mutant peptide in HLA-ET6I system form similar interactions with HLA-E with also comparable intensities. The obtained dynophores also demonstrated that in models containing wild-type peptides with lower affinity for HLA-E than peptide B7, namely peptides G in HLA-EG_wt and B27 in HLA-EB27_wt, the pharmacophore features are more scattered than in the HLA-Ewt model, where they were observed to be more condensed (Figure 5). We also observed that isoleucine on P6 in the HLA-ET6I system forms additional hydrophobic interactions with Trp97 while its branched side chain enables extra contact points with Ile73, Thr70, and Trp97 residues. Moreover, per-residue decomposition of the binding free energy by the MM-GBSA method confirmed that the contribution of the 6I residue of B7_T6I to binding in HLA-ET6I system is indeed greater than the contribution of T6 residue of the B7 peptide (Figure 6).
FIGURE 5. Dynophore models for the wild type peptide models (A) HLA-EG_wt, (B) HLA-EB27_wt, and (C) HLA-Ewt and B7 peptide variant (D) HLA-ET6I calculated for 1,000 uniformly distributed frames between 10 ns and 1 μs of the last 1 μs of MD trajectories. Only interactions that occur in at least 50% of the trajectory are mapped. From the visualization on the right, it is evident that peptides with better binding affinities have more condensed interactions. PI, positive ionizable; HBA, hydrogen bond acceptor; HBD, hydrogen bond donor; H, hydrophobic area.
FIGURE 6. Mean per-residue decomposition of MM-GBSA Binding Free Energies (ΔGb) calculated over 100 equally spaced frames between 600 and 900 ns of the last 1 μs of the MD trajectories for HLA-E/β2m models including variants of the B7 peptides: wild type (wt) (blue), R5V (orange), T6A (grey), T6M (yellow), T6F (pink), T6Y (green), T6W (cyan), T6L (brown), T6I (violet), and capped B7 (cap) (red). Marked whiskers denote the obtained standard deviation.
We further substantiated the dynophore geometric analysis with per-residue MM-GBSA decomposition analysis for the generated peptide variants and compared it to the observations obtained for the three wild type peptides with experimentally reported data. All mutations at the peptides 6th residue, with the exception of T6A, resulted in greater contribution of P6 to the binding (Figures 2, 6; Supplementary Figures S2, S3). However, the point mutations in most of these cases did not result in the overall improved peptide:HLA-E/β2m binding (Supplementary Table S4, Figure 4).
The per-residue analysis further shows that in all simulated cases the C-terminal site of the peptide is more stably bound, as here two residues (P7 and P9) with large contribution to binding could be observed. On the other hand, only one such residue (P2) could be found at the peptide N-terminal site. The central residue (P5) also makes a large contribution to binding (Figures 2, 4, 6; Supplementary Tables S1, S3). Joining our observations from dynophore, MM-GBSA, and RMSF analyses together, we can detect three anchor residues, namely P2, P7, and P9 which is in line with the behavior observed for the wild type peptides. Considering only B7 variants, with exclusion of the peptide variant T6A, P6 could be considered as another anchor residue as binding contributions were substantially improved compared to P6 of the wild type peptides. This confirms that the selection of amino acids that were introduced into P6 can improve the interaction to this part of the peptide binding pocket. It should be noted that P5 also serves as a residue with a significant contribution to peptide binding and is not accommodated in any of the HLA-E binding pockets but protrudes from the binding furrow and is thus putatively involved in recognition by CD94/NKG2A receptor.
In agreement with previous reports (Ruibal et al., 2020), our simulations have confirmed that P6 is indeed intolerant for the incorporation of bulky amino acids, as P6 mutation to W in the HLA-ET6W system resulted in an early detachment of the N-terminal site (residues P1–P4) of the peptide from its binding pocket. Interestingly, the C-terminal site of the peptide (residues P5–P9) remained bound until the end of this simulation (Supplementary Figure S4). Additionally to the model HLA-ET6W, minor detachments of the peptides at the N-terminal side were observed also in models HLA-EG_wt, HLA-ET6L, HLA-ET6Y, HLA-ET6M, and HLA-ET6W. In them only the P1 residue was protruding from the binding pocket, while P2 moves from HLA-E pocket B towards pocket A (Supplementary Table S2). While minor differences in the peptides’ N-terminal side positioning in the binding groove were observed, no major conformational differences in the simulated models were noticed (Supplementary Figures S5–S8). We also detected that in most representative clusters of HLA-E models containing B7 variants, Arg62HLA-E was similarly aligned as in the wild-type models, further confirming the role of capping of the N-terminal amino acid for the orientation of Arg62HLA-E toward the C-terminal peptide site discussed above (Supplementary Figure S9). These results advocate that binding of the peptide might be further improved by securing better interactions at the N-terminal site, thereby preventing the observed detachments, leading to a more active contribution of the N-terminal site to the HLA-E binding.
3.3 Assessment of NKG2A/CD94 Interaction With the Ligand Variants
Finally, we evaluated whether the introduced point mutations in the nonameric peptide can lead to the abrogation of the interaction between HLA-E and inhibitory receptor NKG2A/CD94, resulting in the inhibition of the cytotoxicity of the effector immune cells. To this end, we examined the influence of several designed peptide variants on the occurrence and stability of the hydrogen bond between the Lys135NKG2A and Ser106CD94 and on the salt bridge between Lys135NKG2A and Ser109CD94 formed between receptor proteins. In our previous work, these interactions were shown to play key roles in ensuring inhibitory signal transduction (Figure 7; Supplementary Figure S10) (Prašnikar et al., 2021a).
FIGURE 7. Representation of the full immune complex and dynophores for Lys135NKG2A interactions. (A) Full immune complex involving ligand HLA-E (α1 (magenta), α2 (salmon) and α3 (pink) domains)/β2m (grey)/peptide (yellow) and receptor NKG2A (blue)/CD94 (green) with magnified key residues Lys135NKG2A, Asp106CD94, and Ser109CD94 and illustrated salt bridge (pink) and hydrogen bonds (brown) putatively involved in inhibitory signal transduction (Prašnikar et al., 2021a). (B) Dynophore for Lys135NKG2A of the immune complexes NKG2A/HLA-Ewt, (C) NKG2A/HLA-ER5V, and (D) NKG2A/HLA-ET6M for 1,000 uniformly distributed frames between 10 ns and 1 μs of the equilibrated part of the trajectories. The percentage of trajectories in which an interaction is present is given for each model. More condensed interactions correspond to the more stable interaction. SB, salt bridge; HB, hydrogen bond; BB—the interaction is formed with the oxygen backbone atom of Asp106CD94.
Together, we simulated and analyzed four HLA-E/β2m/peptide/NKG2A/CD94 immune complexes with designed variants of the peptide B7 containing the following mutations: T6M, T6Y, T6W, and T6I. The T6I variant was selected because of its previously discussed favorable HLA-E binding properties. Because the T6W variant partially detached from the HLA-E binding pocket in the HLA-ET6W model, we wanted to further test whether this behavior would lead to the disruption of the proposed signal transduction mechanism. The remaining two simulated peptide variants with mutations T6Y or T6M were selected to expand the “chemical space” studied. Namely, in addition to their general hydrophobic nature, both introduced amino acids, contain functional groups in their side chains that distinguish them from other point mutations. This could lead to different results in terms of the final receptor (NKG2A/CD94) and ligand (HLA-E/β2m/peptide) recognition and outcome.
In addition to the variants of peptide B7 mutated at position 6, we also present here some new data that we extracted from the simulated HLA-E/β2m/peptide/NKG2A/CD94 complexes with peptides G, B7, B7 with R5V mutation, and its apo model, that we generated and simulated in the scope of our previous work where we studied the molecular recognition events necessary for the inhibitory signal transduction (Prašnikar et al., 2021a). According to the literature, peptides B7 and G enable inhibition of cytotoxic activity (Valés-Gómez et al., 1999). In contrast, the R5V mutation of B7 has been reported to completely abrogate the HLA-E-mediated protection against killing by the CD94/NKG2A-expressing NK cells (Michaëlsson et al., 2002).
Measurements of the Asp106CD94-Lys135NKG2A and Ser109CD94-Lys135NKG2A interaction distances as well as their assessment via dynophore models led to the conclusion that T6Y, T6W, and T6I peptide variants are likely to allow the receptor-ligand recognition to take place, thereby preventing the cytotoxic action of the NK cell. On the other hand, the NKG2A/HLA-ET6M model with T6M mutation of the B7 peptide shows the presence of the proposed key interactions at frequencies closer to the negative NKG2A/HLA-ER5V model containing the B7_R5V peptide that abrogates HLA-E-mediated protection. Furthermore, the NKG2A/HLA-ET6M model also has similar average distances between the heavy atoms of the interacting residues as NKG2A/HLA-ER5V. In contrast, the remaining three models with petide variants display comparable features to the models that allow inhibitory signal transduction: NKG2A/HLA-EG_wt and NKG2A/HLA-Ewt (Prašnikar et al., 2021a) with bound peptides G and B7 with respect to the monitored three key interactions (Table 2; Figure 7; Supplementary Figure S10).
TABLE 2. Average distances between key atoms and percentages of trajectory in which interaction (SB, salt bridge; HB, hydrogen bond) is present.
4 Discussion
The MHC class I molecules bind short endogenous peptides during their initial folding in the endoplasmic reticulum, providing the mechanism to monitor abnormal protein production (Jensen et al., 1999). Peptides are generated with the proteasome, transported in the endoplasmic reticulum by an MHC-encoded heterodimeric member of the ATP-binding cassette family of transporters (TAP), the peptide is further cleaved by the endoplasmic reticulum aminopeptidase (ERAP), and loaded in the peptide loading complex to HLA-E (Jensen et al., 1999; Blees et al., 2017; Ogg et al., 2019). HLA-E binds almost exclusively nonameric peptides derived from the signal sequences of MHC Ib molecules, but binding of alternative peptides has also been observed after TAP or ERAP disruption (Ogg et al., 2019).
It has been reported that for strong binding to the HLA-E, it is preferred that the nonameric peptide consists of methionine on position 2 and isoleucine or leucine on position 9 (Michaëlsson et al., 2002). Additionally positions 6, 7 and 3 were also shown to have a marked impact on binding affinity and thereby HLA-E cell-surface stabilization. Together they were labeled as anchor residues (Miller et al., 2003; Hoare et al., 2008). Also, the relevance of position 5 is repeatedly reported, as Arg is preferred on this position (Michaëlsson et al., 2002; Miller et al., 2003; Petrie et al., 2008). In fact, amino acids on peptide positions 5 and 8 presumably serve as a major contact residues to communicate with the CD94/NKG2 receptor on the immune system natural killer (NK) cells and CD8+ T-cells (Miller et al., 2003; Hoare et al., 2008). According to recent work main anchor positions 2 and 9 prefer large hydrophobic residues, while position 7 can accommodate large hydrophobic, small or rigid residues. Lastly, positions 1, 3, and 6 seem to hold preference for smaller residues as this might promote binding of anchor residues at positions 2 and 7 (Figure 8) (Ruibal et al., 2020).
FIGURE 8. Schematic representation of the HLA-E/β2m/peptide/NKG2A/CD94 axis of the immune complex studies. (left) Inhibition of NK cell-mediated cytotoxicity through the harmonious network of interactions between different parts of ligand and receptor, (center) Main study results: 1) the P2, P7, and P9 anchor positions, 2) significant contribution of the P5 residue to binding, 3) securing the P1 anchor with N-terminal capping of the B7 peptide, and 4) introduction of hydrophobic amino acid can improve the binding contribution of P6, nonetheless introduction of bulky amino acid on P6 can result in peptides’ detachment from the binding furrow. (right) Hypothetical scenario in case the designer peptide is introduced and abolishes the interaction between the HLA-E ligand and the NKG2A/CD94 receptor.
By performing an extensive computational study of molecular recognition between wild type peptides G, B27 and B7 with known HLA-E binding properties we identified two peptide positions with the least contribution to binding, namely P1 and P6. We further explained that the negative contribution of P1 to peptide binding is likely a consequence of the N-terminal charge, since capping of P1 resulted in a favorable energy contribution of P1 to peptide binding. Therefore, changing P1s -NH3+ to -CH3 could lead to more optimal interactions between P1 and HLA-E. This retrospective section of our study provided the essential data of the peptides binding, serving as the baseline to assess further peptide alterations.
By generating and simulating several variants of the B7 peptide we demonstrated the preference of position 6 for isoleucine, as the T6I mutation resulted in the highest binding affinity. On the other hand, we confirmed that P6 is indeed incompatible for bulky amino acids as previously suggested (Ruibal et al., 2020). Namely, the T6W mutation resulted in the peptide detachment at the N-terminal site, indicating that tryptophan interferes most with binding of the P2 anchor position and also reduces the contribution of the P3 residue, but apparently has no effect on the interactions of the C-terminal anchors P7 and P9. From this and the calculated per-residue distribution of binding free energy, of all peptides, we can infer that the N-terminal site (P1–P4) is more weakly bound to the binding furrow of the peptide than the C-terminal site (P6–P9) of the peptide. This prompt us with notion that optimizing peptide binding at the N-terminal site would be more prudent to ensure better binding. Peptide processing with ERAP offers a possible explanation for the observed bindings, as this enzyme trims the MHC-I bound precursor peptides at their N-terminal sites, allegedly after peptides bind to its HLA-E binding pocket with their C-terminal residues (Li et al., 2019). However, the stronger binding at the C-terminal site could also be an indication that it is more important for receptor interaction. Thus, mutations in this part might be more likely to lead to impaired interaction with HLA-E receptors.
Although the T6M peptide variant does not result in enhanced binding to HLA-E in the HLA-ET6M compared to the wild type peptide B7, the simulation of the NKG2A/HLA-ET6M model predicts that this variant is likely to abolish recognition and inhibitory signal transduction via the inhibitory receptor CD94/NKG2A when bound to HLA-E. The observed key interactions mimicked the behavior of the model with the R5V peptide variant that abrogates HLA-E-mediated protection. As in our previous work (Prašnikar et al., 2021a) no large-scale conformational changes were observed while simulating HLA-E/β2m/peptide ligand or full HLA-E/β2m/peptide/NKG2A/CD94 immune complexes. This further substantiates the notion how a small change in amino acid sequence can lead to a profound effect on binding and molecular recognition in this important biological system.
Some support of this observation could be found in the experiments performed with the HLA-B leader peptides. The HLA-B leader peptides with methionine on the P2 position (e.g., peptide B7) binds to HLA-E with much higher affinity as those with threonine (e.g., peptide B27) (Braud et al., 1998a) and this was also confirmed in our computational study (see Section 3.1). Peptide binding also secures the HLA-E expression on the cell surface, and consequently HLA-B leader peptides with methionine on P2 result in high and those with threonine on P2 in moderate level of the surface HLA-E ligand (Lauterbach et al., 2015b). In addition, HLA-E surface density is also affected by the *01:01, *01:03 HLA-E polymorphism (Lauterbach et al., 2015b).
Nonetheless, higher peptide affinity for HLA-E does not lead to stronger inhibition of the effector immune cell cytotoxic effect as HLA-G leader sequence (peptide G) is the most effective mediator of the inhibitory signalling (Kaiser et al., 2005; Lauterbach et al., 2015b), despite having relatively low affinity for HLA-E (Braud et al., 1998a) and low surface expression (Lauterbach et al., 2015b). Furthermore, in our previous work, we simulated the HLA-B27 leader peptide in full immune complex in which it displayed features closer to the models that do not allow inhibitory signal transduction via NKG2A/CD94 (Prašnikar et al., 2021a). Thus, the HLA-B27 leader peptide in complex with HLA-E might indeed not be a ligand for the NKG2A/CD94 receptor. Moreover, it could be possible that either the HLA-B leader peptide-dependent HLA-E surface density or the ability of the HLA-B leader peptide to mediate signals through the inhibitory receptor NKG2A/CD94 could influence the (re)education and responsiveness of donor immune effector cells (Watzl et al., 2014; He and Tian, 2017) in hematopoietic stem cell transplantation and lead to different outcomes depending on the polymorphism in the exon 1 of HLA-B (Petersdorf et al., 2020a; Petersdorf et al., 2020b).
In this study we show that MD simulations can be used as tools to predict the peptide binding affinity for HLA-E and the potential of the resulting immune complex to inhibit NK or CD8+ T-cell-mediated cytolysis via the NKG2A/CD94 interaction. The data obtained also lay the necessary structure-based foundation for future development of potential therapeutic peptides, peptidomimetics, or even small molecules that would bind to the HLA-E ligand and then abrogate NKG2A/CD94 recognition. The utility of any dynamic model obtained can be further extended to describe the nature of peptide binding to HLA-E/β2m at the atomistic level, offsetting the disadvantage of the lower throughput of MD simulations compared with cell assays (D’Annessa et al., 2020). It should be noted that our previous studies of molecular recognition between HLA-E ligand and either NKG2A/CD94 or NKG2C/CD94 receptors outlined different key residues and ligand receptor interactions as being vital for signal transduction (Prašnikar et al., 2021a; Prašnikar et al., 2021b). Hence, the results obtained by utilizing the NKG2A/CD94 receptor cannot be straightforwardly transferred to other HLA-E receptors such as NKG2C or T-cell receptor (TCR).
Designer peptides, which are likely the next step in the development towards potential therapeutic application of the outlined approach, could theoretically bind to HLA-E intracellularly or extracellularly. Since MHC class I molecules generally do not bind internalized peptides, but only those formed within the cell (Jensen et al., 1999), designer peptides for intracellular incorporation into the HLA-E complex would likely need to be introduced into the cell in a form of vaccine that allows their intracellular assembly (vector, mRNA vaccines). However, this fact results in several challenges that would have to be overcome, such as selective uptake, TAP permissibility for uptake of designer peptide to the endoplasmic reticulum, and the influence of peptide loading complex machinery on designer peptide binding. In cellular immunology, peptides are commonly added to cell cultures to test the NK and T-cell responses, and this results in antigen specific activation. There are several processes that can mediate this presentation, including cross-presentation (Adiko et al., 2015; Embgenbroich and Burgdorf, 2018; Montealegre and Van Endert, 2019), which allows the presentation of internalized peptides on the MHC-I molecules such as HLA-E.
On the other hand, extracellularly, designer peptides could bind to the HLA-E/β2m as a “rescue peptide” (Chang et al., 2011; Walters et al., 2020) after dissociation of the peptide originally presented by HLA-E/β2m on the cell surface (Jensen et al., 1999; Kraemer et al., 2015; Sturm et al., 2020). Lower pH values could be beneficial for such peptide exchange, as they would facilitate dissociation of the original peptides (Jensen et al., 1999; Sturm et al., 2020). Indeed, mild acid elution is one of the most commonly used methods to secure the cellular MHC peptide ligands (Sturm et al., 2020). In fact, a lower pH was detected at tumor sites (Pillai et al., 2019) and in inflamed tissues (Menkin and Warner, 1937; Punnia-Moorthy, 1987), which is characteristic of the environment of senescent cells (Prašnikar et al., 2020). This circumstance could favorably complement the upregulation of HLA-E on these cells (Pereira et al., 2019; Eugène et al., 2020; Yang et al., 2021) and provide the necessary selectivity in the exchange of the original with designer peptide in case of therapeutic application. Furthermore, if such designer peptide in complex with HLA-E is not a ligand for the inhibitory receptor, NKG2A/CD94+ immune cells would be able to remove the target cells previously protected with HLA-E/β2m/peptide (Figure 8).
To conclude, our extensive molecular dynamics (MD) simulations of the HLA-E/β2m/peptide complex enhanced the atomistic understanding of the peptide binding to its binding pocket on HLA-E and the influence of single mutations to binding. Subsequent simulations of HLA-E/β2m/peptide/NKG2A/CD94 immune complex shed a light on the effect of introduced point mutations on the inhibitory signal transduction that protects the cell from the cytotoxic effects of the immune NK and CD8+ T-cells. Furthermore, performed simulations provide an atomistic interpretation of how a small change in the amino acid sequence can lead to a profound effect on binding and molecular recognition ultimately effecting the signal transduction. Subsequent experimental corroboration would be possible via cytotoxicity assays using NKG2A+ NK cells as effector cells and HLA-E-expressing cells previously incubated with proposed peptides as target cells. In them caspase-3 activity and CD107a expression levels could be monitored to detect the NK cell cytotoxicity and activity. Knowledge provided under the scope of this study can form the basis for the targeted development of designer peptides that would bind to the HLA-E expressed on the target cell surface and abrogate the NKG2A/CD94+ immune cell inhibition, over pronounced in some pathological conditions and in senescent cell accumulation, potentially leading to new cancer immunotherapies as well as new treatment of age-related diseases, by modulating the number of senescent cells.
Data Availability Statement
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.
Author Contributions
Conceptualization, AP and JB; methodology, EP, AP, and JB; validation, EP, AP, and JB; formal analysis, EP, AP, and JB; investigation, EP, AP, and JB; data curation, EP, AP, and JB; supervision, AP and JB, writing—original draft preparation, EP, AP, and JB; visualization, EP, AP, and JB. All authors have read and agreed to the published version of the manuscript.
Funding
This research was funded by the Slovenian Research Agency (ARRS), through research programs P1-0017 (JB) and P1-0012 (AP) research project J1-3019 (JB) and young researcher’s grant 39011 (EP).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
JB, AP, and EP thank the Slovenian Research Agency (grants P1-0017, P1-0012, J1-3019 and young researcher’s program number 39011) for financial support. We acknowledge the Azman high-performance computing (HPC) center at the National Institute of Chemistry in Ljubljana for providing computational support. We thank professor Gerhard Wolber from FU Berlin, Germany for enabling us to calculate dynophore models on the FU computer cluster.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2022.925427/full#supplementary-material
References
Abd Hamid, M., Wang, R.-Z., Yao, X., Fan, P., Li, X., Chang, X.-M., et al. (2019). Enriched HLA-E and CD94/NKG2A interaction limits antitumor CD8+ tumor-infiltrating T lymphocyte responses. Cancer Immunol. Res. 7, 1293–1306. doi:10.1158/2326-6066.CIR-18-0885
Adiko, A. C., Babdor, J., Gutiérrez-Martínez, E., Guermonprez, P., and Saveanu, L. (2015). Intracellular transport routes for MHC I and their relevance for antigen cross-presentation. Front. Immunol. 6, 335. doi:10.3389/fimmu.2015.00335
Berendsen, H. J. C., Postma, J. P. M., van Gunsteren, W. F., DiNola, a., and Haak, J. R. (1984). Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81, 3684–3690. doi:10.1063/1.448118
Bermudez, M., Rakers, C., and Wolber, G. (2015). Structural characteristics of the allosteric binding site represent a key to subtype selective modulators of muscarinic acetylcholine receptors. Mol. Inf. 34, 526–530. doi:10.1002/minf.201500025
Blees, A., Januliene, D., Hofmann, T., Koller, N., Schmidt, C., Trowitzsch, S., et al. (2017). Structure of the human MHC-I peptide-loading complex. Nature 551, 525–528. doi:10.1038/nature24627
Bock, A., Bermudez, M., Krebs, F., Matera, C., Chirinda, B., Sydow, D., et al. (2016). Ligand binding ensembles determine graded agonist efficacies at a G protein-coupled receptor. J. Biol. Chem. 291, 16375–16389. doi:10.1074/jbc.M116.735431
Borišek, J., Saltalamacchia, A., Spinello, A., and Magistrato, A. (2020). Exploiting cryo-EM structural information and all-atom simulations to decrypt the molecular mechanism of splicing modulators. J. Chem. Inf. Model. 60, 2510–2521. doi:10.1021/acs.jcim.9b00635
Borrego, F., Kabat, J., Kim, D.-K., Lieto, L., Maasho, K., Peña, J., et al. (2002). Structure and function of major histocompatibility complex (MHC) class I specific receptors expressed on human natural killer (NK) cells. Mol. Immunol. 38, 637–660. doi:10.1016/s0161-5890(01)00107-9
Borrego, F., Ulbrecht, M., Weiss, E. H., Coligan, J. E., and Brooks, A. G. (1998). Recognition of human histocompatibility leukocyte antigen (HLA)-E complexed with HLA class I signal sequence-derived peptides by CD94/NKG2 confers protection from natural killer cell-mediated lysis. J. Exp. Med. 187, 813–818. doi:10.1084/jem.187.5.813
Borst, L., van der Burg, S. H., and van Hall, T. (2020). The NKG2A-HLA-E axis as a novel checkpoint in the tumor microenvironment. Clin. Cancer Res. 26, 5549–5556. doi:10.1158/1078-0432.CCR-19-2095
Braud, V., Yvonne Jones, E., and McMichael, A. (1997). The human major histocompatibility complex class Ib molecule HLA-E binds signal sequence-derived peptides with primary anchor residues at positions 2 and 9. Eur. J. Immunol. 27, 1164–1169. doi:10.1002/eji.1830270517
Braud, V. M., Allan, D. S., O’Callaghan, C. A., Sӧderstrӧm, K., D’Andrea, A., Ogg, G. S., et al. (1998a). HLA-E binds to natural killer cell receptors CD94/NKG2A, B and C. Nature 391, 795–799. doi:10.1038/35869
Braud, V. M., Allan, D. S., Wilson, D., and McMichael, A. J. (1998b). TAP-and tapasin-dependent HLA-E surface expression correlates with the binding of an MHC class I leader peptide. Curr. Biol. 8, 1–10. doi:10.1016/s0960-9822(98)70014-4
Brooks, A. G., Borrego, F., Posch, P. E., Patamawenu, A., Scorzelli, C. J., Ulbrecht, M., et al. (1999). Specific recognition of HLA-E, but not classical, HLA class I molecules by soluble CD94/NKG2A and NK cells. J. Immunol. 162, 305–313.
Case, D. A., Ben-Shalom, I. Y., Brozell, S. R., Cerutti, D. S., and Cheatham, T. E. (2018). Amber 2018. San Francisco: University of California.
Celik, A. A., Kraemer, T., Huyton, T., Blasczyk, R., and Bade-Dӧding, C. (2016). The diversity of the HLA-E-restricted peptide repertoire explains the immunological impact of the Arg107Gly mismatch. Immunogenetics 68, 29–41. doi:10.1007/s00251-015-0880-z
Chang, C., Dai, L., Tan, Z. W., Choo, J., Bertoletti, A., Grotenbreg, G. M., et al. (2011). Sources of diversity in T cell epitope discovery. Front. Biosci. 16, 3014–3035. doi:10.2741/3895
Chen, L., Reyes-Vargas, E., Dai, H., Escobar, H., Rudd, B., Fairbanks, J., et al. (2014). Expression of the mouse MHC class Ib H2-T11 gene product, a paralog of H2-T23 (Qa-1) with shared peptide-binding specificity. J. Immunol. 193, 1427–1439. doi:10.4049/jimmunol.1302048
D’Annessa, I., Di Leva, F. S., La Teana, A., Novellino, E., Limongelli, V., Di Marino, D., et al. (2020). Bioinformatics and biosimulations as toolbox for peptides and peptidomimetics design: Where are we? Front. Mol. Biosci. 7, 66. doi:10.3389/fmolb.2020.00066
Dolinsky, T. J., Nielsen, J. E., McCammon, J. A., and Baker, N. A. (2004). PDB2PQR: An automated pipeline for the setup of Poisson-Boltzmann electrostatics calculations. Nucleic Acids Res. 32, W665–W667. doi:10.1093/nar/gkh381
Dror, R. O., Dirks, R. M., Grossman, J., Xu, H., and Shaw, D. E. (2012). Biomolecular simulation: A computational microscope for molecular biology. Annu. Rev. Biophys. 41, 429–452. doi:10.1146/annurev-biophys-042910-155245
Embgenbroich, M., and Burgdorf, S. (2018). Current concepts of antigen cross-presentation. Front. Immunol. 9, 1643. doi:10.3389/fimmu.2018.01643
Eugène, J., Jouand, N., Ducoin, K., Dansette, D., Oger, R., Deleine, C., et al. (2020). The inhibitory receptor CD94/nkg2a on CD8+ tumor-infiltrating lymphocytes in colorectal cancer: A promising new druggable immune checkpoint in the context of HLAE/β2m overexpression.. Mod. Pathol. 33, 468–482. doi:10.1038/s41379-019-0322-9
Farhadi, T., and Hashemian, S. M. (2018). Computer-aided design of amino acid-based therapeutics: A review. Drug Des. devel. Ther. 12, 1239–1254. doi:10.2147/DDDT.S159767
Genheden, S., and Ryde, U. (2015). The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities. Expert Opin. Drug Discov. 10, 449–461. doi:10.1517/17460441.2015.1032936
Harvey, M., and De Fabritiis, G. (2009). An implementation of the smooth particle mesh Ewald method on GPU hardware. J. Chem. Theory Comput. 5, 2371–2377. doi:10.1021/ct900275y
He, Y., and Tian, Z. (2017). NK Cell Education via Nonclassical MHC and Non-MHC Ligands. Cell Mol. Immunol. 14, 321–330. doi:10.1038/cmi.2016.26
Hoare, H. L., Sullivan, L. C., Clements, C. S., Ely, L. K., Beddoe, T., Henderson, K. N., et al. (2008). Subtle changes in peptide conformation profoundly affect recognition of the non-classical MHC class I molecule HLA-E by the CD94-NKG2 natural killer cell receptors. J. Mol. Biol. 377, 1297–1303. doi:10.1016/j.jmb.2008.01.098
Humphrey, W., Dalke, A., and Schulten, K. (1996). Vmd: Visual molecular dynamics. J. Mol. Graph. 14, 33–38. doi:10.1016/0263-7855(96)00018-5
Janežič, M., Valjavec, K., Loboda, K. B., Herlah, B., Ogris, I., Kozorog, M., et al. (2021). Dynophore-based approach in virtual screening: A case of human dna topoisomerase IIα. Int. J. Mol. Sci. 22, 13474. doi:10.3390/ijms222413474
Jensen, P. E., Weber, D. A., Thayer, W. R., Westerman, L. E., and Dao, C. T. (1999). Peptide exchange in MHC molecules. Immunol. Rev. 172, 229–238. doi:10.1111/j.1600-065x.1999.tb01368.x
Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., and Klein, M. L. (1983). Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935. doi:10.1063/1.445869
Kaiser, B. K., Barahmand-pour, F., Paulsene, W., Medley, S., Geraghty, D. E., Strong, R. K., et al. (2005). Interactions between NKG2x immunoreceptors and HLA-E ligands display overlapping affinities and thermodynamics. J. Immunol. 174, 2878–2884. doi:10.4049/jimmunol.174.5.2878
Kraemer, T., Celik, A. A., Huyton, T., Kunze-Schumacher, H., Blasczyk, R., Bade-Dӧding, C., et al. (2015). HLA-E: Presentation of a broader peptide repertoire impacts the cellular immune response—implications on HSCT outcome. Stem Cells Int. 2015, 346714. doi:10.1155/2015/346714
Lau, J. L., and Dunn, M. K. (2018). Therapeutic peptides: Historical perspectives, current development trends, and future directions. Bioorg. Med. Chem. 26, 2700–2707. doi:10.1016/j.bmc.2017.06.052
Lauterbach, N., Wieten, L., Popeijus, H. E., Voorter, C. E., and Tilanus, M. G. (2015a). HLA-E regulates NKG2C+ natural killer cell function through presentation of a restricted peptide repertoire. Hum. Immunol. 76, 578–586. doi:10.1016/j.humimm.2015.09.003
Lauterbach, N., Wieten, L., Popeijus, H., Vanderlocht, J., van Zon, P., Voorter, C., et al. (2015b). Peptide-induced HLA-E expression in human PBMCs is dependent on peptide sequence and the HLA-E genotype. Tissue antigens 85, 242–251. doi:10.1111/tan.12525
Li, L., Batliwala, M., and Bouvier, M. (2019). ERAP1 enzyme-mediated trimming and structural analyses of MHC I-bound precursor peptides yield novel insights into antigen processing and presentation. J. Biol. Chem. 294, 18534–18544. doi:10.1074/jbc.RA119.010102
Loncharich, R. J., Brooks, B. R., and Pastor, R. W. (1992). Langevin dynamics of peptides: The frictional dependence of isomerization rates of N-acetylalanyl-N’-methylamide. Biopolymers 32, 523–535. doi:10.1002/bip.360320508
Massova, I., and Kollman, P. A. (2000). Combined molecular mechanical and continuum solvent approach (MM-PBSA/GBSA) to predict ligand binding. Perspect. Drug Discov. Des. 18, 113–135. doi:10.1023/A:1008763014207
Menkin, V., and Warner, C. R. (1937). Studies on inflammation: XIII. Carbohydrate metabolism, local acidosis, and the cytological picture in inflammation. Am. J. Pathol. 13, 25–44.1.
Michaëlsson, J., Teixeira de Matos, C., Achour, A., Lanier, L. L., Kärre, K., Söderström, K., et al. (2002). A signal peptide derived from hsp60 binds HLA-E and interferes with CD94/NKG2A recognition. J. Exp. Med. 196, 1403–1414. doi:10.1084/jem.20020797
Miller, J. D., Weber, D. A., Ibegbu, C., Pohl, J., Altman, J. D., Jensen, P. E., et al. (2003). Analysis of HLA-E peptide-binding specificity and contact residues in bound peptide required for recognition by CD94/NKG2. J. Immunol. 171, 1369–1375. doi:10.4049/jimmunol.171.3.1369
Monaco, E. L., Sibilio, L., Melucci, E., Tremante, E., Suchanek, M., Horejsi, V., et al. (2008). HLA-E: Strong association with beta2-microglobulin and surface expression in the absence of HLA class I signal sequence-derived peptides. J. Immunol. 181, 5442–5450. doi:10.4049/jimmunol.181.8.5442
Montealegre, S., and Van Endert, P. M. (2019). Endocytic recycling of MHC class I molecules in non-professional antigen presenting and dendritic cells. Front. Immunol. 9, 3098. doi:10.3389/fimmu.2018.03098
O’Callaghan, C. A., Tormo, J., Willcox, B. E., Braud, V. M., Jakobsen, B. K., Stuart, D. I., et al. (1998). Structural features impose tight peptide binding specificity in the nonclassical MHC molecule HLA-E. Mol. Cell 1, 531–541. doi:10.1016/s1097-2765(00)80053-2
Ogg, G., Cerundolo, V., and McMichael, A. J. (2019). Capturing the antigen landscape: HLA-E, CD1 and MR1. Curr. Opin. Immunol. 59, 121–129. doi:10.1016/j.coi.2019.07.006
Pereira, B. I., Devine, O. P., Vukmanovic-Stejic, M., Chambers, E. S., Subramanian, P., Patel, N., et al. (2019). Senescent cells evade immune clearance via HLA-E-mediated NK and CD8+ T cell inhibition. Nat. Commun. 10, 2387. doi:10.1038/s41467-019-10335-5
Petersdorf, E. W., Carrington, M., O’hUigin, C., Bengtsson, M., De Santis, D., Dubois, V., et al. (2020a). Role of HLA-B exon 1 in graft-versus-host disease after unrelated haemopoietic cell transplantation: A retrospective cohort study. Lancet. Haematol. 7, e50–e60. doi:10.1016/S2352-3026(19)30208-X
Petersdorf, E. W., Stevenson, P., Bengtsson, M., De Santis, D., Dubois, V., Gooley, T., et al. (2020b). HLA-B leader and survivorship after HLA-mismatched unrelated donor transplantation. Blood 136, 362–369. doi:10.1182/blood.2020005743
Petrie, E. J., Clements, C. S., Lin, J., Sullivan, L. C., Johnson, D., Huyton, T., et al. (2008). CD94-NKG2A recognition of human leukocyte antigen (HLA)-E bound to an HLA class I leader sequence. J. Exp. Med. 205, 725–735. doi:10.1084/jem.20072525
Pillai, S. R., Damaghi, M., Marunaka, Y., Spugnini, E. P., Fais, S., Gillies, R. J., et al. (2019). Causes, consequences, and therapy of tumors acidosis. Cancer Metastasis Rev. 38, 205–222. doi:10.1007/s10555-019-09792-7
Prašnikar, E., Borišek, J., and Perdih, A. (2020). Senescent cells as promising targets to tackle age-related diseases. Ageing Res. Rev. 66, 101251. doi:10.1016/j.arr.2020.101251
Prašnikar, E., Perdih, A., and Borišek, J. (2021a). All-atom simulations reveal a key interaction network in the HLA-E/NKG2A/CD94 immune complex fine-tuned by the nonameric peptide. J. Chem. Inf. Model. 61, 3593–3603. doi:10.1021/acs.jcim.1c00414
Prašnikar, E., Perdih, A., and Borišek, J. (2021b). Nonameric peptide orchestrates signal transduction in the activating HLA-E/NKG2C/CD94 immune complex as revealed by all-atom simulations. Int. J. Mol. Sci. 22, 6670. doi:10.3390/ijms22136670
Punnia-Moorthy, A. (1987). Evaluation of pH changes in inflammation of the subcutaneous air pouch lining in the rat, induced by carrageenan, dextran and Staphylococcus aureus. J. Oral Pathol. 16, 36–44. doi:10.1111/j.1600-0714.1987.tb00674.x
Ruibal, P., Franken, K. L., van Meijgaarden, K. E., van Loon, J. J., van der Steen, D., Heemskerk, M. H., et al. (2020). Peptide binding to HLA-E molecules in humans, nonhuman primates, and mice reveals unique binding peptides but remarkably conserved anchor residues. J. Immunol. 205, 2861–2872. doi:10.4049/jimmunol.2000810
Ryckaert, J.-P., Ciccotti, G., and Berendsen, H. J. (1977). Numerical integration of the cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes. J. Comput. Phys. 23, 327–341. doi:10.1016/0021-9991(77)90098-5
Salari, R., and Chong, L. T. (2010). Desolvation costs of salt bridges across protein binding interfaces: Similarities and differences between implicit and explicit solvent models. J. Phys. Chem. Lett. 1, 2844–2848. doi:10.1021/jz1010863
SchrӧdingerLLC (2015). The PyMOL molecular graphics system, version 2.0. Schrӧdinger, LLC. Available at: https://pymol.org/2/support.html?.
Shenkin, P. S., and McDonald, D. Q. (1994). Cluster analysis of molecular conformations. J. Comput. Chem. 15, 899–916. doi:10.1002/jcc.540150811
Shrivastava, A., Nunn, A., and Tweedle, M. (2009). Designer peptides: Learning from nature. Curr. Pharm. Des. 15, 675–681. doi:10.2174/138161209787315620
Sturm, T., Sautter, B., Wörner, T. P., Stevanović, S., Rammensee, H.-G., Planz, O., et al. (2020). Mild acid elution and MHC immunoaffinity chromatography reveal similar albeit not identical profiles of the HLA class I immunopeptidome. J. Proteome Res. 20, 289–304. doi:10.1021/acs.jproteome.0c00386
Sweet, R. M., and Eisenberg, D. (1983). Correlation of sequence hydrophobicities measures similarity in three-dimensional protein structure. J. Mol. Biol. 171, 479–488. doi:10.1016/0022-2836(83)90041-4
Sydow, D. (2015). Dynophores: Novel dynamic pharmacophores. Master’s Thesis. Berlin (Germany): Humboldt-Universität zu Berlin.
Valés-Gómez, M., Reyburn, H. T., Erskine, R. A., López-Botet, M., and Strominger, J. L. (1999). Kinetics and peptide dependency of the binding of the inhibitory NK receptor CD94/NKG2-A and the activating receptor CD94/NKG2-C to HLA-E. EMBO J. 18, 4250–4260. doi:10.1093/emboj/18.15.4250
Van Der Spoel, D., Lindahl, E., Hess, B., Groenhof, G., Mark, A. E., Berendsen, H. J. C., et al. (2005). Gromacs: Fast, flexible, and free. J. Comput. Chem. 26, 1701–1718. doi:10.1002/jcc.20291
Vivier, E., Tomasello, E., Baratin, M., Walzer, T., and Ugolini, S. (2008). Functions of natural killer cells. Nat. Immunol. 9, 503–510. doi:doi:doi:10.1038/ni1582
Walters, L. C., McMichael, A. J., and Gillespie, G. M. (2020). Detailed and atypical HLA-E peptide binding motifs revealed by a novel peptide exchange binding assay. Eur. J. Immunol. 50, 2075–2091. doi:10.1002/eji.202048719
Watzl, C., Urlaub, D., Fasbender, F., and Claus, M. (2014). Natural Killer Cell Regulation - Beyond the Receptors. F1000 Med. Rep. 6, 87. doi:10.12703/P6-87
Wolber, G., and Langer, T. (2005). LigandScout: 3-D pharmacophores derived from protein-bound ligands and their use as virtual screening filters. J. Chem. Inf. Model. 45, 160–169. doi:10.1021/ci049885e
Keywords: immune system checkpoint, peptide design, HLA-E, NK cell, T-cell
Citation: Prašnikar E, Perdih A and Borišek J (2022) What a Difference an Amino Acid Makes: An All-Atom Simulation Study of Nonameric Peptides in Inhibitory HLA-E/NKG2A/CD94 Immune Complexes. Front. Pharmacol. 13:925427. doi: 10.3389/fphar.2022.925427
Received: 21 April 2022; Accepted: 08 June 2022;
Published: 04 August 2022.
Edited by:
Luciano Saso, Sapienza University of Rome, ItalyReviewed by:
Suzana Katarina Straus, University of British Columbia, CanadaDaniel Fürst, Institute for Clinical Transfusion Medicine and Immunogenetics Ulm, Germany
Simone A. Joosten, Leiden University Medical Center, Netherlands
Copyright © 2022 Prašnikar, Perdih and Borišek. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Andrej Perdih, andrej.perdih@ki.si; Jure Borišek, jure.borisek@ki.si