Cannabinoid CB2 Receptor Functional Variation (Q63R) is Associated with COVID-19 Severity: from Human Study to Molecular Docking

Background: Evidence supports the role of host genetic diversity for the clinical course variation of coronavirus disease 2019 (COVID-19). Variation in the cannabinoid CB2 receptor gene (CNR2) could affect the endocannabinoids regulatory actions on the immune system, resulting in an increased risk of various inammatory diseases. The present study investigated the relationship between the CNR2 rs35761398 (Q63R) functional variation and COVID-19 severity. Results: A total of 200 Iranian COVID-19 patients (100 expired and 100 discharged) were enrolled in the study and genotyped through TaqMan assay. The co-dominant, dominant, recessive, over-dominant, and additive inheritance models were analyzed using SNPStats software. In silico molecular docking was also performed to simulate the effects of Q63R variation on CB2 binding with a ligand and with G-protein. A signicant difference in the Q63R allele and genotype distributions was found between COVID-19 expired and discharged patients in co-dominant (OR: 3.33, 95% CI: 1.25-8.88, p = 0.043), recessive (OR: 2.92, 95% CI: 1.16-7.33, p = 0.017), and additive inheritance (OR: 1.62, 95% CI: 1.06-2.48, p = 0.025) models. The molecular docking results showed that the predicted structure of mutant CB2 (63R type) could not bind to G-protein in the correct position. Conclusions: The data implied the involvement of the CNR2 gene in the severity of COVID-19 in Iranian patients. Identication of genes related to susceptibility and severity of COVID-19 may lead to specic targets for repurposing or drug development. EC; cannabinoid receptor 2, CB2; cannabinoid CB2 receptor gene, CNR2; reverse transcription-polymerase chain reaction, RT-PCR; Hardy-Weinberg equilibrium, HWE; 3-dimensional, 3D; protein data bank, PDB; odds ratios, OR; condence intervals, CI; Akaike information criterion, AIC; root-mean-square deviation, RMSD; angiotensin-converting enzyme 2, ACE2; during respiratory syncytial virus, RSV; G-protein-coupled receptors, GPCRs;


Background
Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) as a newly emerging virus causes mild-tosevere respiratory disease, which has been named coronavirus disease 2019 (COVID-19) (1,2). All people are susceptible to the virus infection, but there is considerable variation in the disease course and outcome among infected individuals (3). While many infected cases do not exhibit any symptoms, others proceed to develop COVID-19; however, severe illness and death occur only in a small minority of patients (4). Although our understanding of the SARS-CoV-2 and COVID-19 is still in its infancy, there is now strong evidence supporting the role of host genetic diversity alongside with other host, viral and environmental factors for the clinical course variation (5)(6)(7)(8)(9)(10)(11)(12). Host genetic diversity could dictate the clinical response to respiratory viruses through susceptibility to viral infection, and propensity to develop harmful pulmonary in ammation (13). Finding a relationship between host genetics and the clinical outcome of SARS-CoV-2 infection may be necessary for identifying high-risk individuals.
The endocannabinoid (EC) system is a biological system composed of endogenous cannabinoids and their respective receptors, CB1 and CB2 (14). The system has been identi ed as a critical endogenous regulator of immune system homeostasis due to its effects on immune cells development, migration, proliferation, and effector functions (15). Cannabinoids have been proposed as a promising immunomodulator to reduce the SARS-CoV-2 immunopathology (16). Variations in the cannabinoid CB2 receptor gene (CNR2) could affect intracellular signaling and reduce the ECs function, which has been associated with an unbalanced immune response and increased risk of various in ammatory diseases (17)(18)(19)(20)(21)(22)(23). The CB2-Q63R polymorphism is a missense mutation of the second and third bases at codon 63 of the CNR2 gene, which leads to a Q/R substitution, causing a different polarization state of the protein (24). This variation has been shown to affect the response of CB2 to cannabinoids and differently modulate the EC-induced inhibition of lymphocyte proliferation (25). While evidence indicated that the mutation does not affect receptor-ligand binding (24), the exact mechanism behind this action is still unknown.
Focusing on the immunopathogenesis of SARS-CoV-2 and EC effects on the immune system, here we describe how variability in the CNR2 gene can conceivably explain variability in COVID-19 clinical phenotype. Besides, in silico molecular docking was performed to simulate the effects of CB2-Q63R variation on receptorligand and receptor-G-protein interaction. The data implied the involvement of the CNR2 gene in the severity of COVID-19 in Iranian patients.

Human study
All patients were COVID-19 con rmed cases and hospitalized in the central hospital for COVID-19 patients in Gorgan city (Sayyad Medical and Educational Center). The details of the demographic, gender, age, and clinical data of all cases are presented in Table 1. The age distribution between expired patients (mean, 62.08 years) was signi cant compared with discharged patients (mean, 54.45 years) (p < 0.05). Moreover, the age distribution between female (mean, 63.12 years) and male (mean, 61.04 years) expired patients was signi cant compared with the discharged subjects (female mean, 55.64 years; male mean, 53.26 years) (p < 0.05). Of all patients enrolled in the study, the most observed symptoms were dyspnea (66.5%), cough (66%), fever or chill (59.9%), sore throat (20.5%), myalgia (16.5%), ageusia and anosmia (14.5%), nausea or vomiting (10.5%), diarrhea (7%), headache (6%), chest pain (5%), and fatigue (3.5%). The allelic frequencies and the genotype distributions in expired and discharged patients are shown in

Molecular docking
The I-TASSER server predicted ve models for the submitted sequences, and we selected the best predicted 3D structure of mutant CB2 (63R) with c-score, i.e., -0.9 with a 0.74 TM-score. C-score shows the quality of the predicted model by I-TASSER, and a model with a higher c-score signi es a higher con dence predicted model. TM-score shows the similarity between the predicted model and the template (6PT0), and the TM-score of more than 0.5 exhibits the correct topology of the predicted model. The best model of homology with 100% con dence and 94% identity were selected. The 3D structure of the predicted 63R models and wild type CB2 (PDB ID 6PT0, Chain R) are shown in Figure 1. Ramachandran plots are shown in Figure 2.
Molecular docking was performed to examine the binding a nity and possible binding pocket of 2-AG with both the built CB2-63R models and wild-type CB2-Q63. The binding energy between 2-AG and CB2-Q63, 2-AG and I-TASSER built model, and 2-AG and Phyre2 built model were -7.5, -6.9, and -7.0 Kcal/mol, respectively. Lower binding energy indicates better interaction, more a nity, and more stability between ligand and protein.
The binding position, hydrogen bonds, alkyl bonds, van der Waals forces, and binding residues are shown in Figure 3.
The HDOCK server was used for protein-protein docking and the possibility of interaction between CB2 and Gprotein. The docking results for wild-type CB2 and G-protein were -452.97 docking score and 1.06 A° rootmean-square deviation (RMSD), and the binding a nity was -13.3 (Kcal/mol). These data showed the correct binding position of G-protein and CB2-Q63 after docking compared with CB2 structure coupled with G-protein (PDB ID: 6PT0). The docking results for the I-TASSER model and G-protein were -350.14 docking score and 80.32 A° ligand RMSD, and for the Phyre2 model and G-protein were -322.83 docking score and 192.66 A°l igand RMSD. These results show that both predicted structures of mutant CB2 cannot bind to G-protein in the correct position. Interaction between both CB2-Q63 and CB2-63R models and G-protein is shown in Figure 4.

Discussion
As a newly emerging virus, all people are susceptible to SARS-CoV-2 infection, though the nature and severity of COVID-19 vary signi cantly among cases. Notably, reported disease burden and case fatality rates differ considerably from one country to another (36). The exact in uence of host genetic makeup in this variation has remained mostly unknown. The importance of host genetics' contribution in differential responses to SARS-CoV-2 is highlighted by a modeling study, revealing that 50% of the variance of the 'predicted COVID-19' phenotype is due to genetic factors (6). For example, studies reported that the variable expression pattern and genetic variation in angiotensin-converting enzyme 2 (ACE2), as a functional receptor for SARS-CoV-2, might be associated with the susceptibility to infection and severity of the disease (37). Notably, a recent genomewide association study revealed that critical illness in COVID-19 is related to host antiviral defense mechanisms (IFNAR2 and OAS genes) and mediators of in ammatory organ damage (DPP9, TYK2, and CCR2) (38). Thus, genes related to immune responses are of particular interest for our understanding of predisposition to severe COVID-19 because of the immunopathology of SARS-CoV-2.
The EC system has received much attention due to its regulatory roles in immune response and its effects on immune-associated disease progression (39). There is evidence supporting the system's speci c involvement in respiratory viral associated immunopathology and in modulating in ammation following infection. We previously provided evidence that the EC system plays an essential role during respiratory syncytial virus (RSV) infection in humans and mice (17). Other studies also showed that the lack of cannabinoid receptors could increase in ammation and tissue damage following in uenza virus infection, and their activation can impair immune responses induced by the virus (40)(41)(42). Such data supported the idea of using cannabinoids as a potential therapeutic approach in COVID-19 patients (16,43). The rationale behind the current study was based on SARS-CoV-2 immunopathology, together with the data available on the immune regulatory role of CB2 signaling.
The CB2-Q63R polymorphism is caused by two missense mutations in the CNR2 gene changing CAA to CGG, which lead to a substitution of an uncharged polar amino acid (glutamine) with a positively charged polar amino acid (arginine) at position 63 located at the rst intracellular loop of the CB2 receptor (24). Studies indicated that the CB2-63R variant has less functions than 63Q in the modulation of immune responses, especially T-cell proliferation (25). While the exact mechanism is unknown, a study reported that the signal intensity caused by 63R activation is relatively weaker than that caused by 63Q activation (44). The present study reports for the rst time, the association between the CB2 receptor and COVID-19 severity. A signi cant difference in the Q63R allele and genotype distributions was found between COVID-19 expired and discharged patients ( Table 2). The co-dominant, recessive, and additive inheritance models showed a signi cant association between Q63R and COVID-19 severity (Table 3). According to the co-dominant model, RR subjects showed a risk for developing severe COVID-19 more than thrice of QQ subjects.
The association between the CB2-Q63R variation and autoimmune conditions such as thrombocytopenic purpura (45), celiac disease (23), juvenile idiopathic arthritis (19), in ammatory bowel disease (46), and rheumatoid arthritis (22) has been reported. The data reported by the current authors have been implied the involvement of the Q63R variation in susceptibility to multiple sclerosis in Iranian patients (18). Interestingly, our previous study showed that the in ammatory response to virus is more inhibited in cases with QQ variants, allowing the virus to replicate and induce severe infection (17). In the case of SARS-CoV-2 infection, while a robust innate immune response is essential to eliminate viral pathogens, a prolonged or dysregulated/exuberant manner can damage the respiratory tract (47). The current results are consistent with previous studies that showed a reduced EC-induced modulation of the immune system in human subjects carrying the RR variant of CB2 compared with those having the QQ variant (2,19,23,45,46).
Importantly, experimental data from previous competition ligand binding assay showed that the binding a nity of 2-AG and CB2-63R is similar to CB2-Q63, but CB2-63R had a signi cantly lower maximum response after binding to 2-AG compared with Q63 type (24). Our molecular docking results con rm that the CB2-63R induces a more inadequate response to ligand binding. The biological effects of cannabinoids are mediated through the activation of G-protein-coupled cannabinoid receptors (48). G-proteins act as adaptors that link Gprotein-coupled receptors (GPCRs) to other signaling and regulatory proteins to operate or modulate intracellular signaling pathways (15). The molecular docking results showed that the predicted structures of mutant CB2 could not bind to G-protein in the correct position, resulting in EC signaling dysfunction. This data is consistent with Wang et al., nding that CB2-63R induces lower signaling transduction than CB2-63Q in human primary T-cells (44). The CB2 receptor is predominantly expressed in the immune and immune-derived cells, and its activation indirectly affects viral infections by altering host immune responses, particularly in ammation, along different signaling pathways (49,50).

Conclusion
Data from the current study points toward host genetic involvement in the severity of COVID-19. Host genetic certainly affects the balance of immune responses during any viral immunopathogenesis, leading to different clinical phenotypes. Results indicate that people with CB2-63R variant are more prone to develop severe COVID-19. Considering the potential of this polymorphism as a biomarker in COVID-19 severity, there is an urgent need to deepen these ndings through further studies. Regarding the sample size limitation and different genetic backgrounds in various populations, other studies using whole-genome sequencing with a large cohort for multiple people would be required. Identi cation of genes related to susceptibility and severity of COVID-19 may lead to speci c targets for repurposing or drug development. We hope that, with great efforts, scienti c support, and information sharing, the overcoming of COVID-19 will come soon.

Human study
A total of 200 Iranian COVID-19 patients were included in this study. The case group consisted of 100 expired cases (50 women and 50 men) with a mean age of 62.08 years, and the control group consisted of 100 discharged cases (50 women and 50 men) with a mean age of 54.45 years. All COVID-19 patients were con rmed by real-time RT-PCR assay targeting the SARS-CoV-2 nucleoprotein (N) and ORF1ab genes (Pishtazteb, Iran). A clinical questionnaire was developed for this study and used to collect data from all patients, including gender, age, and clinical symptoms at admission (Additional le 1). Nasopharyngeal samples were collected from all patients and divided into two groups according to their disease outcome. All of the subjects in this study were from Golestan, Province, and had the same geographical origin, and none were related.
Genomic DNA was extracted from the collected nasopharyngeal samples using a DNA extraction kit following manufacturers' instructions (GeneAll, South Korea). Extracted samples were genotyped for the CNR2 rs35761398 (Q63R) variant using a TaqMan assay with commercial primers/probes (Thermo Fisher, USA). The reaction conditions were as follows: 95 °C for 4 min, followed by 50 cycles of 95 °C for 15 s, and 60 °C for 90 s. Both PCR and post-PCR allelic discrimination was performed on an ABI PRISM 7300 system (Applied Biosystems, USA). Genotypes of ten percent random samples per each group were con rmed by direct PCR sequencing as described previously (18).
The demographic and clinical data were analyzed with SPSS23 software (IBM, Chicago, IL). The Hardy-Weinberg equilibrium (HWE) and differences in allele and genotypic frequencies were calculated using the SNPStats software (a web tool for the analysis of association studies: http://bioinfo.iconcologia.net/SNPstats) (26). Inheritance models such as co-dominant, dominant, recessive, over-dominant, and log-additive were analyzed using the SNPStats software. The OR was adjusted by age in the logistic regression model. The power of the test was calculated using G-power 3.1.9.4 software (Universitat Kiel, Germany). Odds ratios (OR) and 95% con dence intervals (CI) were calculated, and p-values less than 0.05 were considered statistically signi cant.

Molecular docking
To predict the 3-dimensional (3D) structure of mutant CB2 (63R), CB2 protein sequence (NCBI Reference Sequence: NP_001832.1) with a change in position 63 was submitted to I-TASSER server (https://zhanglab.ccmb.med.umich.edu/I-TASSER), which use a hierarchical approach to protein structure prediction (27). The best model is selected from the output based on the con dence score (c-score). The Phyre2 server was also used to predict the 3D structure of 63R based on homology (http://www.sbg.bio.ic.ac.uk/phyre2/html/page.cgi?id=index) (28). The best model based on con dence and identity was selected as a homology model. Both models were submitted to ModRe ner (https://zhanglab.ccmb.med.umich.edu/ModRe ner) for atomic-level, high-resolution protein structure re nement to energy minimization and structure re nement (29). Ramachandran validation was performed using the MolProbity server (http://molprobity.biochem.duke.edu/) for backbone and structure validation of the built 3D models (30

Consent for publication
Not applicable.   Visualization of docking analysis of 2-Arachidonoylglycerol (2-AG) with A) wild type CB2, B) I-TASSER built mutant CB2 model, and C) Phyre2 built mutant CB2 model. The binding position, hydrogen bonds, alkyl bonds, van der Waals forces, pi-sigma, carbon, and binding residues are shown.