Allosteric regulation of binding specificity of HVEM for CD160 and BTLA ligands upon G89F mutation.

Molecular interactions mediated by engagement of the Herpes virus entry mediator (HVEM) with members of TNF and Ig superfamily generate distinct signals in T cell activation pathways that modulate inflammatory and inhibitory responses. HVEM interacts with CD160 and B and T lymphocyte attenuator (BTLA), both members of the immunoglobulin (Ig) superfamily, which share a common binding site that is unique from that of LIGHT, a TNF ligand. BTLA or CD160 engagement with HVEM deliver inhibitory or stimulatory signals to the host immune response in a context dependent fashion, whereas HVEM engagement with LIGHT results in pro-inflammatory responses. We identified a mutation in human HVEM, G89F, which directly interferes with the human LIGHT interaction, but interestingly, also differentially modulates the binding of human BTLA and CD160 via an apparent allosteric mechanism involving recognition surfaces remote from the site of the mutation. Specifically, the G89F mutation enhances binding of CD160, while decreasing that of BTLA to HVEM in cell-based assays. Molecular dynamics simulations for wild-type and G89F mutant HVEM, bound to different sets of ligands, were performed to define the molecular basis of this unexpected allosteric effect. These results were leveraged to design additional human HVEM mutants with altered binding specificities.

Introduction

Herpes virus entry mediator (HVEM) delivers co-stimulatory and co-inhibitory signals (Murphy et al., 2006; Ward-Kavanagh et al., 2016) that modulate the course of the immune response following engagement of the T cell receptor (TCR) and peptide/MHC complex (Bretscher and Cohn, 1970; Lafferty and Cunningham, 1975). HVEM, a member of the tumor necrosis factor receptor superfamily (TNFRSF), interacts with members of both the immunoglobulin (IgSF) and tumor necrosis factor (TNFSF) superfamilies. Depending upon multiple factors, including the identity of bound ligands (Cai and Freeman, 2009; Cai et al., 2008), HVEM can deliver stimulatory or inhibitory signals. A co-inhibitory response is associated with the interaction between HVEM and CD160 (Cai et al., 2008) or BTLA (B- and T-lymphocyte attenuator) (Sedy et al., 2005), both members of the IgSF. In contrast, a co-stimulatory signal is generated upon HVEM binding to LIGHT or lymphotoxin α (LTα) from the TNFSF (Mauri et al., 1998).

HVEM is a type 1 membrane protein, with an ectodomain composed of four cysteine-rich domains (CRD1 to CRD4) (Bodmer et al., 2002). BTLA (Sedy et al., 2005) and CD160 (Cai et al., 2008) compete for binding to CRD1, whereas LIGHT and LTα (Ware, 2008; del Rio et al., 2010) compete with each other for binding to CRD2 and CRD3. Notably, despite the lack of direct steric overlap, cross-superfamily interactions are cooperative as demonstrated by the fact that LIGHT binding promotes CD160 and BTLA binding (Cai and Freeman, 2009; Cai et al., 2008). LIGHT ligands form a homo-trimer at the cell surface and cluster HVEM, forming assemblies with a 3:3 stoichiometry (Mauri et al., 1998). For its part, the HVEM:BTLA complex initiates an inhibitory signal by inducing phosphorylation of the intracellular immunotyrosine inhibitory motif (ITIM) that recruits tyrosine phosphatases SHP1 and SHP2 (Sedy et al., 2005; Gavrieli et al., 2003; Watanabe et al., 2003). In contrast, both activating and inhibitory signals have been attributed to CD160 (Le Bouteiller et al., 2011). HVEM:CD160 suppresses CD4+ T cell proliferation and cytokine production by regulating negative signals (Cai et al., 2008; Vigano et al., 2014). CD160 was shown to costimulate CD8+ T cells by cross linking MHC ligands (HLA-C) and activate NK cell cytrotoxicity and cytokine production (Barakonyi et al., 2004; Nikolova et al., 2002; Le Bouteiller et al., 2002).

Dysfunctional HVEM is linked to multiple pathologies such as autoimmunity, inadequate responses to infections and cancer (Ward-Kavanagh et al., 2016; del Rio et al., 2010; Coenen et al., 2009; Boice et al., 2016). Blocking the HVEM:LIGHT regulated co-stimulatory signal induces the immune system to control inflammatory and autoimmune pathologies (Shaikh et al., 2001) (del Rio et al., 2010). The role of BTLA:HVEM in immune suppression ranges from counterbalancing LIGHT-activated inflammation (Murphy and Murphy, 2010) to inhibiting tumor-specific immune response (Derre et al., 2010). Interaction with CD160 was shown to deliver antigen-dependent co-inhibitory signal in NKT cells during early innate immune activation (Kim et al., 2019). All of these observations underscore a crucial role for HVEM in regulating and balancing the immune response and support HVEM as an important therapeutics target.

We recently reported the development of a residue-based pharmacophore computational approach (Shrestha et al., 2019, 2021; Yap and Fiser, 2016) to discover single and double mutant variants of human HVEM that alter selectivity for its three ligands (Shrestha et al., 2020). We identified a mutation in HVEM, G89F (HVEMG89F) that induces a dramatic change in binding selectivity for ligands (i.e., CD160 and BTLA) at a distant interface. Notably, this distant allosteric effect increases affinity for CD160, while diminishing the affinity for BTLA.

Here, we use molecular dynamics (MD) simulations to understand the molecular basis of this unique allosteric effect. MD is a useful computational approach for understanding the microscopic details underlying the dynamics of biomolecular systems (Frauenfelder et al., 1991; McCammon et al., 1977; Henzler-Wildman et al., 2007). High-dimensional time-dependent data are collected during the course of a simulation, which capture dynamic features, including correlated motions of different regions of the macromolecule (Ichiye and Karplus, 1991). Decoding these high-dimensional data provide insight into the transfer of dynamic information between distal sites (Gorfe et al., 2008; Scarabelli and Grant, 2013). Additional details of molecular structural alterations are provided by the comparison of contact maps over the course of simulations (Gorfe et al., 2008). Monitoring the formation and disruption of inter-residue contacts over time enables the characterization of conformational transitions (Scarabelli and Grant, 2013). Dynamic networks constructed from cross-correlation analysis of residue contacts reveal the paths for allosteric communication (Doshi et al., 2016; Sethi et al., 2009; Rivalta et al., 2012; Gasper et al., 2012).

Our goal is to use MD simulations to characterize the impact of the G89F mutation on the dynamic properties of HVEM, particularly CRD1 domain (HVEM_CRD1). This analysis will aid in defining the molecular details and contributions of allosteric communication to the experimentally observed binding selectivity of HVEM for CD160 and BTLA. First, we focused on understanding how the conformational dynamics of the mutant complex differs from those in the wild type. We subsequently identified remote intermolecular interactions associated with the G89F mutation in HVEM and their role in establishing selectivity for HVEM ligands. Next, we explored the role played by the third ligand LIGHT, in establishing this selectivity. Finally, based on this new insight regarding the molecular mechanisms of HVEM allostery, we designed new mutants that reverted the observed selectivities. These new designs were experimentally evaluated through mutagenesis and cell-based binding assays.

Results

HVEM has two spatially distinct binding surfaces, which enable binding to ligands from either the Ig or TNF superfamilies. TNF ligands (LIGHT and LTα) share the same surface formed by CRD2 and CRD3, while the IgSF ligands, BTLA and CD160, utilize and compete for binding to a surface on CRD1 (Fig. 1). Our recent study explored mutant variants that cause HVEM to selectively bind one or two of its three ligands (Shrestha et al., 2020). That study explored 15 different single and double mutations to produce HVEM variants with six different specificity modes, binding to only one or to only two out of the three possible ligands. One mutation, HVEMG89F was particularly interesting, as it is located on the LIGHT recognition surface, but unexpectedly induced a strong shift in specificity at the remote surface involved in recognition of BTLA and CD160. Surprisingly, CD160 binding increased, while BTLA binding was drastically reduced at the same remote site. To further explore this mutation, the wild-type and G89F variant of HVEM were expressed in HEK293 ​cells as GFP fusion proteins and queried with PE-conjugated ligands in flow cytometry experiments, and the percent of ligand-bound cells was determined (Fig. 2). While CD160 showed an approximately 50% enhancement in binding, BTLA binding to HVEMG89F was greatly reduced. Titration experiments confirmed the stronger binding of CD160 to HVEMG89F-expressing cells relative to wild type, and demonstrated that BTLA binding only becomes detectable at high concentrations. Next, we explored the molecular basis of this apparent allosteric regulation of binding specificity.

Fig. 1

LIGHT:HVEM:BTLA complex. (A) ribbon model of modeled ternary complex of LIGHT (cyan):HVEM (magenta):BTLA (yellow). Red box indicates the site of G89F mutation (B) enlarged figure of the G89F mutation site, with possible interactions between R172 (LIGHT) G89 and Y47 (HVEM). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 2

Results of cell based essays of HVEM binding to its ligands upon G89F mutation as measured by flow cytometry. Normalized % bound for three ligands – LIGHT, CD160 and BTLA. HVEM binding of ligands BTLA (left) CD160 (center) and LIGHT (right) are shown respectively.

Domain flexibility changes in HVEM upon G89F mutation

HVEM is composed of four cystine rich (CRD) domains of which three have been characterized structurally (Liu et al., 2015). To dissect intra and inter-domain conformational plasticity within HVEM, molecular dynamics (MD) simulations were performed for various complexes, including: wild type (HVEMWT:CD160:LIGHT and HVEMWT:BTLA:LIGHT) and mutant (HVEMG89F:CD160:LIGHT and HVEMG89F:BTLA:LIGHT) ternary complexes as well as all of the corresponding binary complexes lacking the LIGHT ligand. Each of these eight complexes were subjected to 250 ns long MD simulations, each with three replicates. The first 25 ns were not considered in the analysis, to allow sufficient time to relax the starting conformations.

Due to their elongated architecture, TNFR superfamily proteins can be highly flexible. First, we tested whether CRD1 displays different domain movement with respect to LIGHT, which binds CRD2 and CRD3. The distance between the center of masses of LIGHT and CRD1 domains was monitored in a combined 678 snapshots extracted from three replicate MD simulations of LIGHT:HVEM:CD160 and in LIGHT:HVEM:BTLA. Indeed, compared to the HVEMWT, the CRD1 domain of HVEMG89F shifts away from the rest of the HVEM structure and LIGHT by about ∼2 ​Å on average (Fig. 3A). We also examined how this shift changes as a function of bound IgSF ligands: BTLA or CD160. For this, we monitored the RMSD of CRD1 with respect to the rest of the HVEM structure (using the central domain, CRD2 as a reference for superposition) with different ligands bound (Fig. 3B and C). Not only did we observe a shift of CRD1 compared to the rest of the HVEM structure in these cases but interestingly, the domains moved in opposing directions when bound to each IgSF ligand. When HVEMG89F is bound to CD160 there is a tendency toward reduced flexibility, with a mean shift of 0.27 ​Å towards smaller RMSD values compared to the wild type complex (Fig. 3B). This shift in the distribution of RMSD values is statistically significant with a p-value of 2.65∗10−10. Meanwhile, when HVEMG89F is bound to BTLA there is a tendency to increased flexibility, with a mean shift of 0.51 ​Å of the RMSD distributions (p ​= ​1.24∗10−32) (Fig. 3C.). These observations suggest that the mutation results in a significant change in domain movements and the direction of these changes depends on the type of bound IgSF ligands.

Fig. 3

Conformational dynamics of HVEM_CRD1. (A) Distribution of distances (Å) from 678 MD snapshots measured between the center of masses of HVEM_CRD1 and LIGHT domains. (B) Distribution of Cα-RMSD values of CRD1 domain of HVEM with (in red) and without (in green) the G89F mutation over the 678 MD snapshots in reference to wild type structure (in green) when bound to CD160 or (C) BTLA ligand. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Conformational alterations in HVEM associated with the G89F mutation

Next, we explored if in addition to the observed differential domain movements there is a concomitant intradomain conformational change as well. From the MD simulations, we calculated cross-correlations of residue movements (Fig. 4) between HVEMWT and HVEMG89F. Irrespective of the bound ligand (CD160 or BTLA) the cross correlation heatmaps between the domains in the wild-type and mutant complexes (upper and lower triangles along the diagonal) are virtually indistinguishable, suggesting there are no significant differences in the intra-domain conformational dynamics associated with the G89F mutation. In addition to looking at the interdomain sections of the heatmap, we also characterized each area of the heatmap (within each domain boundary) by the median values of cross-correlations (Suppl. Fig 1.). It is clear that intra-domain correlations are the highest (i.e. boxed areas along the diagonal) with corresponding averages of CRD1WT vs CRD1G89F: 0.49 vs 0.49, CRD2WT vs CRD2 G89F: 0.54 vs 0.52 etc. suggesting a coherent intradomain correlation of movements both in wild type and mutant HVEM. The only exception visible in the heatmap is the short blue area in CRD2, which corresponds to a loop stacked on the surface of CRD2. When looking at the off-diagonal areas (i.e. correlation between CRD1 vs CRD2, CRD1 vs CRD3, CRD2 vs CRD3), CRD1 and CRD2 are rather cross-corelated in their movements, while CRD3 has a strong anti-correlation with both CRD1 and CRD2 (Fig. 4 and Suppl. Fig. 1.). These cross-correlation maps allow us to further confirm that the CRD1 of HVEM has a stronger tendency to synchronize its movement with CRD2 when it is bound to CD160 than when bound to BTLA. This can be seen from the increasing white and blue colors in the heatmap of BTLA when correlating CRD2 and CRD1 domains and is reflected in the difference of median values of residue cross correlations of the corresponding areas (0.27 (CD160) vs 0.11 (BTLA) in HVEMWT, and 0.17 (CD160) vs 0.11 (BTLA) in HVEMG89F).

Fig. 4

Heatmaps of residue cross-correlations between wild type and mutant trajectories of HVEM from two different complexes – (A) HVEM:CD160 and (B) HVEM:BTLA. Domain structure is shown on both axes. (C) Corresponding domain structure of HVEM in a ribbon model.

Structural ensembles of HVEM

The analysis of the distribution of RMSD values and the inspection of the cross-correlations only provide insight into the flexibility of CRD1 as a function of the bound ligand, but does not reveal what type of conformations are dominating the ensemble. To gain more insight, we structurally clustered (using a density based spatial clustering approach (Schubert et al., 2017)) the conformations from 678 snapshots over the MD simulation in both cases of binary complexes, using a 1.0 ​Å cutoff value (Table 1). When HVEM is bound to CD160, the number of clusters increases from 5 to 12, but the size of the largest cluster remains similar (196 vs 148), and the number of outliers, i.e., conformations that do not cluster, is also similar (459 vs 455). This is in contrast to the situation when HVEM binds to BTLA, where there is a precipitous drop in the size of the largest cluster by 10 fold (from 664 to 65), and a corresponding five-fold increase of outlier conformations (105–544). This observation suggests that in the context of the HVEM G89F mutation the complex formed with BTLA becomes conformationally more variable, and with a large number of alternative ternary states.

Table 1

Clustering the distribution of HVEM_CRD1 conformations from 678 MD snapshots.

complexes	# of clusters	Size of largest clusters	# of outliers
HVEM:CD160W	5	196	459
HVEM:CD160M	12	148	455
HVEM:BTLAW	1	664	105
HVEM:BTLAM	3	65	544

Dissecting the interface contacts for binding specificity

We observed a structural shift of the CRD1 domain with respect to LIGHT ligand and the rest of HVEM domains as well as differences in distribution of RMSD values of CRD1 when bound to BTLA and CD160. We hypothesized that the underlying cause of the changes in binding affinity observed for the G89F mutation is that this amino acid substitution pushes the HVEM_CRD1 domain into a new configuration, establishing new interfaces and new contacts with the ligands. These new contacts appear to be favorable for CD160 binding but unfavorable for binding BTLA. To explore this hypothesis, we analyzed the formation and disruption of inter-molecular contacts of wild-type and mutant HVEM:CD160 and HVEM:BTLA in MD sampled conformational ensembles.

For each of the 678 MD snapshots for each ternary complex, we ran the CSU program (Sobolev et al., 1999) to detect interface contacts and then normalized the frequency of each observation. We identified 36 contacts at the interfaces of HVEMWT:CD160 and HVEMG89F:CD160, with statistically significant differences in frequency (Suppl. Table 1). These contacts are formed by a total of 18 residues from HVEM and CD160, and as expected these are all from the CRD1 region of HVEM, except residue L87 (CRD2). Among the top twenty contacts with the most significant differences between the two interfaces (ranked by p-value (Fig. 5A (left) – red color)), 70% of contacts (13/20) were observed in HVEMG89F compared to its wild-type counterpart, suggesting that an increased number of contacts are established in the HVEMG89F interface with CD160 (Suppl. Table 1). For example, contacts, K43(HVEM):T109(CD160) (observed frequency of occurrence in HVEMWT 8.40% vs 44.40% in HVEMG89F), V74(HVEM):I121(CD160) (28.17% vs 55.01%), and Y61:R115 (6.04% vs 28.61%) are characteristic for HVEMG89F:CD160. These contacts define new interfacial interactions in HVEMG89F:CD160, which are potentially responsible for increasing the binding affinity. In contrast, in HVEMWT:CD160 only seven out of the twenty most significant contacts are present, including L87(HVEM):D67(CD160), and D45(HVEM):G66(CD160), which were largely disrupted in the G89F mutation. However, further inspection revealed that the number of disrupted contacts are not only fewer than the number of new contacts established in the HVEMG89F (13) but the seven contacts lost upon G89F mutation induced conformational change are predominantly formed between backbone atoms and thus may have less impact on binding strength.

Fig. 5

Distribution of residue contact preferences. (A) frequency of observed contacts are plotted between the wild type (horizontal axis) and mutant (vertical axis) of HVEM and CD160 (left panel) or BTLA (right panel). (B) As in (A) but using MD data when LIGHT ligand was also bound to HVEM. Statistically significantly differently preferred contacts are shown in blue and red colors. Red colors are marking the top 20 most significantly different contact preferences. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

HVEMG89F binding to BTLA is drastically reduced compared to wild type (Fig. 2). When we inspected the altered interface we found 69 residue interactions with significantly altered frequencies between HVEMWT:BTLA and HVEMG89F:BTLA (Suppl. Table 2). In contrast to the HVEM:CD160 system, the number of interactions between the two systems of HVEM:BTLA are comparable (29 and 24) (Fig. 5). Once again we examined the top 20 contacts, which happen to be equally distributed between the two systems, 10 vs.10. However, several key contacts with BTLA that are based on charge or hydrogen bond interactions were dramatically reduced or lost, such as E46(HVEM):H127(BTLA) (HVEMWT 94.70% and HVEMG89F 1.00%), D45(HVEM):S128(BTLA) (HVEMWT 93.80% and HVEMG89F 2.06%), S58(HVEM):E125(BTLA) (HVEMWT 100% and HVEMG89F 12.09%) and S58(HVEM):R114(BTLA) (HVEMWT 73.60% and HVEMG89F 0.00%). HVEMG89F does not appear to establish any physio-chemically compatible new contacts that might compensate for loss of these interactions (Suppl. Table 2), which could explain the loss of binding affinity.

Effect of bound LIGHT ligand on the binding selectivity of HVEMG89F to its IgSF ligands

The analysis of contact formation and disruption between HVEMWT and HVEMG89F was extended by including LIGHT in the MD simulation and considering the LIGHT:HVEM:CD160 and LIGHT:HVEM:BTLA ternary complexes. First, we revisited the conformational clustering analysis when LIGHT is included as well (Table 2). The size of the largest cluster decreases less dramatically in LIGHT:HVEMG89F:BTLA than in HVEMG89F:BTLA from 627 to 465 (Table 1), respectively; however, the number of outlier conformations showed a similar trend, in this case about a 4-times increase from 51 to 206, similarly to the case of HVEM:BTLA binary complex (105–544) (Table 1). At the same time LIGHT:HVEM:CD160 showed similar, statistically insignificant changes of dominant and outlier conformations between the wild-type and mutant versions.

Table 2

Clustering the distribution of HVEM_CRD1 conformations from 678 MD snapshots.

Molecules	# of clusters	Size of largest clusters	# of outliers
LIGHT:HVEM:CD160W	3	251	257
LIGHT:HVEM:CD160M	5	288	255
LIGHT:HVEM:BTLAW	1	627	51
LIGHT:HVEM:BTLAM	3	465	206

Next, we focused on the interface contacts of CD160 and BTLA, in the presence of LIGHT, that dominate their interaction with HVEMWT and HVEMG89F. As before, we looked at the top 20 contacts that changed most significantly between HVEMWT and HVEMG89F according to their p-values (Suppl. Table 3). In the case of LIGHT:HVEM:CD160 an even more significant gain of contacts was observed in HVEM upon mutation: 19 of the 20 most significantly different contacts show up in HVEMG89F only, suggesting the presence of LIGHT ligand enhances contacts within the HVEMG89F:CD160 interface (Fig. 5B). Among the most significant preferred contacts formed in the HVEMG89F:CD160 interface are Y61(HVEM):R115(CD160), S58(HVEM):Q124/D67(CD160), and T71(HVEM):N28(CD160), all of which have the potential to establish hydrogen bonds (Fig. 6). Additional polar contacts that are unique for the HVEMG89F:CD160 interface are K43(HVEM):109T (CD160), E69(HVEM):28N(CD160). In addition, a small network of unique hydrophobic interactions L70(HVEM):I27/123L (CD160) is also established. In contrast, Y61(HVEM):G120(CD160) is the only contact uniquely present in the wild-type interface. Unlike in the LIGHT:HVEM:CD160 ternary complex, the highly significant contacts are distributed more uniformly in LIGHT:HVEM:BTLA. However, unlike in the absence of LIGHT, two thirds of the most significant contacts were observed in LIGHT:HVEMWT:BTLA, and are consistent with the experimentally observed loss in binding associated with the mutant, only one third of contacts are unique in LIGHT:HVEMG89F:BTLA. Among these, E46(HVEM):H127(BTLA) was observed in 77% of the HVEMWT trajectories, but only in 22% of the HVEMG89F trajectories (Fig. 5B, Suppl. Table 4). Similarly, K43(HVEM):S128(BTLA), E44(HVEM):T77(BTLA), D45(HVEM):R144(BTLA), and L70(HVEM):I40(BTLA) appear as pivotal contacts observed in HVEMWT and lost in HVEMG89F (Fig. 6).

Fig. 6

Structural models showing interface residues on the CRD1 domain of HVEM that form contacts with ligands (A) CD160 and (B) BTLA. Silver color marks residues participating both in wild type and mutant HVEM contacts. Light blue and orange colors mark residues that are part of the original binding interface but have strong preferences to HVEMWT and HVEMG89F, respectively. Unique contacts (that were not part of the original interface) that form preferably in HVEMWT or HVEMG89F are shown in dark blue and dark red, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

The overlap between the top 20 most significant contacts between HVEM:CD160 and LIGHT:HVEM:CD160 is 55%, while in the case of HVEM:BTLA and LIGHT:HVEM:BTLA the overlap is lower, 30%, perhaps also reflecting the less well defined, or less systematic binding modes of BTLA upon G89F mutation in HVEM.

Experimental probing of unique interface interactions

To validate our findings and gain additional insight into the allosteric mechanisms contributing to the altered binding selectivity of the HVEMG89F mutant, we attempted to revert the observed binding properties through the introduction of additional mutations. We examined a subset of the newly established contacts correlated with the increased binding of CD160 to HVEMG89F, with the expectation that removal of these contacts observed in the MD simulations would reduce binding. Similarly, we targeted specific residue contacts from the MD simulations that were correlated with BTLA binding to HVEMWT, but were lost in HVEMG89F simulations. By mutating these residues, we expected to recover BTLA binding in the background of HVEMG89F. In case of CD160, the top three ranked positions were altered: Y61 was replaced with R, K, H. All of these mutants diminished CD160 binding to the HVEMG89F mutant as expected (Fig. 7). Next, position S58 was mutated to I and L and these also both diminished CD160 binding as expected. The last position explored was L70, which was replaced with E,R or K. L70R and L70K did not express, but L70E showed diminished binding to CD160, as predicted.

Fig. 7

Effect of mutations on HVEM ligand binding in the background of HVEMG89F. (A) for BTLA (B) CD160. Five concentrations (0, 10, 39, 156, 625) were tested in three repeats in binding assay of HVEM double mutants against CD160 and BTLA. All values are normalized to the WT binding at the highest concentration (set to 1.0).

We also explored whether we could introduce mutations into the HVEMG89F background that would rescue binding to BTLA. E52 was replaced with Q, T, Y or S. Of these only E52S expressed, but it showed similar binding to the BTLA as before. We also mutated G68 to E or F, but these mutants were not successfully expressed. The last position explored was S58 to Q or R, both of which showed a moderate increase in BTLA binding, as expected, but not to the original level of HVEMWT. Overall, the tested mutations were consistent with our observations regarding the new allosterically-coupled interfaces present in the HVEMG89F mutant, especially when removing contacts in case of CD160. Increasing binding affinity is typically more challenging, and the partial success in case of BTLA reflects that difficulty.

Interaction with LIGHT

While the most dramatic effects of G89F mutation appear across the CRD1 domain and alter CD160 and BTLA binding, G89 is also a direct interfacial residue with LIGHT. HVEMG89F maintains binding to LIGHT that is comparable to HVEMWT with only a very modest increase observed. The most likely explanation for this observation is that G89 of HVEM and R172 of LIGHT form contacts through their backbone atoms in HVEM:LIGHT interface, and are therefore not as sensitive to side chain permutations (Fig. 1). The only speculative effect when mutating G89 with bulky aromatic residue, phenylalanine (G89F) is that it can create a new stacking interaction with Y47 in HVEM (Fig. 1 and Suppl. Fig. 2).

Discussion

Using MD simulations we identified a possible mechanism by which the G89F mutation on HVEM, which modifies the binding surface for LIGHT, induces a significant increase in CD160 binding and a near complete loss of BTLA binding at a remote site. The MD results suggest that the bulky G89F mutation forces the CRD1 domain of HVEM to a more open conformation. In turn, this change in configuration establishes a modified recognition surface to the CD160 and BTLA ligands (Suppl. Fig. 3). While a number of original interface contacts are preserved, in the case of CD160, an increased number of new contacts are formed with apparent complementarity between hydrogen bond donors and acceptors. In contrast, while BTLA has the same or fewer contacts in the new interface formed in the HVEMG89F complex, those changes seem to establish less favorable interactions.

When LIGHT is added to the analysis, these observations not only hold, but become even more pronounced. When LIGHT is present, the G89F mutation not only pushes the CRD1 domain away from CRD2 but also from LIGHT, and the pattern of contact formation becomes more extreme. In the case of CD160, 19 out of the top 20 most significantly changed contacts show up in the new CD160 interface. Meanwhile, BTLA has a net loss of 50% of its contacts compared to WT among the top 20 most significantly changed interactions.

Our observations about the dramatic allosteric effects suggest that HVEM, with its flexible domain structure may respond in unpredictable ways to certain mutations, and may have a significant impact in the control of disease processes by the immune system (Cheung et al., 2010). Understanding allosteric changes in promiscuous TNFRs can provide engineering opportunities to introduce novel ligand binding selectivity for the realization of new tool reagents or therapeutic lead compounds.

Material and methods

Crystal structures of complexes

Crystal structures of HVEM:BTLA (2AW2) (Compaan et al., 2005) and HVEM:CD160 (6NG3) (Liu et al., 2019) were obtained from the Protein Data Bank (Berman et al., 2000). Missing residues and atoms were added to crystal structures using Modeller (Sali and Blundell, 1993; Fiser and Sali, 2003). We also obtained the ternary complex of LIGHT:HVEM:CD160 (unpublished); since the crystal structure of LIGHT:HVEM:CD160 is unknown, we used Modeller (Sali and Blundell, 1993; Fiser and Sali, 2003) to build a homology model of complex of HVEM:BTLA:LIGHT using templates of two known complexes of HVEM:BTLA (2AW2) and HVEM:LIGHT (4RSU).

Mutation selection and input preparation

We prepared eight systems in total for MD simulations. Four of them are wild type systems: HVEMWT:CD160, HVEMWT:BTLA, LIGHT:HVEMWT:CD160, and LIGHT:HVEMWT:BTLA. The corresponding four mutant versions were obtained by using Modeller (Sali and Blundell, 1993; Fiser and Sali, 2003), HVEMG89F:CD160, HVEMG89F:BTLA, LIGHT:HVEMG89F:CD160, and LIGHT:HVEMG89F:BTLA. Human LIGHT has a homo-trimeric structure, but we included only a homo-dimer and its interacting partner for the MD simulation since only homo-dimer can interact with HVEM.

Molecular dynamic simulation and analysis

Molecular dynamics simulations for each system were carried out for trajectory length of 250 ns using Gromacs v5.0.6 (Mark James Abraham et al., 2015). All MD runs were replicated 3 times with different initial velocities. AMBER99sb force field (Hornak et al., 2006) and the TIP3P water model (William et al., 1983) in a dodecahedron box were used for simulation. Sodium or chloride ions were added in each system for neutralization of charges followed by energy minimization with a force threshold of 1000.0 ​kJ/mol/nm. Subsequently, a 5ns MD simulation for each system was run to equilibrate the NVT ensemble followed by the NPT ensemble. The trajectory was simulated until 250ns after successful completion of two equilibration steps. The V-rescale thermostat (Bussi et al., 2007) and Parrinello-Rahman pressure (Parrinello aAR, 1981) coupling was employed during the simulation. LINCS algorithm (Hess et al., 1997) was used to constraint the bonds. Periodic boundary conditions were utilized and the cut-offs for electrostatic and van der Waals interactions were set to 1 ​nm. Long range electrostatic interactions were preserved using the particle-mesh Ewald method (Tom Darden and Lee, 1993).

After replicating the simulations three times for each system, the simulated trajectory of final 226ns were considered for subsequent analysis after trimming first 24ns to allow the systems to relax. Snapshots from trimmed trajectory were extracted in 1ns interval resulting in 226 snapshots. Afterward, these snapshots from three replications were combined (total of 678 snapshots) for each system for subsequent analysis.

Cross-correlation

The cross-correlation of Cα atom displacements were computed for HVEM:BTLA/CD160 to identify intra and interdomain correlated motions. We employed an equation described in the by Ichiye and Karplus (1991) to compute the cross-correlation between and atoms, as follows:Where, is the displacement of atom from the mean position of the configuration. ranges between 1 for correlated motion and −1 for anti-correlated motions.

Cluster cysteine rich domain

The conformations of HVEM_CRD1 were clustered with respect to HVEM_CRD2 to understand its dynamics. We computed the distance among the CRD1 domains of HVEM after superposing their corresponding CRD2 using Modeller (Sali and Blundell, 1993; Fiser and Sali, 2003). The precomputed matrix is then used as input to the clustering program, DBSCAN (density-based spatial clustering of application with noise) (Schubert et al., 2017). DBSCAN package from scikit-learn was used with 1.0 ​Å as the distance threshold and 5 members as the value of min_samples.

Interface contact analysis

An interface of each complex from MD simulation was assessed using CSU. Interface residues with a pair of atoms at or closer than 4.0 ​Å were considered as interacting. We obtained the percentage of observed contacts among the 678 snapshots. The percentage occurrence of contacts was compared between wild type and mutant trajectories of each system (HVEM:BTLA and HVEM:CD160 and with and without LIGHT. The significance of the difference in contact formation and loss in wild type and mutant configurations were tested using a statistical z-test.

Cell-based flow cytometry HVEM binding assay

HEK 293 Freestyle cells were maintained in Freestyle media and grown in a shaking, humidified CO2 incubator. Cells for transfection were resuspended in fresh media and plated at a density of 1 ​× ​106 ​cells/mL with 1 ​mL total volume per well of a non-treated 24-well sterile tissue-culture plate. Transfections were carried out using linear Polyethylenimine (PEI) in 4-fold molar excess to DNA (2 ​μg PEI to 0.5 ​μg DNA per 1 ​mL transfection). Two days post transfection cells were counted and diluted to 1 ​× ​106 ​cells/mL with 1X phosphate buffered saline pH7.0 (PBS) and 0.2% BSA to be used directly for FACS binding or titration experiments.

For binding experiments, 100K HVEM transfected cells were plated into 96-well v-bottomed plates and incubated with recombinant his-tagged LIGHT, BTLA and CD160 (R&D systems). For initial binding experiments 0.5 ​μg of BTLA or CD160 was used and 0.25 ​μg LIGHT. Binding was carried out at RT using a 96-well bench top plate shaker (900 ​rpm) for 1 ​h. After binding cells were washed 2X with 1X PBS 0.2%BSA by centrifugation at 500×g. After washing, 0.1 ​μg ​PE conjugated anti-HIS secondary antibody (Abcam) was added to the cells in a final volume of 100 ​μL and allowed to bind for 30min as described above. After secondary antibody binding, cells were washed again 2X and immediately analyzed on a BD Accuri flow cytometer equipped with an Intellicyte Hypercyte autosampler. Titration experiments were carried out in an identical manner except that the BTLA/CD160 protein concentration was titrated from 10−8 -10−6 ​M and the LIGHT protein concentration was increased from 10−8.5-10−7 ​M, as indicated. All titration wells and controls received the same concentration of secondary anti-HIS antibody for detection of bound protein.

All flow data was analyzed by gating for live cells based on FSC-A and SSC-A and sub-gating for all GFP positive cells (HVEM expression). From the HVEM positive population, cells were further gated to identify all PE positive cells (ligand bound). Data for the percent GFP positive of cells, the percent PE positive as a percentage of HVEM positive cells and the GeoMean of PE-A for the HVEM positive gate were collected. For titration experiment analysis, the GeoMean dataset for each experiment was normalized to wild-type binding to the highest ligand concentration. Data show the average of three independent experiments with the calculated standard deviation.

CRediT authorship contribution statement

Rojan Shrestha: Conceptualization, Formal analysis, Investigation, Methodology, Writing – original draft. Sarah Garrett-Thomson: Investigation, Validation, Writing – review & editing. Weifeng Liu: Resources. Steven C. Almo: Conceptualization, Supervision, Writing – review & editing. Andras Fiser: Conceptualization, Funding acquisition, Project administration, Supervision, Writing – original draft, Writing – review & editing.

Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Andras Fiser reports financial support was provided by National Institute of Health. Andras Fiser reports a relationship with National Insitite of Health that includes: funding grants.