Single-cell Transcriptome Mapping Identifies Common and Cell-type Specific Genes Affected by Acute Delta9-tetrahydrocannabinol in Humans

Delta-9-tetrahydrocannabinol (THC) is known to modulate immune response in peripheral blood cells. The mechanisms of THC’s effects on gene expression in human immune cells remains poorly understood. Combining a within-subject design with single cell transcriptome mapping, we report that THC acutely alters gene expression in 15,973 blood cells. We identified 294 transcriptome-wide significant genes among eight cell types including 69 common genes and 225 cell-type-specific genes affected by THC administration, including those genes involving in immune response, cytokine production, cell proliferation and apoptosis. We revealed distinct transcriptomic sub-clusters affected by THC in major immune cell types where THC perturbed cell-type-specific intracellular gene expression correlations. Gene set enrichment analysis further supports the findings of THC’s common and cell-type-specific effects on immune response and cell toxicity. This comprehensive single-cell transcriptomic profiling provides important insights into THC’s acute effects on immune function that may have important medical implications.

With the increasing rates of cannabis use for recreational and medical purposes, it is important to address a gap in our understanding of its impact on immune and inflammatory functions 1,2 . Delta-9-tetrahydrocannabinol (THC), the principal psychoactive constituent of cannabis, powerfully modulates immune function in peripheral blood cells 3 , in part, through activating cannabinoid receptor 2 (CBR2) 1,[3][4][5][6] . In vitro studies of cannabis exposure, which contains over 450 compounds, show that it modulates immune function 7-10 , changes cytokine production 8,11,12 , inhibits cell proliferation 2 , and induces apoptosis 13,14 . However, little is known about the mechanisms of in vivo THC exposure on the transcriptomes of distinct types of peripheral blood mononuclear cells (PBMCs) in humans.
Single cell RNA-seq (scRNA-seq) offers an unprecedented resolution to detect drug effects on cell-specific gene expression 15,16 and enables the evaluation of molecular aspects of immune cell heterogeneity 17 . Few studies have applied scRNA-seq to detect differentially expressed genes (DEGs) induced by drug exposure, and none have evaluated the effects of THC in humans. This limitation is due mostly to high inter-individual transcriptomic variability and types of cells that confound the assessment of the impact of environmental factors. Most recently, a scRNA-seq study identified a large number of common and cell-type-specific DEGs for Alzheimer disease, suggesting the improvement of analytical methods to overcome the challenge of high transcriptomic variability 18 . Here, we report the first scRNA-seq study using within-subject combined with linear mixed model (LMM) analysis to detect genes affected by intravenous (IV) THC at single cell resolution. It is conceivable that other routes of Delta-9-tetrahydrocannabinol (THC) alters gene expression shared in multiple cell types of PBMCs. We next applied LMM to detect individual genes affected by THC infusion, with participant included as a random variable to limit confounding effects by individual's genomic background. We identified 294 transcriptome-wide significant genes in eight cell types affected by THC infusion (false discovery rate, FDR < 0.05) (Fig. 2a). DCs and FCGR3A monocytes were excluded from further analyses because no gene reached transcriptome-wide significance in DCs and both cell types had low cell numbers. Among the 294 DEGs, 69 were observed in at least two cell types while 225 were significant in only one cell type ( Fig. 2b; Tables S4-S11). Overall, THC infusion resulted in more upregulation than downregulation in gene expression.
We sought to identify THC-regulated genes common across the six common cell types. We found 28 DEGs in at least three cell types (Fig. 2c). The majority (89%) of the DEGs showed consistent directions of gene expression changes by THC in different cell types; only three genes displayed opposing direction in CD14+ monocytes compared to other cell types (TMSB4X, JUNB, TXNIP). A group of THC-regulated gene functions are in the domains of immune response and inflammatory process. For example, expression of S100A9 and S100A8, which play prominent roles in the regulation of proinflammatory processes and immune response [31][32][33] , decreased in response to THC infusion in five cell types. A major HIV-1 suppressive gene, CCL4, showed increased expression after THC infusion. THC decreased expression of GNLY that is involved in activating antigen presentation 34 . The altered expression of genes involved in humoral immunity were also observed (IGLC2, IGKC). Changes in expression of these genes suggest that THC activated the adaptive immune system shortly after administration, consistent with findings showing an immunomodulatory effect by THC and more complicated than solely immunosuppression 3,8,[35][36][37][38] .
Among the 28 DEGs shared by at least three cell types, a subset of DEGs supports previous findings that acute THC exposure inhibits cell proliferation and induces apoptosis. Genes responsible for cell death were upregulated (BTG1 39 , DDIT4 40 , GZMB) and genes involving in cell growth and differentiation were down-regulated (TMSB10 41 , RPS21 42 , RPL41 43 by THC exposure. The alteration of these genes in distinct cell types may indicate potentially deleterious effects of THC on cell differentiation and survival. We carried out RT quantitative PCR (qPCR) in attempt to validate gene expression differences of common DEGs in at least four cell types using the PBMCs from the same blood samples. Using bulk RNA extracted from the PMBCs, we found that the fold changes of 13 DEGs were highly correlated between bulk PBMCs and CD8+ T-cells (r = 0.76, p = 0.002), and marginally correlated with CD4+ T-cells (r = 0.48, p = 0.09). Figure S3 is a heatmap showing the fold change of gene expression induced by THC by scRNA-seq in six cell types and the fold change of gene expression in bulk RNA by qPCR. This result suggests that THC altering gene expression in 13 common genes is consistent between PBMCs and subtypes of cells. To compare fold changes of DEGs between CD8+ T-cells and PBMCs, we found that expression changes by THC in 9 out of 13 DEGs were concordant. Differential expression of pre-and post-THC in two genes, TMSB4X and MT-ND5, were in opposite direction between bulk RNA in PBMCs and scRNA-seq in CD8+ T-cells. Two genes, ID2 and ZFP36L2, were differentially expressed between pre-THC and post-THC by scRNA-seq method, but not by qPCR. Therefore, this experiment validates 70% of DEGs regulated by acute THC administration in PBMCs.
THC perturbs gene expression and KEGG-based gene network in unique cell type and sub-cell clusters. Given that the majority of THC-regulated DEG are unique to each cell type, we then focused on DEGs and gene-gene networks in each cell type. DEGs were determined both in individual clusters in each cell type and among all clusters in each cell type. Here, we conducted a gene network analysis by leveraging gene-gene relationships cataloged in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database 44 . The rationale is to identify gene networks within the biologically known genes or pathways. Thus, the KEGG-based gene network is defined as the inferred gene-gene relationships by KEGG in each cell type. We constructed gene networks independently in pre-and post-THC samples in each cell type (FDR < 0.05). Hub genes in each cell-type based network were defined as a node (gene) with at least four edges (gene links) in at least one condition (pre-, post-THC, or both). We then performed Gene Ontology (GO) term enrichment analysis for each cell-type-specific KEGG-based gene network (FDR < 0.05).
In CD4+ T-cells, significant DEGs are involved in cytotoxic T-cell activation (IL7R 45 ), histone modification (H3F3B 46 ), and transcriptional regulation (MYC 47 (Fig. 3b), showing that acute THC exposure perturbed gene-gene relationships in CD4+ T-cells and appears to increase the gene-gene statistical correlation involving cell-mediated immunity.
In CD8+ T-cells, we identified 18 unique DEGs that are involved in immune response and inflammation in response to THC (e.g. IL32 48 , SOCS1 49 , and IRF1 50 ). Of note, THC infusion resulted in the differential expression of CXCR4, TSC22D3, DDIT4, BTG1, JUN in cluster 18 of CD8+ T-cells (Fig. 3c), suggesting that cells in cluster 18 may function differently in response to THC as compared to the other CD8+ T-cell sub-group. The KEGG-based gene network included 35 nodes that enriched on 12 GO terms (e.g., immune system process), similar to the CD4+ T-cell network. One hub gene, HLA-A (strongly correlated with HLA-B and HLA-C in pre-and post-THC samples), was observed in the network ( Fig. 3d; Table S14).
In B cells, genes unique to B cells were observed that are involved in B cell maturation (VPREB3 51 ), MHC function (HLA-DQA1, HLA-DQA2), calcium signaling (CALM2 52 ), Toll-like receptors (CD180 53 ), and response to environmental stress via activating MAP kinase MAPK1/ERK2 (DUSP1 54 ). Consistent with in vitro findings, we found that acute THC exposure reduced expression of BCL2 in B cells 45 . The majority of DEG originated from cluster 14 ( Fig. 3e; Table S15), with only 8 DEGs from cluster 12, which were all upregulated by THC. In the KEGG-based gene network, we found 34 nodes with 28 edges (Fig. 3f) enriched on GO pathways including MHC protein complex, immune response, peptide antigen binding. Four hub genes were identified that differed in response to THC (HLA-A and HLA-BQA1) HLA-DPA1 [became more strongly correlated with HLA-DQB1 post-THC], ANX1 [connected to CXCL1, CXCR4, XOL2 in pre-but not post-THC samples]). Notably, two B cell marker genes (BLNK, CD79B) were correlated in post-but not pre-THC samples.
In NK cells, which is defined by one cluster, many DEGs (e.g. DDIT4, CCL4, BTG1 ID2) are involved in immune response and cell proliferation. One DEG unique to NK cells, CD53, was downregulated by THC. Genes in the NK cell network (17 nodes and 12 edges; Fig. 3g) were enriched on 17 pathways including chemokine-mediated signaling pathway, inflammatory response, and chemokine receptor activity. Notably, KLRC1 regulates specific humoral and cell-mediated immunity and is implicated in the recognition of the MHC class I HLA molecules in NK cells and was correlated with HLA-C, which increased in post-THC samples. No hub genes were observed.
These observations from individual gene and KEGG-based gene networks in different cell types suggest a diverse functional response to acute THC exposure in heterogenous immune cells.
Significant pathways for all KEGG-based gene networks in each major cell type are presented in Tables S17-S22.
Common biological pathways including immune response and cell survival/apoptosis are associated with THC administration. We subsequently performed a cell-type based gene set enrichment analysis (GSEA). We found 39 significant KEGG pathways in at least two cell types ( Fig. 4; Tables S24-S31); significant pathways were involved in the domains of immune response, inflammation, and cell survival and GLM tested the association of cell type and marker gene expression in each cell. Significance threshold is set at p < 0.01. Each individual cell is assigned to a cell type based on the predominant proportion of marker genes in a given cell type. Small cell clusters are merged into the closest cell type in the t-SNE plot. The cell mapping deconvolutes cells to eight PBMC cell types: CD4+ T-cells, ILR7+/CD4+ T-cells, CD8+ T-cells, B cells, natural killer cells, CD14+ monocytes, FCGR3A monocytes, and DCs. (f) Percentage of cell numbers from each sample in a given cell type: pre-THC-S1, post-THC-S1, pre-THC-S2, post-THC-S2.  www.nature.com/scientificreports www.nature.com/scientificreports/ apoptosis. Several pathways associated with autoimmune disease (e.g., rheumatoid arthritis) were significant in multiple cell types. The ribosomal pathway, which plays a role cell growth and cellular response to stress, was the most significantly enriched pathway across all five cell types. These results further support the effects of THC on functional domains in immune regulation that have also been associated with immunological disease, although the causality of these relationships is unknown.
THC changes CNR2 co-expression network in B cells. Finally, we were interested in understanding how cannabinoid receptor genes co-expressed with other genes in each cell type following THC administration. As expected, CNR2, encoding for CBR2, was highly expressed in B cells, followed by NK cells, then CD8+ T-cells, and lowest in CD4+ T-cells (Fig. 5a) 55 . Little CNR1 expression was detected in any of the cell types. GPR55 (cannabinoid receptor 3 56 ) showed the highest expression in CD4+ T-cells. We observed gene co-expression between CNR2 and 84 and 74 genes in pre-and post-THC samples, respectively (Tables S31-S32). The co-expressed genes in post-THC B cells were enriched on functional domains of immune processes, biological regulation, cell proliferation, signaling, and response to stimulus (Fig. 5b). Interestingly, little overlapped genes in two CNR2 co-expressed network, suggesting that THC significantly reorganizes CNR2 connected gene-gene relationships in B cells (test of overlapping genes between two group p = 0.302).

Discussion
Our results from a well-controlled human study demonstrate transcriptomic regulation in distinct immune cells following administration of a single dose of THC. In our model, THC increased subjective effects ("high", "stoned") and heart rate consistent with the known effects of cannabis. Thus, THC altered gene expression at a dose that produces effects consistent with the known effects of cannabis. Our study design and analytical approach enabled the identification of common and cell-type-specific DEGs regulated by THC. The cell-type-specific gene and co-expression patterns revealed by scRNA-seq showed little overlap across cell types, which would have been obscured by bulk RNA-seq analysis. Cell type-based DEGs, co-expression networks, and GSEA revealed important THC effects on immune function, cytokine production, signal transduction, and cell apoptosis and survival.
We used within-subject design and mixed regression model to minimize confounding factors on THC's effects and to control high individual transcriptome variability. The approach also reduces potential batch effects across different individuals and experimental batches. Subjective, behavioral and cognitive tests confirm that the dose of THC was relevant to but limited the confounding effects of other compounds in cannabis. These strengths enabled us to reveal differentially expressed genes from highly variable transcriptomes. In addition, we applied GLM method for cell-type identification combining with knowledge-based cell markers. The approach identified biologically meaningful cell-types sub-clusters in major cell types in PBMCs. The sub-clusters present specific DEG patterns altered by THC and different co-expression networks regulated by THC, suggesting that the sub-clusters may be functionally distinct cell types. The biological verification of sub-clusters is warranted in a future study.
Although the majority DEGs were cell-type-specific, we observed top significant differentially expressed genes shared across multiple cell types. Importantly, we validated the top 18 common genes affected by THC in at least four cell types in bulk RNA from PBMCs by RT q-PCR in the same samples. The functions of common DEGs were relevant to the biological pathways shared in multiple cell types as well. These common DEGs present almost identical direction by THC regulation, suggesting that THC's immune response is involved in important genes to maintain immune response across different cell types. We identified more up-regulated genes than down-regulated genes by THC, which is expected after acute THC administration. Importantly, the large number of DEGs shows cell-type-specific fashion and a single low dose of THC perturbates cell-type-specific co-expression network. The complexity of THC's cell-type-specific effect raises a challenge of therapeutic treatment using cannabis for medical diseases.
Limitations. The results of this exploratory, pilot study should be considered in the context of its limitations and strengths. First, the findings need to be validated in a larger sample. Second, while the IV THC paradigm offers excellent experimental control, it has limited ecologically validity. Therefore, the effects of acute administration of smoked or vaporized cannabis should be studied. Third, while the effects of a single dose of THC were studied, given the known biphasic effects of cannabinoids, the dose-related effects of cannabinoids on gene expression need to be investigated. Furthermore, while we demonstrated the changes in gene expression at one time point after the administration of THC, the full temporal profile of changes in gene expression need to be studied. Fourth, while the effects of a single dose of THC were demonstrated, the effects of chronic exposure to medicinal and recreational cannabinoid need to be studied. Fifth, while the effects of THC in healthy subjects was demonstrated, how cannabinoids people with altered immune function (e.g., HIV) or chronic inflammatory conditions (e.g., rheumatoid arthritis) needs to be studied. Sixth, the next logical step is to investigate the precise mechanisms underlying THC altered gene expression observed in the exploratory study. For example, CHIP-qPCR analysis to detect histone markers may further enhance our understanding of the mechanism of the identified immune/inflammation-related genes function 57 . Finally, while the effects of THC on transcriptomic profiling and differential gene expression has been demonstrated, the effects of the other constituents of cannabis needs to be studied. monocytes. (i) KEGG-based gene network in CD14+ monocytes. Genes with nominal p < 0.001 in CD4+ T-cells were selected to construct a KEGG-based network. We observed 18 nodes and 13 edges. One hub gene, HLA-DQB1, was identified in both pre-THC and post-THC samples. (2020) 10:3450 | https://doi.org/10.1038/s41598-020-59827-1 www.nature.com/scientificreports www.nature.com/scientificreports/ Nevertheless, the study provides a foundation for future studies of cell-type-specific immunologic effects of cannabis and cannabinoid constituents. These findings highlight the complexity of cannabinoid effects on immune function and nuance of immune cell types that may have medical relevance.

Methods
The study was conducted at the Neurobiological Studies Unit (VA Connecticut Healthcare System, West Haven, CT) with the approval of the Institutional Review Boards at VA Connecticut and Yale University, and the Protocol Review Committee of the Department of Psychiatry, Yale University. The study was amended to include prospective measures addressing safety. All methods were performed according to the relevant guidances and regulations by Yale University and VA Connecticut Healthcare System. Participants. Two healthy participants were recruited from the community by advertisements. Both participants were male, 21-year old, and of European American descent. Subjects were informed about the potential for psychosis, anxiety, and panic. After obtaining informed consent, subjects underwent a structured psychiatric interview for DSM-IIIR 58 and were carefully screened for any DSM-IV Axis I or Axis II lifetime psychiatric or substance abuse disorder (excluding nicotine) and family history of major Axis I disorder. The history provided by subjects was confirmed by a telephone interview conducted with an individual (spouse or family member) identified by the subject prior to screening. In order to avoid exposing cannabis-naïve individuals to a potentially addictive substance, only subjects who had been exposed to cannabis but did not meet lifetime criteria for a cannabis use disorder were included. Past month cannabis use was quantified using a time-line-follow-back approach. Finally, subjects underwent a general physical and neurologic examination, EKG, and laboratory tests (serum electrolytes, liver function tests, complete blood count with differential and urine toxicology). Subjects were instructed to refrain from caffeinated beverages, alcohol, and illicit drugs from 2 weeks prior to testing until study completion. Urine toxicology was conducted on the morning of each test day to rule out recent illicit drug use.
Procedure. Subjects received 0.03 mg/kg of THC, the principal active ingredient of cannabis 23,24 . This dose equivalent to 2.1 mg in a 70 kg individual has been shown in previous studies to produce effects consistent with the known effects of cannabis, in a safe manner [59][60][61] . THC was administered on its own, without the >450 other chemical constituents of cannabis because THC is the principal active constituent of cannabis and the other chemical constituents could render the results challenging to interpret. The intravenous route of administration was chosen to reduce inter and intraindividual variability in plasma THC levels with the inhaled route 62 . Timeline of behavioral assessment and blood draw is presented in Fig. S5. Subjects were attended to by a research psychiatrist, a research nurse, and a research coordinator. Clear 'stopping rules' were determined a priori and rescue medication (lorazepam) was available if necessary. Medical condition and psychiatric status of participants were monitored closely during and after THC challenging. Subjective and clinical ratings were repeatedly assessed. www.nature.com/scientificreports www.nature.com/scientificreports/ www.nature.com/scientificreports www.nature.com/scientificreports/ Medical and behavioral assessment. Vital signs, Cannabis intoxication, Psychotic symptoms, Perceptual alteration, and Cognitive test battery were measured prior, during, and post THC infusion as illustrated in Fig. S5.
Physiological and psychoactive assessments of two participants are presented in Table S33.
Single cell RNA sequencing in 10X genomics platform. PBMCs from pre-THC (N = 2) and post-THC samples (N = 2) from fresh whole blood were isolated at the same time using a standard protocol. The cells were washed twice with phosphate-buffered saline containing 0.04% bovine serum albumin, and the final cell concentration was adjusted to 1000 cells/mL for library preparation. 5000 cells/sample were prepared for single cell capture.
Sequencing data processing and quality control. We used the 10X Genomics Chromium Single Cell 3′ v2.0 platform to prepare individually barcoded single cell RNA-Seq libraries following the manufacturer's protocol. Data normalization. Only genes with at least one UMI count detected in at least one cell were retained for analysis Single cells were excluded when >10% of reads mapped to mitochondrial RNA to ensure that all of the single cells originated from nucleated cells. The cells with fewer than 370 expressed genes or possible doublet cells (>100,000 reads) were also discarded. Applying these three criteria resulted in retention of 21,430 genes and 15,973 single cells for downstream analysis. Data normalization was performed using a standard protocol in Seurat 28 . An exploratory analysis showed cells clustered by each participant. Batch effects were then removed using an empirical Bayesian framework (ComBat function in R package sva) with the individual labeled as the random variable. The normalized data was dimensionally reduced in two dimensions using t-distributed stochastic neighbor embedding (t-SNE) with a perplexity parameter of 20 and 3000 iterations after an initial principle component analysis. The distance matrix of the single cells was computed with the output of t-SNE and converted into a graph. Cluster_louvain function was used by default options (https://igraph.org/r/doc/cluster_louvain.html). Then, the graph was clustered with the cluster_ louvain function in R package igraph, which implements the multi-level optimization of modularity algorithm for finding community structure.
Cell type identification. A gene set of the cell type markers in PBMC was manually curated from the literature. A binary cell type matrix with cell types as columns and genes as rows was generated with value 1 representing a marker in a cell type and value 0 denoting non-markers for a cell type. The generalized linear model (GLM) was constructed by deciding on one vector, representing one cell type as a response, of the binary matrix and corresponding gene expressions as explanatory variables in one cell. In details, = × Y y ( ) gc G C is the binary cell type matrix of G marker genes of C cell types and = × X x ( ) gs G S is the gene expression matrix of G marker genes of S single cells. For the g th gene, y 0 gc = or 1 denotes the cell type-specific gene indicator of c th cell type and x gs denotes the gene expression of s th single cell, where = … g G . One statistical value (t value) and p value is obtained for each cell and cell type; t cs and p cs denote the test statistic value and P-value, respectively for the c th cell type of s th single cell.
We used the cutoff, p < p 0 and t value >0 to assign one cell to one cell type, where p 0 is a predefined cut-off value (p 0 = 0.01 in this study). Then, the proportion of each assigned cell type was calculated for every cluster. The dominant proportions are used to assign one cluster to one cell type. where P c k denotes the proportion of the c th cell type in the k th cluster. S k is the set of cell indices for clustering k, with N S k k = being the total number of cells in the k th cluster. Finally, the cell type for each cluster was confirmed manually by cell type marker gene expressions mapping on the 2D t-SNE plot. The cell clusters mapping to the same cell type were merged for downstream analysis.
Differential gene expression. We applied a linear mixed regression model to identify genes associated with THC infusion. To control transcriptomic variation between two participants, we used participant as a random effect. Transcriptome-wide significance was set at false discovery rate (FDR) < 0.05.
KEGG-based gene network analysis. Here, gene network was defined as significant gene-gene correlation in each given cell type leveraging the relationship among genes in KEGG database. Using this supervised gene network analysis approach, we tested whether gene links in KEGG were significant in each cell type using linear regression model. Significant KEGG-based network was set at FDR < 0.05. The analysis was performed separately in pre-THC and post-THC samples to test if THC alters gene-gene correlation.
Gene Ontology (GO) enrichment analysis. Genes in each co-expression network in major cell types were tested the enrichment on GO terms using The DAVID Gene Functional Classification Tool 63 . Significant level was set at FDR < 0.05.
Cell type-based gene set enrichment analysis. We separately analyzed gene set enrichment in each cell type using KEGG database. Differentially expressed genes in a given cell type with nominal p < 0.001 were selected. Significant pathway was set at FDR < 0.05.
Co-expression of cannabinoid receptor genes with other genes Four cannabinoid receptor genes, CNR1, CNR2, CPR18, and GPR55 were tested differential expression across seven major cell types using linear mixed regression model. Correlation of CNR2 with other genes in B cells was performed in pre-and post-THC samples independently. Significant correlation was set at FDR < 0.05. We tested for statistical significance of the overlap between the two groups (pre-and post-THC) of genes (http://nemates.org/MA/progs/overlap_stats.cgi).

Data availability
Single cell transcriptome data have been deposited in Gene Expression Omnibus and are available under project number GSE130228.