Molecular basis of the potential interaction of SARS-CoV-2 spike protein to CD147 in COVID-19 associated-lymphopenia.

Lymphopenia is considered one of the most characteristic clinical features of the coronavirus disease 2019 (COVID-19). SARS-CoV-2 infects host cells via the interaction of its spike protein with the human angiotensin-converting enzyme 2 (hACE2) receptor. Since T lymphocytes display a very low expression level of hACE2, a novel receptor might be involved in the entry of SARS-CoV-2 into T cells. The transmembrane glycoprotein CD147 is highly expressed by activated T lymphocytes, and was recently proposed as a probable route for SARS-CoV-2 invasion. To understand the molecular basis of the potential interaction of SARS-CoV-2 to CD147, we have investigated the binding of the viral spike protein to this receptor in-silico. The results showed that this binding is dominated by electrostatic interactions involving residues Arg403, Asn481, and the backbone of Gly502. The overall binding arrangement shows the CD147 C-terminal domain interacting with the spike external subdomain in the grove between the short antiparallel β strands, β1' and β2', and the small helix α1'. This proposed interaction was further confirmed using MD simulation and binding free energy calculation. These data contribute to a better understanding of the mechanism of infection of SARS-CoV-2 to T lymphocytes and could provide valuable insights for the rational design of adjuvant treatment for COVID-19. Communicated by Ramaswamy H. Sarma.

Introduction

In December 2019, an outbreak of a pneumonia-causing virus was first reported in Wuhan, China (Lu, Stratton, et al., 2020; Zhou et al., 2020), and has been spreading globally ever since (Sohrabi et al., 2020). The virus was found to phylogenetically belong to the severe acute respiratory syndrome coronavirus (SARS-CoV) (Lu, Zhao, et al., 2020). The virus was named SARS-CoV-2 (C. S. G. of the International, 2020; Gorbalenya, 2020), and the disease named coronavirus disease-19 (COVID-19). The entry of the virus into host cells is facilitated by binding of its transmembrane spike (S) protein with angiotensin-converting enzyme 2 (ACE-2) receptor (Hoffmann et al., 2020). The spike protein is a class I viral fusion transmembrane glycoprotein that plays a dual role in viral entry and fusion with host cell membrane (Bosch et al., 2003; Li, 2016). Cryo-electron microscopy (cryo-EM) revealed the structure of S protein as trimeric spikes composed of three S1 and S2 heterodimer subunits, protruding out of the viral envelope (Gui et al., 2017; Wrapp et al., 2020). The S1 subunit facilitates viral binding to cell receptor through the receptor-binding domain (RBD), whereas the S2 subunit encloses a hydrophobic fusion domain, allowing viral fusion (Song et al., 2018). Cleaving the S1-S2 boundary is primed by host cell proteases such as transmembrane protease serine 2 (TMPRSS2), and is required for viral activation and entry (Hoffmann et al., 2020). TMPRSS2 and ACE-2 are highly expressed in the epithelial tissue of multiple organs and thus significantly contribute to viral pathogenicity, clinical manifestations, and mortality (Sungnak et al., 2020).

One of the clinical indicators of the severity of SARS-CoV-2 is low lymphocytic count (lymphopenia) (Tan et al., 2020), especially in the CD8+, CD4+, and CD3+ T lymphocyte populations (Zeng, 2020; Zheng et al., 2020). Lymphopenia is a risk factor in viral infections, and was reported to predict the severity of the disease in COVID-19 patients (Gui et al., 2017). Although several studies examined the causes of lymphopenia during viral infections (Lalueza et al., 2017, 2019), the mechanisms underlying this feature are not fully understood in COVID-19 patients. Proposed mechanisms for the observed lymphopenia include direct invasion of lymphocytes (Xu et al., 2020), hyperlactic acidemia causing suppression of lymphocyte proliferation (Fischer et al., 2007), as well as an inflammatory cytokine storm resulting in lymphocyte apoptosis (Liao et al., 2002), similar to that reported earlier with MERS-CoV infection (Chu et al., 2016).

Angiotensin-converting enzyme 2, or ACE-2 receptor, has been shown to provide an entry point to SARS-CoV-2 in different types of human cells (Li et al., 2020; Zhou et al., 2020). However, T lymphocytes were shown to be consistently negative for ACE-2 receptors (Hamming et al., 2004), suggesting an alternate route for viral entry. CD147 was recently proposed as a probable route for SARS-CoV-2 invasion into Vero E6 host cells (Wang, 2020). CD147 (also known as Basigin or EMMPRIN) is a transmembrane glycoprotein that was originally identified in 1992 as a T lymphocyte activation-associated antigen (Kasinrerk et al., 1992). Further research confirmed its expression on activated T lymphocytes and its relatively weaker expression in resting T lymphocytes, and that its expression increases rapidly upon T lymphocyte activation (Koch et al., 1999). CD147 induces matrix metalloproteinase (MMP) production during physiological and pathological events including inflammatory response, wound healing, and tissue homeostasis (Cruzat et al., 2018; Jin et al., 2019; Kato et al., 2009; Ungern-Sternberg et al., 2018). Interestingly, CD147 receptor is expressed on red blood cells (RBCs), providing a route for malaria entry (Zhang et al., 2018). It is worth mentioning that PfRH5 complex on the cell membrane of Plasmodium falciparum, the causative agent of malaria, interacts with CD147 on human RBCs, a binding which was necessary for invasion by all strains when tested in vitro (Crosnier et al., 2011). Anti-CD147 antibodies were shown to prevent entry of P. falciparum and SARS-CoV-2 into host cells (Bian, 2020; Zhang et al., 2018).

To understand the mechanism of interaction of the SARS-CoV-2 spike protein with the CD147 receptor, we have performed a four-stage in-silico study. First, the binding site of both proteins had to be accurately predicted, taking into consideration the few number of their available crystal structures. Second, molecular docking was performed using three different protein-protein interaction servers. Then, the best docking pose was analyzed in MD simulation to check its stability. Finally, the binding free energy was estimated, using snapshots of the MD simulation, to get in-depth knowledge of the binding determinants.

Materials and methods

The study starts with determination of the most probable binding sites for spike and CD147. This was specially challenging for the latter because only few crystal structures were reported for this membrane glycoprotein. Later, the molecular docking was performed using three different web servers and the best pose was selected using the PDBePISA website. Then, the best complex was subjected to a 100 ns MD simulation in NPT ensemble. Finally, free energy calculation was performed to gain insights into the major contributing forces in binding and the critical residues involved.

Binding site determination

The SARS-CoV-2 receptor binding domain (RBD) was obtained from its crystal structure complexed with human angiotensin-converting enzyme 2 (hACE2) (PDB ID: 6LZG), while that of CD147 was isolated from its crystal structure complexed with malaria invasion protein, reticulocyte-binding protein homologue (RH5) attained from the parasite Plasmodium falciparum (PDB ID: 4U0Q). The crystal structures of both targeted proteins, SARS-CoV-2 and CD147, were closely examined for integrity and the active residues taking part in the interaction were identified. The surface of CD147 was examined using the Computed Atlas of Surface Topography of Proteins (CASTp) server. The default probe radius used by CASTp is 1.4 Å. We opted to use a smaller probe of 0.8 Å to get more exhaustive coverage of all the potential binding sites on the CD147 surface.

Molecular docking

Prior to the docking simulation, both proteins, spike and CD147, were checked for errors and prepared in the Protein Preparation Wizard, within MOE software. The preparation steps include removal of any bound ligands, water of crystallization, and salts. Also, all the hydrogens were added and the protein termini were capped. Both proteins were minimized individually in MOE to a gradient of 0.001 kcal/mol.Å2. First, all heavy atoms were kept fixed, and the hydrogen atoms were allowed to optimize their positions. Second, the entire backbone of the protein was kept fixed and the side chains were minimized using the same parameters. Finally, the whole molecule including the backbone was minimized (Helal et al., 2011; Helal & Avery, 2012). Three different web-docking servers, HADDOCK v2.2, ZDOCK v3.0.2, and HawkDock 2018, were utilized to computationally predict the protein-protein interactions of the two unbound protein structures. The docking trial on each server was repeated twice using the default settings; in the first run no restraints was given to the server so that it can freely determine the binding region (interface) for each protein. This step was performed to find out to which extent the software would bias toward a specific docking pose. In the second run on each server, specific active residues were defined for the software to be considered as restraints (see the discussion part). These residues were selected based on our understanding of the crystal structures of both SARS-CoV-2-RBD/hACE2 and CD147 RH5 as well as the prediction by CASTp server.

Molecular dynamics simulation

Molecular dynamics (MD) simulation was performed using GROMACS 2019 software package and the all-atom OPLS force field. The protein complex was solvated in a cubic box of SPC water model of 100 × 100 × 100 Å size and Periodic boundary conditions were implemented. Ionizable residues were assigned standard ionization states at PH 7. The system was then neutralized by the addition of five sodium ions. The total number of atoms in the system, including the water and ions, was 117,013. Short-range interactions were treated with a 10 Å cutoff for both Lennard-Jones and Coulomb potentials, while the Particle Mesh Ewald (PME) algorithm was used for long-range electrostatics. The MD simulation was performed in four steps. During the first three steps, a force constant of 1000 kJ mol−1 nm−2 was used to restrain all the heavy atoms to preserve the original fold of the proteins. First, the solvated protein system was minimized using a 5000 steps of steepest descent algorithm to resolve any steric clashes or inappropriate geometry. Then, to ensure a reasonable starting structure, the system was equilibrated under constant Number of particles, Volume, and Temperature (NVT) ensemble for 100 ps using Berendsen thermostat (Golo & Shaĭtan, 2002). The second round of equilibration was performed under constant pressure (NPT ensemble) using the Parrinello-Rahman barostat for an additional 100 ps (Tuble et al., 2004). Finally, the position restraints were released and the system was simulated in the production run under an NPT ensemble (Nosé‐ Hoover thermostat and Parrinello‐Rahman barostat) for 100 ns using a time step of 2 fs. Simulation results were analyzed using the Visual molecular dynamics (VMD) software, ver. 1.9.3 (Humphrey et al., 1996). Pymol graphical software ver. 2.3 was used for the generation of the figures of protein-protein interactions.

Binding free energy calculation

Binding free energy was calculated using the Molecular Mechanics Poisson Boltzmann Surface Area (MM/PBSA) method with the g_mmpbsa script. This method calculates the binding energy using molecular mechanics potential energy (electrostatic and van der Waals interactions) and free energy of solvation (polar and nonpolar solvation energies). A total number of 100 snapshots, extracted at regular intervals from the whole MD trajectory, were used for calculating each energy term. The average for each term was then calculated. In addition, a per-residue binding free energy decomposition at each snapshot was performed and the average binding energy for each residue was calculated and depicted. The solvent-accessible surface area (SASA) model was used to compute non-polar solvation energy with a surface tension of 0.02267 (kJ/mol2) and a probe radius of 1.4 Å.

Results and discussion

The recently reported crystal structure of SARS-Cov-2 spike protein complex with ACE2 (PDB ID: 6LZG) reveals that the virus utilizes the external subdomain of the spike Receptor Binding Domain (RBD) to recognize the human ACE2 receptor (Wang, Zhang, et al., 2020). This external subdomain is composed mainly of a flexible loop, S438-Y508, connecting two anti-parallel small β strands. Also, a part of this loop, the segment F486-G502, represents the most solvent accessible part of the spike RBD and makes direct contact with ACE2 in the crystal structure. In addition, in a recent in-silico study, it was proposed that a small circular region of this loop, C480-C488, could be responsible for the interaction with the human GRP78 (Ibrahim et al., 2020). Collectively, these findings have convinced us to use the solvent-accessible F486-G502 segment as the interaction site of the spike protein in our docking simulation. Regarding the CD147 receptor, the only crystal structure available of this protein, bound to another biological molecule, is for its complex with the PfRH5 protein (PDB ID: 4U0Q) (Wright et al., 2014). In this structure, the PfRH5 binds to the cleft between the two domains of CD147 assisted by a network of H-bond and salt bridge interactions. Notably, the glycoprotein CD147 is a relatively small protein and few pockets could be noticed on its surface. To get an accurate prediction of the potential binding site on CD147, we have performed a comprehensive search for its surface pockets using the Computed Atlas of Surface Topography of Proteins (CASTp) server (Tian et al., 2018). This tool enables an exhaustive and quantitative characterization of the protein topographic features, providing a complete list of the available pockets with the amino acid residues lining these pockets. Figure 1(A) shows the pockets with significant volume (>15 Å3) on CD147 surface. Among them, pockets 1, 3, 4, and 6 are the most solvent accessible and constitute a nearly continuous grove extending at the interface of the two domains of CD147. Furthermore, according to the crystal structure 4U0Q, the PfRH5 protein uses pockets 1 (124.67 Å3) and 6 (72.57 Å3) for interaction with CD147. In addition to the CASTp server, we have also used the CPROT tool (Bonvinlab) which is an algorithm for the prediction of protein-protein interface merging six calculation methods into a consensus result (de Vries & Bonvin, 2011). Interestingly, the CPROT results were very similar to those of CASTp with most residues of pockets 1 and 6 identified (Figure 1(B)). Therefore, we have selected the residues of these pockets, 1 and 6, as the active residues in our docking study, namely Ser78, Asp79, Gln81, Trp82, Gly103, Pro104, Pro105, Arg106, Val 107, Lys108, Lys127, Ser128, Glu129, Ser193, and Asp194. (Table S1, Supporting Information).

Figure 1.

(A) Pockets detected on the surface of the CD147 receptor using the CASTp server. The protein is shown as a grey cartoon and the pockets are displayed as colored spheres. (B) Interacting residues of the CD147 as predicted by the CPROT tool. Predicted binding residues and surrounding residues are shown as red and green lines, respectively. Residues not participating in the interaction are shown as blue lines.

Molecular docking was performed using three different protein-protein docking servers, namely HADDOCK v2.2, ZDOCK v3.0.2, and HawkDock 2018 (Dominguez et al., 2003; Pierce et al., 2014; Weng et al., 2019). We specified the active, i.e. interacting, residues as the segment F486-G502 of the spike RBD and the residues of pockets 1 and 6 of the CD147 as illustrated in the previous section. In order to evaluate the docking results of the three software servers, the top ranking pose from each was examined using the PDBePISA online tool, developed by The European Bioinformatics Institute, EMBL-EBI, for the exploration of macromolecular interfaces (Krissinel & Henrick, 2007). Moreover, a Molecular Mechanics energies combined with Generalized Born and Surface Area (MM/GBSA) calculation, as implemented within the Hawkdock server, was used to estimate the binding affinities (Hou et al., 2011; Zhong et al., 2020). The results of the examination of the top complex generated by each software are illustrated in Table 1. Interestingly, all the three servers returned very similar arrangements of the interacting proteins with RMSD of 2.23 for ZDOCK/HADDOCK and 0.24 for HawkDOCK/HADDOCK. Moreover, for each server, performing the docking run without specifying any constraints returned quite a similar docking pose, to the one with constraints, within the top ten out of more than one hundred solutions (data not shown). Based on the RMSD values and visual examination, it was noticed that the poses proposed by HADDOCK and HawkDOCK are very close and almost converging. We decided to use the pose generated by the HADDOCK server as it shows a higher interface surface area of 1,069 Å2 and the best score of −68.00 in the MM/GBSA calculations (Table 1).

Table 1.

Parameters of the interface of the interaction between the spike RBD and CD147 as predicted by the PDBePISA server and the Hawkdock MM/GBSA calculations.

Docking Server	Spike		CD147		Interface
	Interface Residues	aInterface Surface, Å2	Interface Residues	Interface Surface, Å2	bInterface Surface, Å2	cΔG P-Value	MM/GBSA kcal/mol
HADDOCK	32	10,339	34	10,130	1,069	0.297	−68.00
ZDOCK	32	10,292	39	10,500	1,154	0.367	−66.77
HawkDock	36	10,259	29	10,284	1,042	0.181	−59.26

Total solvent accessible surface area in square angstroms for each protein.

Interface area, calculated as difference in total accessible surface areas of isolated and interfacing structures divided by two.

Δ indicates the P-value of the observed solvation free energy gain. The P-value measures the probability of getting a lower than observed ΔiG, when the interface atoms are picked randomly from the protein surface. indicates interfaces with surprising (higher than would-be-average for given structures) hydrophobicity, implying that the interface surface can be interaction-specific.

In addition, thorough investigations were carried out using PDBePISA tool to determine the key interacting residues in the selected docked complex generated by the HADDOCK server. Only interfacing residues within 4 Å van dar Waals distance cutoff between the SARS-CoV-2 RBD and CD147 were examined by the PDBePISA algorithm (Table S3, Supporting Information). It can be easily noted that the docked SARS-CoV-2 RBD adopted a similar orientation to that observed in its co-crystal structure with hACE2, which highlights the binding role of the region that extends between the residues C480 and Q505 (Figure 2). The spike external subdomain forms a surface to which parts of CD147 domains dock. This surface is dominated by the interaction of the C-terminal domain of CD147 which binds to the spike external subdomain in the grove between the short antiparallel β strands β1’ and β2’ and the small helix α1’ (Figure 2(A)). This grove makes contacts with nearly 50% of the C and F strands of CD147. On the other domain of CD147, parts of the A and G strands forms a complementary surface to the cyclic highly solvent-accessible segment C480-C488 of the spike RBD.

Figure 2.

Proposed binding mode between the RBD of the SARS-CoV-2 spike and CD147 as predicted by the HADDOCK server v2.2. (A) Side-view of the interaction between the two proteins. Key structural elements participating in the interaction are labeled. (B) Top-view of the proposed interaction. (C) Surface representation of the proposed complex. Spike RBD is shown in cyan color, while CD147 is displayed as an orange surface.

In the predicted complex, a row of hydrophilic residues along the interface formed strong and intricate polar contacts (Figure 3). The most important interacting residues from both sides, together with the distances between each pair, are listed in Table 2. Strong polar contacts were found to involve the spike RBD residues Gln493 and Ser494 with the side chain of Thr135 in CD147. Also, the side chain nitrogen of the spike Gln493 makes a 3.00 Å H bond with the hydroxy group of Ser190 in CD147. In addition, the backbone carbonyl of spike Glu484 formed a strong H bond (2.5 Å) with the backbone NH of Gln100 in CD147. Moreover, the side chain NH2 of Asn481 on the cyclic segment of the spike β1’/β2’ loop formed a bidentate H-bonds with the backbone of Phe27 and side chain OH of Thr28 in CD147. Another set of adjacent hydrophilic residues located on the spike RBD external domain (Gln498, Thr500, Asn501, Gly502 and Tyr505) were found to contribute a network of H-bonds with CD147 amino acids Trp137, Thr140, Ser163 and Thr188, which further augments the two proteins interaction. Furthermore, two strong salt bridges were observed to stabilize the binding interface at both ends of the spike F486-G502 segment. The first one was noticed between Arg403 and Asp136 of the CD147 and the other was a stronger interaction (2.7 Å) emerging from spike Glu484 towards CD147 Lys191. The latter was acting as a centrally located anchor point that might have a key role in stabilizing the complex. It is worth noting that the same residue, Lys191, is an important player in the interaction of CD147 with the P. falciparum protein (PfRH5) (Wright et al., 2014). Notably, fewer hydrophobic interactions were noticed at the interface as shown in Table 2. The most important of these interactions was the one between the spike Tyr489 and His102 on the linker between the two domains of CD147.

Figure 3.

Polar interactions at the interface as predicted by the HADDOCK server v2.2. Key interacting residues are shown as sticks, while the proteins backbone is depicted as a transparent cartoon.

Table 2.

The interacting residues from both proteins and the distances in Å between each interacting pair.

No.	Residue from Spike	Residue from CD147	Distance (Å)
Ionic Bonds
1	Arg403	Asp136	5.00
2	Glu484	Lys191	2.72
H-Bonds
1	Arg403	Thr135	3.10
2	Asn481	Thr28	3.50
3	Glu484	Gln100	3.28
4	Gly496	Asp136	2.62
5	Gln498	Thr188	2.60
6	Thr500	Tyr140	2.80
7	Asn501	Trp137	2.46
8	Gly502	Trp137	3.21
9	Tyr505	Ser163	2.42
Hydrophobic Interactions
1	Tyr489	His102	5.50
2	Leu455	Pro133	4.00
3	Val483	Thr28	3.30

In addition, to further validate the proposed binding mode, we have docked the spike glycoprotein of SARS-CoV, the causative agent of the 2003 epidemic (PDB ID: 5XLR) to the target protein CD147. This protein is known to bind with the human cell CD147 to facilitate entry of the virus responsible for the 2003 outbreak. We used HADDOCK server and the same parameters exploited for docking of the SARS-CoV-2 spike. To our delight, this protein adopted a very similar binding pose to that of the SARS-CoV-2 predicted above, with a slight shift towards the C-domain of the CD147 (Figure S2, Supp. Info.). Also, CD147 protein made polar contacts to the SARS-CoV protein using the same residues, namely, Gly100, His102, and Lys191. This gives added validity to our proposed binding pose of the SARS-CoV-2 spike and CD147.

The complex was then subjected to a Molecular Dynamics (MD) simulation study in order to check the stability of the proposed binding mode and to study the time-dependent behavior of this system. In absence of 3D structure information, MD simulations is especially useful for studying protein-protein complexes as it can help the molecules “explore” their conformation space more efficiently than using the static “image” provided by the molecular mechanics energy minimization (Karplus & Petsko, 1990; Zhao et al., 2013). The system was solvated in a cubic SPC water box and simulated for 100 ns in GROMACS software package under NPT ensemble using Nosé-Hoover thermostat (Van Der Spoel et al., 2005). The final frame was extracted, minimized to a gradient of 0.001 kcal.Å−2.mol−1 in MOE software, and was further analyzed for key structural changes (Vilar et al., 2008). To our delight, the overall fold of both interacting proteins was preserved as well as most of the binding interface with an RMSD of just 1.25 Å. The spike RBD has moved slightly towards the CD147 N-terminal domain, while most of the polar interactions were retained or replaced by similar interactions with nearby residues (Figure 4(A, B)). The backbone RMSD fluctuation of the whole complex as well as the interface residues were monitored during the MD simulation to check the stability of the system (Figure 4(C)). Both proteins were stable with some RMSD fluctuations that stabilized after 40 ns. To our delight, the residues of CD147 at the interface attained a significant stability over time and reached 1 Å after 60 ns. Distances of the most important polar contacts at the interface were calculated throughout the simulation and represented as a heatmap in Figure 4(D). It was noticed that the salt bridge between Arg403 of the spike and Asp136 of CD147 was stable during the simulation with only minor fluctuation for a short time between 20 and 36 ns. This ionic bond could be critical for anchoring the β-strand B of the CD147 C-terminal domain to the binding grove created by the spike RBD. Similarly, the H bond between Gln498 and Thr188 at the end of stand F of CD147, although weak, could contribute to some extent to the interaction of this strand to the surface of the spike. On the other side, at the beginning of the spike F486-G502 segment, the H bond of Asn481 to Thr28 was replaced by a stronger and more stable one with Thr44. The most prominent difference was the loss of the characteristic salt bridge between Glu484 and Lys191. In the starting complex, this interaction was centrally located below the loop connecting the two domains of CD147 and we initially thought it would be essential for stabilizing the interaction. As a compensation, both residues Glu484 and Lys191, upon liberation from the ionic bond, were able to form H bonds with Gln100 and Tyr449, respectively. Finally, we noticed that Tyr140 adopted a swing movement to bind with the backbone NH of Gly502 instead of Thr500, thus, stabilizing the strand D of CD147 to the edge of the spike F486-G502 segment.

Figure 4.

Molecular dynamics simulation of the proposed spike RBD-CD147 complex. (A) Alignment of the complexes before and after the MD simulation. (B) Polar interactions at the interface. Key interacting residues are shown as sticks, while the proteins backbone is depicted as a transparent cartoon. (C) Backbone RMSD fluctuation of the whole complex, spike residues at the interface, and CD147 residues at the interface during the MD simulation. (D) Heat map representing the change in the distances of the critical polar interactions at the interface during the whole MD simulation.

The residues proposed to be critical for the Spike-CD147 interaction were mutated to alanine to validate their important role in the interaction. These residues are: Arg403 and Asn481 in spike, and Tyr140 in CD147. The mutated proteins were docked using HADDOCK server with the same parameters used for docking of the wild-type proteins. As expected, the number of residues at the interface and the interface surface area were lower than those of the wild-type proteins (Table S4, Supp. Data). Also, the MM/GBSA score, as predicted by the Hawkdock server, was −60.64 which is lower than the scores of the wild-type protein complexes generated by any of the docking servers we used. Overall, the mutant proteins bind in a similar fashion to the wild-type one with an RMSD of 1.332 Å. However, the critical salt bridges and H bonds formed by the mutated residues were replaced by weaker contacts involving Phe486, Tyr489, Gln493, and Asn501 from spike with Val160, Asp136, Ser109, and Gln100 from CD147, respectively (Figure S3, Supp. Data).

In addition, the complex with the mutated residues was energy minimized and then subjected to MD simulation for 30 ns following the same protocol used for the wild-type complex. The purpose of this simulation is to compare the time-dependent behavior of both proteins and validate our proposed binding mode. It was noticed that the H bonds Tyr489-Asp136 and Gln493-Ser109 were maintained throughout the simulation. However, the complex was significantly less stable compared to the wild-type complex. This was evident from the higher fluctuation of the backbone RMSD of the whole complex and, especially, the interface residues (Figure S4, Supp. Data). These findings highlight the role of the mutated polar residues, Arg403 and Asn481 in spike, and Tyr140 in CD147, for the interaction of these proteins and give more confidence to our proposed binding mode.

To get more insight into the nature of interaction between the two proteins and detailed information about individual ligand contribution, we have performed binding free energy calculation (Cavasotto, 2020). We have used the MD-based method Molecular Mechanics Poisson Boltzmann Surface Area (MM/PBSA) as implemented in the g_mmpbsa tool. This technique can achieve similar accuracy to the Free Energy Perturbation (FEP) methods but at a smaller computational cost (Kumari et al., 2014). The calculation was performed with the single trajectory approach using 100 frames saved at regular intervals of the complex MD simulation. According to the MM/PBSA method, the binding free energy of the interacting proteins is calculated using the following equations: where G is the total free energy of the complex and G is the total free energy of the isolated proteins. ⟨EMM⟩ is the vacuum molecular mechanics potential energy. TS refers to the entropic contribution to the free energy where T and S denote the temperature and entropy, respectively. The last term ⟨Gsolvation⟩ is the free energy of solvation and is composed of the polar part Gpolar and the nonpolar part Gnon-polar; the former could be estimated by solving the Poisson − Boltzmann (PB) equation and the latter could be calculated from the solvent-accessible surface area (SASA). The term γ is a coefficient related to surface tension of the solvent and b is fitting constant (Kumari et al., 2014). The results of the binding free energy calculation are shown in Table 3. In accordance with our initial prediction, the binding energy is dominated by the electrostatic term with a contribution of 66% compared to 30% of the van der Waals interaction. Therefore, the attractive electrostatic term greatly exceeded the repulsive term and was the main driving force for the binding. The MM/PBSA method also returned a Per-residue binding-free energy decomposition (Figure S5, Supp. Data). This helps in identifying the key residues involved in the interaction. In our case, a contribution of more than 5 kJ/mol to the binding-free energy of the complex was considered significant. As expected, most of the residues participated in the binding of the initial and MD-refined complexes showed a significant contribution. Particularly, the calculation highlighted the role of the spike Arg403, Tyr449, Asn481, Val483, and Tyr505 in the attachment to the CD147 receptor. The overall proposed binding mode of the host cell CD147 and the SARS-CoV-2 spike is illustrated in Figure 5.

Table 3.

Breakdown of the binding free energy as calculated by the g_mmpbsa tool.

Energy Term	Value, KJ/mol	Percent of Contribution
van der Waal	−243.107	30
Electrostatic	−534.870	66
Solvation (polar)	366.279	0
Solvation (SASA)	−31.926	4
Binding Energy	−443.623	100

Figure 5.

Proposed mechanism of the SARS-CoV-2 fusion to the lymphocytes as mediated by the interaction between the spike RBD and the human CD147.

Involvement of CD147 in SARS-CoV-2-assosciated lymphopenia

Lymphopenia was found to be associated with the severity of SARS-CoV-2 infection in many clinical reports (Guan et al., 2020; Liu, 2020; Liu et al., 2020; Wang, Yang, et al., 2020; Zhang et al., 2020). A pooled analysis of 2282 cases found that the occurrence of lymphopenia corresponds to a nearly 3-fold increased risk of severe disease (Zhao et al., 2020), supporting its potential usage as a clinical sign for serious morbidity (Tan et al., 2020). Although the cause of lymphopenia during viral infection still needs further investigations, several studies have postulated potential mechanisms for the evident lymphopenia. These include direct viral invasion into T lymphocytes resulting in cell death (Tan et al., 2020), indirect upregulation of P53 apoptotic pathways in T lymphocytes by excessive secretion of pro-inflammatory chemokines and cytokines (IL10, IL7, IL6, IL2, TNFα, CXCL10/IP-10, CCL3/MIP-1A CCL2/MCP-1, and CCL4/MIP1B) (Chen et al., 2020; Liao et al., 2002; Xiong et al., 2020), and hyperlactic-acidemia-induced apoptosis (Díaz et al., 2018). Notably, T lymphocytes barely express ACE-2 receptor (Hamming et al., 2004) - the currently accepted point of entry of SARS-CoV-2 (Li et al., 2020). CD147 may thus present an alternative point for viral spike protein entry. Indeed, CD147 was reported to mediate the entry of human immunodeficiency virus (HIV) into T lymphocytes (Pushkarsky et al., 2001). It was also reported to mediate entry of cytomegalovirus (CMV) into endothelial and epithelial cells (Vanarsdall et al., 2018), and Plasmodium falciparum into red blood cells (Zhang et al., 2018).

CD147, also termed as EMMPRIN or basigin, is a pleiotropic molecule with various effects (Grass & Toole, 2016; Iacono et al., 2007). Given the emerging use of bioinformatics in biological research (Karim et al., 2020; Mehmood et al., 2014), it is worth mentioning that, through pathway analyses, CD147 was found to be a significant cell-surface protein contributing to chemoresistance and cell survival in cancer-stem-cell-like cells (Kang et al., 2013), and the invasion of various cancers (Li et al., 2015; Liang et al., 2016; Zheng & Gong, 2017). Additionally, it was found to be highly expressed in activated regulatory T lymphocytes, enabling it to be a potential marker for distinguishing cytokine-producing activated T lymphocytes from resting T lymphocytes (Landskron & Taskén, 2013; Solstad et al., 2011). The structure of transmembrane CD147 receptor consists of three domains: extracellular, transmembrane, and a short cytoplasmic tail (Muramatsu, 2016). Those domains facilitate its interaction with several proteins, such as monocarboxylate transporters (MCTs) and cyclophilins (Cyps), especially cyclophilin A, to modulate specific functions (Kirk et al., 2000; Yurchenko et al., 2010). CypA is a ubiquitously expressed chaperon protein that has a high affinity to bind to CD147 receptor for enhancing viral invasion and inflammation (Chen et al., 2005; Yurchenko et al., 2010); a study in 2013 reported that blocking of CypA inhibits SARS-CoV infection and replication (Tanaka et al., 2013).

CD147 acts as the major chaperone protein of MCT 1, MCT 3 and MCT 4 proton transporters, that control the export of lactic acid molecules outside the cell (Kirk et al., 2000; Philp et al., 2003). Activated T lymphocytes have been reported to undergo increased glycolysis with increased lactate production (Cham & Gajewski, 2005), as they need high energy supplies for proliferation and producing cytokines (Frauwirth & Thompson, 2004). The expression of MCTs for lactate export helps T lymphocytes to cope with the increased lactic acid levels (Merezhinskaya et al., 2004), since the accumulation of lactic acid induces T cells apoptosis (Díaz et al., 2018). These data provide an insight into the possible impact of SARS-CoV-2 infection on the balance of CD147-MCT interaction, and direct a need for further research on this role.

Conclusion

Angiotensin-converting enzyme 2, ACE-2 receptor, has been proved to be the entry point of SARS-CoV-2 into the lungs and other types of human cells. However, T lymphocytes are rarely expressing this receptor. Several recent studies have suggested that CD147 is a probable route for SARS-CoV-2 invasion into certain host cells. These findings encouraged us to design the current study in which we employed docking and MD simulation to investigate this possibility. The results were also compared with the molecular modeling data of the mutant proteins interaction as well as the interaction of CD147 with the spike of the SARS-CoV, the causative agent of the 2003 outbreak.

Our model shows interaction between CD147 C-terminal domain with SARS-CoV-2 spike external subdomain in the grove between the short antiparallel β strands, β1’ and β2’, and the small helix α1’. We confirmed this proposition using MD simulation and binding free energy calculation. This model provides an appreciation for targeting CD147 for adjuvant drug therapy strategies, to improve the prognosis of COVID-19 pateint. Our data is supported by the results of using Meplazumab, an anti-CD147 humanized antibody, for viral clearance and lymphocyte count normalization (Bian, 2020). In addition, of proposing azithromycin as a therapeutic strategy to treat P. falciparum via interfacing in PfRH5/CD147 interaction (Muralidharan & Striepen, 2015; Wilson et al., 2015). The mechanisms proposed in our study encourage research on CD147 involvement in SARS-CoV-2 in direct invasion to T lymphocytes and associated lymphopenia.

Click here for additional data file.

Click here for additional data file.