Single-nucleus transcriptomics reveals time-dependent and cell-type-specific effects of psilocybin on gene expression

Single-nucleus transcriptomics reveals time-dependent and cell-type-specific effects of psilocybin on gene expression | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Single-nucleus transcriptomics reveals time-dependent and cell-type-specific effects of psilocybin on gene expression View ORCID Profile Clara Liao , Ethan O’Farrell , View ORCID Profile Yaman Qalieh , View ORCID Profile Neil K. Savalia , View ORCID Profile Matthew J. Girgenti , Kenneth Y. Kwan , View ORCID Profile Alex C. Kwan doi: https://doi.org/10.1101/2025.01.04.631335 Clara Liao 1 Meinig School of Biomedical Engineering, Cornell University , Ithaca, NY 14853, USA 2 Interdepartmental Neuroscience Program, Yale University School of Medicine , New Haven, CT 06511, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Clara Liao Ethan O’Farrell 3 Michigan Neuroscience Institute, University of Michigan , Ann Arbor, MI 48109, USA 4 Department of Human Genetics, University of Michigan , Ann Arbor, MI 48109, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Yaman Qalieh 3 Michigan Neuroscience Institute, University of Michigan , Ann Arbor, MI 48109, USA 4 Department of Human Genetics, University of Michigan , Ann Arbor, MI 48109, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Yaman Qalieh Neil K. Savalia 1 Meinig School of Biomedical Engineering, Cornell University , Ithaca, NY 14853, USA 2 Interdepartmental Neuroscience Program, Yale University School of Medicine , New Haven, CT 06511, USA 5 Medical Scientist Training Program, Yale University School of Medicine , New Haven, Connecticut, 06511, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Neil K. Savalia Matthew J. Girgenti 6 Department of Psychiatry, Yale School of Medicine , New Haven, CT 06511, USA 7 Clinical Neuroscience Division, National Center for Posttraumatic Stress Disorder, Veterans Affairs Connecticut Healthcare System , West Haven, CT 06516, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Matthew J. Girgenti Kenneth Y. Kwan 3 Michigan Neuroscience Institute, University of Michigan , Ann Arbor, MI 48109, USA 4 Department of Human Genetics, University of Michigan , Ann Arbor, MI 48109, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: kykwan{at}umich.edu alex.kwan{at}cornell.edu Alex C. Kwan 1 Meinig School of Biomedical Engineering, Cornell University , Ithaca, NY 14853, USA 8 Department of Psychiatry , Weill Cornell Medicine, New York, NY 10065, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Alex C. Kwan For correspondence: kykwan{at}umich.edu alex.kwan{at}cornell.edu Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF ABSTRACT There is growing interest to investigate classic psychedelics as potential therapeutics for mental illnesses. Previous studies have demonstrated that one dose of psilocybin leads to persisting neural and behavioral changes. The durability of psilocybin’s effects suggests that there are likely alterations of gene expression at the transcriptional level. In this study, we performed single-nucleus RNA sequencing of the dorsal medial frontal cortex of male and female mice. Samples were collected at 1, 2, 4, 24, or 72 hours after psilocybin or ketamine administration and from control animals. At baseline, major excitatory and GABAergic cell types selectively express particular serotonin receptor transcripts. The psilocybin-evoked differentially expressed genes in excitatory neurons were involved in synaptic plasticity, which were distinct from genes enriched in GABAergic neurons that contribute to mitochondrial function and cellular metabolism. The effect of psilocybin on gene expression was time-dependent, including an early phase at 1-2 hours followed by a late phase at 72 hours of transcriptional response after administration. Ketamine administration produced transcriptional changes that show a high degree of correlation to those induced by psilocybin. Collectively, the results reveal that psilocybin produces time-dependent and cell-type specific changes in gene expression in the medial frontal cortex, which may underpin the drug’s long-term effects on neural circuits and behavior. INTRODUCTION Psilocybin is a serotonergic psychedelic that has shown promise as a potential therapeutic for mental health conditions. Data from several clinical trials indicate that psilocybin with psychological support can relieve symptoms in patients with major depressive disorder or treatment-resistant depression 1 , 2 , 3 , 4 , 5 . Although the neurobiological basis for the therapeutic effects is not understood, there is growing evidence that the antidepressant effect may relate to the drug-evoked structural neural plasticity in the prefrontal cortex 6 . The plasticity-promoting effect of psilocybin and its active metabolite psilocin was demonstrated with cultured cortical neurons, where the application of psilocin induced dendritic outgrowth and spinogenesis 7 . Moreover, in live mice, psilocybin administration elevates the density of dendritic spines in the mouse medial frontal cortex, an effect that persists for over a month 8 , 9 . Therefore, studies in humans and model systems suggest psilocybin can induce enduring effects on neuronal connectivity and behavior. Long-lasting modifications to synaptic connections and behavior require changes in gene expression 10 , 11 . Indeed, previous studies of mRNA transcript and protein levels identified differential expression of genes for synaptic function, glutamate receptor signaling, and immune response in the neocortex after a single dose of a classic psychedelic such as 2,5-dimethoxy-4-iodoamphetamine (DOI) 12 , 13 , 14 , 15 , lysergic acid diethylamide (LSD) 12 , 16 , 5-methoxy- N , N -dimethyltryptamine (5-MeO-DMT) 17 , or psilocybin 15 , 18 , 19 . Although these results demonstrated transcriptional alterations after drug exposure, the early studies relied on methods such as quantitative PCR and microarrays, which were limited by their focus on preselected candidate genes. A genome-wide transcriptomic investigation can yield novel insights, as exemplified by recent studies of transcriptional responses to drugs such as cocaine 20 , 21 , 22 , heroin 23 , fentanyl 24 , DOI 14 , MDMA 25 , ketamine and imipramine 26 , as well as in mouse models of chronic stress 27 , 28 and early life stress 29 , 30 . In addition to the goal of sequencing the transcriptome, understanding the temporal dynamics of the transcriptional response is also important. In humans, the oral administration of psilocybin induces subjective effects on perception and emotion within 30 minutes 31 , 32 . Peak subjective effects for some individuals arrive as early as 60 minutes, coinciding with peak plasma psilocin levels 33 and pronounced changes in brain functional connectivity 34 . The subjective effects subside in 3 – 6 hours 32 , because of rapid clearance with an elimination half-life of 3 hours 35 .Long after the subjective effects wear off, therapeutic effects for major depression can be detected with an onset at 2 – 8 days 1 , 4 , and persist for at least 6 weeks after the initial administration 4 . Although reverse translation of this time course to rodents is inexact, psilocybin produces likewise acute effects in mice within a couple hours, including head-twitch responses 8 , 36 , ultrasonic vocalizations 37 , and cognitive flexibility 38 . These are followed in subsequent days by changes in behavioral assays such as learned helplessness 8 , 9 , forced swim 39 , tail suspension 9 , sucrose preference 40 , 41 , fear extinction 42 , and social reward learning 25 . Therefore, efforts to understand psilocybin-evoked changes in gene expression should consider time as a key parameter. With these goals in mind, we designed a study to perform single-nucleus RNA sequencing (snRNA-seq) to determine the transcriptomic changes in the dorsal medial frontal cortex of male and female mice after psilocybin administration. To capture both acute and long-term changes, we collected tissues at 1, 2, 4, 24, and 72 hours after drug exposure, as well as from control animals. We delineated differentially expressed genes (DEGs) due to psilocybin administration in excitatory and GABAergic cell types. For excitatory neurons, our results show that psilocybin leads to an overrepresentation of genes linked to synaptic plasticity. By contrast, GABAergic interneurons showed an enrichment in biological processes associated with mitochondrial function and cellular metabolism. Temporally, we identified two waves of elevated transcriptional responses: an initial phase at 1-2 hours and a late phase at 72 hours after psilocybin administration. Collectively, the results reveal cell-type-specific and time-dependent changes to the transcriptional landscape in the mouse medial frontal cortex after a single dose of psilocybin. RESULTS Single-nucleus RNA sequencing of the mouse medial frontal cortex at various timepoints after a single dose of psilocybin In this study, we focus on the mouse dorsal medial frontal cortex, which encompasses the anterior cingulate cortex (ACAd) and the premotor cortex (medial MOs). This region is important for the drug action of psilocybin and ketamine, showing robust c-Fos response 43 and structural neural plasticity 8 , 9 following psilocybin administration. We administered C57/BL6J mice with psilocybin (1 mg/kg, i.p.), and then collected tissues corresponding to the dorsal medial frontal cortex of both hemispheres via microdissection at 1, 2, 4, 24, or 72 hours after drug administration (n = 20 animals; 2 male and 2 female for each time point; Figure 1a ). Additionally, we collected control tissues via the same method from animals that received no injection (n = 9 animals; 3 male and 6 female). The tissues were processed through nuclear dissociation, nuclei sorting, and barcoding for snRNA-seq via the 10x Genomics Chromium Single Cell 3’ v3 platform. We started with 145,593 cells that were sequenced at a mean depth of 9023±366 reads per cell for the psilocybin group, and 67,600 cells that were sequenced at a mean depth of 8423±302 reads per cell for the controls. After removing cells that did not meet quality control criteria (see Methods ), our data set consisted of 122,685 cells for the psilocybin group (range = 18,944 – 31,443 single nuclei per time point; n = 17 animals) and 56,693 cells for the control group (n = 8 animals). Download figure Open in new tab Fig. 1: Single-nucleus RNA sequencing of mouse medial frontal cortex after psilocybin administration. (a) Experimental workflow. (b) UMAP representation of the snRNA-seq data set including 122,685 nuclei from psilocybin-treated mice (range = 18,944 – 31,443 nuclei per time point; n = 17 animals) and 56,693 nuclei from the control mice (n = 8 animals). Color denotes the cell-type cluster. (c) Heatmap for the expression of cell-type-specific genes for the cell-type clusters. Dendrogram was generated by performing hierarchical clustering using expression profiles across all conditions. CPM, count per million. L6-CT, layer 6 corticothalamic. L5-PT, layer 5 pyramidal tract. L4/5-IT, layer 4/5 intratelencephalic. L4-IT, layer 4 intratelencephalic. L2/3-IT, layer 2/3 intratelencephalic. L5/6-IT, layer 5/6 intratelencephalic. CIN-VIP, cortical vasointestinal peptide-expressing interneurons. CIN-PV, cortical parvalbumin-expressing interneurons. CIN-SST, cortical somatostatin-expressing interneurons. OL, oligodendrocytes. AC, astrocytes. OPC, oligodendrocyte progenitor cells. We segmented the cells into clusters by performing dimensionality reduction and then Leiden clustering. Among all clusters, we focused on 12 clusters identified based on expression profiles of established marker genes and annotated them with corresponding cell types ( Figure 1b, c ; see Methods ). The majority of cells were classified into one of these clusters, including 108,060 cells (88.13%) in the psilocybin group and 49,963 cells (88.08%) in the control group. The 12 clusters included 6 classes of excitatory neurons – layer 6 corticothalamic neurons (L6-CT; number of cells per time point or control condition: 3952±750, 13.2% of all cells that met quality control criteria), layer 5 pyramidal tract neurons (L5-PT; 1420±241, 4.8%), layer 4/5 intratelencephalic neurons (L4/5-IT; 3582±761 cells, 12.0%), layer 4 intratelencephalic neurons (L4-IT; 2469±450, 8.3%), layer 2/3 intratelencephalic neurons (L2/3-IT; 4402±856, 15.0%), and layer 5/6 intratelencephalic neurons (L5/6-IT; 2090±374, 7.0%), alongside 3 classes of GABAergic neurons – vasoactive intestinal peptide-expressing interneurons (CIN-VIP; 853±134, 2.9%), parvalbumin-expressing interneurons (CIN-PV; 1425±234, 4.8%), and somatostatin-expressing interneurons (CIN-SST; 1039±150, 3.5%), as well as 3 classes of non-neuronal cells – oligodendrocytes (OL; 2179±626, 7.3%), astrocytes (AC; 1680±323, 5.6%), and oligodendrocyte progenitor cells (OPC; 1156±220, 3.9%) ( Supplementary Table 1 ). The same cell-type clusters were detected in samples collected at different time points after psilocybin administration and in controls ( Figure S1 ). Cell-type specific expression of serotonin receptor transcripts in the mouse medial frontal cortex Psilocin binds not only to the 5-HT 2A receptor, but also to other serotonin receptor subtypes expressed in the neocortex, including the 5-HT 1A and 5-HT 2C receptors 44 , 45 . We visualized the presence of the Htr1a , Htr2a , and Htr2c transcripts which encode the 5-HT 1A , 5-HT 2A , and 5-HT 2C receptors respectively in the control samples in our snRNA-seq data set ( Figure 2a – c ). The distribution of the three transcripts clearly differed, prompting us to quantify how the transcript levels vary across the 12 cell-type clusters ( Figure 2d – f ). To further support this analysis and incorporate additional identified cell types, we also extracted publicly available single-cell sequencing data from the Allen Institute 46 , focusing on cells residing in the mouse frontal cortex, including the ACA, ALM, ORB, and PL-ILA regions ( Figure 2g ). There was a strong agreement between the two data sets. Many IT and PT excitatory neurons express a high level of Htr2a as well as moderate amounts of Htr1a and Htr2c . Interestingly, there is also an abundance of Htr2a transcripts in the deep-lying Layer 6b neurons, which are known for their excitatory effects on apical dendrites and higher-order thalamus 47 . Unlike the excitatory neurons, the expression pattern is more selective in the GABAergic cell subtypes. The three major interneuron subtypes – CIN-PV, CIN-SST, and CIN-VIP – predominantly express only one of the receptor transcripts: Htr2a , Htr1a , and Htr2c respectively. Lamp5 -expressing interneurons, which include the Ndnf -expressing neurogliaform cells 48 , express Htr2a , with relatively minimal expression of Htr1a and Htr2c . Finally, although serotonin receptors have been reported in microglia 49 , 50 and astrocytes 51 , the current data suggest that the Htr transcript levels are relatively low in the non-neuronal cell types. Altogether, these analyses reveal a highly selective pattern for how serotonin receptors are expressed in cortical cell types. Because psilocin is a non-selective serotonergic agonist that binds to a range of receptor subtypes, we anticipate that the drug can directly influence the transcriptional programs of most major excitatory and GABAergic cell types. Download figure Open in new tab Fig. 2: Cell-type specific expression of serotonin receptor transcripts in the mouse medial frontal cortex (a) UMAP representation of the snRNA-seq data set. Blue circle, cell in control mice that expresses Htr1a . Gray circle, cell in control mice that does not express Htr1a . (b, c) Similar to (a) for Htr2a and Htr2c respectively. (d) The expression level of Htr1a for the 12 cell-type clusters. (e, f) Similar to (d) for Htr2a and Htr2c respectively. (g) Single-cell transcript counts summed from the mouse ACA, ALM, ORB, and PL-ILA regions, extracted from the SMART-Seq data set from Allen Institute. Slc17a7 and Gad1 are markers of excitatory and GABAergic neurons, respectively. Pvalb , Cck , Vip , Sst , and Ndnf are markers of subtypes of GABAergic neurons. Psilocybin-induced transcriptional changes in frontal cortical excitatory neurons are involved in synaptic plasticity To examine how psilocybin alters the broad transcriptomic landscape, we implemented memento 52 , an end-to-end framework that uses a hierarchical model for statistical analyses of multidimensional sequencing data. This analysis showed that certain cell types had similar patterns of transcriptionally responsive genes evoked by psilocybin ( Figure S2 ). Informed by this, we grouped L4/5-IT, L4-IT, L2/3-IT, L5/6-IT cell types together as IT neurons, and CIN-VIP, CIN-PV, CIN-SST cell types together as cortical interneurons (CINs). The IT neurons and CINs were compared to L5-PT neurons (PT neurons). This aggregation allowed us to proceed next to identify DEGs with more cells in each group at each time point (PT: 1420±241 cells, range = 971–1415 for psilocybin and 2568 for control; IT: 12633±2430 cells, 7354–13102 and 24123; CIN: 3318±516 cells, 2308–3415 and 5773). To identify specific DEGs, we used pseudobulk analysis where sequence reads were summed for each cell group and individual replicate. Differential expression analysis was performed using edgeR 53 . Pseudobulk analysis has been shown to yield reliable quantifications with fewer false positives 54 , 55 . We first determined DEGs in the excitatory pyramidal neurons ( Figure 3a ). For PT neurons, we detected a high number of DEGs (1698 with FDR < 0.05) at 1 hour after psilocybin administration ( Figure 3b ). Psilocybin caused both upregulation and downregulation at this 1-hour time point (57% and 43% respectively among DEGs). Subsequently, the drug’s effect on transcriptional activity waned precipitously with only 108, 103, and 14 DEGs detected at 2, 4, and 24 hours after psilocybin administration. Intriguingly, the number of DEGs climbed back to 822 at 72 hours after drug exposure (35% upregulated and 65% downregulated among DEGs). A similar time course for DEGs was observed in IT neurons ( Figure 3b ). Psilocybin evoked 3401 DEGs (FDR < 0.05) in IT neurons at the 1-hour time point (52% upregulated and 48% downregulated), which dropped to 887, 192, and 6 DEGs at the 2, 4, and 24-hour time points, before rising to 2549 DEGs at the 72-hour time point (43% upregulated and 57% downregulated). Download figure Open in new tab Fig. 3: Psilocybin-induced transcriptional changes in frontal cortical excitatory neurons are involved in synaptic plasticity (a) Schematic of the cortical microcircuit highlighting the PT and IT excitatory pyramidal neurons. (b) Differentially expressed genes (DEGs) identified by pseudobulk analysis for each time point for frontal cortical PT neurons and IT neurons. Each circle represents an individual gene. Red, upregulated with FDR < 0.05. Blue, downregulated with FDR < 0.05. (c) Top 15 terms from gene ontology enrichment analysis based on the DEGs, ranked based on the mean adjusted p-values across time points in both PT and IT neurons. n.p., not present. Grey, gene ontology enrichment analysis was not performed for time points with <100 DEGs with FDR < 0.05. (d) Expression levels of specific transcripts in L5-PT, L4/5-IT, L4-IT, L2/3-IT, L5/6-IT cell types following psilocybin administration, normalized to controls. To gain insights into the functions of the psilocybin-induced transcriptional changes, we performed gene ontology (GO) enrichment analysis for the DEGs identified in PT and IT neurons at different time points. Figure 3c shows the top 15 overrepresented GO terms when we considered the prevalence of each term averaged across all time points and for both excitatory cell types. The top terms involve the regulation of glucocorticoid receptor signaling in the brain and numerous processes related to synaptic plasticity and neuronal connectivity, including the regulation of neuronal synaptic plasticity, regulation of neuron projection development, axon guidance, neuron projection guidance, axon extension involved in axon guidance, neuron projection extension involved in neuron projection guidance, axogenesis, and regulation of short-term neuronal synaptic plasticity. The engagement of mechanisms included in these terms was time-dependent (main effect of time: P = 1.6 x 10 -12 , main effect of cell type: P = 0.03, two-way ANOVA with post-hoc pairwise statistics reported; Supplementary Table 2 ). We visualized select upregulated genes that are part of the enriched pathways such as Camk1d , Grik3 , Grin2b , Grin2a , Shisa9 , and Shisa7 (mean log 2 (fold-change) = +0.351, +0.303, +0.393, +0.312, +0.296, +0.446, respectively; Figure 3d ). Camk1d is a member of the calcium/calmodulin-dependent protein kinase 1 family that regulates dendritogenesis 56 and NMDA receptor-mediated calcium elevation 57 . Grik3 , Grin2b , and Grin2a encode subunits of ionotropic glutamate receptors 58 . Shisa9 and Shisa7 encode AMPA receptor-associated proteins that modulate AMPA channel kinetics and receptor trafficking 59 , 60 . These transcriptionally upregulated genes are expected to contribute to the regulation of excitatory synaptic neurotransmission. Of note, the transcriptional changes were coordinated in the excitatory neurons, because mostly similar time courses were observed in L5-PT, L4/5-IT, L4-IT, L2/3-IT, and L5/6-IT cell types (main effect of time for Camk1 d: P = 3 x 10 -14 ; Grik3 : P = 0.1; Grin2b : P = 6 x 10 -16 ; Grin2a : P = 1 x 10 -14 ; Shisa9 : P = 2 x 10 -9 ; Shisa7 : P= 5 x 10 -10 ; main effect of cell type for Camk1 d: P = 0.002; Grik3 : P = 0.8; Grin2b : P = 0.1; Grin2a : P = 0.2; Shisa9 : P = 0.2; Shisa7 : P= 0.2; two-way ANOVA with time and cell type as factors; Figure 3d ). Altogether, these results indicate that a main transcriptional effect of psilocybin on frontal cortical excitatory neurons is to regulate synaptic function and neuronal connectivity. Psilocybin-induced transcriptional changes in frontal cortical GABAergic neurons contribute to mitochondrial function and metabolic processes We next focused on GABAergic neurons and performed pseudobulk analysis on the CINs ( Figure 4a ). The largest change in the number of DEGs in the CINs occurred immediately following psilocybin administration, with 2222 DEGs detected at the 1-hour time point (FDR < 0.05; 61% upregulated and 39% downregulated; Figure 4b ). The number of DEGs declined to 3 at 2 hours following administration, before climbing slowly to 8 and 8 at 4 and 24 hours and then surging to 1383 at 72 hours after drug exposure (37% upregulated and 63% downregulated). GO enrichment analysis revealed that DEGs in CINs were overrepresented for pathways related to mitochondrial function and metabolic processes ( Figure 4c ). The engagement of these pathways was strongest at the most acute 1-hour time point (main effect of time: P = <2 x 10 -16 , two-sample t-test; Supplementary Table 2 ). Download figure Open in new tab Fig. 4: Psilocybin-induced transcriptional changes in frontal cortical GABAergic neurons contribute to mitochondrial function and metabolic processes (a) Schematic of the cortical microcircuit highlighting the PV, SST, and VIP cell types that are grouped as cortical interneurons (CINs). (b) Differentially expressed genes (DEGs) identified by pseudobulk analysis for each time point for the CINs. Each circle represents an individual gene. Red, upregulated with FDR < 0.05. Blue, downregulated with FDR < 0.05. (c) Top 15 terms from gene ontology enrichment analysis based on the DEGs, ranked based on the mean adjusted p-values across time points in CINs. n.p., not present. Grey, gene ontology enrichment analysis was not performed for time points with <100 DEGs with FDR < 0.05. (d) Expression levels of specific transcripts in PV, SST, and VIP cell types following psilocybin administration, normalized to controls. As exemplars, we plotted genes most frequently associated with the top pathways, including Cox5b , Cox6a1 , Cox6a2, Eno2 , Ndufb5 , and Ndufs2 ( Figure 4d ). These genes had sharply increased expression in CINs 1 hour following psilocybin administration (mean log 2 (fold-change) = 0.87, 0.53, 0.33, 0.72, 0.64, and 0.72, respectively) that dropped by the 2-hour timepoint (mean log 2 (fold-change) = -0.09, 0.15, 0.02, 0.03, 0.18, and 0.09, respectively). Cox5b , Cox6a1, and Cox6a2 encode subunits of the terminal enzyme of the mitochondrial electron transport chain, cytochrome c oxidase, critical for supporting the metabolic demands of neurons through oxidative phosphorylation 61 . Of note, Cox6a2 is known to express selectively in PV interneurons only to support their morphological maturation and high firing rates 62 . Disruption of Cox6a2 has been linked to mitochondrial dysfunction in the brain 63 , 64 . Eno2 encodes neuron-specific enolase, a key enzyme in glycolysis and neuronal metabolism 65 with neuroprotective effects 66 . Ndufb5 and Ndufs2 encode NADH dehydrogenases that are subunits of Complex 1 in the mitochondrial electron transport chain, essential for mitochondrial respiration, ATP production, and oxidative stress protection 67 , 68 . Transcriptional changes were largely coordinated among these genes across PV, SST, and VIP cell types (main effect of time for Cox5b : P = 2 x 10 -8 ; Cox6a1 : P = 0.008; Cox6a2 : P = 0.5 ; Eno2 : P = 1 x 10 -5 ; Ndufb5 : P = 9 x 10 -4 ; Ndufs2 : P = 2 x 10 -5 ; main effect of cell type for Cox5b : P = 0.005; Cox6a1 : P = 0.1; Cox6a2 : P = 0.02 ; Eno2 : P = 0.7; Ndufb5 : P = 0.1; Ndufs2 : P = 0.04; two-way ANOVA with time and cell type as factors; Figure 4d ). Overall, the analyses shows that the psilocybin-induced transcriptional changes in frontal cortical GABAergic neurons are involved in acutely engaging mitochondrial function and metabolic processes. Early and late phases of transcriptional changes evoked by a single dose of psilocybin We next investigated the temporal dynamics in the psilocybin-induced emergence of DEGs. Figure 5a shows the number of DEGs for the PT neurons, IT neurons, and CINs, which varied significantly as a function of time (main effect of time: P = 4.3 x 10 -4 , main effect of cell group: P = 0.03, two-factor ANOVA with time and cell type as factor). The influence of psilocybin on the molecular program exhibited two phases with an early response at 1 – 2 hours followed by a late response detected at 72 hours after drug administration. Download figure Open in new tab Fig. 5: Early and late phases of transcriptional changes evoked by a single dose of psilocybin (a) The number of DEGs (FDR < 0.05) identified by pseudobulk analysis for each time point for PT, IT, and CIN. (b) Time course of expression changes for DEGs for frontal cortical PT neurons belonging to each of the groups, identified via hierarchical clustering based on the temporal profile of expression changes (log 2 fold change from control group at each of the 5 time points). Each line is an individual DEG. The lines are colored based on their proximity to the centroid of the group. (c) The heatmaps show the DEGs sorted based on hierarchical clustering, for frontal cortical PT neurons. Each row is an individual DEG. (d, e) Similar to (c) for frontal cortical IT neurons and CINs. We surmise that the two phases of transcriptional changes could involve either the same or different DEGs. To disambiguate these possibilities, we considered all DEGs in PT neurons, IT neurons, and CINs, and analyzed them by performing hierarchical clustering to sort them based on their temporal profile of expression changes, which revealed 4 main patterns ( Figure 5b , S3 ). For PT and IT neurons, most genes showed early upregulation (Group 1; PT neurons: 39% of 2229 DEGs, IT neurons: 35% of 4934 DEGs) or early downregulation (Group 3; PT neurons: 34%, IT neurons: 38%; Figure 5c , d ). There were fewer genes that were altered during both early and late phases (Group 2; PT neurons: 10%, IT neurons: 10%; Group 4; PT neurons: 17%, IT neurons: 16%). Likewise, most DEGs in the GABAergic interneurons exhibited early changes (33% of 3013 DEGs in Group 1 and 34% in Group 3; Figure 5e ), rather than showing a bimodal alteration across time (17% in Group 2 and 16% in Group 4). The early and late phases share only a minor fraction of the DEGs, with many unique DEGs engaged only in the acute early time point following psilocybin administration. GO enrichment analysis of each temporally profiled group revealed differences in functional enrichment ( Figure S4 ). PT neurons with early upregulation (Group 1) or early downregulation (Group 3), as well as IT neurons with early upregulation (Group 1) were the groups overrepresented by DEGs relating to synaptic function. Meanwhile, CINs that showed upregulation during early phases (Group 1) were overrepresented by DEGs relating to metabolic processes. These results provide insights into the distinct phases of transcriptional responses evoked by a single dose of psilocybin. Similar sets of differentially expressed genes evoked by psilocybin and ketamine Given that psilocybin is being investigated as a treatment for major depressive disorder, how do these transcriptional changes in the frontal cortical cell types compare to those induced by a more established rapid-acting antidepressant with a distinct molecular target? To answer this question, we generated additional data from samples collected after the administration of ketamine, which were obtained as part of the same experiment. Briefly, we administered C57/BL6J mice with ketamine (10 mg/kg, i.p.), and dissected tissues corresponding to the dorsal medial frontal cortex of both hemispheres at 1, 2, 4, 24, or 72 hours after drug administration (n = 20 animals; 2 male and 2 female for each time point; Figure 6a ). We chose 10 mg/kg for ketamine because this dose induces structural neural plasticity in the mouse frontal cortex 69 , 70 , 71 . Control data were the same no-drug samples as those used for analyzing the psilocybin data. The ketamine samples were processed using the same experimental pipeline, yielding 139,628 cells that were sequenced at a mean depth of 8687±256 reads per cell. After removing cells that did not meet quality control criteria, we had 119,815 cells from the ketamine samples (range = 21,542 - 28,008 single nuclei per time point; n = 17 animals). The 12 cell-type clusters captured the majority of cells from ketamine-treated animals: L6-CT (number of cells per time point or control condition: 3840±751), L5-PT (1418±252), L4/5-IT (3643±721), L4-IT (2565±426), L2/3-IT (4290±881), L5/6-IT (2062±374), CIN-VIP (839±136), CIN-PV (1379±240), CIN-SST (1003±156), oligodendrocytes (2386±547), astrocytes (1530±308), and OPCs (1105±218) ( Figure 6b , Supplementary Table 1 ). The same cell-type clusters were detected across samples collected at different time points after ketamine administration ( Figure S5 ). We identified differentially expressed genes in PT neurons, IT neurons, and CINs at different times points after the single dose of ketamine ( Figure S6, S7, S8, S9 ). Download figure Open in new tab Fig. 6: Similar sets of differentially expressed genes evoked by psilocybin and ketamine (a) Experimental workflow. (b) UMAP representation of the snRNA-seq data set including 119,815 nuclei from psilocybin-treated mice (range = 21,542 - 28,008 nuclei per time point; n = 17 animals) and 56,693 nuclei from the control mice (n = 8 animals). Color denotes the cell-type cluster. (c) Correlation between psilocybin- and ketamine-evoked DEGs in the PT neurons (top row), IT neurons (middle row), and CINs (bottom row) at the different time points after drug administration. Dark gray circle, individual DEGs (FDR < 0.05 in either ketamine or psilocybin condition). Light gray circle, all other genes. Colored circles highlight specific example genes. To compare the transcriptional changes evoked by psilocybin and ketamine, we generated scatter plots to correlate the drug-evoked expression changes in all genes. This analysis revealed a highly significant positive correlation in the psilocybin- and ketamine-evoked sets of altered genes in PT neurons, IT neurons, and CINs (all with r >= 0.6; Figure 6c ). Essentially, genes that were upregulated by psilocybin tended to be also upregulated by ketamine, whereas genes that were downregulated by psilocybin were similarly downregulated by ketamine. These results show similar transcriptional changes in the mouse medial frontal cortex evoked by a single dose of ketamine or psilocybin. DISCUSSION The present study applied single-nucleus transcriptomics to determine how the classic psychedelic psilocybin alters gene expression in the mouse medial frontal cortex. There are three main findings. One, excitatory and GABAergic neurons exhibit distinct sets of differentially expressed genes with roles in regulating synaptic plasticity and controlling metabolic processes respectively. Two, psilocybin’s impact is time-dependent with early and late phases of transcriptional responses. Three, psilocybin and ketamine evoke highly similar patterns of gene expression changes. Our results therefore unveil the dynamics and diversity of transcriptional responses in frontal cortical cell types after a single dose of psilocybin. Analyses of our snRNA-seq data set and the Allen Institute’s database show that most excitatory and GABAergic neurons express serotonin receptor subtypes and are thus likely direct targets of classic psychedelics. The psilocybin-evoked transcriptional changes in excitatory neurons indicate a coordinated response across PT and IT cell types on synaptic plasticity pathways. This spotlight on regulation of synaptic function is consistent with the numerous reports of psychedelic-induced structural neural remodeling in cortical pyramidal neurons 7 , 8 , 37 , 72 , 73 . In GABAergic neurons, we found that psilocybin-induced transcriptional changes were dominated by pathways related to mitochondrial and metabolic processes. This may be a direct effect, because 5-HT 2A receptor activation can influence mitochondrial biogenesis and function in cortical neurons by engaging the SIRT1-PGC-1α axis 74 . Psilocybin may engage similar mechanisms to regulate cellular metabolism in CINs. Alternatively, this may be an indirect consequence due to drug-evoked spiking activity in GABAergic neurons, as that would likely trigger an acute upregulation of oxidative stress protection, especially to support the high-energy repetitive firing in PV interneurons 62 . The bolus injection of psilocybin causes rapid elevation of drug concentration in the central nervous system followed by full clearance within a few hours. The pharmacokinetics hints at a time course of transcriptional changes, and here we uncovered an early phase at 1 – 2 hours followed by a late phase at 72 hours with high numbers of differentially expressed genes. Our findings echo the framework that there are often multiple waves of transcription following a perturbation 10 , 75 . For example, there are early and late response genes in visual cortical cells in dark-reared animals after their first experience of visual stimuli 76 . The timing of the activity-dependent transcription varies depending on the neuronal spiking pattern 77 , suggesting that different phases of transcriptional responses may relate to psilocybin’s distinct effects on the firing of different frontal cortical cell types. Furthermore, a previous study suggested that the early and sustained antidepressant effects of ketamine may depend on distinct neural mechanisms 71 . It is tempting to speculate that the multiple phases of transcriptional responses observed in this study may reflect therapeutic mechanisms at different time scales. There are limitations for this study. Transcriptomics is a rapidly growing field with numerous approaches for sequencing and analysis. Here, we chose to work with single nuclei instead of single cells, because we were concerned that cortical neurons with extended axons and dendrites may be damaged in protocols needed to dissociate single cells, which can lead to aberrant readout of gene expression 78 , 79 . Sequencing with single nuclei may better reflect recent changes in de novo transcription, and has been shown to yield excellent sensitivity and reliable classification of cell types 80 . Our results should be readily comparable to the increasing number of snRNA-seq data sets, including those collected from individuals with major depression 81 , 82 or substance use phenotypes 83 , and from animal models following exposure to drugs such as cocaine 21 , 22 and DOI 14 . A caveat is that we compared the effects of psilocybin at different time points against controls that were acquired at the “zero time point”, without drug administration. An ideal control would be vehicle-injected animals collected at the same 5 time points, but we were constrained by the costs associated with sequencing. Anticipating the importance of having high-quality data from controls for our analysis, we increased the size of the drug-free samples to 9 animals instead of the 4 animals used for each time point for each drug condition. Like psilocybin, ketamine is a compound that elicits acute psychological changes followed by longer-term antidepressant effects. Ketamine is primarily a NMDA receptor antagonist. Although psilocybin and ketamine have distinct molecular targets, they seemingly converge to produce similar effects on neural structural plasticity and behavior. In this study, a head-to-head comparison between psilocybin and ketamine reveals striking parallels in the transcriptional responses at the medial frontal cortex, suggesting shared molecular mechanisms that may underpin their behavioral effects. The overlapping patterns of transcriptional regulation raise an intriguing possibility of a transcriptional fingerprint that may be unique to rapid-acting antidepressants. We anticipate the growing number of studies revealing the molecular and cellular impact of psilocybin and ketamine will eventually reveal mechanisms of action that can accelerate the development of more effective and safer drugs for treating mental illnesses. METHODS Animals Animal care and experimental procedures were approved by the Animal Care & Use Committee (IACUC) at Yale University. C57BL/6J (Stock No. 000664) mice were purchased from Jackson Laboratory and habituated in our animal facility at Yale University for 1 week. Animals were housed in same-sex groups of 2-5 mice per cage in a temperature-controlled room, operating in normal 12 hr light - 12 hr dark cycle (8:00 AM to 8:00 PM for light). Food and water were available ad libitum. Animals were randomly assigned to different experimental groups. 20 mice (10 males and 10 females) were assigned to the psilocybin-treated group and 20 mice (10 males and 10 females) were assigned to the ketamine-treated group. 8-week-old animals were injected with psilocybin (1 mg/kg, i.p, prepared from working solution, which was made fresh monthly from powder; Usona Institute) or ketamine (10 mg/kg, i.p; Zoetis), returned to their home cage, and used 1, 2, 4, 24, or 72 hours after injection. 9 mice (3 males and 6 females) were used with no drug administration as controls. Tissue collection and generation of single-cell suspensions For tissue collection, decapitation immediately followed cervical dislocation and brains were immediately dissected from the skull. The dorsal medial frontal cortex, consisting of ACAd and medial MOs, was microdissected in cold RNAse-free 1x PBS, flash frozen in ethanol and dry ice, and stored at -80°C. Samples were shipped overnight on dry ice to the University of Michigan for processing. Nuclei were isolated and processed for single-nucleus sequencing using RNase-free materials and reagents. Fresh buffers were made immediately before nuclei processing. Buffers were made in accordance with recommended protocol (CG000393 Rev A, “Nuclei isolation from adult mouse brain tissue for single cell RNA sequencing”, 10x Genomics). Lysis buffer consisted of 10 mM Tris-HCl ph7.4 (Millipore-Sigma T2194), 10 mM NaCl (Millipore-Sigma 59222C), 3 mM MgCl 2 (Millipore-Sigma, M1028), 0.1% Nonidet P40 (Millipore-Sigma 74385), and nuclease-free water (Invitrogen AM9932). Wash buffer consisted of 1% BSA solution (Thermo Fisher Scientific AM2616), 0.2U/µl RNase Inhibitor (Millipore-Sigma 3335399001), and 1x PBS (Gibco, 10010023). Tissue was dounce homogenized in 1 mL lysis buffer and centrifuged at 500g for 5 minutes at 4°C. Supernatant was removed and samples were resuspended in 1mL of wash buffer. Samples were strained through 40 µm cell strainers (Bel-Art H13680-0040) and centrifuged again at 500 g for 5 minutes at 4°C. Supernatant was removed and samples were resuspended in 1mL of wash buffer and passed through 40 µm cell strainers once more. Nuclei were visually checked with Trypan Blue Stain (Thermo Fisher Scientific T10282) to ensure minimal debris and efficient lysis. Nuclei were stained with 0.25 µL 7AAD (Invitrogen A310) and immediately flow sorted and counted on a MoFlo Astrios (Beckman Coulter) at the North Campus Research Complex Flow Core at University of Michigan to reduce ambient RNA. Samples were centrifuged and resuspended in wash buffer at a concentration of 1,000 cells/µL. Single cell library preparation and sequencing Library prep and next-generation sequencing were carried out in the Advanced Genomics Core at the University of Michigan. Resulting libraries were sequenced on the Illumina NovaSeq S4 300 cycle, at a target depth of 25,000 reads/cell and a target of 10,000 cells/sample. Single cell quality control, clustering, and cell type identification Single nucleus RNA sequencing at the University of Michigan Advanced Genomics Core returned demultiplexed FASTQ files from all 49 samples. The 10x Genomics Mouse mm10 (GENCODE vM23/Ensembl98) reference version 2020-A was downloaded as a reference genome. FASTQ files were processed with 10x Genomics cellranger count version 7.1.0 and aligned to this reference, including intronic reads. Preprocessing and further analysis relied on anndata 84 (version 0.8.0) and scanpy 85 (version 1.10.3). Filtered feature-barcode matrices for all samples and all drug treatments were imported to AnnData objects and concatenated using anndata.concat() on the intersection of genes. One sample was not included in the study due to low cell numbers, 4 were excluded for oversequencing, and 2 were dropped for large batch effects. A total of 42 samples that met all quality control measures were further processed. All genes detected in fewer than 3 unique nuclei or having fewer than 5 total transcript counts were removed. All nuclei with fewer than 3 unique genes detected or having fewer than 5 total transcript counts were removed. Multiplet nuclei were identified and removed with scanpy.pp.scrublet() including a batch key identifying individual samples 86 . Nuclei with low transcript counts on a per-sample basis (20% of counts or more than 3 mean absolute deviations above the median) were removed. Ribosomal genes based on the Ribosomal Protein Gene Database 87 ( http://ribosome.med.miyazaki-u.ac.jp/ ) were removed from analysis. The ubiquitously overexpressed Malat1 gene was removed from analysis. All nuclei libraries were randomly downsampled without replacement to 1000 transcripts with scanpy.pp.downsample_counts(), and nuclei with less than 1000 transcripts were removed. Each cell’s reads were normalized to counts per million (CPM) with scanpy.pp.normalize_total(), and a natural logarithm was taken with scanpy.pp.log1p() 88 . This yielded ln(CPM+1), which was used for all downstream expression analysis except pseudobulk. We selected the top 5000 highly variable genes using scanpy.pp.highly_variable_genes() with flavor=‘seurat’ and these were used to compute the top 236 principal components, which encompassed 50% of the overall variation in the highly variable genes. The nearest neighbors graph was constructed with scanpy’s batch corrected K nearest neighbors algorithm scanpy.pp.bbknn(). This was input to Uniform Manifold Approximation Projection 89 (UMAP) and Leiden clustering with resolution 1.0. Cell types were hand-annotated based on the marker gene results from scanpy.tl.rank_genes_groups() with the Wilcoxon rank-sum test. Cell type identification was based on marker gene expression congruence with earlier work 90 . Example marker genes for L2/3-IT: Rasgrf2 , Calb1 ; L6-CT: Foxp2 , Hs3st4 ; L4/5-IT: Il1rapl2 ; L4-IT: Rorb ; L5/6-IT: Galnt14 , Sulf1 ; L5-PT: Coro6 , Parm1 ; CIN-PV: Pvalb ; CIN-SST: Sst ; CIN-VIP: Vip ; oligodendrocytes: Mbp , Plp1 ; astrocytes: Slc1a2 , Clu ; OPC: Olig2 . Cell-type specific analysis We performed more in-depth analysis on the response relationship of the following cell types: excitatory neurons (L2/3-IT, L6-CT, L4/5-IT, L4-IT, L5/6-IT, L5-PT), cortical interneurons (CIN-PV, CIN-SST, CIN-VIP), and glial cells (oligodendrocytes, astrocytes, OPC). We use hierarchical clustering to classify these cell types based on transcriptional profiles over all cells with scanpy.tl.dendrogram() and we recomputed the top marker genes based on these cell types only with scanpy.tl.rank_genes_groups(). We plotted the expression of the top 3 marker genes for each cell type with scanpy.pl.heatmap(). 5-HT receptor expression To examine cell-type specific expression patterns of psilocybin’s main 5-HT receptor targets, we extracted Htr1a , Htr2a , and Htr2c transcript levels from the control samples in our snRNA-seq data set (i.e., control samples from no-drug condition). Cell type-specific patterns were visualized on a per-cell basis on the data set’s UMAP representation. Group-wise variation among distinct cell subtypes (see below) was visualized with violin plots. UMAPs were created using matplotlib and the scanpy UMAP embeddings of the dataset, while violin plots were constructed using scanpy.pl.violin(). The scale of expression was unified across all genes by identifying the maximum expression level. Transcript levels from our snRNA-seq dataset were compared to publicly available single-cell sequencing data set from the Allen Institute for Brain Science ( https://portal.brain-map.org/atlases-and-data/rnaseq/mouse-whole-cortex-and-hippocampus-smart-seq ). Briefly, this independent dataset used SMARTseq-v4 to analyze 76,381 single cell transcriptomes across the mouse cortex and hippocampus, and cluster these cells into transcriptional cell subtypes 46 , 91 . We selected the 17,461 cells sampled in frontal cortex regions (ACA, ALM, ORB, and PL-ILA), including 84 micro_PVM, 118 oligo, 277 astro, 1403 L2/3 IT, 3118 L4/5 IT, 1541 L5 IT, 640 L6 IT, 471 L5 PT, 2159 L6 CT, 1009 L5/6 NP, 579 L6b, 1174 Pvalb, 382 Sncg, 1708 Vip, 126 Sst Chodl, 1573 Sst, and 1099 Lamp5 as identified and annotated already by the Allen Institute. Transcript levels for 5-HT receptors ( Htr1a , Htr1b , Htr1d , Htr2a , Htr2b , Htr2c , Htr3a , Htr3b , Htr4 , Htr5a , Htr5b , Htr6 , Htr7 ), as well as cell type markers for excitatory ( Slc17a7 ) or inhibitory ( Gad1 , Pvalb , Cck , Vip , Sst , Ndnf ) neurons were extracted for each cell subtype. Expression levels were depicted in a dot plot showing log 2 (CPM+1) per cell type, with dot size scaled by the percentage of cells of that type with non-zero transcript reads. Processing and visualization of Allen Institute data were done using MATLAB scripts. Differential gene expression using pseudobulk We performed primary cell-type specific differential expression analysis using pseudobulk techniques to identify genes of interest while minimizing the risk of false discoveries. For pseudobulk, raw un-normalized transcript counts were sum aggregated within cell types 92 or within a group of cell types when appropriate. Differential expression testing was performed with edgeR 53 . Genes with less than 10 total counts were excluded using edgerR::filterByExpr(). Differential expression was measured with log 2 (fold-change) and p-values were adjusted false discovery rate (FDR) correction. Genes were identified as significant with an FDR below 5%. ggplot2 and kwanlibr were used to create volcano plots. Genes with log 2 (fold-change) or -log 10 (p-value) beyond the bounds of the plot area were clamped to be displayed on the axis limits. We used two-side hypothesis test on H 0 : log 2 FC = 0 and corrected p-values with false discovery rate correction. We used an FDR threshold of 0.05 to identify significantly differentially expressed genes. DEG similarity between cell types using memento We conducted an additional cell-type specific differential expression analysis of control vs drug-treated cells using an alternative independent tool: Memento 52 . Memento uses an efficient bootstrapping method alongside method of moments estimation to compare the distributional parameters of genes’ differential expression, while making no assumptions about the underlying transcript count distribution. Memento analysis estimated both drug dependent and time dependent changes by using sex as a covariate and time after treatment as a factor of interest. Sequencing saturation was taken from the cellranger count 10x_metrics.csv output and multiplied by 0.25 to obtain the capture rate, as instructed by memento documentation. We used the memento ht_1d_moments function with 5000 bootstraps. Differential expression was measured with log2 fold change and p-values were adjusted False Discovery Rate correction. We used two-side hypothesis test on H 0 : log 2 FC = 0 . We used an FDR threshold of 0.05 to identify significantly differentially expressed genes. We generated heatmaps to show the similarity in differential expression between cell types within a certain drug treatment and time after treatment. For each pair of cell types, we retrieved their memento differential expression results compared to control from the corresponding timepoint. For each cell type, we found the 1000 differentially expressed genes with the most significant FDR. We computed cosine similarity on the log 2 (fold-change) of these vectors, filling in missing values from the join with zeros for compatibility between the vectors. We performed this analysis at each treatment time point (1, 2, 4, 24, and 72 h) and for each drug (psilocybin, ketamine). Grouping of cell types for downstream analysis To support our choice to group IT neurons together, CIN neurons together, and analyze other cell types individually, we assessed many potential clustering of cell types. The Calinski-Harabasz index measures how tightly members within the same cluster are grouped compared to how spread out different clusters are. Higher values indicate better clustering. This metric utilizes within-cluster sum of squares (WCSS) and between-cluster sum of squares (BCSS). We computed the relative cosine distances between cell types based on the cosine similarity of their differential expression under psilocybin compared to control cells. However, both BCSS and WCSS rely on absolute embeddings of individual items. We therefore use mCHI, a modified Calinski Harabasz index here. When we define a cluster of cell types, the cluster does not have a centroid for the above reason. Therefore, instead of computing the BCSS, which measures the sum-squared distances from cluster centroids to the data centroid, we compute a modified BCSS (mBCSS). For any pair of clusters A and B, we define the distance between them as the average distance between their bipartite members. mBCSS is the square of this average distance, summed over all pairs of clusters. Let C be the clustering. Let P C be the set of all combinatorial pairs (a,b) of clusters of the clustering C. Let d(i,j) be the cosine distance between 2 cell types summed across all 5 treatment timepoints. Then Similarly, we cannot compute the distance between a cluster member and its centroid as WCSS requires, so instead we compute modified WCSS (mWCSS). mWCSS is the sum of the squared distance between two cell types, for all pairs of cell types in the same cluster and for all clusters. Let C be the clustering. Let P a be the set of all combinatorial pairs (i,j) of cell types in cluster a. Let d(i,j) be the cosine distance between 2 cell types summed across all 5 treatment timepoints. Then In computing the modified Calinsky-Harabasz index, we directly replace BCSS and WCSS with mBCSS and mWCSS respectively, and maintain the standard degrees of freedom. Then for a clustering C, For the clustering of [L6-CT, L5-PT, L4/5-IT, L4-IT, L2/3-IT, L5/6-IT] and [CIN-VIP, CIN-PV, CIN-SST], mCHI is 0.127. For the clustering of [L5-PT], [L6-CT, L4/5-IT, L4-IT, L2/3-IT, L5/6-IT], [CIN-VIP, CIN-PV, CIN-SST], mCHI is 0.203. For the clustering of [L6-CT], [L5-PT], [L4/5-IT, L4-IT, L2/3-IT, L5/6-IT], [CIN-VIP, CIN-PV, CIN-SST], mCHI is 0.278. Temporal profiles of DEGs To identify genes responsible for the early and late waves of differential expression under psilocybin, we clustered differentially expressed genes based on their temporal expression profile. Due to our interest in specific genes, data for this analysis came from pseudobulk differential expression to minimize the risk of false discoveries. For each of PT, IT, and CIN, we identified genes that were significantly differentially expressed at any time point and retrieved their profile of differentially expressed log 2 (fold-change) across each time point, yielding a vector. We accumulated this list across the 3 cell types, maintaining multiple profiles if a single gene was differentially expressed in multiple cell types. This yielded 10,176 unique differentially expressed profiles. We used hierarchical clustering, specifically bisecting K-means from sklearn to cluster the genes by their temporal profiles. We tested between 2 and 20 clusters to tune for the optimal number of clusters. Inertia, or within-cluster sum of squares, continuously decreases as expected as we increase the number of clusters, but the incremental drops in inertia from adding one more cluster becomes smaller as the number of clusters increases. This drop in efficiency led us to prioritize fewer clusters. Additionally, we saw that after 4 clusters, adding any additional clusters did not yield any new distinctive temporal profile shapes, or that new temporal profiles shared significant overlap with already existing clusters. Furthermore, adding additional clusters seemed to only split up already small clusters instead of dominant clusters, which means outlier genes were driving the clustering at higher numbers of clusters. This led us to choose 4 clusters for our analysis. After grouping genes into one of 4 clusters with bisecting K-means, we computed the Euclidean distance of each gene g to its cluster centroid, d(g). We also defined centroid proximity as how close a gene is to its centroid, inversely to distance, as prox(g) = exp(-2 d(g)). Since distance is strictly positive real numbers, this proximity metric yields values where smaller distances to centroid have proximity close to 1, while farther distances to centroid have proximity close to 0. Taking the square yields a symmetric distribution of proximity. This distribution was plotted for each cell type using seaborn.kdeplot() with default bandwidth and without common normalization. Line plots of genes for each temporal profile group were generated independently for each cell type. Darker lines have higher proximity to centroid, while lighter lines have lower proximity. Additionally, we plotted the mean differential expression of genes in each group. The width of the error band was computed using the 95th percentile interval over all genes within the group using seaborn. Heatmaps show the expression progression profiles of these genes. Red represents positive log 2 (fold-change), blue represents negative log 2 (fold-change), and white represents no change. Each row in the heatmap represents a gene. Comparison of differential expression between psilocybin and ketamine To compare the transcriptional reaction of genes between psilocybin and ketamine treatment, we directly compared differential expression of psilocybin cells to ketamine cells from the pseudobulk analysis. For each time point and the 3 major cell types we analyzed, we made a scatter plot of log 2 (fold-change) response to psilocybin on the x-axis and log 2 (fold-change) response to ketamine on the y-axis. Genes that were significantly differentially expressed (FDR < 0.05) in either psilocybin or ketamine drug treatment were colored dark grey, and some genes of interest were highlighted. We computed and plotted the linear regression of ketamine log 2 (fold-change) and psilocybin log 2 (fold-change) using scipy.stats.linregress() using all genes. Similarly, we also computed the Pearson correlation coefficient using all genes with scipy.stats.pearsonr(). Gene ontology analysis Gene ontology (GO) analysis was performed using Enrichr 93 , 94 , 95 and Gene Ontology Biological Processes 2023 database to identify biological processes that were overrepresented for genes that were differentially expressed. For each time point with a minimum of 100 DEGs (FDR < 0.05), the DEGs were analyzed using Enrichr. The adjusted p-values were computed by Enrichr using the Benjamini-Hochberg procedure for multiple hypotheses testing correction. To select the top GO terms for excitatory cells, we considered PT and IT neurons together, averaging the -log 10 (adjusted p-value) from GO analysis across PT and IT neurons and across time points to rank the GO terms. We only included GO terms that were present from GO analysis for at least two time points. To select the top GO terms for cortical interneurons, we averaged adjusted -log 10 (adjusted p-value) from GO analysis across timepoints to rank the GO terms. For statistical analysis, we performed ANOVA using the adjusted p-values. If a term was not present at a certain time point, the value was imputed to 1 for the purpose of the statistical analysis. To identify representative genes, we first determined the 3 most frequent associated genes for each of the top 30 GO terms across all time points. From this refined list, we then hand-selected genes of interest. DATA AVAILABILITY Raw sequencing data associated with this study have been submitted to the NIH Sequence Read Archive under BioProject ID PRJNA1204073 ( http://www.ncbi.nlm.nih.gov/bioproject/1204073 ). CODE AVAILABILITY Computer code to process and analyze the snRNA-seq data is provided as Snakemake scripts, Jupyter notebooks, and R scripts on GitHub ( https://github.com/thekwanlab/psilocybin_cortex_snRNAseq ). CONTRIBUTIONS C.L, K.Y.K., and A.C.K planned the study. C.L. conducted the experiment. E.O. and Y.Q. performed the analysis. C.L. assisted in the analysis. N.K.S. and A.C.K. analyzed the data from Allen Institute. M.J.G. contributed to the planning of the study. C.L., K.Y.K., and A.C.K. drafted the manuscript. All authors reviewed the manuscript before submission. DISCLOSURES A.C.K. has been a scientific advisor or consultant for Boehringer Ingelheim, Empyrean Neuroscience, Freedom Biosciences, and Psylo. A.C.K. has received research support from Intra-Cellular Therapies. The other authors report no financial relationships with commercial interests. SUPPLEMENTAL TABLES Table S1. The number of cells in each experimental group. Table S2. Statistical analyses associated with Figure 3c , 3d , 4c , 4d , and 5a . Table S3. Details for gene ontology (GO) enrichment analyses including pathway and ID. ACKNOWLEDGEMENTS We thank Rosemarie Terwilliger and Mandy Lam for coordination assistance and guidance on sample preparation. Psilocybin was generously provided by Usona Institute’s Investigational Drug & Material Supply Program; the Usona Institute IDMSP is supported by Alexander Sherwood, Robert Kargbo, and Kristi Kaylo in Madison, WI. This work was supported by NIH grants R01MH121848 (A.C.K.), R01MH128217 (A.C.K.), R01MH137047 (A.C.K.), R01NS125313 (K.Y.K.), R01MH126287 (K.Y.K.), One Mind – COMPASS Rising Star Award (A.C.K.), NIH training grant T32NS041228 (C.L.), NIH fellowships F30MH129085 (N.K.S.), and NCI P30CA046592 (Cancer Center Shared Resource: Single Cell and Spatial Analysis Shared Resource). Footnotes http://www.ncbi.nlm.nih.gov/bioproject/1204073 https://github.com/thekwanlab/psilocybin_cortex_snRNAseq REFERENCES 1. ↵ Goodwin GM , et al. Single-Dose Psilocybin for a Treatment-Resistant Episode of Major Depression . N Engl J Med 387 , 1637 – 1648 ( 2022 ). OpenUrl CrossRef PubMed 2. ↵ Davis AK , et al. Effects of Psilocybin-Assisted Therapy on Major Depressive Disorder: A Randomized Clinical Trial . JAMA Psychiatry 78 , 481 – 489 ( 2021 ). OpenUrl CrossRef PubMed 3. ↵ Carhart-Harris R , et al. Trial of Psilocybin versus Escitalopram for Depression . N Engl J Med 384 , 1402 – 1411 ( 2021 ). OpenUrl CrossRef PubMed 4. ↵ Raison CL , et al. Single-Dose Psilocybin Treatment for Major Depressive Disorder: A Randomized Clinical Trial . JAMA 330 , 843 – 853 ( 2023 ). OpenUrl CrossRef PubMed 5. ↵ von Rotz R , et al. Single-dose psilocybin-assisted therapy in major depressive disorder: A placebo-controlled, double-blind, randomised clinical trial . EClinicalMedicine 56 , 101809 ( 2023 ). 6. ↵ Liao C , Dua AN , Wojtasiewicz C , Liston C , Kwan AC . Structural neural plasticity evoked by rapid-acting antidepressant interventions . Nat Rev Neurosci , in press ( 2024 ). 7. ↵ Ly C , et al. Psychedelics Promote Structural and Functional Neural Plasticity . Cell Rep 23 , 3170 – 3182 ( 2018 ). OpenUrl CrossRef PubMed 8. ↵ Shao LX , et al. Psilocybin induces rapid and persistent growth of dendritic spines in frontal cortex in vivo . Neuron 109 , 2535 – 2544 e2534 ( 2021 ). OpenUrl CrossRef PubMed 9. ↵ Shao LX , et al. Pyramidal cell types and 5-HT2A receptors are essential for psilocybin’s lasting drug action . biorxiv , ( 2024 ). 10. ↵ Yap EL , Greenberg ME . Activity-Regulated Transcription: Bridging the Gap between Neural Activity and Behavior . Neuron 100 , 330 – 348 ( 2018 ). OpenUrl CrossRef PubMed 11. ↵ Ma H , et al. Excitation–transcription coupling, neuronal gene expression and synaptic plasticity . Nature Reviews Neuroscience 24 , 672 – 692 ( 2023 ). OpenUrl CrossRef PubMed 12. ↵ Gonzalez-Maeso J , et al. Transcriptome fingerprints distinguish hallucinogenic and nonhallucinogenic 5-hydroxytryptamine 2A receptor agonist effects in mouse somatosensory cortex . J Neurosci 23 , 8836 – 8843 ( 2003 ). OpenUrl Abstract / FREE Full Text 13. ↵ Desouza LA , et al. The Hallucinogenic Serotonin2A Receptor Agonist, 2,5-Dimethoxy-4-Iodoamphetamine, Promotes cAMP Response Element Binding Protein-Dependent Gene Expression of Specific Plasticity-Associated Genes in the Rodent Neocortex . Front Mol Neurosci 14 , 790213 ( 2021 ). OpenUrl CrossRef PubMed 14. ↵ Muir J , et al. Isolation of psychedelic-responsive neurons underlying anxiolytic behavioral states . Science 386 , 802 – 810 ( 2024 ). OpenUrl CrossRef PubMed 15. ↵ Gonzalez-Maeso J , et al. Hallucinogens recruit specific cortical 5-HT(2A) receptor-mediated signaling pathways to affect behavior . Neuron 53 , 439 – 452 ( 2007 ). OpenUrl CrossRef PubMed Web of Science 16. ↵ Nichols CD , Sanders-Bush E . A single dose of lysergic acid diethylamide influences gene expression patterns within the mammalian brain . Neuropsychopharmacology 26 , 634 – 642 ( 2002 ). OpenUrl CrossRef PubMed Web of Science 17. ↵ Nogueira M , et al. Serotonergic psychedelic 5-MeO-DMT alters plasticity-related gene expression and generates anxiolytic effects in stressed mice . Mol Psychiatry , ( 2024 ). 18. ↵ Jefsen OH , Elfving B , Wegener G , Muller HK . Transcriptional regulation in the rat prefrontal cortex and hippocampus after a single administration of psilocybin . J Psychopharmacol 35 , 483 – 493 ( 2021 ). OpenUrl CrossRef PubMed 19. ↵ Donovan LL , et al. Effects of a single dose of psilocybin on behaviour, brain 5-HT2A receptor occupancy and gene expression in the pig . Eur Neuropsychopharmacol 42 , 1 – 11 ( 2021 ). OpenUrl CrossRef PubMed 20. ↵ Walker DM , et al. Cocaine Self-administration Alters Transcriptome-wide Responses in the Brain’s Reward Circuitry . Biol Psychiatry 84 , 867 – 880 ( 2018 ). OpenUrl CrossRef PubMed 21. ↵ Savell KE , et al. A dopamine-induced gene expression signature regulates neuronal function and cocaine response . Sci Adv 6 , eaba4221 ( 2020 ). OpenUrl FREE Full Text 22. ↵ Phillips RA , 3rd , et al. Distinct subpopulations of D1 medium spiny neurons exhibit unique transcriptional responsiveness to cocaine . Mol Cell Neurosci 125 , 103849 ( 2023 ). OpenUrl CrossRef PubMed 23. ↵ Browne CJ, et al. Transcriptional signatures of heroin intake and relapse throughout the brain reward circuitry in male mice . Sci Adv 9 , eadg8558 ( 2023 ). OpenUrl CrossRef PubMed 24. ↵ Engeln M , et al. Transcriptome profiling of the ventral pallidum reveals a role for pallido-thalamic neurons in cocaine reward . Mol Psychiatry 27 , 3980 – 3991 ( 2022 ). OpenUrl CrossRef PubMed 25. ↵ Nardou R , et al. Psychedelics reopen the social reward learning critical period . Nature 618 , 790 – 798 ( 2023 ). OpenUrl CrossRef PubMed 26. ↵ Bagot RC , et al. Ketamine and Imipramine Reverse Transcriptional Signatures of Susceptibility and Induce Resilience-Specific Gene Expression Profiles . Biol Psychiatry 81 , 285 – 295 ( 2017 ). OpenUrl CrossRef PubMed 27. ↵ Walker DM , et al. Sex-Specific Transcriptional Changes in Response to Adolescent Social Stress in the Brain’s Reward Circuitry . Biol Psychiatry 91 , 118 – 128 ( 2022 ). OpenUrl CrossRef PubMed 28. ↵ Bagot RC , et al. Circuit-wide Transcriptional Profiling Reveals Brain Region-Specific Gene Networks Regulating Depression Susceptibility . Neuron 90 , 969 – 983 ( 2016 ). OpenUrl CrossRef PubMed 29. ↵ Pena CJ , et al. Early life stress confers lifelong stress susceptibility in mice via ventral tegmental area OTX2 . Science 356 , 1185 – 1188 ( 2017 ). OpenUrl Abstract / FREE Full Text 30. ↵ Pena CJ , et al. Early life stress alters transcriptomic patterning across reward circuitry in male and female mice . Nat Commun 10 , 5098 ( 2019 ). OpenUrl CrossRef PubMed 31. ↵ Vollenweider FX , Vollenweider-Scherpenhuyzen MF , Babler A , Vogel H , Hell D . Psilocybin induces schizophrenia-like psychosis in humans via a serotonin-2 agonist action . Neuroreport 9 , 3897 – 3902 ( 1998 ). OpenUrl CrossRef PubMed Web of Science 32. ↵ Madsen MK , et al. Psychedelic effects of psilocybin correlate with serotonin 2A receptor occupancy and plasma psilocin levels . Neuropsychopharmacology 44 , 1328 – 1334 ( 2019 ). OpenUrl CrossRef PubMed 33. ↵ Hasler F , Bourquin D , Brenneisen R , Bar T , Vollenweider FX . Determination of psilocin and 4-hydroxyindole-3-acetic acid in plasma by HPLC-ECD and pharmacokinetic profiles of oral and intravenous psilocybin in man . Pharm Acta Helv 72 , 175 – 184 ( 1997 ). OpenUrl CrossRef PubMed 34. ↵ Preller KH , et al. Psilocybin Induces Time-Dependent Changes in Global Functional Connectivity . Biol Psychiatry 88 , 197 – 207 ( 2020 ). OpenUrl CrossRef PubMed 35. ↵ Brown RT , et al. Pharmacokinetics of Escalating Doses of Oral Psilocybin in Healthy Adults . Clin Pharmacokinet 56 , 1543 – 1554 ( 2017 ). OpenUrl CrossRef PubMed 36. ↵ Halberstadt AL , Koedood L , Powell SB , Geyer MA . Differential contributions of serotonin receptors to the behavioral effects of indoleamine hallucinogens in mice . J Psychopharmacol 25 , 1548 – 1561 ( 2011 ). OpenUrl CrossRef PubMed Web of Science 37. ↵ Jefferson SJ , et al. 5-MeO-DMT modifies innate behaviors and promotes structural neural plasticity in mice. Neuropsychopharmacology , Online ahead of print ( 2023 ). 38. ↵ Torrado Pacheco A , Olson RJ , Garza G , Moghaddam B . Acute psilocybin enhances cognitive flexibility in rats . Neuropsychopharmacology 48 , 1011 – 1020 ( 2023 ). OpenUrl CrossRef PubMed 39. ↵ Hibicke M , Landry AN , Kramer HM , Talman ZK , Nichols CD . Psychedelics, but Not Ketamine, Produce Persistent Antidepressant-like Effects in a Rodent Experimental System for the Study of Depression . ACS Chem Neurosci 11 , 864 – 871 ( 2020 ). OpenUrl CrossRef PubMed 40. ↵ Cameron LP , et al. 5-HT2ARs Mediate Therapeutic Behavioral Effects of Psychedelic Tryptamines . ACS Chem Neurosci 14 , 351 – 358 ( 2023 ). OpenUrl CrossRef PubMed 41. ↵ Hesselgrave N , Troppoli TA , Wulff AB , Cole AB , Thompson SM . Harnessing psilocybin: antidepressant-like behavioral and synaptic actions of psilocybin are independent of 5-HT2R activation in mice . Proc Natl Acad Sci U S A 118 , ( 2021 ). 42. ↵ Woodburn SC , Levitt CM , Koester AM , Kwan AC . Psilocybin Facilitates Fear Extinction: Importance of Dose, Context, and Serotonin Receptors . ACS Chem Neurosci 15 , 3034 – 3043 ( 2024 ). OpenUrl CrossRef PubMed 43. ↵ Davoudian PA , Shao LX , Kwan AC . Shared and Distinct Brain Regions Targeted for Immediate Early Gene Expression by Ketamine and Psilocybin . ACS Chem Neurosci 14 , 468 – 480 ( 2023 ). OpenUrl CrossRef PubMed 44. ↵ Nichols DE. Psychedelics . Pharmacol Rev 68 , 264 – 355 ( 2016 ). OpenUrl Abstract / FREE Full Text 45. ↵ Savalia NK , Shao LX , Kwan AC . A Dendrite-Focused Framework for Understanding the Actions of Ketamine and Psychedelics . Trends Neurosci 44 , 260 – 275 ( 2021 ). OpenUrl CrossRef PubMed 46. ↵ Yao Z , et al. A taxonomy of transcriptomic cell types across the isocortex and hippocampal formation . Cell 184 , 3222 – 3241.e3226 ( 2021 ). OpenUrl CrossRef PubMed 47. ↵ Zolnik TA , et al. Layer 6b controls brain state via apical dendrites and the higher-order thalamocortical system . Neuron 112 , 805 – 820 e804 ( 2024 ). OpenUrl CrossRef PubMed 48. ↵ Schuman B , Machold RP , Hashikawa Y , Fuzik J , Fishell GJ , Rudy B . Four Unique Interneuron Populations Reside in Neocortical Layer 1 . J Neurosci 39 , 125 – 139 ( 2019 ). OpenUrl Abstract / FREE Full Text 49. ↵ Krabbe G , Matyash V , Pannasch U , Mamer L , Boddeke HW , Kettenmann H . Activation of serotonin receptors promotes microglial injury-induced motility but attenuates phagocytic activity . Brain Behav Immun 26 , 419 – 428 ( 2012 ). OpenUrl CrossRef PubMed 50. ↵ Kolodziejczak M , et al. Serotonin Modulates Developmental Microglia via 5-HT2B Receptors: Potential Implication during Synaptic Refinement of Retinogeniculate Projections . ACS Chem Neurosci 6 , 1219 – 1230 ( 2015 ). OpenUrl CrossRef PubMed 51. ↵ Carson MJ , Thomas EA , Danielson PE , Sutcliffe JG . The 5-HT5A serotonin receptor is expressed predominantly by astrocytes in which it inhibits cAMP accumulation: A mechanism for neuronal suppression of reactive astrocytes . Glia 17 , 317 – 326 ( 1996 ). OpenUrl CrossRef PubMed 52. ↵ Kim MC , et al. Method of moments framework for differential expression analysis of single-cell RNA sequencing data . Cell 187 , 6393 – 6410 e6316 ( 2024 ). OpenUrl CrossRef PubMed 53. ↵ Chen Y , Lun AT , Smyth GK . From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline . F1000Res 5 , 1438 ( 2016 ). OpenUrl 54. ↵ Squair JW , et al. Confronting false discoveries in single-cell differential expression . Nat Commun 12 , 5692 ( 2021 ). OpenUrl CrossRef PubMed 55. ↵ Murphy AE , Skene NG . A balanced measure shows superior performance of pseudobulk methods in single-cell RNA-sequencing analysis . Nat Commun 13 , 7851 ( 2022 ). OpenUrl CrossRef PubMed 56. ↵ Takemoto-Kimura S , et al. Calmodulin kinases: essential regulators in health and disease . J Neurochem 141 , 808 – 818 ( 2017 ). OpenUrl CrossRef PubMed 57. ↵ Schmitt JM , Guire ES , Saneyoshi T , Soderling TR . Calmodulin-dependent kinase kinase/calmodulin kinase I activity gates extracellular-regulated kinase-dependent long-term potentiation . J Neurosci 25 , 1281 – 1290 ( 2005 ). OpenUrl Abstract / FREE Full Text 58. ↵ Myers SJ , Yuan H , Kang JQ , Tan FCK , Traynelis SF , Low CM . Distinct roles of GRIN2A and GRIN2B variants in neurological conditions . F1000Res 8 , ( 2019 ). 59. ↵ von Engelhardt J , et al. CKAMP44: a brain-specific protein attenuating short-term synaptic plasticity in the dentate gyrus . Science 327 , 1518 – 1522 ( 2010 ). OpenUrl Abstract / FREE Full Text 60. ↵ Schmitz LJM , et al. The AMPA receptor-associated protein Shisa7 regulates hippocampal synaptic function and contextual memory . Elife 6 , ( 2017 ). 61. ↵ Arnold S . The power of life--cytochrome c oxidase takes center stage in metabolic control, cell signalling and survival . Mitochondrion 12 , 46 – 56 ( 2012 ). OpenUrl CrossRef PubMed Web of Science 62. ↵ Sanz-Morello B , et al. Complex IV subunit isoform COX6A2 protects fast-spiking interneurons from oxidative stress and supports their function . EMBO J 39 , e105759 ( 2020 ). OpenUrl CrossRef PubMed 63. ↵ Holper L , Lan MJ , Brown PJ , Sublette EM , Burke A , Mann JJ . Brain cytochrome-c-oxidase as a marker of mitochondrial function: A pilot study in major depression using NIRS . Depress Anxiety 36 , 766 – 779 ( 2019 ). OpenUrl CrossRef PubMed 64. ↵ Khadimallah I , et al. Mitochondrial, exosomal miR137-COX6A2 and gamma synchrony as biomarkers of parvalbumin interneurons, psychopathology, and neurocognition in schizophrenia . Mol Psychiatry 27 , 1192 – 1204 ( 2022 ). OpenUrl CrossRef PubMed 65. ↵ Lee BC , Tsai JC , Huang YH , Wang CC , Lee HC , Tsai HJ . The 419th Aspartic Acid of Neural Membrane Protein Enolase 2 Is a Key Residue Involved in the Axonal Growth of Motor Neurons Mediated by Interaction between Enolase 2 Receptor and Extracellular Pgk1 Ligand . Int J Mol Sci 25 , ( 2024 ). 66. ↵ Haque A , Polcyn R , Matzelle D , Banik NL . New Insights into the Role of Neuron-Specific Enolase in Neuro-Inflammation, Neurodegeneration, and Neuroprotection . Brain Sci 8 , ( 2018 ). 67. ↵ Bandara AB , Drake JC , James CC , Smyth JW , Brown DA . Complex I protein NDUFS2 is vital for growth, ROS generation, membrane integrity, apoptosis, and mitochondrial energetics . Mitochondrion 58 , 160 – 168 ( 2021 ). OpenUrl CrossRef PubMed 68. ↵ Dunham-Snary KJ , et al. Ndufs2, a Core Subunit of Mitochondrial Complex I, Is Essential for Acute Oxygen-Sensing and Hypoxic Pulmonary Vasoconstriction . Circ Res 124 , 1727 – 1746 ( 2019 ). OpenUrl CrossRef PubMed 69. ↵ Phoumthipphavong V , Barthas F , Hassett S , Kwan AC . Longitudinal Effects of Ketamine on Dendritic Architecture In Vivo in the Mouse Medial Frontal Cortex . eNeuro 3 , ENEURO.0133-0115.2016 ( 2016 ). 70. ↵ Ng LHL , Huang Y , Han L , Chang RC , Chan YS , Lai CSW . Ketamine and selective activation of parvalbumin interneurons inhibit stress-induced dendritic spine elimination . Transl Psychiatry 8 , 272 ( 2018 ). 71. ↵ Moda-Sava RN , et al. Sustained rescue of prefrontal circuit dysfunction by antidepressant-induced spine formation . Science 364 , eaat8078 ( 2019 ). OpenUrl Abstract / FREE Full Text 72. ↵ de la Fuente Revenga M , et al. Prolonged epigenomic and synaptic plasticity alterations following single exposure to a psychedelic in mice . Cell Rep 37 , 109836 ( 2021 ). 73. ↵ Cameron LP , et al. A non-hallucinogenic psychedelic analogue with therapeutic potential . Nature 589 , 474 – 479 ( 2021 ). OpenUrl CrossRef PubMed 74. ↵ Fanibunda SE , et al. Serotonin regulates mitochondrial biogenesis and function in rodent cortical neurons via the 5-HT(2A) receptor and SIRT1-PGC-1alpha axis . Proc Natl Acad Sci U S A 116 , 11028 – 11037 ( 2019 ). OpenUrl Abstract / FREE Full Text 75. ↵ Tullai JW , Schaffer ME , Mullenbrock S , Sholder G , Kasif S , Cooper GM . Immediate-early and delayed primary response genes are distinct in function and genomic architecture . J Biol Chem 282 , 23981 – 23995 ( 2007 ). OpenUrl Abstract / FREE Full Text 76. ↵ Hrvatin S , et al. Single-cell analysis of experience-dependent transcriptomic states in the mouse visual cortex . Nat Neurosci 21 , 120 – 129 ( 2018 ). OpenUrl CrossRef PubMed 77. ↵ Tyssowski KM , et al. Different Neuronal Activity Patterns Induce Different Gene Expression Programs . Neuron 98 , 530 – 546 e511 ( 2018 ). OpenUrl CrossRef PubMed 78. ↵ van den Brink SC , et al. Single-cell sequencing reveals dissociation-induced gene expression in tissue subpopulations . Nat Methods 14 , 935 – 936 ( 2017 ). OpenUrl CrossRef PubMed 79. ↵ Lacar B , et al. Nuclear RNA-seq of single neurons reveals molecular signatures of activation . Nat Commun 7 , 11022 ( 2016 ). 80. ↵ Ding J , et al. Systematic comparison of single-cell and single-nucleus RNA-sequencing methods . Nat Biotechnol 38 , 737 – 746 ( 2020 ). OpenUrl CrossRef PubMed 81. ↵ Nagy C , et al. Single-nucleus transcriptomics of the prefrontal cortex in major depressive disorder implicates oligodendrocyte precursor cells and excitatory neurons . Nature Neuroscience 23 , 771 – 781 ( 2020 ). OpenUrl CrossRef PubMed 82. ↵ Zeng L , et al. A Single-Nucleus Transcriptome-Wide Association Study Implicates Novel Genes in Depression Pathogenesis . Biol Psychiatry 96 , 34 – 43 ( 2024 ). OpenUrl CrossRef PubMed 83. ↵ Tran MN , et al. Single-nucleus transcriptome analysis reveals cell-type-specific molecular signatures across reward circuitry in the human brain . Neuron 109 , 3088 – 3103 e3085 ( 2021 ). OpenUrl CrossRef PubMed 84. ↵ Virshup I , Rybakov S , Theis FJ , Angerer P , Wolf FA. anndata: Annotated data . biorxiv , ( 2021 ). 85. ↵ Wolf FA , Angerer P , Theis FJ . SCANPY: large-scale single-cell gene expression data analysis . Genome Biol 19 , 15 ( 2018 ). 86. ↵ Wolock SL , Lopez R , Klein AM . Scrublet: Computational Identification of Cell Doublets in Single-Cell Transcriptomic Data . Cell Syst 8 , 281 – 291 e289 ( 2019 ). OpenUrl CrossRef PubMed 87. ↵ Nakao A , Yoshihama M , Kenmochi N . RPG: the Ribosomal Protein Gene database . Nucleic Acids Res 32 , D168 – 170 ( 2004 ). OpenUrl CrossRef PubMed Web of Science 88. ↵ Herring CA , et al. Human prefrontal cortex gene regulatory dynamics from gestation to adulthood at single-cell resolution . Cell 185 , 4428 – 4447 e4428 ( 2022 ). OpenUrl CrossRef PubMed 89. ↵ McInnes L , Healy J , Melville J . UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction . biorxiv , ( 2020 ). 90. ↵ Fischer S , Gillis J . How many markers are needed to robustly determine a cell’s type? iScience 24 , 103292 ( 2021 ). 91. ↵ Tasic B , et al. Shared and distinct transcriptomic cell types across neocortical areas . Nature 563 , 72 – 78 ( 2018 ). OpenUrl CrossRef PubMed 92. ↵ Hafemeister C , Halbritter F . Single-cell RNA-seq differential expression tests within a sample should use pseudo-bulk data of pseudo-replicates . biorxiv , ( 2023 ). 93. ↵ Chen EY , et al. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool . BMC Bioinformatics 14 , 128 ( 2013 ). 94. ↵ Kuleshov MV , et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update . Nucleic Acids Res 44 , W90 – 97 ( 2016 ). OpenUrl CrossRef PubMed 95. ↵ Xie Z , et al. Gene Set Knowledge Discovery with Enrichr . Curr Protoc 1 , e90 ( 2021 ). OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted January 04, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Single-nucleus transcriptomics reveals time-dependent and cell-type-specific effects of psilocybin on gene expression Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Single-nucleus transcriptomics reveals time-dependent and cell-type-specific effects of psilocybin on gene expression Clara Liao , Ethan O’Farrell , Yaman Qalieh , Neil K. Savalia , Matthew J. Girgenti , Kenneth Y. Kwan , Alex C. Kwan bioRxiv 2025.01.04.631335; doi: https://doi.org/10.1101/2025.01.04.631335 Share This Article: Copy Citation Tools Single-nucleus transcriptomics reveals time-dependent and cell-type-specific effects of psilocybin on gene expression Clara Liao , Ethan O’Farrell , Yaman Qalieh , Neil K. Savalia , Matthew J. Girgenti , Kenneth Y. Kwan , Alex C. Kwan bioRxiv 2025.01.04.631335; doi: https://doi.org/10.1101/2025.01.04.631335 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Neuroscience Subject Areas All Articles Animal Behavior and Cognition (7629) Biochemistry (17660) Bioengineering (13881) Bioinformatics (41909) Biophysics (21435) Cancer Biology (18576) Cell Biology (25479) Clinical Trials (138) Developmental Biology (13366) Ecology (19887) Epidemiology (2067) Evolutionary Biology (24301) Genetics (15598) Genomics (22482) Immunology (17726) Microbiology (40359) Molecular Biology (17162) Neuroscience (88529) Paleontology (666) Pathology (2830) Pharmacology and Toxicology (4820) Physiology (7636) Plant Biology (15125) Scientific Communication and Education (2044) Synthetic Biology (4290) Systems Biology (9817) Zoology (2269)