Previously-published RNA sequencing datasets were retrieved from the Gene Expression Omnibus (GEO) repository (Edgar et al., 2002; Clough and Barrett, 2016). Briefly, GEO datasets were scanned for adequate output using specific search terms. Table 1 describes the datasets used in the current study, with corresponding accession links given in Supplementary Table S1. We focused on 13 specific genes from the eCB system (Supplementary Table S2). Supplementary material M1 summarizes all DESeq2 results of the present study.

Differential gene expression analysis was performed in the freeware R Studio version 4.0.3 (R Core Team, 2013; RStudio Team, 2020). Following the workflow described by Sanchis et al. (2021), we first retrieved count data, using the GEOquery package. If necessary, count data were then normalized using the DESeq2 package. Differential expression analyzes were performed using the DESeq2 package. Gene annotations were performed with either the org.Hs.eg.bd (Homo sapiens), org.Mm.eg.db (Mus musculus) or AnnotationDbi packages, as appropriate. Volcano plots were drawn with the EnhancedVolcano package. Gene ontology (over-representation analysis, ORA) was performed with the GOplot package using Biological Processes. Gene set enrichment analysis (GSEA) was performed with clusterProfiler and enrichplot packages. Scatterplot graphs were drawn using the ggplot2 package. All of the above-mentioned packages are available under open access at the Comprehensive R Archive Network (CRAN) repository.¹

Phylogenetic analyzes were performed in Mega X (Kumar et al., 2018) following alignment of FASTA sequences. Evolutional analysis was performed using the maximum likelihood method. Initial trees were obtained by applying Neighbor-Join (NJ) and BioNJ algorithms (Gascuel, 1997) to a matrix of pairwise distances estimated using the Jones-Taylor-Thornton (JTT)-based model (Jones et al., 1992). The final drawn trees are those with the highest log likelihoods. One hundred bootstraps (r = 100) were used to infer confidence in the output trees. BLASTP was used to scan different species with Homo sapiens as bait (reference sequence), available at the National Center for Biotechnology Information (NCBI). ²

To search for conserved domains/residues, protein sequence (Homo sapiens) alignments was performed in Clustal Omega, available at the European Molecular Biology Laboratory – European Bioinformatics Institute (EMBL-EBI). ³

Structural prediction of transmembrane protein topology (TMHMM, Transmembrane prediction using Hidden Markov Models) was performed from protein FASTA sequences using the platform provided by the Department of Health Technology at the Technical University of Denmark,⁴ which was made publicly available in 2001 (Krogh et al., 2001). Homology models were created using Phyre2 (Kelley and Sternberg, 2009; Kelley et al., 2015). For docking of anandamide into TRPV1, the Cryo-EM structure of human TRPV1 in complex with the analgesic drug SB-366791 (PDB: 8GFA) was used. The A chain of the multimer was used with ligand and phospholipids removed from the .pdb file. The anandamide ligand (CHEMBL15848) was retrieved from the ChEMBL webserver as a .csv file and converted to a .pdb file using the Online SMILES Translator and Structure File Generator and then to a .pdbqt file with merged non-polar hydrogens with AutoDockTools version 1.5.7. The ligand was docked into the receptor using AutoDock Vina 1.1.2 (Trott and Olson, 2010) with a docking grid 20 Å in length in the x, y, and z directions centered on the “tunnel” anandamide-binding pocket previously identified (Morales et al., 2022), flexible N438, Y555, Y554, Y487 and D708 side-chains and an exhaustiveness of 20. Structural analysis was performed by superposition in PyMol (Schrödinger and DeLano, 2020). A professional academic license was acquired for PyMol.

Protein–protein interaction networks (Szklarczyk et al., 2019) were predicted using the STRING (Search Tool for the Retrieval of Interacting Genes) database server, ⁵ which contains both known and predicted protein interactions. Queries were pasted with either a single protein name or multiple names (as appropriate). When relevant, outer networks were searched for by expanding the nodes on display, using the built-in database tool, to reflect the immediate, intermediate and extended protein interactomes. For protein clustering within the networks, k mean clustering (based on centroids) was applied (Dalmaijer et al., 2022) directly in the server.

We first examined dataset GSE80655 (Ramaker et al., 2017), which contains RNA sequencing of 281 postmortem samples extracted from patients with either major depressive disorder (MDD, n = 69), bipolar disorder (BP, n = 71), schizophrenia (SCZ, n = 71) or controls (n = 70), in an attempt to decipher the potential role of the eCB system in these pathologies (Table 1). Out of 13 transcripts of interest (CNR1, CNR2, GPR55, GPR110, TRPV1, MGLL, DAGLA, PLCB1, NAPEPLD, FAAH, ABHD6, PLA2G4A and ALOX8, Supplementary Table S2), none were considered significantly up- or down-regulated in MDD (Figure 1A), respectively defined as log₂(fold change) >1 or < −1 and -log₁₀(p) ≥ 1.3 (p = 0.05). Similar results were observed with BP (Figure 1B) and SCZ (Figure 1C). These results indicate that our 13 transcripts of interest are not dysregulated under these conditions.

Next, we investigated whether we could detect transcriptional changes in our 13 transcripts of interest following acute or chronic THC exposure. To that end, we used GSE189821, GSE116813 and GSE70689 (Table 1). In GSE189821, 119 samples were analyzed from adolescent mice chronically exposed to THC (3 weeks at 10 mg/kg/day) or not (vehicle) and following a 14-day washout period, mimicking human behavior. In accordance with previous observations (Zuo et al., 2022), we found only minor transcriptomic dysregulations of our 13 eCB-related transcripts of interest following THC exposure in these mice (males and females pooled). This might be due to the long washout period that was applied. However, we noted 31 differentially expressed genes, with 13 transcripts down-regulated and 18 transcripts up-regulated, outside of our list of interest, when pooling all different brain regions together (Figure 2A). In the ventral tegmental area, GPR55 was found to be significantly upregulated after THC exposure, the only one from our list of interest and out of all brain regions examined. In GSE116813, 12 hippocampal samples were sequenced from mice exposed to chronic THC (21 days at 8 mg/kg/day) or not (Jouroukhin et al., 2019) and after 21 days of washout. After subsetting the dataset to only include wild-type mice, we observed 134 differentially expressed transcripts, 38 up-regulated and 96 down-regulated (Figure 2B). As observed with the previous dataset, none of these differentially-expressed transcripts were linked to our list of interest. The next dataset, GSE70689, was used to observe acute transcriptional changes in mice microglia exposed to in vitro THC (10 μM for 2 h) (Juknat et al., 2013). Here, we performed subsetting of microglia exposed to either THC or vehicle, excluding samples exposed to cannabidiol. Using such a method, 60 transcripts were found to be significantly differentially expressed, 2 up-regulated and 58 downregulated (Figure 2C). None of our transcripts of interest were found to be significantly dysregulated in such a dataset. A further analysis revealed that none of the previously differentially expressed genes (outside our 13 transcripts of interest) overlapped in the three datasets (Venn diagram in Figure 2D). In addition to these observations, we also ran ORA Gene Ontology (GO) analyzes for each dataset. GSE189821 and GSE116813 returned no GO enrichment, while GSE70689 revealed many GO enrichments. These included transcripts involved in immune response, such as T cell activation, lymphocyte migration, cytokine-mediated signaling and ERK (extracellular signal-regulated kinases) cascades (Figure 2D). These surprising results indicate that our 13 transcripts of interest are only minimally involved in the response to in vivo THC exposure, at least in mice. However, in vitro exposure to THC likely activates microglia, in an attempt to prevent THC-induced inflammation. Representations of the three datasets (GSE189821, GSE116813 and GSE70689) and our 13 transcripts of interest can be found on Figure 2E. Note that GPR55 is the only transcript of interest that was significantly upregulated after chronic in vivo THC exposure in mice, and only in the VTA (Figure 2E). This suggests that our 13 transcripts of interest were not affected by THC exposure, at least under these specific conditions, except for GPR55. GSEA analyzes did not reveal any biological processes of interest (not shown).

Our next step was to investigate the eCB system with regard to temporal lobe epilepsy, in light of extensive evidence of the links between the two (Alger, 2004; Kwan Cheung et al., 2019; Bilbao and Spanagel, 2022; Martínez-Aguirre et al., 2022; Sugaya and Kano, 2022). Three datasets were thus investigated: GSE77578 (Srivastava et al., 2018), GSE190451 (Chen et al., 2023) and GSE74150 (Damasceno et al., 2020). In GSE77578, 56 hippocampal samples were compared from epileptic or control mice. Here, we only compared control mice to epileptic mice (temporal lobe epilepsy), excluding mice treated with an antiepileptic drug. We found 129 differentially expressed transcripts in such a model of temporal lobe epilepsy, with 43 transcripts up-regulated and 86 down-regulated (Figure 3A). In GSE190451, 6 neocortical samples were compared between patients with temporal lobe epilepsy (n = 3) and aged-matched non-epileptic controls (n = 3), an analysis that differs slightly from the published article, in which the authors analyzed a total of 8 samples (Chen et al., 2023). We found a total of 1873 differentially expressed transcripts in temporal lobe epilepsy, 833 up-regulated and 1,040 down-regulated (Figure 3B). In GSE74150, an animal model of genetic epilepsy was used (Wistar audiogenic rat). Here, we found 67 differentially expressed transcripts, 53 down-regulated and 14 up-regulated (Figure 3C). Up- and down-regulated transcripts in the three datasets presented some overlap (Figure 3D, Supplementary Table S3), showing that transcripts involved in, or responding to, epilepsy seem to be found across species and models. Several gene ontology enrichments (ORA) were detected in the first two datasets, ranging from regulation of neurophysiology to neuronal death. Interestingly, 37 biological process terms were found to be common across the two first datasets (Supplementary Table S4), such as synapse organization, neuron death, myelination and plasma membrane organization (Figure 3C). GSEA analyzes did not reveal any biological processes of interest (not shown). Out of our 13 transcripts of interest, TRPV1 was found to be significantly down-regulated in the neocortex of patients with temporal lobe epilepsy (Figure 3D). These results confirm, once again, that TRPV1 plays an important role in the pathophysiology of epilepsy, as observed previously (Jia et al., 2015; Saffarzadeh et al., 2015; Kong et al., 2019; Wang et al., 2019; Lazarini-Lopes et al., 2022).

Using the STRING database (Szklarczyk et al., 2019), containing both known and predicted protein interactions, we first fed into the database single protein names (from our list of 13 transcripts/proteins). Analysis of the immediate interactomes revealed great disparities in terms of protein overlap in pair-wise comparisons (Figure 4A). Interestingly, great overlap was observed with CB1R, CB2R and ADGRF1 (GPR110). We also noted great protein overlap in the interactomes of DAGLA, NAPE-PLD and ABHD6. As we expanded the interactome searches, similar conclusions were reached (Figure 4A). Surprisingly, TRPV1 and ALOX15B presented the lowest overlap in all levels, whether immediate, intermediated or extended interactomes (Figure 4A). We next scanned the database using our 13 proteins of interest by bulk feeding all 13 proteins into the search equation. When examining the immediate interactome of these 13 proteins, only the same 13 proteins were returned (Figure 4B), providing no additional information of interest. However, when the interactome search was expanded to an intermediate level, new proteins appeared, including GNAI1, GNAI2, GNAQ, and TRMP8 (Figure 4B). When the search was expanded to another level (extended interactome), new protein interactions arose, such as CALM3, CALML3, CALML5, CALML6 and TRPA1 (Figure 4B). Our final step was to apply k means clustering using a fixed number of clusters (k = 3). Very surprisingly, 4 proteins clustered outside of the main cluster containing the proteins of interest: TRPV1, PLC1B, ALOX15B and PLA2G4A (Figure 4C). We also noted that TRPM8 and GATB clustered with the majority of our proteins of interest, via FAAH, acting as a link within the extended interactome (Figure 4C).

Our next step was to try to determine how CB1R (the protein produced from the CNR1 gene) evolved within the human proteome. Using GeneCards (Stelzer et al., 2016) with Homo sapiens CB1R as input, 18 paralogues were returned. As expected, the closest paralogue was CNR2 (Figure 5A). CNR1 and CNR2 formed a distinct clade, with the next most closely related G-coupled protein receptors (GPCRs) being lysophosphatidic acid receptors (LPARs) and sphingosine-1-phosphate receptors (S1PRs). The furthest clade included melanocortin receptors 1–4 (MCRs) and GPR119. Inspection of the sequences revealed three conserved regions (Figure 5B), which could be divided into residues conserved within the clade containing CNR, LPAR and S1PR families (Figure 5B, blue) and residues only conserved within the CNR clade (Figure 5B, orange). Three conserved motifs (CMs) were identified: CM1, CM2 and CM3 are located within transmembrane (TM) 2, 3 and 7, respectively (Figure 5C). We then examined the location and function of these residues in the Homo sapiens CB1R structure bound to an anandamide analog (Krishna Kumar et al., 2023) and in the lysophosphatidic acid receptor structure bound to its native lysophosphatidic acid ligand (PDB: 7TD1) (Liu et al., 2022). The residues conserved within all three CNR, LPAR and S1PR families are largely found in the lower part of the receptor and include five of the six residues identified as being important in receptor activation (Figure 5D), as observed before (Hua et al., 2020). The ligand binding pocket in both receptors is formed by a large pocket of hydrophobic residues. Most of the residues conserved only within the CNR clade are located in this area and include some hydrophobic residues and H178 which forms a hydrogen bond with the terminal hydroxyl of the ligand (Figures 5D,E).

Out of 13 deuterostome species spanning from mammals to echinoderms (plus insects as protostomes), CB1R orthologues were found in placental mammals, marsupials, monotremes, birds, reptiles, amphibians, teleosts, elasmobranchs and Agnatha (Supplementary Table S5), but absent in cephalochordates, tunicates, echinoderms and insects. Furthermore, CB1R orthologues were also not found in plants, bacteria, Cnidaria and Porifera, except in Klebsiella (gram-negative bacteria) which presented 100% homology with the bird CB1R sequence (not shown), thus likely accounting for a database error in GenBank. When the same experiment was conducted with CB2R, we observed the presence of orthologues in placental mammals, marsupials, monotremes, birds, reptiles, amphibians and teleosts, the latter presenting two orthologues (Supplementary Table S5). Within the BLASTP database, unknown CBRs were found in cephalochordates and tunicates. These results suggest that the classical endocannabinoid receptors (CB1R and CB2R) are confined to a subgroup of deuterostomes, containing urochordates, cephalochordates and chordates, but absent from other deuterostomes, other animals, plants and prokaryotes, thus partly corroborating previous observations (McPartland et al., 2006; Elphick, 2012). Furthermore, it is suggested that a duplication event has occurred in teleosts since they diverged from other chordates. A further phylogenetic analysis revealed a clear division between CB1R and CB2R into two distinct clades (Figure 6A). We noted that most taxa presented copies of both proteins. There are two explanations for the lack of CNR2 paralogues in agnathans and elasmobranchs: (i) a duplication may have occurred since the divergence of other chordates and elasmobranchs or (ii) CNR2 paralogues were lost by both agnathans and elasmobranchs, the first scenario being by far the more likely explanation. Finally, we noted that the CBRs of Cephalochordates and Tunicates are very different, making it difficult to distinguish these two proteins from other related GPCRs. To address this, we performed multiple alignments with CB1R and CB2R protein sequences from the above-mentioned species. Our first alignment included cephalochordates and tunicates and led to great disruptions of the CMs (Figure 6B). The second alignment, which was performed without cephalochordates and tunicates, produced more defined, but yet not optimal, CMs (not shown). Our final alignment was performed without cephalochordates, tunicates and Danio rerio. This final approach produced robust identification of the three CMs (Figure 6C), thus corroborating the above results. Subsequently, we created homology models of the agnathan CB1R and cephalochordate unclassified CBR. Both the agnathan (Petromyzon marinus: ModelArchive ma-yvfgy) and cephalochordate (Branchiostoma floridae: ModelArchive ma-mi5fb) receptors contain the six residues important for receptor activation. However, the cephalochordate receptor lacks the majority of the residues involved in the AEA-binding pocket, including a histidine equivalent (H178) needed for hydrogen-bond formation with the terminal hydroxyl (Figures 6D,E).

Our final analyzes focused on TRPV1, as this is the second eCB receptor in the central nervous system (together with CB1R). Using similar approaches, we first determined the phylogeny of the six Homo sapiens paralogues, using TRPV1 as the reference sequence. We observed that TRPV1 and TRPV2 diverged recently, which arose from a duplication event (Figure 7A). Likewise, divergence of TPRV3, TRPV5 and TRPV6 arose earlier, with TPRV5 and TRPV6 being the most recent split (Figure 7A). Two small CMs were noted throughout those 6 different proteins, including Y375-G376-P377, L382-Y383-D384-L385, L547-G548-W549 and R557-G558 (Figure 7B). Using Homo sapiens TRPV1 sequence as bait, we found that TPRV5 or TRPV6 were present across the most diverse range of taxa (Supplementary Table S6). These receptors are present in Bilatera but absent from Porifera, the latter presenting weak protein homology. These observations suggest that the TRPV family is specific to the Bilatera. Furthermore, TPRV1 and TRPV2 resulted from later duplication events in the common ancestor of jawed vertebrates and amniotes, respectively (Figure 7C). In addition, our phylogenetic analysis highlighted that eCB-binding receptors (TRPV1 and TRPV2) form one clade, which includes taxa from elasmobranchs onwards, with a more recent duplication forming TRPV2, found in amniotes only (Figure 7C), as suggested above (Supplementary Table S6). Alignments of TPRV1 and TRPV2 from these species correctly identified conserved residues involved in channel gating and selectivity, as well as ligand binding [anandamide, vanilloid, double-knot toxin (DkTx), Figure 7D]. Conserved residues were often found throughout species, indicative of the importance of these residues for TRPVs function. Finally, out of these 28 functional residues in TRPV1, only 8 were conserved in TRPV5, whether fully conserved or present with conservative substitutions (Figure 7E). Taken together, these results pinpoint a divergence of eCB receptors between jawed and jawless vertebrates, around 450 million years ago during the late Ordovician. These observations are consistent with previous studies (Elphick et al., 2003; McPartland, 2004; McPartland et al., 2006).

We examined the structures of TRPV1 and TRPV4 to identify residues important for ligand binding. TRPV1 can be activated by endogenous ligands such as AEA and exogenous ligands such as capsaicin, the component of chili peppers which creates a sensation of burning (Lam et al., 2005). Structural studies have identified the capsaicin binding site, known as the vanilloid-binding pocket (Nadezhdin et al., 2021; Kwon et al., 2022). There are no experimentally determined structures of the TRPV1 receptor bound to AEA. However, using a model of human TRPV1 based on the experimental rat structure, molecular dynamics and docking studies have identified two possible binding sites for anandamide (Muller et al., 2020; Morales et al., 2022), represented in Figure 8A. One site was the vanilloid-binding pocket (VBP) but the preferred site was a tunnel (AEA tunnel) involving interactions with Y554, Y555, Y487, D708, and N438. We docked AEA into the AEA tunnel using the recently experimentally determined structure of human TRPV1 (PDB: 8GFA) (Neuberger et al., 2023). One of nine docking solutions had a similar orientation to that shown in Morales et alia (Morales et al., 2022), with Y555, N438 and Y554 forming hydrogen bonds with the amide group hydrogen and the terminal hydroxyl group oxygen (Figure 8B, ModelArchive ma-dol4y). The five interacting residues are well conserved across species in TRPV1 and TRPV2. AEA itself does not directly bind to TPRV4, however this receptor is activated by an AEA metabolite 5′,6′-epoxyeicosatrienoic acid (5′,6′-EET) (Watanabe et al., 2003). Homo sapiens TRPV4 (PDB: 8TD1) (Nadezhdin et al., 2021) is missing equivalents to TRPV1 Y847 and Y555 which stabilize AEA with hydrogen bonds in the AEA tunnel. Simulations have suggested that residues K535, F549, Q550, Y591 and R594 of TRPV4 form the 5′,6′-EET binding pocket (Berna-Erro et al., 2017), as displayed in Figure 8C. Of these, Y591 is equivalent to Y554 in TRPV1 which forms a hydrogen bond with the terminal hydroxyl of AEA, but the remaining pocket is very different. We also examined differences in the vanilloid-binding pocket by superimposing capsaicin from the rat (Rattus norvegicus) TRPV2 complex structure (PDB: 8SLY) (Gochman et al., 2023) onto the Homo sapiens TRPV1 receptor (PDB: 8GFA). Whilst some important hydrophobic residues are conserved in the equivalent pocket in TRPV4, there are no equivalent residues for the hydrogen-bonding Y511 and T550 which were previously observed to be involved in propagating ligand-activated conformational changes (Figure 8D) (Kwon et al., 2022).