Abstract
During adolescence, cannabis experimentation is common, and its association with interindividual variations in brain maturation well studied. Cellular and molecular underpinnings of these system-level relationships are, however, unclear. We thus conducted a three-step study. First, we exposed adolescent male mice to Δ-9-tetrahydrocannabinol (THC) or a synthetic cannabinoid WIN 55,212-2 (WIN) and assessed differentially expressed genes (DEGs), spine numbers, and dendritic complexity in their frontal cortex. Second, in human (male) adolescents, we examined group differences in cortical thickness in 34 brain regions, using magnetic resonance imaging, between those who experimented with cannabis before age 16 (n = 140) and those who did not (n = 327). Finally, we correlated spatially these group differences with gene expression of human homologs of mouse-identified DEGs. The spatial expression of 13 THC-related human homologs of DEGs correlated with cannabis-related variations in cortical thickness, and virtual histology revealed coexpression patterns of these 13 genes with cell-specific markers of astrocytes, microglia, and a type of pyramidal cells enriched in dendrite-regulating genes. Similarly, the spatial expression of 18 WIN-related human homologs of DEGs correlated with group differences in cortical thickness and showed coexpression patterns with the same three cell types. Gene ontology analysis indicated that 37 THC-related human homologs are enriched in neuron projection development, while 33 WIN-related homologs are enriched in processes associated with learning and memory. In mice, we observed spine loss and lower dendritic complexity in pyramidal cells of THC-exposed animals (vs controls). Experimentation with cannabis during adolescence may influence cortical thickness by impacting glutamatergic synapses and dendritic arborization.
Significance Statement
Cells and molecular underpinning of cannabis-related variations in cortical thickness observed in adolescents are poorly understood. We assessed differences in gene expression in the frontal cortex of adolescent mice exposed (or not) to Δ-9-tetrahydrocannabinol (THC) or WIN 55,212-2 (WIN). THC preferentially targeted glial cells, while WIN affected pyramidal neurons; both modified nuclear-coded subunits in mitochondrial respiratory complexes and excitatory synapse genes. In humans, we evaluated a spatial correlation between group differences in cortical thickness and the expression of human homologs of the cannabinoid-sensitive genes in mice. These genes coexpressed with those specific to astrocytes, microglia, and a type of pyramidal cells enriched in dendrite-regulating genes. Experimentation with cannabis during adolescence may influence cortical thickness via synaptic processes and dendritic arborization.
Introduction
Cannabis is the most widely used psychoactive drug, with an annual prevalence of consumption of 2.5% of the world's population (World Health Organization, 2022). Consumption peaks in adolescence (Government of Canada, 2018; Johnston et al., 2019), when anatomical remodeling of circuits in the cerebral cortex supports cognitive and emotional maturation into adult life (Tau and Peterson, 2010; Larsen and Luna, 2018; Chini and Hanganu-Opatz, 2021). Δ-9-tetrahydrocannabinol (THC) is the psychoactive component in cannabis and produces its effects via the endocannabinoid system, which regulates developmental plasticity in the adolescent brain (French et al., 2015; Lu and Mackie, 2021).
Cannabis and thickness of the cerebral cortex
A growing body of evidence confirms the relationship between adolescent cannabis exposure and cortical thickness in humans (Chye et al., 2020; Albaugh et al., 2021; Lichenstein et al., 2022; Owens et al., 2022; E. A. Robinson et al., 2023). During this period, structural maturation in the cerebral cortex is characterized by a linear decrease in its thickness between the ages of 5 and 20 years (Ducharme et al., 2016). Taking a neurodevelopmental perspective, three recent studies revealed associations between cannabis exposure during middle to late adolescence and age-related cortical thinning in the frontal lobe (Albaugh et al., 2021, 2023; Owens et al., 2022).
The structural effects of cannabinoids, notably at the cellular level, have also been studied in animal models. Such approaches are particularly useful because their methodological design provides mechanistic insight on how cannabinoids influence both brain structure and function. Repeated exposure of male adolescent rats to THC was shown to induce premature pruning of dendritic spines on layer 3 (L3) pyramidal neurons of the prelimbic cortex (vs vehicle-treated) and a protracted atrophy of their dendritic trees (Miller et al., 2019). The transcriptome of the L3 pyramidal neurons, which was studied in parallel with these structural variations, indicated alterations in gene networks associated with cell morphogenesis, cytoskeleton, and dendritic development (Miller et al., 2019). Regarding cognition, adolescent male (CD1) mice that received daily THC injections for 21 d showed lower (vs all other treatment groups) performance on object recognition memory in the days immediately following the end of the exposure and persisted during adolescence (M. Murphy et al., 2017). In addition, daily injections of THC for 21 d reduced (vs vehicle) acquisition and consolidation of spatial working memory during a delayed alternating T maze task performed by adolescent mice (C57BL6/J) of both sexes (H. T. Chen and Mackie, 2020). In rodents, cannabinoid-related deficits in working memory have been related to synaptic variations that were mediated via glial targets (Han et al., 2012; R. Chen et al., 2013).
Objective of this study
Considering the high prevalence of adolescent cannabis use in North America (Hammond et al., 2021; Lim et al., 2022) and Europe (EMCDDA, 2017), the increasing concentrations of THC in commercial cannabis products (ElSohly et al., 2021; Freeman et al., 2021), and the emergence of synthetic cannabinoids (Vardakou et al., 2010), it is imperative to improve our understanding of the effects of cannabinoids on brain maturation and cognition (Morin et al., 2019). For this reason, we sought to identify the cellular and molecular underpinnings of cannabinoid actions, which could shed light on the relationship between cannabis experimentation and the cerebral cortex during adolescence. To achieve this aim, we used a two-pronged strategy. First, we used adolescent male mice to identify genes and cell types in the frontal cortex that were susceptible to THC or the synthetic cannabinoid WIN, as assessed by their transcriptomic response to these compounds. Second, we used transcriptomic information to determine how human homologs of murine genes that were affected by the treatment are related to inter-regional differences in cortical thickness as a function of cannabis experimentation in human adolescents. In a species-specific way, we further assessed how genes modified by the two ligands are related to different molecular processes as identified by gene ontology enrichment analysis.
Materials and Methods
Mouse experimental design and methods
Animals
Adult male CD56/Black6 mice weighing 17–21 g were purchased from Charles River Laboratories and housed in a control environment on a 12 h light/dark cycle with ad libitum access to food and water. All experimental methods and animal care procedures were approved by the Animal Care Committee of the CHU Ste Justine Research Centre (CIPBAR protocol 2019-4097), in accordance with the guiding principles as enunciated by the Canadian Council on Animal Care.
Drug preparation
WIN 55,212-2 and THC were purchased from Cayman Chemical. Stock solutions of the drugs were prepared in ethanol.
Treatments
Adolescent (male) mice were injected twice a day with THC (10 mg/kg, s.c.) from postnatal day (P) 35–38, plus a single injection on P39 after which they were killed. The daily dose used in this study is a variation on a dosing scheme that has been frequently used to assess adolescent exposure to THC in mice, where the drug is given for 12 d from the beginning to end of adolescence with dosing escalation every 4 d for a total 160 mg of THC (Rubino et al., 2009, 2015; Zamberletti et al., 2015; Renard et al., 2017; Qvist et al., 2022). We modified this scheme by constraining THC administration to 4.5 d during the second half of mouse adolescence. The daily dose we administered was the one used in the original model during later stages of adolescence, resulting in a cumulative dose of 90 mg THC (n = 6). WIN 55,212-2 was administered following the same scheme of two daily injections, at a dose of 3 mg/kg each (n = 4). The 1:3 dosing ratio between the two treatments takes into consideration affinity of WIN for CB1/CB2 receptors (Showalter et al., 1996; Wallace et al., 2009; AlKhelb et al., 2023; Harding et al., 2024) and its lower brain penetration as compared with THC (Barna et al., 2009; Torrens et al., 2020). Control mice received sesame oil (n = 6), which was the vehicle for both cannabinoids.
Cortical dissection and cell suspensions for scRNASEq
On the day of killing, male mice aged P39 were anesthetized with ketamine/xylazine (4:1) and transcardially perfused with ice-cold oxygenated aCSF solution of the following composition (in mM): 62.5 NaCl, 70 sucrose, 2.5 KCl, 25 NaHCO3, 1.25 NaH2PO4, 7 MgCl2, 1 CaCl2, 20 HEPES, and 20 glucose. Brains were then quickly removed, transferred to ice-cold oxygenated aCSF solution, and fixed to a stage allowing to cut 300 μm slices between 1.5 and 3.8 mm anterior from the Bregma with a VT1200 S, Leica Vibratome. Slices were individually transferred to a cutting surface in ice-cold aCSF, where the frontal cortex was microdissected. Dissected tissue was then dissociated into cell suspension (Muñoz-Manchado et al., 2018) using the Papain dissociation system (Worthington) as per manufacturer's instructions, with the addition of AP-V (0.5 mM), kynurenic acid (0.4 mM), and 5% D-trehalose (5%) as described in Saxena et al. (2012). Post trituration, cells were counted using iNCYTO C-Chip hemocytometers and resuspended at 150–200 cells/μl in Drop-seq buffer (in mM): ×1 HBSS, 2 MgSO4, 2.13 MgCl2, 1.26 CaCl2, 10 glucose, 5% trehalose, and 0.01% BSA.
Single-cell capture and sequencing
Resuspended cells were run through a custom Drop-seq setup allowing the capture of single cells and barcoded beads (ChemGenes) within nanoliter droplets as described in the original Drop-seq protocol (http://mccarrolllab.com/dropseq/; Macosko et al., 2015). The setup closely resembled the original one (Macosko et al., 2015), with minor modifications, as previously described (Lavertu-Jolin et al., 2023). In particular, Medfusion 3500 syringe pumps (Smiths Medical) instead of KD Scientific Legato 100 pumps were used, as well as a different RNase inhibitor (SUPERase In; catalog #AM2694; Thermo Fisher Scientific). Single-cell suspensions (150–200 cells/μl) were run alongside barcoded beads (140 beads μl/1) in 1.7 ml of cell lysis buffer (Tris 200 mM, EDTA 20 mM, 20% Sarkosyl, 20% Ficoll, DTT 50 mM; also containing 50 μl of RNase inhibitor), to produce beads attached to single-cell transcriptomes. These beads were subsequently washed, reverse transcribed, treated with exonuclease, and counted using a Fuchs-Rosenthal hemocytometer (INCYTO C-Chip, catalog #82030-474), as described in the original Drop-seq protocol (Macosko et al., 2015). Bead counting was performed using water and bead counting solution (10% polyethylene glycol, 2.5 M NaCl) in a 1:1 ratio as described in Gierahn et al. (2017). For PCR amplification, 4,000 beads (∼200 single-cell transcriptomes attached to microparticles) were used as input for each PCR. Individual PCR reactions were then pooled to achieve the desired number of single-cell transcriptomes attached to microparticles (STAMPS). In this experiment we used 2,000 STAMPS/reaction, which were cleaned up with AMPure XP (Beckman Coulter Life Sciences) beads at a ratio of 0.6× prior to sequencing. Correct size distribution and concentration of complementary DNA was determined using a Bioanalyzer High Sensitivity DNA assay (catalog #5067-4626; Agilent Technologies). Library preparation was performed using the Nextera XT DNA Library Preparation Kit (Illumina) with 600 pg input according to the Drop-seq protocol using Illumina i7 primers (N701–N706) together with a custom P5 primer (GCCTGTCCGCGGAAGCAGTGGTATCAACGCAGAGTAC) in the library PCR amplification (see the Drop-seq protocol; Macosko et al., 2015). Libraries were quality controlled for size and concentration using the Bioanalyzer High Sensitivity DNA assay and quantified using the Universal Library Quantification Kit (KAPA) before sequencing on a NextSeq 550 system (Illumina) at the Institute for Research in Immunology and Cancer (IRIC) in Montreal. Sequencing was performed on a NextSeq 500 system with the read settings R1 = 20 bp, R2 = 63 bp, and index = 8 bp. Cells were sequenced to a median depth of 16,000 aligned reads/cell, representing a median of 1,540 genes/cell.
Confocal imaging, analysis of spines, and of dendritic architecture
Male Thy-1 GFP mice from Jackson Laboratories (JAX: 007788) were treated with THC or vehicle, and at the end of treatment they were perfused with 4% paraformaldehyde (PFA). Brains were then removed, cryoprotected in 30% (w/v) sucrose (Sigma)/PBS, and stored at −80°C until they were sectioned coronally at 30 μm on a Leica VT100 cryostat. Floating sections stored in PBS before mounting in IMMU-MOUNT. Images from prelimbic cortex layer 2/3 were acquired with a Leica confocal microscope (SP8) using a 63× oil objective and a z-step size of 0.2 μm. All confocal parameters were maintained constant throughout acquisition within an experiment. Images were processed with Imaris (Oxford Instruments) software.
Spine analysis
Spines were automatically identified by the software within a radius of 100 μm from the soma. At least 5–6 neurons per mouse were considered. Spines were counted on primary and secondary apical branches and on 4–7 basal dendrites per neuron. Each of five independent experiments included control and treated mice.
Dendritic architecture
Fully labeled PFC-pyramidal neurons having triangular cell soma with unbroken apical and basal dendrites were considered. Dendritic complexity was quantified using Sholl analysis in Imaris. Briefly, apical and basal dendrites were drawn semiautomatically, separately, across all stacks, and Sholl rings were placed 20 μm apart, centered on the cell soma. Morphometric parameters, including dendritic length, branch points and branch depth, as well as the number of intersections with Sholl rings, were calculated. Each of five independent experiments included control and treated mice and four neurons were selected per mouse. The investigator was blind to treatment conditions during image acquisition and analysis. Results were analyzed with GRAPHPAD10.
Immunohistochemistry
Mice that received either THC, WIN, or vehicle were anesthetized and then perfused intracardially with PBS followed by 4% (w/v) paraformaldehyde (PFA) in PBS. Intact brains were extracted and postfixed in 4% PFA/PBS for 3 h. The tissue was then cryoprotected in 30% (w/v) sucrose (Sigma)/PBS, sectioned coronally at 40 μm on a Leica VT100 cryostat, and stored as floating sections in PBS. Brain sections were then blocked in 10% normal goat serum (NGS, Invitrogen, 10000C) in PBS containing 1% (v/v) Triton X-100 for 2 h at room temperature. Primary antibodies were diluted in 5% goat serum in PBS containing 0.1% (v/v) Triton X-100 and incubated with brain sections for 24 h at 4°C. The primary antibodies used in this study were the following: mouse monoclonal anti-NDUFB8 (1:200, Abcam, AB110242), rabbit monoclonal anti-NeuN (1:500, Millipore, MABN140), and rabbit monoclonal anti-TMEM119 (1:500, Abcam, AB209064). Slices were then washed in PBS (3 × 10′), incubated in the appropriate Alexa-conjugated antibodies in 5% NGS, 0.1% Triton in PBS for 2 h at room temperature, washed again in PBS (3 × 10′), and mounted in IMMU-MOUNT before imaging. All images were acquired using Leica confocal microscopes (SP8). We imaged prelimbic cortex layer 2/3 using a 63× oil objective. At least six confocal images per mouse were acquired with z-step size of 0.5 μm. All the confocal parameters were maintained constant throughout acquisition within an experiment. The number of NDUFB8 particles was manually counted through all the cell volume using ImageJ-Fiji software. The investigator was blind to treatment conditions during imaging and analysis. Each experiment included control and treated mice. Results were analyzed with GRAPHPAD10.
Monitoring CB1/CB2 receptor signaling with BRET (bioluminescence resonance energy transfer)
Ligand preparation
Ligands were initially dissolved or diluted in ethanol to achieve a stock concentration of 10 mM. The dissolved ligands were then serially diluted in logarithmic steps of 100.5, resulting in final treatment concentrations ranging from 10−4.5 M down to 10−9.5 M.
Cell culture and transfection
HEK293 SL cells were a kind gift from Dr. Laporte, McGill University. They were cultured in 100 mm Petri dishes at 37°C and 5% CO2 in the Dulbecco's modified Eagle's medium (DMEM) supplemented with 10% fetal bovine serum, 2 mm l-glutamine and 100 unit/ml penicillin-streptomycin.
Transient transfections of vectors encoding BRET biosensors in combination with complementary signaling partners were performed in 96-well culture plates (PerkinElmer; 35,000 cells/well), using the polyethylenimine transfection reagent (Polysciences) at a PEI/DNA ratio of 3:1.
BRET assays
Gαi1, Gαi2, Gαi3, Gao activation assay: HEK 293 cells were cotransfected with either human CB1 or CB2 receptors and the BRET biosensor for Gαi/o-coupled subunits: Rap1Gap-RlucII/rGFP10-CAAX (Avet et al., 2022). The ratio of DNA in the transfection solution was 1/0.65/0.065 for receptor/Rap1Gap-RlucII/rGFP10-CAAX (μg per 96-well plate surface). The corresponding untagged Gα subunits were cotransfected at a 1/1 ratio with Rap1Gap-RlucII and their responses separately monitored in different cells.
β-Arrestin recruitment was tested in the presence of GRK2 (Benredjem et al., 2019). To this end, HEK 293 cells were cotransfected with either human CB1 or CB2 receptors and the BRET biosensor for βarr recruitment: βarr2-RlucII/rGFP-CAAX at the following DNA ratios (μg per 96-well plate surface), 1/0.2/2.3 receptor/βarr2-RlucII/rGFP10-CAAX (Namkung et al., 2016). GRK2 was cotransfected at a 1:1 ratio with the receptor.
Forty-eight hours after the transfection, culture media was removed, and the cells were washed with Tyrode's solution (in mM: 140 NaCl, 2.7 KCl, 1 CaCl2, 12 NaHCO3, 5.6 D-glucose, 0.49 MgCl2, 0.37 NaH2PO4, 25 HEPES, pH 7.4). Following the washing step, 85 μl of Tyrode's solution was added to each well for 30–60 min at room temperature.
Prolume Purple (MeO-e-CTZ; Nanolight Technology) was added to the cells to a final concentration of 5 μM for 3 min. At this time, the plate was treated with the respective ligands from the lowest concentration to the highest concentration. Plates were read 5 min after on the Mithras LB 940 microplate reader (Berthold Technologies), 0.5 s per well, with filters set at 400 nm (RlucII) and 515 nm (GFP10). BRET ratios were calculated as GFP10/RlucII emissions.
Dose–response curves were analyzed with GRAPHPAD10, comparing THC and WIN curves by two-way ANOVA.
Computational analysis of single-cell transcriptomic data
Demultiplexing of raw Illumina sequencing data was performed with bcl2fastq v.2.17 at the Institute for Research in Immunology and Cancer (IRIC) in Montreal. Fastq files were quality controlled and processed using the Drop-seq computational pipeline, similarly to an earlier report (Wünnemann et al., 2020). Drop-seq tools v.1.12 (Drop-seq alignment cookbook v.1.2, January 2016) were used to produce digital gene expression matrices for each sample (Macosko et al., 2015). These matrices were processed in R (3.5.2) using Seurat (v.2.2.0; Butler et al., 2018). Cells with <200 genes expressed and genes detected in less than three cells were excluded from further analysis. Cell-cycle scores for included cells were calculated using the CellCycleScoring function from Seurat with cell–cycle gene lists from Nestorowa et al. (2016) (https://satijalab.org/seurat/articles/cell_cycle_vignette.html). Expression values were normalized and scaled to 10,000 transcripts with Seurat's Normalize function. Normalized expression values were scaled using the ScaleData function from Seurat, regressing out the number of unique molecular identified count matrices, Mouse ID, percent mitochondrial genes, as well as S phase and G2/M phases scores from cell–cycle prediction. Cells from the 16 mice were then aligned with canonical correlation analysis (CCA; Satija et al., 2015; Butler et al., 2018). To do so, we used the RunMultiCCA function on the top 1,000 highly variable genes overlapping between the samples. Subsequently, we ran the AlignSubspace function in Seurat, using the treatment as a grouping variable with an appropriate number of CCA dimensions. These dimensions were derived from the inflection point of the elbow plot generated by Seurat ElbowPlot function, resulting in 15 dimensions. The output from this analysis was then analyzed with RunTSNE followed by the FindClusters function, with a resolution parameter of 1.2. The process resulted in the identification of a total of 21 clusters within the 28,360 analyzed cells, and these were subsequently plotted as a UMAP (Uniform Manifold Approximation and Projection; Extended Data Fig. 1-1).
Identification of cell types
Cell type identification was done as previously described (Lavertu-Jolin et al., 2023) using the top markers in each cluster. Of the 21 clusters identified, two of them (clusters 20 and 21) had <200 cells and were not considered for further analysis. Among those considered, seven clusters contained Vglut1 (Vesicular glutamate transporter 1) and were considered to be pyramidal neurons. Clusters containing Gad1 and/or Gad2 (Glutamate decarboxylase 1 and/or 2) genes were considered to be GABAergic interneurons (seven clusters). Immature neurons expressing Gad1/2 were also included in this group. Cells that did not express Vglut1 or Gad1/2 genes were considered non-neuronal and added up to a total of seven clusters. In-depth analysis of cell subtypes was done as before (Lavertu-Jolin et al., 2023) by comparing the top 30 markers in each cluster with the DropViz database, a template for murine cortical transcriptomes (Saunders et al., 2018). Identity was deemed confirmed if at least 75% of the top 30 markers could be identified as being at least twofold enriched in the corresponding clusters of the DropViz database as compared with the rest of the cortex (http://dropviz.org/).
Differential expression analysis of mouse single-cell transcriptomes
Differential gene expression analysis was carried out within each cluster comparing cells from vehicle-injected mice to cells from mice that were either injected with THC or WIN 55,212-2. Statistical significance was established using Wilcoxon rank-sum test filtering out genes with a log2 fold change of <0.25 and adjusted p values of >0.05 as previously described (Wünnemann et al., 2020). Chi-square tests of independence were used to compare cellular distribution (i.e., cell type groups) of differentially expressed genes (DEGs) by each treatment (Extended Data Figs. 1-2, 1-3).
Human study participants and methods
The Saguenay Youth Study (SYS), a long-term prospective community-based cohort from a population with a genetic founder effect (European ancestry), recruited participants between 2003 and 2012 in the region of the Saguenay–Lac-Saint-Jean region in Quebec, Canada. The SYS protocol included T1-weighted magnetic resonance imaging (MRI). Additional information on recruitment and assessments for SYS (http://saguenay-youth-study.org/) is available previously (Pausova et al., 2017). For this report, we analyzed data obtained in male adolescents only in order to match the experimental findings acquired in male mice. We retrieved data from 475 adolescent males (mean age = 14.9 ± 1.8; see Table 1 for demographics). It should be noted that our previous analysis of this dataset reported no group differences in cortical thickness in SYS females as a function of cannabis exposure before the age of 16 (French et al., 2015). Assent and consent from all participants and their parents, respectively, as well as ethic research boards approvals, were obtained for this cohort.
Demographics for Saguenay Youth Study adolescents as a function of cannabis experimentation before age 16
Table 1-1
Analytic pipeline between inter-regional variations in cortical thickness as a function of cannabis exposure before the age of 16 in Saguenay Youth Study adolescent males and gene expression for THC-related and WIN-related differentially expressed genes in peri-adolescent male mice. Download Table 1-1, TIF file.
Table 1-2
Lists of human homolog ligand-related differentially expressed genes in mice. Download Table 1-2, XLSX file.
Table 1-3
Magnitude of the age-adjusted group difference between non-exposed adolescent males and cannabis-exposed adolescent males for 34 left cortical regions of the Desikan-Killiany atlas. Download Table 1-3, XLSX file.
Table 1-4
Distribution of group differences in thickness. The blue dashed lines represent the average cannabis-related thickness across 34 regions while the black full lines mark 0. Negative values indicate reduced cortical thickness in not-exposed adolescent males against cannabis-exposed adolescent males before adjusting (A) after adjusting (B) for age. Download Table 1-4, TIF file.
MRI data acquisition and processing
T1-weighted anatomical images were acquired using a Phillips 1.0-T MRI system with 3D RF-spoiled gradient echo scan with 140–160 slices at 1-mm isotropic resolution (TR = 25 ms; TE = 5 ms; flip angle = 30°). The cerebral cortex of SYS participants was parcellated into 68 cortical regions with the Desikan–Killiany atlas using the fully automated and validated segmentation program FreeSurfer version 5.3 (Fischl et al., 2002; Desikan et al., 2006). Additional information on the MRI acquisition and processing was described previously (Pausova et al., 2007).
Statistical analysis
Overall approach
For data obtained with scRNAseq in each group of mice, we first determined the number of neuronal (i.e., in cells expressing Vglut1 or Gad1/2) and non-neuronal DEGs (i.e., all but cells expressing Vglut1 or Gad1/2) according to the mRNA profile clusters (see below, Gene classification). We then intersected these lists (human homologs) with a list of 2,511 genes that show a consistent inter-regional profile of their expression in the human cerebral cortex, across lifespan and across individuals (see below, Virtual histology). For each of these gene lists, we determined the spatial relationship between gene expression and the magnitude of group differences in cortical thickness (see below, Gene expression–group differences spatial relationships). For genes that showed this spatial relationship (at p value ≤ 0.05), we then determined their coexpression patterns with cell-specific gene panels (see below, Virtual histology). Finally, we sought to characterize the functions of these gene lists by performing a gene ontology enrichment analysis (see below, Gene ontology enrichment analysis). A scheme of the analytic pipeline is presented in Extended Data Table 1-1. This figure also shows the number of homolog genes that were differentially expressed for each ligand. All analyses were performed in R (4.1.0) Project for Statistical Computing.
Gene classification
From the list of DEGs in the THC (vs vehicle) and/or WIN (vs vehicle) mice, we first identified genes with a human homolog using the “homologene” R package and then intersected these with the list of 2,511 consistent genes (see below, Virtual histology, for details). This was performed for both neuronal and non-neuronal gene lists obtained in the THC and WIN experiments, respectively (Extended Data Table 1-2).
Cannabis exposure
We started by grouping the adolescent males according to their history of cannabis exposure before the age of 16. Of the 475 adolescent males, 327 adolescents had not experimented with cannabis before age of 16 (mean age = 14.4 ± 1.7), and 140 adolescents did so at least once in their life (mean age = 16.0 ± 1.45; mean age of exposure = 13.8 ± 1.6). We grouped participants according to their exposure to cannabis before age 16 to be consistent with published studies that evaluated the relationship between cannabis use before age 16 with structural variations in the cerebral cortex (French et al., 2015; Scott et al., 2019).
Group difference in cortical thickness
Given the age difference between the two groups of adolescents (see above), we started by extracting residuals of the relationship between thickness and age for each cortical region. Using these residuals, we calculated the average (age-adjusted) cortical thickness for each of the 34 regions of the left hemisphere for each of the two groups. We evaluated only the cortical regions of the left hemisphere due to the availability of gene expression from the Allen Institute for Brain Science (Hawrylycz et al., 2012; see below, Virtual histology, for details). Next, we calculated the regional difference in group means between cannabis-exposed and not-exposed adolescents, thus creating an inter-regional profile of group differences in age-adjusted thickness. Since we are interested in the magnitude of the group differences related to cannabis exposure rather than the relative difference (most regions showed lower age-adjusted thickness in adolescents who did vs did not experiment with cannabis), we applied an absolute value function to the group difference values. Age-adjusted group differences are presented in Extended Data Tables 1-3 and 1-4.
Gene expression–group differences spatial relationships
For each gene from lists shown in Extended Data Table 1-2, we evaluated the inter-regional (spatial) relationship across the 34 regions between log2 intensity gene expression [from the Allen Institute for Brain Science (Hawrylycz et al., 2012); see below, Virtual histology, for details], and the age-adjusted inter-regional profile of group differences in thickness observed in SYS males across the same 34 regions. Genes that showed a nominally significant relationship with age-adjusted group differences across the 34 left cortical regions were used for virtual histology. Since the sole purpose of this step is to select genes for further analysis, no correction for multiple comparison was applied.
Virtual histology
The virtual histology procedure starts with the expression data for 20,737 genes, as assessed at the Allen Institute for Brain Science. These gene expression data were mapped on the 34 left cortical regions of the Desikan–Killiany atlas as described by French and Paus (2015). A two-stage filtering process was performed to select genes that show a consistent expression profile across the cerebral cortex. First, we selected genes that are expressed consistently across individuals if their inter-regional expression profile had a donor-to-median correlation >0.446. We then selected genes that are consistently expressed across the lifespan using the BrainSpan atlas (correlation >0.52), which led to a list of 2,511 “consistent” genes as described by Shin et al. (2018). Using a single-cell RNA sequencing and clustering analysis approach, cell-specific gene panels were developed for nine cell types as described by Zeisel et al. (2015). By intersecting those lists of cell-specific genes with the 2,511 consistent genes as a reference panel, the following number of genes for each cell type were obtained: astrocytes (n = 54), CA1 pyramidal neurons (n = 103), endothelial (n = 57), ependymal (n = 84), interneurons (n = 100), microglia (n = 48), mural (n = 25), oligodendrocytes (n = 60), and S1 pyramidal neurons (n = 73). To assess the spatial correlation between the inter-regional profile of a brain-related phenotype (e.g., cortical thickness) and cell-specific gene panels, we evaluated the distribution of correlation coefficients for each cell type using a false discovery rate (FDR) procedure and a resampling-based approach (n = 10,000) as described by Shin et al. (2018). Significance was assessed by determining whether the mean correlation coefficient was significantly different from 0.
Virtual histology in human adolescents
As explained above (Virtual histology), virtual histology facilitates neurobiological interpretation of the observed MRI-based phenotypes (Shin et al., 2018). To determine which (if any) cell-specific genes coexpress with DEGs that covary with cannabis-related group differences in cortical thickness (as revealed in the spatial relationships step described above), we grouped the age-adjusted nominally significant neuronal and non-neuronal DEGs and performed a virtual histology analysis in a ligand-specific (i.e., THC, WIN) way. Results of virtual histology were corrected for multiple comparisons using a FDR procedure (Benjamini and Hochberg, 1995) based on the number of genes correlated with gene panels. A visual scheme of the different steps of the analysis is shown in Extended Data Figure 6-1. All results are shown in Extended Data Tables 6-2–6-10.
Gene ontology enrichment analysis
In order to characterize the functional role of the ligand-specific neuronal, non-neuronal, and combined lists of DEGs, a gene ontology enrichment analysis for biological processes, cellular compartments, and molecular functions was performed using the “ClusterProfiler” package in R (Wu et al., 2021). This enrichment analysis was performed on lists of “consistent genes” for THC and WIN, and a Benjamini–Hochberg procedure was applied to each analysis to correct for multiple comparisons (Benjamini and Hochberg, 1995). Enrichment analysis of genes modified by THC and WIN in murine frontal cortex was carried out on non-respiratory genes, using the same procedure as above.
Results
Previous imaging studies reported variations in cortical thickness in human adolescents who experimented with cannabis (French et al., 2015; Chye et al., 2020; Albaugh et al., 2021, 2023; Owens et al., 2022). Since refinement of cortical circuits is a trans-species feature of adolescent cortical maturation (Rakic et al., 1994; Huttenlocher and Dabholkar, 1997; Jacobs et al., 1997; Petanjek et al., 2011; Pattwell et al., 2016), we reasoned that identification of cannabinoid-induced changes in the cerebral cortex of adolescent mice could reveal cellular and molecular processes underlying variations in the cerebral cortex observed in human adolescents who experimented with cannabis. We started to address this question by analyzing how repeated exposure to THC or to a synthetic cannabinoid influenced the transcriptome of different cell types in the frontal cortex of adolescent male mice.
Mice results: effects of cannabinoids on single-cell transcriptomes of murine frontal cortex
We analyzed the transcriptome of 28,360 cells recovered from the frontal cortex of 16 mice that received either vehicle, THC, or WIN from P35 to P39. RNASeq analysis revealed clusters that were represented as a UMAP and shown in Figure 1A. Using the DropViz database (Saunders et al., 2018) as a template for murine cortical transcriptomes, we then identified the cell types enriched in each cluster (Fig. 1A). Within each of them, we compared transcriptomes obtained from vehicle-treated mice with those recovered from mice treated with each cannabinoid. The comparison revealed 223 genes modified by THC and 172 modified by WIN. The list of genes modified by each treatment, along with the magnitude of treatment-induced changes in expression and FDR-corrected p values, is available in Extended Data Figures 1-2 and 1-3, and the information is summarized in Figure 1B. As shown therein, both ligands had distinct cellular effects (chi-square p value <2.2 × 10−16), THC predominantly modifying gene expression in non-neuronal cells while WIN primarily targeted Vglut1-expressing cells in layers 2/3 (L2/L3) and 5a (L5a). Both ligands induced intermediate changes in neurons expressing Gad1/2.
Distribution of clustered fronto-cortical cells based on mRNA profiles. The enrichment of genetic markers for the indicated clusters can be found in Extended Data Figure 1-1. A, Clusters expressing Vglut1 or Gad1/2 were identified as pyramidal (Pyr) neurons or interneurons, respectively. Cells that did not express these genes were considered as “non-neuronal” cells. Phenotypic profiles of cells in each category indicated in the figure. B, Frequency of DEGs by THC and WIN 55,212-2 across the different clusters. CGE, caudal ganglionic eminence; Gad1/2, glutamate decarboxylase; MGE, medial ganglionic eminence; Vglut1, type I vesicular glutamate transporter. The lists of DEGs with the statistics for THC and WIN can be found in Extended Data Figures 1-2 and 1-3, respectively.
Figure 1-1
Enrichment of genetic markers for indicated cluster. Download Figure 1-1, XLSX file.
Figure 1-2
Differentially expressed genes by THC in mice for each Cluster. Download Figure 1-2, XLSX file.
Figure 1-3
Differentially expressed genes by WIN in mice for each Cluster. Download Figure 1-3, XLSX file.
Cannabinoids and respiratory chain genes
The majority of genes modulated by THC and by WIN were specific to each ligand, with the exception of 42 genes that were influenced by both treatments. Among these genes common to both ligands, almost one-third (31%) coded for subunits in mitochondrial respiratory complexes. An additional set of genes coding for respiratory genes was modified distinctively by each cannabinoid (listed in Fig. 2A), making genes of the in mitochondrial respiratory chain the category most susceptible to modulation by either drug. All respiratory chain genes were downregulated by the treatment, the majority of them (90% for THC, 92% for WIN) corresponding to nuclear-coded subunits of proton-mobilizing complexes (CI, CIII, and CIV) and the ATP synthase (CV; Scarpulla, 2008). Variations in respiratory genes were most evident in L2/3 pyramidal neurons (clusters 1 and 6), where WIN effects predominated over those of THC (Fig. 2B). The converse was true for glial cells, where microglial (cluster 7) effects of THC were prevalent (Fig. 2C). GABA neurons were the least affected, with no evident difference between treatments (Fig. 2D). Expression effects induced by THC and WIN on respiratory genes and corresponding FDR-corrected p values are shown in Extended Data Figures 2-1 and 2-2, respectively.
THC and WIN 55,212-2 modified expression of respiratory genes in mouse frontal cortex. A, Genes commonly and distinctively modified by THC and/or WIN 55,212-2 are indicated. The expression changes induced by THC and WIN on respiratory genes can be found in Extended Data Figures 2-1 and 2-2, respectively. Comparison of genes differentially expressed by THC and WIN across the different respiratory chain complexes in (B) Vglut1+ pyramidal neurons, (C) non-neuronal cells, and (D) Gad1/2+ neurons. CI, CIII, and CIV, proton-mobilizing complexes; CV, ATP synthase, Ass. F, assembly factor; cluster numbers indicated at bottom of histograms. The cell-specific effects of THC and WIN on the expression of NDUFB8 subunits in microglia and layer 2/3 neurons of the prelimbic cortex can be found in Extended Data Figure 2-3. E, Mean cellular expression of Cnr1 and of genes coding for indicated downstream effectors in pyramidal (C1, C6, C8) and microglial (C7) clusters were obtained using the Average Expression function in Seurat. Inset shows expression ratios for βarr:Gαi/o subunits calculated using all Gαi/o versus βarr1/βarr2 effectors (top) or using cell-prevalent Gαi/o versus βarr effectors (bottom). Note: Expression profiles for Cnr1 and downstream effectors across all other clusters are provided in Extended Data Figure 2-5 for comparison. Only microglia displayed trace amounts of Cnr2. F, Concentration response curves for WIN and THC were obtained for the indicated effectors using a heterologous expression system (HEK293 cells). All responses were normalized to the endogenous ligand 2-AG which was the standard. Pharmacodynamic parameters can be found in Extended Data Figure 2-4. Data correspond to 3–4 independent experiments with four replicates per experiment. Curves generated by WIN and THC at each effector were compared by two-way ANOVA which showed an effect of drug as indicated in the figures, an effect of concentration (p ≤ 0.001), and an interaction (p ≤ 0.001), except for CB1-Gαi2 (p = 0.8700) and CB2-βarr2 (p = 0.4520), where no interaction was observed. C1–18 refer to clusters shown in Figure 1. C1, pyramidal neurons of layers 2/3; C6, pyramidal neurons of layers 2/3; C8, pyramidal neurons of layer 5a; C3, astrocytes; C7, microglia; C11, polydendrocytes; C10, medial ganglionic eminence (MGE) interneurons expressing parvalbumin; C12, MGE interneurons expressing somatostatin; C13, caudal ganglionic eminence (CGE) interneurons expressing vasoactive intestinal polypeptide; C16, CGE interneurons expressing neuropeptide Y; C18, neurons expressing glutamic acid decarboxylase (Gad).
Figure 2-1
Expression changes induced by THC on respiratory genes. The log of fold change is indicated by the grey tone (0: white; 1: dark grey). Download Figure 2-1, XLSX file.
Figure 2-2
Expression changes induced by WIN on respiratory genes. The log of fold change is indicated by the grey tone (0: white; 1: dark grey). Download Figure 2-2, XLSX file.
Figure 2-3
THC and WIN display cell specific effects on the expression of NDUFB8 subunits by microglia and Layer 2/3 neurons of the prelimbic cortex. An experimenter blind to treatment conditions counted in 3D the puncta corresponding to immunoreactivity (IR) of the NDUFB8 subunit I in Complex I (green, white arrows). (A) NDUFB8-IR within microglial silhouettes. Microglia were labelled with an antibody that recognizes transmembrane protein TMEM119 (red). (B) NDUFB8-IR within somatic volume of neurons revealed by NeuN (red). (C) Histogram shows the average number of puncta ± SEM counted in neurons and microglia of CTL and treated mice as indicated. The number of mice per condition is indicated at the bottom of each histogram. The number of microglia assessed per condition were CTL: 385 microglia, WIN: 209 microglia and THC: 260 microglia. The number of neurons assessed per condition were CTL: 1135 neurons, WIN: 300 neurons and THC: 425 neurons. Data were analyzed by Two-way ANOVA, which indicated and overall effect of cell type F(1, 34) = 18.24; p < 0.0001), no overall effect of treatment F(2, 34) = 2.044; p = 0.1451) and an interaction F(2, 34) = 5.037; p = 0.0121). Post-hoc comparisons using Tukey’s test showed a significant downregulation of microglial NDUFB8 by THC (p = 0.0279) but not WIN (p = 0.9996), and a trend for downregulation of neuronal NDUFB8 by WIN (p = 0.0857) but not THC (p = 0.9718). Download Figure 2-3, TIF file.
Figure 2-4
Pharmacodynamic parameters describing THC and WIN 55,212-2 responses at CB1 and CB2 receptors. Download Figure 2-4, XLSX file.
Figure 2-5
Mean cellular expression of Cnr1 and genes coding for indicated signaling effectors in all clusters identified. Note: Only microglia (C7) displayed trace amounts of Cnr2. Download Figure 2-5, TIF file.
Cell- and ligand-specific modulation of respiratory subunits was verified in TMEM119-positive microglia and in NeuN-labeled neurons within superficial layers 2/3 of the prelimbic cortex. In these cells, we quantified immunolabeled NDUFB8 particles in complex I. Comparison of the effects of THC and WIN on the two cell types was done by two-way ANOVA, which corroborated cell-specific effects for both ligands as indicated by an effect of cell type (F(1,34) = 18.24; p < 0.0001), and an interaction between cell type and treatment (F(2,34) = 5,037; p = 0.0121). Subsequent comparisons with Tukey's test revealed significant microglial downregulation of NDUFB8 subunits by THC (p = 0.0279), but not by WIN (p = 0.9996). Conversely, a trend toward downregulation of neuronal NDUFB8 was evident for WIN (p = 0.0857), but not for THC (p = 0.9718; Extended Data Fig. 2-3).
WIN and THC are respective full and partial agonists (Ryalls et al., 2024) at Gαi/o-coupled CB1 and CB2 receptors, which mediate neuronal and glial effects of cannabinoids (Azad et al., 2008; R. Chen et al., 2013; Hebert-Chatelain et al., 2016). While the modulation of genes that code for respiratory subunits in pyramidal neurons was consistent with the difference in ligand efficacy at CB1/CB2 receptors (i.e., a stronger transcriptomic response to WIN vs THC), the rank order of microglial responses was reversed (i.e., a stronger response to THC vs WIN). To gain insight into this apparent paradox, we considered cellular expression of signaling (Gαi/o subunits) and regulatory (β-arrestin: βarr) effectors immediately downstream of CB1/CB2 receptors. As shown in Figure 2E, βarr:Gαi/o ratios were approximately 5–7.5-fold higher in microglia than those in pyramidal neurons. Assessment of ligand efficacy at the different Gαi/o and βarr effectors confirmed that WIN and THC, respectively, behave as full and partial agonists, but their responses across these effectors were not proportional (Fig. 2F; Extended Data Fig. 2-4, curve parameters). THC concentrations that produced maximal activation of microglial-predominant Gαi2/Gαo transducers failed to recruit βarr2, which was also prevalent in these cells. In contrast, WIN concentrations that effectively activated Gαi2/Gαo were similar to those inducing βarr2 recruitment. Therefore, in microglia where βarrs were more prevalent than in neurons, WIN could limit its own Gαi2/GαoB signaling while THC could not, transcriptomic responses being consistent with ligand efficacy to recruit βarr. In pyramidal neurons, predominant Gα subunits exceeded βarr expression by a greater margin than in microglia. Therefore, in these cells, WIN could signal via the predominant Gαo/Gαi1 targets before triggering βarr desensitization. Hence, in superficial pyramidal neurons, gene regulation was consistent with ligand efficacy to induce Gαo/Gαi1 activation.
THC treatment and nonrespiratory genes
Non-respiratory genes modified by THC were analyzed for enrichment of ontology terms, revealing a prevalence of categories associated with synapse (GO:0098984), its organization (GO:0050808), presynaptic (GO:0098793; GO:0030133; GO:0070382; GO:0070382) and post-synaptic (GO:0099572; GO:0014069) components, and excitatory/glutamatergic synapses (GO:0098978; GO:0032279); 33% of the genes associated with these terms were also enriched in categories related to neuronal development (GO:0007409) and differentiation (GO:0021953), while additional terms functionally associated with synaptic transmission such as amine transport (GO:0015837), calcium detection (GO:0005513), calmodulin (GO:0005516), and glutamate receptor binding (GO:0035254) were also over-represented among DEGs induced by THC (Fig. 3A). The genes populating these categories, the cell types where they were modified, and the magnitude of expression changes along with corresponding statistics appear in Extended Data Figure 3-1. Synapse-related genes modified by THC are summarized in Figure 3B, which shows that the predominant effect of treatment in mature neurons was the upregulation of glutamatergic synapse genes. Conversely, in glia and mitotic neural cells, downregulation was prevalent, affecting synaptic and cytoskeletal genes.
THC induces transcriptomic and structural changes to glutamatergic synapses in the frontal cortex of adolescent male mice. A, Gene ontology enrichment results for 203 non-respiratory genes differentially expressed by THC. Statistics can be found in Extended Data Figure 3-1. Categories associated with synapse-related terms are highlighted. B, Cellular distribution of THC DEGs associated with synapse function. In bold, genes (mentioned in text) that are functionally or structurally associated to glutamatergic synapses. C, Effect of THC on spine density of apical and basal dendrites in pyramidal neurons of layer 2/3 of prelimbic cortex. Histograms show mean ± SEM of five mice/condition. Five to six neurons were analyzed per mouse; spines on primary and secondary apical branches and on 4–7 basal dendrites/neuron were considered. Results corresponding to Mann–Whitney comparisons are indicated in the figure.
Figure 3-1
Enrichment results for non-respiratory genes for differentially expressed genes by THC. Download Figure 3-1, XLSX file.
In mature neurons, THC upregulated genes coding for AMPA subunits (Gria2), for proteins that organize AMPA and NMDA receptors at the synapse (Adam22; Hsia et al., 2019), for proteins that can regulate membrane expression of these receptors (Ptpn5; Kamceva et al., 2016), and genes whose products contribute to normal plasticity of glutamatergic synapses (Igf1, Ryr2; Dyer et al., 2016; Lobos et al., 2021). Other gene products that THC upregulated in mature neurons were previously associated with excitatory synapse dysfunction in Alzheimer's disease (Apoe, Glul; S. R. Robinson, 2000; Y. Chen et al., 2010). Genes coding for adhesion proteins at the neuron's membrane (Flrt2) and for cytoskeleton-organizing proteins that support connectivity, dendrite complexity, and structural plasticity of glutamatergic synapses (Dst, Ctnna2) were also upregulated in these cells (Abe et al., 2004; Bhanot et al., 2011; Valdés-Undurraga et al., 2023) along with genes that code for cytoskeletal proteins (Tub1a, Tubb5). The Ywhah gene that codes for the 14-3-3 protein ε, which promotes membrane expression of NMDA2 subunits (Chung et al., 2015), was exceptionally downregulated in these neurons.
In oligodendrocytes, THC modulated genes that code for extracellular anchored proteins (Cntn2, Nrn1), which stabilize dendritic spines and promote connectivity via trans-synaptic organization of AMPARs (Boudkkazi et al., 2023). These cells and astrocytes also displayed upregulation of genes coding for extracellular matrix components that participate in adolescent maturation of cortical synapses (Bcan, Vcan, Hapln2; Härtig et al., 1992; Mauney et al., 2013). Astrocytes additionally displayed upregulation of genes that code for microtubule-organizing proteins (Dclk1; K. E. Murphy et al., 2023) or microtubule constituents (Tubb2a). In contrast, the Stnm2 gene, which codes for a regulator of microtubule stability and cytoskeletal dynamics (Tortoriello et al., 2014), was downregulated across the three major glial subtypes, along with genes for calcium-sensing proteins (Syt1, Calm1). Genes coding for proteins normally expressed in presynaptic terminals (e.g., Sv2a, SNARE components Syt1, Syt11, Snap25) were downregulated by THC in all glial subtypes, with microglia displaying the highest number of SNARE genes downregulated by treatment. In these cells THC also downregulated the Hspa8 gene, whose product may function as a regulator of axonal protein synthesis (Tanaka et al., 2024), as well as genes that code for AMPA subunits (Gria2) and for signaling integrators (Nrgn, Arpp21) implicated in spine morphogenesis (Rehfeld et al., 2018). Genes that code for the lysosomal protein Laptm5 and for the C1qa complement that participates in microglial pruning of glutamatergic synapses (Stevens et al., 2007; Cong et al., 2020) were two microglia-intrinsic genes upregulated by THC. In addition to synaptic terms, Laptm5, C1qa, and SNARE-coding genes (Syt1, Syt11, Stxbp1) were also enriched in terms associated with tagging of inactive synapses for removal (phosphatidylserine binding GO:0001786) and with phagocytosis of neuronal components (regulation of cell death GO:1901214).
Transcriptomic results suggest a scenario where phenotypic changes to glial cells and mature neurons converge to modify glutamatergic synapses. Structural remodeling of glutamatergic synapses is a hallmark of adolescent maturation of rodent cortex (Juraska et al., 2020; Pöpplau et al., 2024). Hence, to determine if the same THC treatment that modified excitatory synapse genes also modified the structure of these synapses, we evaluated spine density on L2/3 pyramidal neurons of the prelimbic cortex (Fig. 3C). Mann–Whitney comparison of spine numbers in apical (p = 0.0159) and basal (p = 0.0079) dendrites showed ∼60% decrease in number of spines in THC-treated mice.
In layer 2/3 pyramidal neurons of adolescent prefrontal cortex, spine density is closely interrelated with dendritic complexity (Pöpplau et al., 2024). Hence, we used Sholl analysis to monitor dendritic morphology in these neurons (Fig. 4). Mann–Whitney comparison of dendritic length (p = 0.0079), number of branch points (p = 0.0079), and branching depth (p = 0.0079) indicated a reduction in overall dendritic complexity in THC-treated mice. Two-way ANOVA analysis of the number of dendritic intersections with concentric circles at increasing distance from the cell soma further revealed that THC reduced the complexity of both apical (p < 0.01) and basal (p < 0.001) dendritic trees.
THC decreases the complexity of dendritic arborization of L2/3 pyramidal neurons in prefrontal cortex of adolescent male mice. A, Characteristic morphology of GFP-labeled L2/3 PYR neurons in prelimbic cortex of adolescent Thy-1 male mice that received vehicle (CTL) or THC. Histograms display (B) total dendrite length, (C) number of branch points, and (D) branch depth. They correspond to mean ± SEM of five mice/condition. Four neurons displaying apical and basal dendrites were reconstituted per mouse. Results corresponding to Mann–Whitney comparisons are indicated in the figure. Intersections of (E) apical dendrites within a 360 μm radius from soma and (F) basal dendrites within 200 μm radius from soma from the soma of prelimbic L2/3 pyramidal neurons in CTL and THC-treated mice were compared by two-way ANOVA. Apical: showed and effect of distance from soma F(17,676) = 25.73, p < 0.001; an effect of treatment F(1,676) = 58.03, p < 0.001 and an interaction F(17,676) = 2.144, p = 0.005. Basal showed and effect of distance from soma F(9,676) = 74.46, p < 0.001; an effect of treatment F(1,380) = 187.4, p < 0.001 and an interaction F(9,380) = 12.65, p < 0.001. Post hoc comparisons with Sidak's test: *p < 0.01; **p < 0.001; ***p < 0.0001.
WIN treatment and nonrespiratory genes
Ontology analysis of non-respiratory genes modified by WIN also revealed a prevalence of terms associated with glutamatergic synapse (GO:0032279; GO:0098978), including postsynaptic components (GO:0014069; GO:0043197) and cellular processes implicated in learning, memory (GO:0007611), and cognition (GO:0050890; Fig. 5A).
WIN 55,212-2 induces transcriptomic changes to glutamatergic synapses in the frontal cortex of adolescent male mice. A, Gene ontology enrichment results for 149 non-respiratory genes differentially expressed by WIN. Statistics can be found in Extended Data Figure 5-1. Categories associated with synapse-related terms are highlighted. B, Cellular distribution of WIN DEGs associated with synapse function. In bold, genes (mentioned in text) that are functionally or structurally associated to glutamatergic synapses.
Figure 5-1
Enrichment results for non-respiratory genes for differentially expressed genes by WIN. Download Figure 5-1, XLSX file.
WIN primarily targeted L2/3 pyramidal neurons, where it induced gene upregulation (Fig. 5B; Extended Data Fig. 5-1). Therein WIN upregulated Prkg1, which codes for cGMP-dependent kinase 1. This signaling effector is a component of the signaling pathway involved in WIN-mediated disruption of long-term potentiation (LTP) and impaired memory consolidation (Maroso et al., 2016). Other genes upregulated in L2/3 pyramidal neurons were also related to impaired cognition or altered plasticity. Thus, Apoe (Y. Chen et al., 2010) and Synj-1 (which codes for the PiP2 phosphatase Synaptojanin-1; Miranda et al., 2018) were previously associated with memory impairment and cognitive dysfunction in Alzheimer's disease, while the chromatin-remodeling enzyme that is coded by Chd4 is known to repress long-term neuronal plasticity by suppressing transcriptional availability of the genome (Goodman et al., 2020). Grim5 and Sorcs3, whose products were also upregulated by WIN, code for proteins that participate in long-term depression (LTD) and memory extinction. Thus, mGluR5 receptors that are coded by Grim5 underlie memory impairment in animal models of Alzheimer's disease (Hamilton et al., 2016), and their activation induces LTD via production of endocannabinoids through a mechanism that is also shared by WIN (Izumi and Zorumski, 2012). Moreover, through modulation of AMPA receptor trafficking at the postsynaptic density, upregulated Sorcs3 products may also participate in LTD and memory extinction (Breiderhoff et al., 2013).
WIN induced additional changes in the expression of genes that code for glutamate-binding membrane proteins, pointing to the functional reorganization of synapses on L2/3 pyramidal neurons. In particular, WIN induced upregulation of Gria1/Gria3 genes that code for AMPA1/AMPA3 subunits and the ectopic expression of glial glutamate transporter genes Slc1a3 (EAAT1) and Slc1a2 (EAAT2) within this neuronal population. Downregulation of Ywhah, a gene associated with membrane expression of NMDA subunits (Chung et al., 2015), is also consistent with functional reorganization of excitatory synapses on these pyramidal neurons of superficial layers. Upregulation of genes coding for the syntenin-interacting protein Rheb and for the pseudo-kinase CamkV additionally suggest possible structural changes to the synapse, as these proteins respectively participate in spine morphogenesis (Sugiura et al., 2015) and stability (Liang et al., 2016).
In other neurons, WIN modified the expression of genes that code for the Na+/K+ ATPase (Atp1a1, Atp1b1, Atp1a2) and for K+ channels (Kcnd3, Kcnj3), resulting in the enrichment of transmembrane transport terms (e.g., GO:0099587, GO:0022853) and suggesting changes to membrane excitability (Fig. 5; Extended Data Fig. 5-1).
Human results: exposure to Δ-9-tetrahydrocannabinol
From mice to humans
A total of 127 unique genes were differentially expressed in mouse neuronal cells (i.e., cells expressing Vglut1 or Gad1/2) following exposure to THC, and 99/127 had a human homolog; 20/99 of the (human) genes intersected with the list of 2,511 “consistent genes” (i.e., with consistent inter-regional expressed profiles in the human cerebral cortex (Shin et al., 2018). For non-neuronal cells, 117 mouse genes were differentially expressed, and 100/117 had a human homolog; 21/100 of the (human) genes intersected with the list of 2,511 consistent genes.
Gene expression–thickness relationships
When examining the spatial correlation between gene expression and group differences in cortical thickness across the cerebral cortex, a total of 7/20 neuronal DEGs by THC and 8/21 non-neuronal DEGs by THC were associated (positively or negatively) with cannabis-related group differences in thickness after adjusting for age. Of the 13 (unique) genes, the following 9 showed a positive age-adjusted relationship between their spatial expression and the cannabis-related group differences in age-adjusted thickness: NTSR2 (r = 0.45; p = 0.007), TUBB2A (r = 0.54; p = 9.92 × 10−4), NFIB (r = 0.41; p = 0.015), UCHL1 (r = 0.45; p = 0.007), APOE (r = 0.51; p = 0.002), ENC1 (r = 0.59; p = 2.51 × 10−4), YWHAH (r = 0.51; p = 0.002), LINGO1 (r = 0.58; p = 3.66 × 10−4), and BCAN (r = 0.36; p = 0.039). On the other hand, the following four nominally significant genes showed a negative (spatial) relationship with the cannabis-related group differences in age-adjusted thickness: NDUFC1 (r = −0.36; p = 0.037), PCP4 (r = −0.44; p = 0.009), SV2A (r = −0.39; p = 0.021), and ATP2B2 (r = −0.37; p = 0.030). Statistics for all THC-related consistent neuronal genes are presented in Extended Data Figure 6-2.
Spatial correlation results for inconsistent THC genes are presented in Extended Data Figure 6-3. Thirty genes were nominally significant in their association with age-adjusted group differences profiles, from which the following 15 showed a positive relationship: LPL (r = 0.42; p = 0.014), GLUL (r = 0.36; p = 0.038), SOX2 (r = 0.35; p = 0.042), TOX3 (r = 0.35; p = 0.041), CST3 (r = 0.44; p = 0.010), CDCA7 (r = 0.34; p = 0.050), SLC30A9 (r = 0.39; p = 0.024), HEPACAM (r = 0.36; p = 0.037), NNAT (r = 0.38; p = 0.028), PRODH (r = 0.45; p = 0.008), TUBB (r = 0.54; p = 0.001), ADI1 (r = 0.35; p = 0.042), IVD (r = 0.57; p = 0.001), ATP6V1H (r = 0.43; p = 0.012), and UQCRFS1 (r = 0.35; p = 0.041). On the other hand, the following 15 genes showed a negative age-adjusted spatial relationship: ZBTB43 (r = −0.39; p = 0.023), LENG8 (r = −0.42; p = 0.013), TIA1 (r = −0.45; p = 0.008), EZH2 (r = −0.46; p = 0.008), AHI1 (r = −0.51; p = 0.002), ATP1B1 (r = −0.43; p = 0.011), VCAN (r = −0.44; p = 0.010), MDH1 (r = −0.47; p = 0.005), DCLK1 (r = −0.36; p = 0.035), DDX5 (r = −0.34; p = 0.050), NFKBID (r = −0.41; p = 0.016), NDUFA4 (r = −0.50; p = 0.003), DTD1 (r = −0.40; p = 0.021), SLC25A4 (r = −0.36; p = 0.037), and PLEKHH1 (r = −0.37; p = 0.030).
Virtual histology
Figure 6A shows virtual histology results for the 13 DEGs by THC that (1) passed the two-stage consistency filtering and (2) are nominally associated with group differences in age-adjusted thickness (statistics provided in Extended Data Fig. 6-4). As shown in this figure, the nine genes with positive (spatial) relations between their expression and the magnitude of age-adjusted cannabis-related group differences also show a positive coexpression—across the 34 cortical regions—with astrocytes, CA1 pyramidal, and microglia genes (except YMHAH for the microglia panel), with Pearson's correlation coefficients ranging between 0 and 0.40. That is, regions with high expression of astrocytes, CA1 pyramidal, and microglia genes show large cannabis-related group differences in cortical thickness.
Spatial coexpression between nine cell-specific gene panels and (A) 13 nominally significant THC-related genes that passed the two-stage consistency filtering (FDR-corrected for 13 × 9 observations) and (B) 18 nominally significant WIN-related DEG that passed our two-stage consistency filtering (FDR-corrected for 18 × 9 observations). A visual schematic of the analytic pipeline can be found in Extended Data Figure 6-1. The spatial correlation coefficients between group differences in cortical thickness and gene expression for THC-related consistent and inconsistent DEGs can be found in Extended Data Figures 6-2 and 6-3, respectively. Statistics for virtual histology for THC-related consistent and inconsistent DEGs can be found in Extended Data Figures 6-4 and 6-5, respectively. The spatial correlation coefficients between group differences in cortical thickness and gene expression for WIN-related consistent and inconsistent DEGs can be found in Extended Data Figures 6-6 and 6-7, respectively. Statistics for virtual histology for WIN-related consistent and inconsistent DEGs can be found in Extended Data Figures 6-8 and 6-9, respectively. Inter-regional coexpression of THC and WIN nominally significant genes can be found in Extended Data Figure 6-10.
Figure 6-1
Virtual histology procedure. Reprinted from Navarri et al., 2023 (PMID: 36963244). Download Figure 6-1, TIF file.
Figure 6-2
Correlation between inter-regional profiles of group differences in cortical thickness (age-adjusted) and gene expression for 37 THC-related consistent differentially expressed. Download Figure 6-2, XLSX file.
Figure 6-3
Correlation between inter-regional profiles of group differences in cortical thickness (age-adjusted) and gene expression for 135 (out of the 149 genes) THC-related DEG that did not pass our two-stage consistency filtering and were used for gene ontology analysis. The following 14 genes were not in Human Allen Brain Database: BMERB1, COLGALT1, HSP90AA1, TRANK1, PHYHIPL, POLR2 M, LOC102724788, NTM, COX14, ZFHX3, ARPP21, GPCPD1, SCHIP1, CELF4. Download Figure 6-3, XLSX file.
Figure 6-4
Inter-regional co-expression of 13 consistent nominally significant THC-related DEG with nine cell-specific genes (corrected for 117 comparisons using a False Discovery Rate procedure). The group difference sign refers to the directionality of the spatial coexpression relationship between absolute age-adjusted cortical thickness (between cannabis-exposed adolescents males and non-exposed adolescent males from the Saguenay Youth Study cohort) and gene expression from the Allen Brain Atlas. Download Figure 6-4, XLSX file.
Figure 6-5
Inter-regional co-expression of 30 inconsistent nominally significant THC-related DEG with nine cell-specific genes (corrected for 270 comparisons using a False Discovery Rate procedure). The group difference sign refers to the directionality of the spatial coexpression relationship between absolute age-adjusted cortical thickness (between cannabis-exposed adolescents males and non-exposed adolescent males from the Saguenay Youth Study cohort) and gene expression from the Allen Brain Atlas. Download Figure 6-5, XLSX file.
Figure 6-6
Pearson correlation coefficients between inter-regional profiles of group differences in cortical thickness (age-adjusted) and gene expression for 33 WIN-related DEG that passed our two-stage consistency filtering and were used for gene ontology analysis. Download Figure 6-6, XLSX file.
Figure 6-7
Correlation between inter-regional profiles of group differences in cortical thickness (age-adjusted) and gene expression for 96 (out of the 108 genes) WIN-related DEG that passed our two-stage consistency filtering and were used for gene ontology analysis. The following 12 genes were not in Human Allen Brain Database: ATP5PF, UQCR10, TRAPPC13, PKM, PCDHA2, BMERB1, NAA15, ATP5PO, UQCC2, ATP5MF, SCHIP1, and ETNPPL. Download Figure 6-7, XLSX file.
Figure 6-8
Inter-regional co-expression of 18 consistent nominally significant WIN-related DEG with nine cell-specific genes (corrected for 162 comparisons using a False Discovery Rate procedure). The group difference sign refers to the directionality of the spatial coexpression relationship between absolute age-adjusted cortical thickness (between cannabis-exposed adolescents males and non-exposed adolescent males from the Saguenay Youth Study cohort) and gene expression from the Allen Brain Atlas. Download Figure 6-8, XLSX file.
Figure 6-9
Inter-regional co-expression of 21 inconsistent nominally significant WIN-related DEG with nine cell-specific genes (corrected for 162 comparisons using a False Discovery Rate procedure). The group difference sign refers to the directionality of the spatial coexpression relationship between absolute age-adjusted cortical thickness (between cannabis-exposed adolescents males and non-exposed adolescent males from the Saguenay Youth Study cohort) and gene expression from the Allen Brain Atlas. Download Figure 6-9, XLSX file.
Figure 6-10
(A) Inter-regional co-expression of 30 nominally significant THC-related DEGs that did not pass the two-stage consistency filtering with cell-specific genes. FDR-corrected for 270 comparisons (30 genes × 9 cell panels). (B) Inter-regional co-expression of 21 nominally significant WIN-related DEGs that did not pass the two-stage consistency filtering with cell-specific genes. FDR-corrected for 189 comparisons (21 genes × 9 cell panels). Download Figure 6-10, TIF file.
Since we observed (predominantly) lower cortical thickness across all 34 regions in male adolescents who had experimented with cannabis (vs those who had not), lower expression of these genes (and the coexpressed panels of cell-specific genes) would be observed in the cerebral cortex of the male adolescents who experimented with cannabis. These results may be interpreted as a reduction in the dendritic arborization of CA1 pyramidal neurons and of the glial cells that contribute to the homeostatic and metabolic activity of this type of pyramidal neurons.
Conversely, the four genes with negative (spatial) relation between their expression and the age-adjusted group differences in thickness show a negative coexpression with astrocyte, CA1 pyramidal, and microglia panels (except NDUFC1 and PCP4 for the microglia panel), with r correlation coefficients ranging between −0.40 and 0. That is, regions with low expression of astrocytes, CA1 pyramidal, and microglia genes show large cannabis-related group differences in cortical thickness. Thus, higher gene expression of these cell-specific genes might be observed in the thinner regions of the cerebral cortex in male adolescents who experimented with cannabis.
BCAN inter-regional profile negatively correlates with the S1 pyramidal panel (r = −0.18, pFDR = 0.001). Virtual histology results for the 30 inconsistent genes are shown in Extended Data Figures 6-5 and 6-10A.
Gene ontology enrichment analysis
Figure 7A shows gene enrichment analysis of all 37 DEGs identified in the THC experiment that passed the two-stage filtering process (Extended Data Figs. 7-1, 7-2). This analysis revealed involvement of neurodevelopmental processes (regulation of neuron projection development; pBH = 6.893 × 10−4), synaptic activity (modulation of chemical synaptic transmission; pBH = 3.472 × 10−3, synaptic vesicle priming: pBH = 1.007 × 10−3) and synaptic plasticity (regulation of synaptic plasticity; pBH = 4.296 × 10−3). DEGs that did not pass the two-stage filtering process are enriched in mitochondria-related processes (ATP metabolic process; pBH = 3.154 × 10−6; Extended Data Figs. 7-3, 7-4).
Gene ontology enrichment results for the genes that were differentially expressed in mice and passed our two-stage filtering process in humans when exposed to a cannabinoid receptor agonist. Enrichment for biological process, cellular compartment, and molecular function for 37 THC-related consistent genes for Δ-9-tetrahydrocannabinol (A) and 33 consistent genes for WIN 55,212-2 (B). All statistics for the consistent THC-related genes, both pooled and separated into neuronal and non-neuronal categories, are available in Extended Data Figures 7-1 and 7-2, respectively. Statistics for GO analysis for inconsistent THC-related DEGs, both pooled and separated into neuronal and non-neuronal categories, are available in Extended Data Figures 7-3 and 7-4, respectively. Statistics for the consistent WIN-related genes, both pooled and separated into neuronal and non-neuronal categories, are available in Extended Data Figures 7-5 and 7-6, respectively. Statistics for GO analysis for inconsistent WIN-related DEGs, both pooled and separated into neuronal and non-neuronal categories, are available in Extended Data Figures 7-7 and 7-8, respectively.
Figure 7-1
Gene ontology results for lists of consistent THC-related differentially expressed genes for biological process (BP), cellular compartments (CC) and molecular function (MF). Gene ratio indicates the enrichment for a given term and p-values underwent a Benjamini-Hochberg procedure. Download Figure 7-1, XLSX file.
Figure 7-2
Gene ontology results for lists of THC-related differentially expressed consistent genes for biological process (BP), cellular compartments (CC) and molecular function (MF). Gene ratio indicates the enrichment for a given term and p-values are corrected following a Benjamini-Hochberg procedure. Top row refers to consistent neuronal genes (n = 20), middle row to consistent non-neuronal genes (n = 21), and the bottom row to all consistent genes (n = 37). Download Figure 7-2, TIF file.
Figure 7-3
Gene ontology results for lists of inconsistent THC-related differentially expressed genes for biological process (BP), cellular compartments (CC) and molecular function (MF). Gene ratio indicates the enrichment for a given term and p-values underwent a Benjamini-Hochberg procedure. Download Figure 7-3, XLSX file.
Figure 7-4
Gene ontology results for lists of inconsistent THC-related differentially expressed genes for biological process (BP), cellular compartments (CC) and molecular function (MF). Gene ratio indicates the enrichment for a given term and p-values are corrected following a Benjamini-Hochberg procedure. Top row refers to consistent neuronal genes (n = 79), middle row to consistent non-neuronal genes (n = 79), and the bottom row to all consistent genes (n = 149). Download Figure 7-4, TIF file.
Figure 7-5
Gene ontology results for lists of consistent WIN-related differentially expressed genes for biological process (BP), cellular compartments (CC) and molecular function (MF). Gene ratio indicates the enrichment for a given term and p-values underwent a Benjamini-Hochberg procedure. Download Figure 7-5, XLSX file.
Figure 7-6
Gene ontology results for lists of WIN-related differentially expressed consistent genes for biological process (BP), cellular compartments (CC) and molecular function (MF). Gene ratio indicates the enrichment for a given term and p-values corrected following a Benjamini-Hochberg procedure. Top row refers to consistent neuronal genes (n = 31), middle row to consistent non-neuronal genes (n = 2), and the bottom row to all consistent genes (n = 33). Download Figure 7-6, TIF file.
Figure 7-7
Gene ontology results for lists of inconsistent WIN-related differentially expressed genes for biological process (BP), cellular compartments (CC) and molecular function (MF). Gene ratio indicates the enrichment for a given term and p-values underwent a Benjamini-Hochberg procedure. Download Figure 7-7, XLSX file.
Figure 7-8
Gene ontology results for lists of WIN-related differentially expressed inconsistent genes for biological process (BP), cellular compartments (CC) and molecular function (MF). Gene ratio indicates the enrichment for a given term and p-values are corrected following a Benjamini-Hochberg procedure. Top row refers to consistent neuronal genes (n = 95), middle row to consistent non-neuronal genes (n = 16), and the bottom row to all consistent genes (n = 108). Biological process and molecular function analysis could not be computed for the inconsistent non-neuronal WIN-related DEG given the gene list. Download Figure 7-8, TIF file.
Exposure to WIN
Considering the emergence of synthetic cannabinoids, we then asked whether WIN 55,212-2, a CB1R/CB2R ligand with an efficacy profile similar to that of recreationally used synthetic cannabinoids (Soethoudt et al., 2017; Ryalls et al., 2024), could also be associated with systems-level measures of the maturation of the cerebral cortex during adolescence. To this end, we carried out a similar analysis for DEGs by WIN in mice.
From mice to humans
A total of 151 unique neuronal genes were differentially expressed following WIN exposure in mice, among which 126/151 had a human homolog; 31/126 of the (human) genes intersected with the list of 2,511 consistent genes. For non-neuronal genes, 27 mouse genes were differentially expressed and 18/27 had a human homolog; only 2 genes with a human homolog intersected with the list of 2,511 consistent genes. In total, 33 DEGs following WIN exposure passed the two-stage filtering process and were used to enrichment analysis.
Gene expression–thickness relationships
When examining the spatial correlation between gene expression and group differences in cortical thickness across the cerebral cortex, a total of 16/31 neuronal genes and 2/2 non-neuronal genes (18 unique genes) associated with WIN exposure were nominally associated (positively or negatively) with cannabis-related group differences in thickness after adjusting for age. Of these 18 genes, the following 14 showed a positive relationship between their spatial expression and the cannabis-related group differences in age-adjusted thickness: TUBB2A (r = 0.54; p = 0.001), NAV1 (r = 0.40; p = 0.020), UCHL1 (r = 0.45; p = 0.007), PCDH17 (r = 0.36; p= 0.037), APOE (r = 0.51; p = 0.023), CRYM (r = 0.43; p = 0.011), PGM2L1 (r = 0.49; p = 0.004), PDE1A (r = 0.36; p = 0.040), CAMKV (r = 0.44; p = 0.009), PRKCG (r = 0.55; p = 0.001), YWHAH (r = 0.51; p = 0.002), GRIA1 (r = 0.46; p = 0.007), DGKB (r = 0.36; p = 0.040), and NKAIN4 (r = 0.37; p = 0.030). On the other hand, the following four nominally significant genes showed a negative relationship with the cannabis-related group differences in age-adjusted thickness: NDUFC1 (r = −0.36; p = 0.040), ATP1A1 (r = −0.39; p = 0.020), PCP4 (r = −0.44; p = 0.009), and ACSL3 (r = −0.48; p = 0.004). Statistics for all genes is presented in Extended Data Figure 6-6.
Spatial correlation results for inconsistent WIN genes are presented in Extended Data Figure 6-7. Twenty-one genes were nominally significant in their association with age-adjusted group differences profiles, from which the following 11 showed a positive age-adjusted spatial relationship: ATP1A2 (r = 0.44; p = 0.010), SLC1A2 (r = 0.55; p = 0.001), LOC728392 (r = 0.37; p = 0.033), KIF1B (r = 0.37; p = 0.030), SLC1A3 (r = 0.49; p = 0.003), CST3 (r = 0.44; p = 0.010), PHACTR3 (r = 0.44; p = 0.010), TTYH1 (r = 0.55; p = 0.002), TUBB (r = 0.54; p = 0.001), ZCCHC18 (r = 0.36; p = 0.035), and SCG3 (r = 0.41; p = 0.06). On the other hand, the following 10 genes showed a negative relationship: AHI1 (r = −0.51; p = 0.002), ATP1B1 (r = −0.43; p = 0.011), TOP1 (r = −0.34; p = 0.046), EGR1 (r = −0.44; p = 0.010), ATP6V1G1 (r = −0.44; p = 0.02), GNG13 (r = −0.38; p = 0.029), COX5A (r = −0.35; p = 0.041), DDX46 (r = −0.45; p = 0.008), PPM1A (r = −0.36; p = 0.037), and NDUFA4 (r = −0.50; p = 0.003).
Virtual histology
Figure 6B shows results of virtual histology for the 18 DEGs by WIN that (1) passed the two-stage consistency filtering and (2) are nominally associated with group differences in age-adjusted thickness (statistics provided in Extended Data Fig. 6-8). As shown in this figure, the 14 genes with positive (spatial) relations between their expression and the magnitude of age-adjusted cannabis-related group differences also show a positive coexpression—across the 34 cortical regions—with astrocytes, CA1 pyramidal, and microglia genes (with the exceptions of DGKB for the astrocyte panel and DGKB, YWHAH, and PDE1A for the microglia panel), with Pearson's correlation coefficients ranging between 0 and 0.40. That is, regions with high expression of astrocytes, CA1 pyramidal, and microglia genes show large cannabis-related group differences in cortical thickness.
Conversely, the four genes with negative (spatial) relation between their expression and the age-adjusted group differences in thickness show a negative coexpression with astrocyte, CA1 pyramidal, and microglia panels with r correlation coefficients ranging between −0.40 and 0 (with the exception of NDUFC1 and PCP4 for the microglia panel). That is, regions with low expression of astrocyte, CA1 pyramidal, and microglia genes show large cannabis-related group differences in cortical thickness.
ATP1A1 inter-regional profile positively correlates with the S1 pyramidal panel (r = 0.16; pFDR = 0.047) and NAV1 inter-regional profile negatively correlates with the S1 pyramidal panel (r = 0.17; pFDR = 0.031). Virtual histology results for the 21 inconsistent genes are visualized in Extended Data Figures 6-9 and 6-10B.
Gene ontology enrichment analysis
Figure 7B shows gene enrichment analysis of all 33 DEGs identified in the WIN experiment that passed the two-stage filtering process (Extended Data Figs. 7-5, 7-6). This analysis revealed involvement of cognitive processes (learning or memory; pBH = 1.372 × 10−4, cognition; pBH = 2.039 × 10−4). DEGs that did not pass the two-stage filtering process, on the other hand, are involved in mitochondria-related processes (ATP metabolic process; pBH = 1.088 × 10−14; Extended Data Figs. 7-7, 7-8).
Discussion
Capitalizing on the identification of cortical genes that in mice were sensitive to THC and to WIN 55,212-2, we investigated the cellular and molecular underpinnings of group differences in cortical thickness as a function of cannabis exposure (before 16 years of age) in human adolescents. We found that cannabis-related group differences in cortical thickness correlated with the regional expression of human genes whose homologs were modulated by THC or WIN in the murine frontal cortex and established that genetic markers for pyramidal neurons that are enriched for dendrite-related genes (CA1 type), astrocytes, and microglia codistributed with cannabinoid-sensitive human homologs that correlated with regional differences in cortical thickness.
Virtual histology analysis revealed that DEGs by THC and WIN, which showed a positive correlation with the inter-regional profile of group differences in cortical thickness, correlated positively with genes specific to astrocytes, microglia, and the type of pyramidal cells (CA1) that are enriched for dendrite-related genes (Shin et al., 2018). The converse was true for the genes that correlated negatively with group differences in cortical thickness. Such consistency supports the premise that cannabis-induced differences in cortical thickness are related to the expression of cannabinoid-sensitive genes specifically associated to the cellular components in the tripartite pyramidal synapse and to microglia, which remodel these synapses over the course of circuit maturation in the adolescent cerebral cortex (Mallya et al., 2019; Cong et al., 2020; Pöpplau et al., 2024). We have previously shown that astrocytes, microglia, and (CA1) pyramidal neurons are associated with inter-regional differences in thickness during the course of normal cortical maturation (Shin et al., 2018). The same cell types were also associated with abnormal cortical thinning observed in psychiatric disorders (Patel et al., 2021) and during aging (Vidal-Pineiro et al., 2020). Now we provide evidence that the same cell types are also associated with cortical effects of cannabis in the adolescent brain. It is important to note, however, that the THC- and WIN-sensitive genes were distinct. Thus, although cellular relationships with cortical thickness are similar for both cannabinoids, the actual murine genes modified by both ligands were distinct, making it difficult to conclude that their influence on glutamatergic synapses involved similar mechanisms. The effects of these treatments in murine cortex provide insight as to how THC and WIN may affect adolescent glutamatergic synapses via different cortical cell types.
Glutamatergic synapse genes that were sensitive to THC were primarily found in glia of mouse frontal cortex, with pyramidal cells and GABA neurons contributing additional DEGs. In the adolescent cortex, inactive glutamatergic synapses are normally removed through microglial engulfment of pre- and postsynaptic elements (Mallya et al., 2019; Pöpplau et al., 2024), explaining why genes that are normally translated within the neuronal components of excitatory synapses (e.g., Syt1, Syt11, Hspa8, Gria2; Hafner et al., 2019) could be detected and modulated by THC in microglial cells. Together with downregulation of these neuronal genes, THC upregulated microglia-intrinsic genes, such as Latmp5 and C1qa, whose products actively participate in phagocytic removal of inactive excitatory synapses (Mallya et al., 2019; Cong et al., 2020; Pöpplau et al., 2024). These same pruning and synaptic genes that were enriched among synaptic terms were also associated with terms that evoke phagocytic removal of neuronal components (neuronal death) and tagging of neuronal cargo for microglial phagocytosis (binding of phosphatidylserine; Scott-Hewitt et al., 2020; Zhuang et al., 2023), further supporting a scenario where THC would influence the pruning of maturing glutamatergic synapses. This interpretation is also consistent with the reduction in spine density we observed in L2/3 pyramidal neurons of the prelimbic cortex of THC-treated mice and with findings from a previous study implicating glutamatergic synapses in the cortical effects of THC in adolescent female rats (Rubino et al., 2015).
THC-sensitive genes in cortical astrocytes were also associated with the glutamatergic synapse, particularly genes coding for calcium sensing proteins. These are of interest since Ca2+ is the intracellular signal that allows astrocytes to communicate with other cells in the synapse (Eraso-Pichot et al., 2023). The observed downregulation of genes that code for astrocytic synaptotagmin 1 (Syt1) and calmodulin-1 (Calm1) could lead to changes in glutamate release (Zhang et al., 2004; Eraso-Pichot et al., 2023) and uptake (Ashpole et al., 2013), disrupting astrocytic buffering of synaptic glutamate. The notion that THC may have disturbed synaptic buffering of glutamate is further supported by the observation that its administration induced the ectopic upregulation of the glutamine synthase gene (Glul) within pyramidal neurons. The results we obtained in the frontal cortex are also consistent with previous studies showing that activation of astrocytic CB1Rs by THC disrupted extracellular glutamate levels in hippocampal synapses (Han et al., 2012; R. Chen et al., 2013). Changes to perisynaptic hippocampal glutamate modified membrane expression of synaptic AMPA/NMDA receptors, disrupted LTP, induced LTD, and were accompanied by impairment in spatial memory (Han et al., 2012; R. Chen et al., 2013). The reduced spine density and dendritic complexity we observed on superficial pyramidal neurons of the prelimbic cortex together with the reduced expression of genes that code for neuronal AMPA subunits (Gria2) and for proteins that enable expression of AMPA/NMDA at the postsynaptic membrane (e.g., Ywhah, Ptpn5, Igf1; Dyer et al., 2016; Kamceva et al., 2016; Hsia et al., 2019) are similarly consistent with THC-driven weakening of glutamatergic transmission in adolescent frontal cortex.
WIN-sensitive genes in murine frontal cortex were also enriched with glutamatergic synapse terms. In particular, transcriptomic changes to genes expressed in L2/3 pyramidal neurons were consistent with reorganization of excitatory synapses, including genes that code for AMPAR subunits and for glutamate transporters. Moreover, the upregulation of Prkg and Grim5 genes, whose protein products are integral components of endocannabinoid signaling pathways, suggested a reduction in glutamatergic transmission in these cells. Indeed, signaling supported by Prkg and Grim5 can be engaged by WIN to weaken glutamatergic transmission and produce impairments in learning and memory (Izumi and Zorumski, 2012). Reduced efficiency of glutamatergic transmission is also consistent with the observed upregulation of Sorcs3, a gene whose product supports CA1-dependent LTD (Breiderhoff et al., 2013), and with downregulation of the gene that codes for protein 14-3-3 ε, which maintains functional NMDARs at the membrane (Chung et al., 2015).
THC and WIN additionally downregulated genes coding for respiratory subunits in multiple mitochondrial complexes, suggesting that both ligands inhibited biogenesis of the respiratory chain. Mitochondrial biogenesis is a conserved cellular process that is controlled by PGC-1α (L. Chen et al., 2022) and under the influence of cAMP/CREB (L. Chen et al., 2022). This pathway can be inhibited by Gαi/o proteins that are downstream effectors of CB1 receptors. Genes for CB1R and its canonical transducers were detected in L2/3 pyramidal neurons and microglia, along with genes coding for βarrs. We showed that THC and WIN displayed distinct efficacies at modulating G-protein and βarr responses and proposed that these differences together with cell-specific levels of βarr expression could account for ligand-specific regulation of respiratory genes in L2/3 pyramidal neurons and microglial cells. It is unlikely, however, that this mechanism may explain why other glial transcriptomes were sensitive to THC but not WIN, since other glial cells lacked the βarr enrichment present in microglia. Immunolabeling allowed to corroborate carryover of cell-specific respiratory responses to the protein level. While microglial expression of NDUFB8 subunits was only reduced by THC, its neuronal expression displayed a tendency to specific downregulation by WIN. Lack of significance in WIN responses can be explained by labeling limitations as NeuN cannot distinguish between interneurons and pyramidal cells (Miller et al., 2019). The functional consequences of respiratory gene downregulation require functional characterization. However, broad inspection did not reveal gross structural changes to immunolabeled neurons or microglia. This is in contrast with a recent study where longer THC treatment induced microglial apoptosis (Hasegawa et al., 2023).
There are several limitations to our translational study. First, we performed analysis in males only, as cannabis-related group differences in cortical thickness were not observed in female participants of our original study (French et al., 2015). Moreover, since experimental models do not fully reflect the human experience, the magnitude and timeframe of cannabinoid exposure in preclinical models are not necessarily equivalent to cannabinoid consumption in the human cohorts. Therefore, the findings of the study should be translated with caution.
Nonetheless, it is conceivable that exposure to THC present in cannabis products or to synthetic cannabinoids before the age of 16 would modify molecular processes within astrocytes, microglia, and pyramidal neurons to alter the functional and structural maturation of glutamatergic synapses in human cortex. We speculate that these synaptic effects modify spine density and dendritic arborization of pyramidal neurons (i.e., CA1 type, whose cell-specific genes are involved in the regulation of dendritic growth; Shin et al., 2018), thus leading to differences in cortical thickness between adolescent males that experimented or not with cannabinoids. Our findings in the mouse model support strongly this possibility. Indeed, transcriptomic changes induced by both ligands converged to indicate impaired functional plasticity of the glutamatergic synapse and, most importantly, THC induced robust differences (vs controls) in the number of spines and dendritic complexity of cortical glutamatergic neurons.
Footnotes
This study was supported by grants from the Canadian Institutes of Health Research (NET203021, MOP86678, PJT153010, CA7170130, MOP324876) and the Fonds d'Accélération des Collaborations en Santé (FACS-008). We thank the members of the Saguenay Youth Study cohort and their families for their participation in the project.
↵*X.N., D.N.R., and I.C. are the joint first authors.
↵**G.P. and T.P. are the joint senior authors.
The authors declare no competing financial interests.
- Correspondence should be addressed to Graciela Piñeyro at graciela.pineyro.filpo{at}umontreal.ca or Tomáš Paus at tomas.paus{at}umontreal.ca.
