Exploring CCRL2 chemerin binding using accelerated molecular dynamics.

Chemokine (C-C motif) receptor-like 2 (CCRL2), is a seven transmembrane receptor closely related to the chemokine receptors CCR1, CCR2, CCR3, and CCR5. Nevertheless, CCRL2 is unable to activate conventional G-protein dependent signaling and to induce cell directional migration. The only commonly accepted CCRL2 ligand is the nonchemokine chemotactic protein chemerin (RARRES2). The chemerin binding to CCLR2 does induce leukocyte chemotaxis, yet, genetic targeting of CCRL2 was shown to modulate the inflammatory response in different experimental models. This mechanism was shown to be crucial for lung dendritic cell migration, neutrophil recruitment, and Natural Killer cell-dependent immune surveillance in lung cancer. To gain more insight in the interactions involved in the CCRL2-chemerin, the binding complexes were generated by protein-protein docking, then submitted to accelerated molecular dynamics. The obtained trajectories were inspected by principal component analyses followed by kernel density estimation to identify the ligand-receptor regions most frequently involved in the binding. To conclude, the reported analyses led to the identification of the putative hot-spot residues involved in CCRL2-chemerin binding.

INTRODUCTION

The chemokine (C‐C motif) receptor‐like 2 (CCRL2) is a seven transmembrane domain receptor that maps within a CC chemokine receptor cluster in human chromosome 9 and mouse chromosome 3 and therefore, CCRL2 is most related to chemokine receptors CCR2, CCR3, and CCR5. Since CCRL2 is unable to activate signal transduction through G proteins or β‐arrestin, it is considered to be related to the family of Atypical Chemokine Receptors. ,

So far, chemerin (encoded by the RARRES2 gene) is the only accepted CCRL2 ligand. Chemerin is produced by mammalian cells as a 163 amino acid (aa) pro‐precursor. The N‐terminal processing of 20 aa leads to the secretion of a precursor form of chemerin. The further C‐terminal cleavage results in both active and deactivated chemerin forms, according to the extent of processing. For example, proteases like plasmin, elastase, and cathepsin G activate chemerin and generate various chemerin isoforms with different affinity to CMKLR1, the active chemerin receptor. Further cleavage of bioactive chemerin by chymase produces an inactive form of chemerin. Thus, the C‐terminal proteolytic processing acts as a regulatory mechanism to control the concentration of active chemerin.

In addition to CCRL2 and CMKLR1, chemerin also binds G protein‐coupled receptor 1 (GPR1); these three receptors have distinct patterns of expression and biological functions.

CMKLR1 is mostly expressed by innate immune cells, such as dendritic cells, macrophages, and Natural Killer (NK) cells , CCRL2 is expressed by a large variety of leukocytes subsets and by barrier cells, such as vascular and lymphatic endothelium and some epithelium. , GPR1 is predominantly expressed in the central nervous system and skin. Among these three receptors, CCRL2 is the only one devoid of the ability to activate an intracellular signaling cascade.

CCRL2 was shown to regulate inflammation‐related diseases such as experimental autoimmune encephalitis, hypersensitivity, and inflammatory arthritis , , , and the recruitment of the NK cells in pathological conditions. , We decided to further characterize CCRL2‐chemerin interaction, by a protein–protein docking followed by accelerated molecular dynamic (aMD) simulation of the proposed binding conformations.

Instead of a classical molecular dynamic approach, which is limited by kinetic trapping effects and limited sampling of the conformational space, it was established a protocol of aMD. This method reduces the energy barrier between different low‐energy states by applying a potential energy boost, increasing the transition probabilities between two different conformations. We obtained a total of 5.5 μM second trajectories, which were analyzed by principal component analyses (PCA) to reduce the complexity of the data and to identify the greatest variance. , Then, the residues more often involved in binding interactions were highlight as hot‐spot residues of the CCRL2 chemerin complex.

MATERIALS AND METHODS

Homology model of CCRL2 conformational states

CCRL2 primary aminoacidic sequence was derived from Uniprot (UniProtKB: O00421). Since human CCRL2 exist in two isoforms, we modeled the shorter isoform 1 (344 aa, CRAM‐B, O00421‐1). The template choice for CCRL2 active state was performed by aligning all the available chemokine receptor structures, in the active state, resolved by crystallography. The best model was US28 (UniProtKB: P69332; PDB: 5WB1; root mean square deviation (RMSD) 0.776 Å).

The resulting CCRL2 structure was obtained with the I‐Tasser server, and as restrain was inserted, the protein topology derived from UniProtKB. Loop refinement and energy minimization was carried out using ModRefiner. The quality of generated model was validated with respect to backbone and side chain geometry. To validate protein backbone quality, the MolProbidity tool was adopted.

Ab initio modeling of chemerin

For chemerin (UniProtKB: Q99969), we lack a proper resolved homologue structure; therefore, we utilized the RAPTOR‐X prediction server to build the model. Loop refinement and energy minimization was carried out using ModRefiner. The quality of generated models was validated with respect to backbone and side chain geometry. To validate protein backbone quality, the MolProbidity tool was adopted.

Structural comparison of modeled proteins

The optimized CCRL2 and chemerin models were compared to their respective AlphaFold (https://alphafold.ebi.ac.uk/, entry: O00421 CCRL2_HUMAN) conformations using Matchmaker function of UCSF Chimera and RMSD values were obtained.

Protein–protein docking

Protein–protein docking was performed on HADDOCK server using the 3D models of CCRL2 and chemerin proteins. We set up as “active” residues (residues expected to be involved in the interaction between the two molecules) the N‐terminal residues 1–100 of chemerin and the extracellular loops of the CCRL2 receptor (from UniProtKB, 25–67; 120–128; 190–122; 284–310) following literature data ; the other residues were defined as “passive” (residues accessible to the solvent closed of the active residues). This docking protocol consisted of three stages: rigid body (it0), semi‐flexible refinement (it1), and explicit solvent refinement (water). The docking experiments were carried out by using default parameters. Just MD steps in the TAD and cooling stage were increased to 2000 for it0. HADDOCK produced 1000 models in the first step, then refined to 200 best model in the following steps. The final models were automatically clustered based on the fraction of common contacts that measures the similarity of the intermolecular contacts. At the end, we obtained 12 clusters, and for each of them, it was selected as representative structure the conformation with the best score value.

Accelerated molecular dynamics

The membrane embedded complexes structures and the Amber parameters were obtained by CHARMM‐GUI server through Membrane builder module.

The structural information to have a reliable placement of the bilayer membrane at CCRL2 were obtained by Uniprot (UniProtKB: O00421) at the subcellular localization section. Indeed, in this section were defined both the nonmembrane region (topological domain), and the extent of the membrane‐spanning regions (transmembrane) of CCRL2.

All the CCRL2‐chemerin complexes were inserted in a membrane of 20% cholesterol and 80% POPC (1‐palmitoyl‐2‐oleoyl‐sn‐glycero‐3‐phosphocholine); the system was solvated with TIP3P water model and ionized up to a concentration of 0.15 M NaCl still using the CHARMM‐GUI, additional Cl− ions were added to neutralize the systems. All the system (proteins + membrane + solvent) consists of 67 447 atoms. Each system was then submitted to aMD, carried out on Cineca supercomputer using Amber20. The whole system was minimized (5000 cycle) using restraints for CCRL2 and membrane (10 and 2.5, respectively); then, the CHARMM‐steps equilibration protocol with progressive removal of position restraints was applied to the membrane and protein atoms (http://www.charmm-gui.org/demo/amber_ff/2). This equilibration protocol was carried out by Amber and consists of two NVT (constant number of particles (N), volume (V), and temperature (T)) steps to heat the system to 303.15 K employing as thermostat Langevin dynamics (collision frequency 1 ps) and four NPT (constant number of particles (N), pressure (P), and temperature (T)) steps (125 ps each) with SHAKE algorithm and the particle mesh Ewald (PME) (with a cutoff of 9 Å). The required average dihedral energy and average total potential energy were computed during 5 ns classical molecular dynamics for each studied complex. The aMD production (500 ns) was conducted at 315 K with constant pressure (1 bar) and periodic boundary condition, Shake (ntc = 2) and PME with cut of 10 Å were set, each simulation was repeated three times.

Total Accessible Surface Area and Buried Surface Area (BSA) were computed by Pisa server (http://pdbe.org/pisa/). The property maps were calculated by Coco server.

Trajectories analyses

Trajectories analyses were carried out by mdtraj. The PCA analyses was carried out with scikit‐learn using the decomposition module. Scipy library was used to calculate Gaussian Kernel density estimation (KDE). Graphics were done with Matplotlib.

RESULTS AND DISCUSSION

Modeling of CCRL2 and chemerin

To identify the putative residues involved in the CCRL2 chemerin binding, it was followed a protocol based on protein–protein docking and aMDs.

Since the structures of CCRL2 and chemerin were not experimentally resolved, it was decided to follow a homology modeling approach to obtain CCRL2, and an ab‐initio computations for chemerin.

CCRL2, similarly to other chemokine receptors had two conformational states: active and inactive. Engagement by the ligands turns GPCRs in the active state; therefore, it was decided to model only the active state of the receptor. In details, the CCRL2 model was based on the most similar chemokine receptor (see methods) and was devoid of N‐terminal tail.

On the other hand, chemerin was modeled ab‐initio due to the lack of highly conserved homologous proteins.

To be mentioned, meanwhile the designed computations were accomplished, both the structures of CCRL2 and chemerin became available at the AlphaFold database (alphafold.ebi.ac.uk). A comparison of the AlphaFold and our models was carried out by measuring the RMSD. For CCRL2, it was calculated a Cα RMSD of 1.02 Å and the great amount of this distance was related with the extracellular loop 2 (ECL2, residues 169–192) and the TM6 helix (Figure S1). TM6 was embedded in the membrane, far from the chemerin binding site. Therefore, we assumed that it would only have a marginal effect on the ligand binding. For the ECL2, it is challenging to reliable predict a long loop (23 residues) and also AlphaFold listed this loop as at low confidence (per‐residue confidence score between 70 and 50). Furthermore, the implementation of aMD instead classical MD reduced the bias related with the different loop conformations. Indeed, aMD offered a great advantage in modeling conformational change and to simulate infrequent events required for protein conformational change without previous knowledge of conformational states.

For chemerin, the superimposition of our model and the AlphaFold proposed led to Cα RMSD of 1.12 Å. The less fitting domain was the C terminal helix 2 (Figure S2). This region was reported to be not involved in chemerin binding to the CCRL2. In general, it was observed a good superimposition between the AlphaFold and our in‐house models.

CCRL2‐chemerin protein–protein docking

The docking computations were carried out by Haddock providing 12 different clusters. A representative binding conformation for each cluster, named complexes 1–12, was selected by Haddock post docking quality assessment tools (Figure S3). Further refinement of the selected models were not carried out. Similarly, to be as much unbiased as possible were not taken into account the docking score and the energy of the complexes.

Each selected complex was embedded into the membrane bilayer and submitted to aMD (500 ns). The obtained trajectories were first analyzed by RMSD to evaluate the system stability. For complex 9, during the simulation time the C‐terminal helix moved up to the binding site. Given that, it is accepted that the Chemerin C‐terminal binds CMKLR this complex was not further considered. Also, for complex 10, it was observed a dramatic change of chemerin conformation, with the C‐terminal that moved far from the N‐terminal and also this complex was not further studied.

Selection of CCRL2‐chemerin binding models

All the other trajectories were analyzed by PCA and the obtained matrixes were investigated by Gaussian kernel (KDE) to create a probability distribution functions in subspaces spanned by principal components 1 and 2 (PC1 vs. PC2). The approach that combined a dimension reduction step (PCA) with subsequent clustering (KDE) to analyze MD trajectories data were shown to be capable of reducing the noise and to generate more compact and well separated clusters of conformations.

Thus for each trajectory, the KDE plots allowed the identification of the higher populated conformational basins, and for each of them, it was extracted the representative conformations. This approach provided 23 highly frequent conformations assumed to be as the most relevant (Tables S1 and S3–S22).

For all these conformations, the BSA was computed. The conformations with BSA lower than 600 Å2 was rejected as consequence of the small size of the interface between CCRL2 and chemerin (Table S2). The 22 remaining conformations were analyzed by visual inspection focusing on the salt bridge interactions. This type of contact has a considerable contribution to the specificity of interaction of proteins with other biomolecules. Indeed, the energetic penalty due of dehydration of polar groups was paid off by favorable energy of salt bridge formation limiting the number of conformations of a molecule or complex, thus playing a crucial role in determining specificity. ,

By visual inspection, the studied conformations have been grouped into two general chemerin binding modes and it was also possible to identify the CCRL2 and Chemerin regions more often involved in the binding. For CCRL2, the two extracellular loops ECL2 (residues 169–192) and ECL3 (residues 264–270), and the residue lining the entrance of the receptor channel. For Chemerin, the three regions mainly involved in the binding with the cognate receptor were the α1 helix, the β1 sheet, and the loop between β2 and β3 helixes (β2β3‐loop residues 49–73).

The first binding mode (defined BM1) was shared by 12 of the 22 inspected conformations. This binding mode was featured by the contacts between chemerin β2β3‐loop with ECL3 (6 conformations of 12) and with ECL2 (6 conformations of 12). Furthermore, the chemerin α1 helix contacted the entrance of the channel (9 conformations of 12).

For BM2, shared by seven of the studied conformations, the chemerin β2β3‐loop contacted both the CCRL2 ECL2 and ECL3 (seven conformations), the α1 helix interacted with the CCRL2 ECL2 (seven conformations), and the β1 sheet had contacts with both the ECL3 and the residues lining the entrance of the receptor channel (four conformations). The three remaining conformations were featured by the significative involvement of the Chemerin C‐terminal domain in the binding to CCRL2. Since it was reported that the C‐terminal was only involved in the binding of the CMKLR1, these three conformations were rejected.

Worthily, the main differences between the two binding modes, BM1 and BM2, was a 180° rotation of the chemerin conformation. For the BM1, the chemerin α1 helix was located behind the β sheets, in contrast to the BM2 where the α1 helix was located in front of the β sheets (Figure 1).

FIGURE 1

Docking proposed chemerin binding modes

Proposed interaction models for CCRL2‐chemerin binding

To gain more insight, the residues involved in the binding was analyzed the types and the frequencies of the observed interactions. Actually, for the BM1 it was observed two patterns of interactions. For the first one, we had that the chemerin β2β3 loop established contacts with the residues of CCRL2 ECL2. The residues of the chemerin β2β3 loop were mostly polar and the most frequently observed interactions were salt bridges and H‐bonds. Indeed, we found a conserved array of polar contacts (6 conformation of 12) Lys60chem with Asp271CCRL2, Lys61chem with Glu265CCRL2, Glu63chem with Lys197CCRL2, and Lys72chem with Asp176CCRL2. It was also observed hydrophobic interaction between Val66chem and Phe188CCRL2 (Figure 2 and Figure S4).

FIGURE 2

BM1 first proposed pattern of interactions

FIGURE 4

Proposed interactions for BM2

The second pattern of interactions, for the conformation falling within BM1, consisted of the chemerin α1 helix residue Glu1, and Arg4 involved in salt bridge with CCRL2 Lys30 and Glu175, respectively. Also, chemerin Arg5 had polar contact with Glu26 or Asp29 of CCRL2. Worthily, also the residues of the chemerin β1 sheet were involved in interactions with the CCRL2 ECL2 and a polar contact between Glu26chem and Arg185CCRL2 was observed. Another polar interaction was seen between the chemerin β2β3 loop Lys61 and Glu192 of the CCRL2 ECL2 (Figure 3 and Figure S5).

FIGURE 3

BM1 second proposed pattern of interactions

Thus, the analyses of the BM1 conformations highlighted two main positions named as first and second pattern of interactions. Despite during the simulations time, we did not observe the shifting of one position to the other and we speculated that the chemerin β2β3‐loop might interact with the CCRL2 TM6‐TM7 loop, moving the latter far from the CCRL2 entrance channel enabling the chemerin α1 helix to move toward this channel.

For the BM2, we had that the chemerin β2β3‐loop formed extensive polar interactions and hydrophobic contacts. Indeed, the chemerin residues Lys60, Lys65, Arg67, and Lys72 established salt bridge with Glu175 of ECL2, Asp32 and Glu26 of TM1, and Asp271 of ECL3, respectively (five conformations of seven).

Worthily, it seemed that interactions between chemerin β2β3‐loop and the CCRL2 ECL2, forced the latter farm from the receptor entrance channel creating a space filled by β1 sheet residues (QETSV) doing a salt bridge between Glu322chem and Arg161ECL2 and hydrophobic contact between Gln321chem and Phe159EL2 (Figures 4 and S6).

CONCLUSION

The accomplished computations led us to gain more insight in the chemerin binding to CCRL2. A total of 5.5 μs simulations turned back with two binding modes for chemerin, both BMs suggesting a crucial role for the chemerin α1 helix, the β1 sheet and for the β2β3‐loop. It was also postulated that the CCRL2 chemerin complex formation might be dependent by the shift of the CCRL2 ECL2 far from the receptor entrance channel, driven by chemerin approach, lastly facilitating the binding. Furthermore, the analyses of the trajectories produced a short list of hotspot residues that might be crucial in favoring the complex formation and the chemotactic activity. Indeed, we identify for chemerin the α1 helix Glu1, Arg4, and Arg5, at the β2β3‐loop three lysine residues (60, 61, and 65), and for the β1 sheet Gln25 and Glu26. Also, for CCRL2, two regions were highlighted: the ECL2 and the ECL3. For ECL3, a crucial role seemed to be played by Glu175, Asp176, and Asp271 residues. The reported data represent the earliest attempt to shed light to the CCRL2 chemerin interaction. Although these results still need to be experimentally validated, they might help in better clarify CCRL2‐chemerin interaction. Furthermore, the proposed models might pave the way for medicinal chemistry efforts in search for modulators of CCRL2 chemerin interaction and help to better clarify the physio‐pathological role of both the CCRL2 and the chemerin and their potential value as target for therapeutic intervention.

CONFLICT OF INTEREST

The authors declare no competing interests.

Appendix S1 Supporting Information.

Click here for additional data file.