Mol. Cells

Gene Expression Profiling of the Habenula in Rats Exposed to Chronic Restraint Stress

Hyeijung Yoo, Hyun Jung Kim, Soo Hyun Yang, Gi Hoon Son, Jeong-An Gim

Additional article information


Chronic stress contributes to the risk of developing depression; the habenula, a nucleus in epithalamus, is associated with many neuropsychiatric disorders. Using genome-wide gene expression analysis, we analyzed the transcriptome of the habenula in rats exposed to chronic restraint stress for 14 days. We identified 379 differentially expressed genes (DEGs) that were affected by chronic stress. These genes were enriched in neuroactive ligand-receptor interaction, the cAMP (cyclic adenosine monophosphate) signaling pathway, circadian entrainment, and synaptic signaling from the Kyoto Encyclopedia of Genes and Genomes pathway analysis and responded to corticosteroids, positive regulation of lipid transport, anterograde trans-synaptic signaling, and chemical synapse transmission from the Gene Ontology analysis. Based on protein–protein interaction network analysis of the DEGs, we identified neuroactive ligand-receptor interactions, circadian entrainment, and cholinergic synapse-related subclusters. Additionally, cell type and habenular regional expression of DEGs, evaluated using a recently published single-cell RNA sequencing study (GSE137478), strongly suggest that DEGs related to neuroactive ligand-receptor interaction and trans-synaptic signaling are highly enriched in medial habenular neurons. Taken together, our findings provide a valuable set of molecular targets that may play important roles in mediating the habenular response to stress and the onset of chronic stress-induced depressive behaviors.

Keywords: chronic restraint stress, depression, gene enrichment analysis, gene expression profiling, habenula


Recent findings indicate that the habenula is involved in the onset of psychiatric disorders such as depression, bipolar disorder, schizophrenia, and substance use disorders (Fakhoury, 2017; Hikosaka, 2010; McLaughlin et al., 2017). Despite significant advances in the evaluation of the neuroanatomical mechanisms underlying the relationship between the habenula and pathophysiology of psychiatric diseases, especially in depressive disorder and drug addiction (Lee et al., 2019; Yang et al., 2018), the molecular details of many of these changes remain unknown.

In mammals, the habenula is subdivided into the medial (MHb) and lateral (LHb) habenula, whose anatomical border is molecularly marked by varying patterns of gene expression (Aizawa et al., 2012). Both subregions have distinct roles in the stress response and onset of depressive phenotypes. For example, rats subjected to chronic restraint stress (CRS) exhibit a downregulation of MHb-specific choline acetyltransferase (Chat) and calcium-dependent secretion activator 2 (Caps2 or Cadps2) in the habenula, which has also been observed in human suicide cases of a major depressive disorder (MDD) (Han et al., 2017; Yoo et al., 2021). In addition, proteomic analysis of the LHb in rats with congenital learned helplessness revealed enhanced expression of calcium/calmodulin-dependent protein kinase II beta (CaMK2B) and astroglial potassium channel (Kir4.1), which mediate LHb hyperactivity, subsequently contributing to depressive symptoms (Cui et al., 2018; Li et al., 2013b). Furthermore, p11 (S100A10), a multifunctional protein that interacts with 5-HT (5-hydroxytryptamine; serotonin) receptors, is upregulated in mouse LHb by CRS, and p11 knockdown alleviates stress-induced depression-like behaviors (Seo et al., 2018). These molecules have been identified in animal models of depression, as altered gene expression potentially evokes depression-like behaviors in rodents. Therefore, we performed a whole transcriptome analysis of the habenula after CRS to investigate the pathogenesis of depression at the molecular level.

There have been some microarray-based transcriptomic analyses of gene expression in the habenula of animal models. For example, microarray data combined with high-throughput in situ hybridization were used to define a unique habenula profile for neurotransmitters, receptors, ion channels, and regulatory factors that are critical for habenular development using Brn3a null mouse embryos (Quina et al., 2009). In Wistar rats, chronic mild stress-induced changes in LHb gene expression were comparatively analyzed in responders and non-responders to antidepressant medication, which confirmed the importance of LHb in depression and drug responsiveness (Christensen et al., 2013). In another study, the habenular transcripts had a differentially altered expression when C57BL/6 mice were exposed to intravenous self-administration of nicotine, leading to drug-seeking behavior (Lee et al., 2015).

CRS has been widely used as a rodent model for depression, representing key features of the disorder for the study of a wide range of topics, from behavioral phenotypes to molecular signatures (Chiba et al., 2012; Xu et al., 2017; Yi et al., 2020). We have determined that CRS results in lower body weight gain, higher plasma corticosterone levels, enlarged adrenal glands, and depression-like behaviors (Han et al., 2017), consistent with previous findings (Andrus et al., 2012; Buynitsky and Mostofsky, 2009). We hypothesized that repeated exposure to stress, which is a key risk factor in the etiology of depression, may drive transcriptional changes in the habenula. Next, we examined transcriptomic profiles in the rat habenula to identify CRS-induced differentially expressed genes (DEGs). We then performed pathway enrichment and protein–protein interaction (PPI) network analyses on the DEGs to identify functionally relevant biological pathways. In addition, the expression of DEGs in enriched pathways was explored using single-cell RNA sequencing data.



All experimental procedures involving animals were approved by the Institutional Animal Care and Use Committee (IACUC) of Korea University and were performed according to the guidelines of the university (No. KOREA-2016-0059). Adult male Sprague-Dawley rats were obtained from Orient Bio (Korea). The rats were individually housed in cages under a 12-h light/dark cycle (lights on at 8 a.m.) and given ad libitum access to food and water.

Chronic restraint stress model of depression

The CRS model was established as previously reported (Andrus et al., 2012). The animals were randomly assigned to either of the two groups: the CRS group (n = 4), in which the animals were exposed to restraint stress for 2 h per day (at Zeitgeber time 2-4) for 2 weeks, and a control non-stress group (n = 4), in which the animals were kept in their home cages without being exposed to any stressors. Restraint stress was induced by placing the animals in breathable plastic film envelopes (DC 200, DecapiCones; Braintree Scientific, USA). On the last day of the stress period, the rats in the experimental and control groups were sacrificed following the last stress exposure. The CRS model was validated in our previous study (Han et al., 2017). Accordingly, we confirmed increased plasma corticosterone levels and adrenal weight in the CRS rats used in the present study (data not shown).

Microarray study

Tissue preparation

The habenula was isolated from the brain immediately after decapitation and stored at –80°C in RNAlater RNA stabilization solution (AM7020; Ambion, USA). Total RNA was extracted using the RNeasy Mini Kit (74104; Qiagen, USA), and a double-amplified cDNA was synthesized using the Ambion WT Expression Kit (4411973; Ambion) starting with 100 ng of mRNA and labeled using the GeneChip WT Terminal Labeling and Controls Kit (901525; Affymetrix, USA), according to the manufacturer’s instructions.

Microarray hybridization and scanning

Hybridization cocktails containing biotin-labeled DNA fragments for the target and hybridization controls were prepared. The target was then hybridized to the GeneChip Rat Gene 1.0 ST array (901172; Affymetrix) in a GeneChip Hybridization Oven 640 (Affymetrix) for 16 h. Immediately after hybridization, the array was washed and stained with streptavidin phycoerythrin conjugate using GeneChip Hybridization, Wash and Stain Kit (900720; Affymetrix) on the GeneChip Fluidics Station 450 according to Affymetrix recommendations. GeneChip arrays were scanned using the GeneChip Scanner 3000 7G (Affymetrix). Raw scanned images (*.DAT) and processed DAT files (*.CEL) were acquired using the Affymetrix GeneChip Command Console software (AGCC).

Microarray data processing and analysis

Raw data CEL files from Affymetrix were imported into GeneSpring Analysis software version GX12.1 (Agilent Technologies, USA). For the analysis of expression changes between stress and control samples, the files were pre-processed using the probe logarithmic intensity error, and replicas from each of the stress or control group arrays were binned using the grouping feature of the program. Normalized mean values for the two experimental groups were then generated for the experimental interpretation. Signal intensities from the chronic stress group were compared to those from the corresponding non-stress group to generate fold changes. All microarray data were minimum information about a microarray experiment (MIAME)-compliant.

Pathway enrichment and protein–protein interaction network analysis

Using selected 379 DEGs, Gene Ontology (GO) pathway analysis was performed using the gprofiler2 package (Raudvere et al., 2019) on Rattus norvegicus using ordered query by fold change, and other settings were remained as default. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed using a cluster profiler (Yu et al., 2012) on Rattus norvegicus with default settings. Statistical significance was set at adjusted P < 0.05. The PPI network was constructed using Cytoscape 3.8.0, according to the Search Tool for the Retrieval of Interacting Genes (STRING) database (Szklarczyk et al., 2019). PPIs with a confidence ≥ 0.4 were selected and visualized. The number of degrees for each node was calculated using the cytoHubba plugin (Chin et al., 2014). Subclusters of the PPI network were analyzed using the Molecular Complex Detection (MCODE) plugin, with a degree cutoff of 3 and K-core of 2 (Bader and Hogue, 2003). Pie charts for pathway enrichment analysis were generated using enhancedGraphics plugin (Morris et al., 2014).

Single-cell RNA sequencing data analysis

Habenular single-cell RNA sequencing data (Gene Expression Omnibus accession No. GSE137478) were obtained from the Gene Expression Omnibus database ( These datasets were generated on an Illumina NextSeq 500 (v2.5) (Illumina, USA) and aligned using the 10× Genomics Cell Ranger pipeline (V2). Data analysis was performed as described previously (Hashikawa et al., 2020). The genes expressed in fewer than three cells and cells with low gene expression or a high mitochondrial gene percentage (total unique molecular identifiers [UMI] < 700 or > 15,000, or mitochondrial gene percentage > 20%) were filtered out. Doublet cells were explored and removed computationally using DoubletDecon package (DePasquale et al., 2019). Normalization, scaling, and clustering using shared nearest neighbors were done using Seurat V4 package (Hao et al., 2021). Principal component analysis was performed to reduce the dimensions, and the first 30 principal components (PCs) were used for clustering (resolution = 0.8). Uniform Manifold Approximation and Projection (UMAP) was used for the visualization of the clusters on PC1 to PC30. Then the clusters were assigned to cell types and habenular regions using canonical markers (on Supplementary Fig. S1) and known habenular (Tac2, Slc17a7 for MHb marker, Pcdh10, Gabra1 for LHb marker) and peri-habenular regional markers (Slitrk6, Ramp3) (Aizawa et al., 2012; Hashikawa et al., 2020; Saunders et al., 2018; Wallace et al., 2020; Zeisel et al., 2018). Finally, the list of CRS DEGs was compared with the genes expressed in each cell type and the habenular region.

Quantitative real-time polymerase chain reaction (qRT-PCR)

qRT-PCR was performed as described previously (Han et al., 2017; Kim et al., 2020). Briefly, the habenula was isolated from the animal brain immediately after decapitation and placed in TRIzol solution (Ambion). Total RNA samples (2 μg) were reverse-transcribed into cDNA using Moloney murine leukemia virus reverse transcriptase (M-MLV RT; Promega, USA) and oligo (dT) primers (Novagen, USA). qRT-PCR was performed with 0.5 μg of the RT product in the presence of specific primer sets (Supplementary Table S1). PCR amplification with iQ SYBR Green Supermix was performed in triplicate using the CFX96 Touch-Time System (Bio-Rad, USA). The final products of qRT-PCR were electrophoresed on 2% agarose gels and visualized with SafeView Nucleic Acid Stain (G108; Applied Biological Materials, Canada). The cycle numbers (Ct) of the critical point at which the fluorescent signal exceeded the background were determined by qRT-PCR, and the expression values for each gene were normalized to the expression of GAPDH, the endogenous control within each sample. Relative quantification was used to calculate fold change using the comparative Ct method (ΔΔCt).


Statistical analysis was conducted using R ver.3.6.3 and IBM SPSS Statistics for Windows software (ver. 24; IBM, USA). Microarray data were analyzed using the Student’s t-test in R 3.6.3. The Student’s t-test for qRT-PCR results and linear regression for the comparison of qRT-PCR and microarray results were done using IBM SPSS Statistics software (ver. 24). Statistical significance was set at P < 0.05, and the values are expressed as mean ± SEM.

Data availability

The datasets generated for this study can be found in the NCBI GEO database (GSE183386, and the single-cell RNA sequencing data datasets analyzed for this study is GSE137478.


Generating the animal depression model and differentially expressed gene screening

To discover habenular DEGs in response to aversive stimuli, rats were subjected to restraint stress for 14 consecutive days, and the habenula tissue was harvested. A comparison between the habenular expression profiles of the CRS and control (non-stressed) groups revealed that 379 genes were differentially expressed in the CRS group (Supplementary Table S2). Among these genes, 210 were downregulated and 169 genes were upregulated. We determined the significant DEGs by the |fold change| > 1.2 with an unpaired t-test with P < 0.05. Volcano plots and heatmaps were generated based on the obtained data (Fig. 1).

Figure F1
Genome-wide gene expression analysis showed 379 DEGs, selected by P < 0.05 and |fold change| > 1.2 criteria. (A) A volcano plot revealed upregulated (red dots) and downregulated (blue dots) ...

Pathway enrichment analysis

To investigate the biological implications of DEGs, we performed KEGG and GO pathway enrichment analyses. The KEGG pathway analysis identified 10 significant pathways (adjusted P < 0.05), including neuroactive ligand-receptor interaction, cyclic adenosine monophosphate (cAMP) signaling pathway, circadian entrainment, cholinergic and serotonergic synapse, and oxytocin signaling pathway (Fig. 2A). Of note, the neuroactive ligand-receptor interaction pathway included DEGs associated with tachykinin-signaling related genes, such as tachykinin receptor 1 (Tacr1), tachykinin 1 and 2 (Tac1, 2), and neuropeptide signaling-related genes, such as neuropeptide Y receptor 1 (Npy1r), somatostatin receptor 2, 4 (Sstr2, 4), oxytocin (Oxt), vasoactive intestinal peptide (Vip), and endothelin-1 (Edn1) (Fig. 3A). Using the GO analysis, 119 GO terms, including 104 biological processes (BPs), eight molecular functions (MFs), and seven cellular components (CCs), were obtained. The most significantly enriched biological pathways in BPs were responses to corticosteroids, positive regulation of lipid transport, anterograde trans-synaptic signaling, trans-synaptic signaling, and chemical synaptic signaling (Figs. 2B and 3B). Channel activity and transmembrane transporter activities in MFs (Fig. 2C), components of the plasma membrane, neuron projection, cell body, ion channel complex, and synapse in CCs (Fig. 2D) were included in the significantly enriched GO terms. Both the KEGG and GO analyses indicated synapse and synaptic signaling-related pathways, including cholinergic genes (Chat, Slc5a7, Chrna3, Chrnb4); serotonergic receptors (Htr3a, Htr5b, Htr4); membrane-trafficking proteins such as synaptotagmin 6, 10 (Syt6, 10); and Cadps2. Pathway enrichment analysis results, including adjusted P values and related DEGs, are listed in Supplementary Table S3.

Figure F2
A total of 379 DEGs were used to analyze the functionally enriched pathways. Only cases with an adjusted P < 0.05 were selected. (A) KEGG pathway and (B-D) GO analysis ...
Figure F3
The network indicates the genes affiliated with the top five significantly enriched pathways in order of their P value in (A) KEGG and (B) GO enrichment analyses. The size of ...

Protein–protein interaction network analysis

Next, we analyzed PPIs to decipher the interrelationship between DEGs and related biological functions underlying the habenular response to CRS. Using the selected DEGs, the PPI network consisted of 359 nodes and 543 edges that were constructed based on the STRING database (Fig. 4A). A total of 21 nodes had more than 10 degrees. Subcluster analysis was performed to compare the functional clusters of interacting proteins based on pathway enrichment analysis. The two most significant subclusters were related to neuroactive ligand-receptor interactions (scores of 10 and 6.571). Other subclusters were related to circadian oscillation (score 5), including period circadian regulator 1-3 (Per 1-3), albumin D-box binding protein (Ddp), nuclear receptor subfamily 1 group D (Nr1d1), and cholinergic synapse (score 3) (Fig. 4B).

Figure F4
(A) Enriched pathway and clustered modules are circled in the PPI network. (B) Results of the clustering analysis related to the enriched pathways. Neuroactive ligand-receptor, circadian oscillatory gene, and cholinergic ...

Categorization of DEGs into cell types and regional expression

To gain insight into cell type-specific and regional expression of the DEGs in the habenula, we used the publicly available single-cell RNA sequencing datasets of mouse habenula (GSE137478) (Hashikawa et al., 2020). Using this published dataset, we identified 26 cell clusters and assigned cell types using canonical markers and known habenular regional markers (Hashikawa et al., 2020; Saunders et al., 2018; Wallace et al., 2020; Zeisel et al., 2018) (Supplementary Fig. S1). On plotting DEGs related to enriched pathways from the KEGG and GO pathway analyses among different cell types and region-specific neuronal clusters, the genes constituting neuroactive ligand-receptor and trans-synaptic signaling-related genes were found to be mainly expressed in the MHb neurons (Fig. 5). For example, genes related to tachykinin signaling (Tac2 and Tacr1) and several neuronal acetylcholine receptor subunits (Chrna3, Chrnb3, and Chrnb4) among the CRS DEGs were specifically expressed in MHb neurons. Moreover, Slc5a7, encoding high-affinity choline transporter 1, was also highly expressed in MHb neurons. Therefore, these results strongly suggest that the MHb neurons are the main neuronal target of chronic stress in the habenula, which may mediate behavioral alterations. Graph indicating cell type of all DEGs related to neuroactive ligand-receptor and trans-synaptic signaling pathway (Supplementary Fig. S2) also shows the expression of DEGs in MHb. The full list of expression levels of DEGs in each cell cluster is shown in Supplementary Table S4.

Figure F5
Violin plot indicating the expression levels and number of cells expressing genes in habenular cell-type clusters from GSE 137478 that are related to (A) neuroactive ligand-receptor pathway and (B) trans-synaptic ...

qRT-PCR validation

We validated the expression profiles of key DEGs that constitute significant pathways using qRT-PCR. Microarray and qRT-PCR results are shown in Fig. 6A, and the correlation analysis indicated a significant correlation between the two datasets (R = 0.937, P < 0.001; Supplementary Fig. S3). The mRNA expression profiles of Tac1, Tac2, and Tacr1 from the neuroactive ligand-receptor interaction pathway; Per1 and Per2 from circadian oscillation; Syt6 for synaptic transmission; and Chrnb4, Slc5a7, and Chat for cholinergic synapses are shown in Fig. 6B-6E. The qRT-PCR results showed reduced expression of Tac1 and Tac2 mRNA, mediating tachykinin neuropeptide signals (Fig. 6B). In contrast, the expression of Per1 and Per2, which are related to the circadian rhythm, was increased in stressed rats (Fig. 6C). The expression of Syt6 and cholinergic genes (Chrnb4,Slc5a7, and Chat) was reduced by CRS (Figs. 6D and 6E). Although the expression of Tacr1 was not significantly reduced, the expression showed a tendency to decrease after CRS. Further details of the qRT-PCR are provided in Supplementary Table S5.

Figure F6
(A) Changes in the expression levels from microarray results validated by qRT-PCR for DEGs. (B) Genes related to neuroactive ligand-receptor pathway. (C) Genes related to circadian entrainment. (D and E) ...


This study revealed that chronic stress mediates significant alterations in the gene expression of the rat habenula, allowing the investigation of distinct groups of genes in stress-mediated depression research. Pathway enrichment analyses indicated that several signaling pathways related to neuronal function were highly enriched in genes related to mood-related disorders. Altered gene expression in the habenula of chronic stress-mediated depression can be categorized into neuroactive ligand-receptor interaction, cAMP signaling pathway, circadian entrainment, and synaptic signaling. Using recently published habenular single-cell RNA-seq data (Hashikawa et al., 2020), we also identified habenular regions and putative cell types in which CRS-evoked DEGs were abundantly expressed.

There is convincing evidence that neuroactive ligand-receptor interactions and neuropeptides are crucial factors involved in mood-related disorders. Indeed, the ablation of the dorsal part of the MHb, which expresses Tac1, results in reduced voluntary wheel-running activity and shorter immobility time in the tail suspension test in mice, suggesting that tachykinin signaling in the MHb is associated with the regulation of mood-related behaviors (Hsu et al., 2014; 2016). Substance P levels are also significantly elevated in the cerebrospinal fluid (CSF) of depressed patients, and selective TACR1 antagonists have shown antidepressant-like effects in several rodent studies (Dableh et al., 2005; Geracioti et al., 2006; Kramer et al., 1998). In addition, the injection of TACR3 antagonist into the brain regions, which blocks tachykinin peptide produced from Tac2, has shown antidepressant-like effects (Zelikowsky et al., 2018). A recent study reported that the expression of Cadps2, which plays a critical role in the neuropeptide secretion of MHb, was decreased in animal models of depression, and downregulation of Cadps2 in mice induced despair-like behaviors (Yoo et al., 2021). In contrast to previous reports that tachykinin antagonism has an antidepressant-like effect, the expression of Tac1 and Tac2 genes was reduced in CRS rats (Fig. 5), which may be a compensatory response to the upregulated tachykinin signaling, and the exact mechanism underlying attenuated habenular tachykinin signaling-induced depressive phenotypes needs further investigation.

Several other neuropeptides among the selected DEGs have also been closely related to neuropsychiatric symptoms. For example, the activation of neuropeptide Y receptor 1 (Npy1r) prevents anxiety and depressive-like symptoms in a rodent model of post-traumatic stress disorder (Nwokafor et al., 2020; Redrobe et al., 2002). It should be noted that the hyperactivation of the LHb neurons induces depression-like phenotypes through the activation of GABAergic neurons in the rostromedial tegmental nucleus, resulting in the suppression of dopaminergic neurons in the midbrain (Yang et al., 2018). A recent study reported that NPY administration into the LHb reduced the spontaneous firing rate of LHb neurons, leading to hypoactivation of the LHb. These results suggest that NPY signaling partially mediates the protective roles in mood-related disorders (Cheon et al., 2019). Thus, the upregulation of habenular Npy1r in CRS rats might be involved in suppressing the hyperactivation of LHb neurons (Fig. 3A, Supplementary Table S2). Consistent with our results, significantly reduced somatostatin expression levels have been reported in patients with both neuropsychiatric (e.g., MDD, bipolar disorder, schizophrenia) and neurodegenerative (e.g., Alzheimer’s and Parkinson’s disease) disorders (Lin and Sibille, 2013). In addition, the analysis of CSF and plasma in depressed patients has shown a negative correlation between oxytocin levels and depression symptom severity (Cochran et al., 2013; Ozsoy et al., 2009; Scantamburlo et al., 2007). However, post-mortem studies have shown that oxytocin mRNA and the number of oxytocin-producing neurons both increase in the hypothalamus of patients with MDD, which is presumably due to a compensatory effect for insufficient levels of oxytocin (Meynen et al., 2007; Purba et al., 1996). Similarly, post-mortem diagnosis of MDD patients shows increased oxytocin receptor mRNA levels in the dorsolateral prefrontal cortex (Lee et al., 2018). Neuropeptides and their receptors are widely distributed in the brain, especially in the regions related to emotional regulation, and neuropeptide signaling modulates monoamine systems, particularly during pathological and stressful conditions (Hökfelt et al., 2000). Considering that the habenula is emerging as a crucial brain region for the regulation of emotions, further studies on neuropeptide signaling in the habenula would be interesting.

The circadian disturbance is a prominent symptom of depression (American Psychiatric Association, 2013). Circadian entrainment-related genes show the rhythmicity of gene expression in regulating the sleep-wake cycle, and circadian patterns of these genes are disrupted in MDD (Lamont et al., 2007; Li et al., 2013a). In this study, although gene expression was not assessed based on the circadian rhythm, Per1-3, Dbp, and Nr1d1 expression levels were altered in the habenula of stressed animals. Some studies have revealed a relationship between habenula and sleep, and LHb is known to contribute to circadian timekeeping and sleep fragmentation by the firing of habenular cholinergic cells (Baño-Otálora and Piggins, 2017; Ge et al., 2021). In the LHb of the chronic mild stress rat model of depression, Per2 gene expression is increased and its altered expression is likely associated with the induction of a depression-like state (Christiansen et al., 2016). Furthermore, silencing the Per2 gene in the LHb of normal rats induces a shorter climbing time in the forced swimming test (Li et al., 2021).

Genes related to synaptic signaling were also altered in CRS rats. DEGs were associated with the cholinergic synapses, serotonergic synapses in the KEGG pathway, and trans-synaptic and chemical synaptic signaling in the GO pathway, which indicated that synaptic signaling functions are significant in mood disorders. Cholinergic synapse-related genes were clustered as a module in the PPI network, and their expression was reduced in the stress group in the qRT-PCR results (Fig. 4B). Acetylcholine precursors and cholinergic agonists have been reported to cause depressive symptoms, and antimuscarinic drugs have antidepressant effects on the prefrontal cortex (Dulawa and Janowsky, 2019; Ghosal et al., 2018; Navarria et al., 2015). In contrast to the systemic antagonism of cholinergic signaling leading to antidepressant-like phenotypes, the downregulation of Chat in the MHb evokes anhedonia-like behavior in rats (Han et al., 2017). The LHb is reciprocally interconnected with a set of serotonin (5-HT) neurons in the dorsal and median raphe nuclei by a complex network of parallel topographically organized direct and indirect pathways. DEGs from the habenula of CRS rats included three serotonin receptors: Htr3a, Htr4, and Htr5b (Supplementary Table S2). In particular, Htr5b, which is enriched in the habenula, was downregulated in CRS rats (Fig. 3A).

In PsyGeNET, which is the database of previously reported psychiatric disease-related genes (Gutiérrez-Sacristán et al., 2015), 62 genes of selected DEGs were related to various psychiatric disorders, and 39 genes, including neuropeptide and synaptic signaling-related genes, were related to depressive disorders (Supplementary Table S6). In addition, synaptic signaling-related DEGs include well-known drug targets of neuropsychiatric disorders (Liao et al., 2021; Wishart et al., 2006); for example, Chrnb4, Chrna3, Htr3a, and Htr4 were downregulated and Adra1d and Gabrd were upregulated in CRS rats (Supplementary Table S2, Supplementary Fig. S4). These are target molecules of antidepressants and antianxiety agents such as amitriptyline, imipramine, fluoxetine, bupropion, and benzodiazepines. This notion implies that habenular responses to chronic stress can be associated with alterations in the responsiveness to medications, as well as the onset of depressive disorder.

The present study identified CRS-associated habenular DEGs and proposed significantly enriched functional pathways. Some of the pathways are well known to be involved in mood disorders, and others need further analysis to confirm their disease association. Some CRS-associated molecular targets strongly suggest neuronal signaling, particularly for the MHb regions, which are significantly affected by chronic stress, and imply a possible role for the MHb neurons in mediating the onset of stress-evoked depressive phenotypes. Although there were few studies on gene expressions of LHb, genome-wide gene expression of MHb has not been studied. Since MHb is one of the important structures in psychiatric disorders, we attempted to analyze the gene expression profile of habenula as a whole and informatically sort genes from regional expression data. The habenula is an essential structure linked to neuropsychiatric disease, and the pathways and genes implicated in this study are candidates for a better understanding of its role in mood disorders.

Supplemental Materials

Note: Supplementary information is available on the Molecules and Cells website (

Article information

Mol. Cells.May 31, 2022; 45(5): 306-316.
Published online 2022-04-25. doi:  10.14348/molcells.2022.2257
1Department of Anatomy, College of Medicine, Korea University, Seoul 02841, Korea
2Department of Legal Medicine, College of Medicine, Korea University, Seoul 02841, Korea
3Medical Science Research Center, College of Medicine, Korea University, Seoul 02841, Korea
4Graduate School of Medical Science and Engineering, Korea Advanced Institute of Science and Technology (KAIST), Daejeon 34141, Korea
5Department of Biomedical Sciences, BrainKorea21 Four, College of Medicine, Korea University, Seoul 02841, Korea
*Correspondence: (HK); (HWL)
Received October 21, 2021; Accepted February 7, 2022.
Articles from Mol. Cells are provided here courtesy of Mol. Cells


  • Aizawa, H., Kobayashi, M., Tanaka, S., Fukai, T., Okamoto, H. (2012). Molecular characterization of the subnuclei in rat habenula. J. Comp. Neurol.. 520, 4051-4066.
  • , (2013). . Desk Reference to the Diagnostic Criteria from DSM-5, , ed. (Washington, DC:American Psychiatric Publishing), pp. .
  • Andrus, B.M., Blizinsky, K., Vedell, P.T., Dennis, K., Shukla, P.K., Schaffer, D.J., Radulovic, J., Churchill, G.A., Redei, E.E. (2012). Gene expression patterns in the hippocampus and amygdala of endogenous depression and chronic stress models. Mol. Psychiatry. 17, 49-61.
  • Bader, G.D., Hogue, C.W. (2003). An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics. 4, 2.
  • Baño-Otálora, B., Piggins, H.D. (2017). Contributions of the lateral habenula to circadian timekeeping. Pharmacol. Biochem. Behav.. 162, 46-54.
  • Buynitsky, T., Mostofsky, D.I. (2009). Restraint stress in biobehavioral research: recent developments. Neurosci. Biobehav. Rev.. 33, 1089-1098.
  • Cheon, M., Park, H., Rhim, H., Chung, C. (2019). Actions of neuropeptide Y on synaptic transmission in the lateral habenula. Neuroscience. 410, 183-190.
  • Chiba, S., Numakawa, T., Ninomiya, M., Richards, M.C., Wakabayashi, C., Kunugi, H. (2012). Chronic restraint stress causes anxiety- and depression-like behaviors, downregulates glucocorticoid receptor expression, and attenuates glutamate release induced by brain-derived neurotrophic factor in the prefrontal cortex. Prog. Neuropsychopharmacol. Biol. Psychiatry. 39, 112-119.
  • Chin, C.H., Chen, S.H., Wu, H.H., Ho, C.W., Ko, M.T., Lin, C.Y. (2014). cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst. Biol.. 8, S11.
  • Christensen, T., Jensen, L., Bouzinova, E.V., Wiborg, O. (2013). Molecular profiling of the lateral habenula in a rat model of depression. PLoS One. 8, .
  • Christiansen, S.L., Bouzinova, E.V., Fahrenkrug, J., Wiborg, O. (2016). Altered expression pattern of clock genes in a rat model of depression. Int. J. Neuropsychopharmacol.. 19, .
  • Cochran, D.M., Fallon, D., Hill, M., Frazier, J.A. (2013). The role of oxytocin in psychiatric disorders: a review of biological and therapeutic research findings. Harv. Rev. Psychiatry. 21, 219-247.
  • Cui, Y., Yang, Y., Ni, Z., Dong, Y., Cai, G., Foncelle, A., Ma, S., Sang, K., Tang, S., Li, Y. (2018). Astroglial Kir4.1 in the lateral habenula drives neuronal bursts in depression. Nature. 554, 323-327.
  • Dableh, L.J., Yashpal, K., Rochford, J., Henry, J.L. (2005). Antidepressant-like effects of neurokinin receptor antagonists in the forced swim test in the rat. Eur. J. Pharmacol.. 507, 99-105.
  • DePasquale, E.A.K., Schnell, D.J., Van Camp, P.J., Valiente-Alandí, Í., Blaxall, B.C., Grimes, H.L., Singh, H., Salomonis, N. (2019). DoubletDecon: deconvoluting doublets from single-cell RNA-sequencing data. Cell Rep.. 29, 1718-1727.e8.
  • Dulawa, S.C., Janowsky, D.S. (2019). Cholinergic regulation of mood: from basic and clinical studies to emerging therapeutics. Mol. Psychiatry. 24, 694-709.
  • Fakhoury, M. (2017). The habenula in psychiatric disorders: more than three decades of translational investigation. Neurosci. Biobehav. Rev.. 83, 721-735.
  • Ge, F., Mu, P., Guo, R., Cai, L., Liu, Z., Dong, Y., Huang, Y.H. (2021). Chronic sleep fragmentation enhances habenula cholinergic neural activity. Mol. Psychiatry. 26, 941-954.
  • Geracioti, T.D., Carpenter, L.L., Owens, M.J., Baker, D.G., Ekhator, N.N., Horn, P.S., Strawn, J.R., Sanacora, G., Kinkead, B., Price, L.H. (2006). Elevated cerebrospinal fluid substance p concentrations in posttraumatic stress disorder and major depression. Am. J. Psychiatry. 163, 637-643.
  • Ghosal, S., Bang, E., Yue, W., Hare, B.D., Lepack, A.E., Girgenti, M.J., Duman, R.S. (2018). Activity-dependent brain-derived neurotrophic factor release is required for the rapid antidepressant actions of scopolamine. Biol. Psychiatry. 83, 29-37.
  • Gutiérrez-Sacristán, A., Grosdidier, S., Valverde, O., Torrens, M., Bravo, À., Piñero, J., Sanz, F., Furlong, L.I. (2015). PsyGeNET: a knowledge platform on psychiatric disorders and their genes. Bioinformatics. 31, 3075-3077.
  • Han, S., Yang, S.H., Kim, J.Y., Mo, S., Yang, E., Song, K.M., Ham, B.J., Mechawar, N., Turecki, G., Lee, H.W. (2017). Down-regulation of cholinergic signaling in the habenula induces anhedonia-like behavior. Sci. Rep.. 7, 900.
  • Hao, Y., Hao, S., Andersen-Nissen, E., Mauck, W.M., Zheng, S., Butler, A., Lee, M.J., Wilk, A.J., Darby, C., Zager, M. (2021). Integrated analysis of multimodal single-cell data. Cell. 184, 3573-3587.e29.
  • Hashikawa, Y., Hashikawa, K., Rossi, M.A., Basiri, M.L., Liu, Y., Johnston, N.L., Ahmad, O.R., Stuber, G.D. (2020). Transcriptional and spatial resolution of cell types in the mammalian habenula. Neuron. 106, 743-758.e5.
  • Hikosaka, O. (2010). The habenula: from stress evasion to value-based decision-making. Nat. Rev. Neurosci.. 11, 503-513.
  • Hökfelt, T., Broberger, C., Xu, Z.Q., Sergeyev, V., Ubink, R., Diez, M. (2000). Neuropeptides--an overview. Neuropharmacology. 39, 1337-1356.
  • Hsu, Y.W., Morton, G., Guy, E.G., Wang, S.D., Turner, E.E. (2016). Dorsal medial habenula regulation of mood-related behaviors and primary reinforcement by tachykinin-expressing habenula neurons. Neuro. 3, .
  • Hsu, Y.W., Wang, S.D., Wang, S., Morton, G., Zariwala, H.A., de la Iglesia, H.O., Turner, E.E. (2014). Role of the dorsal medial habenula in the regulation of voluntary activity, motor function, hedonic state, and primary reinforcement. J. Neurosci.. 34, 11366-11384.
  • Kim, S.M., Cho, S.Y., Kim, M.W., Roh, S.R., Shin, H.S., Suh, Y.H., Geum, D., Lee, M.A. (2020). Genome-wide analysis identifies NURR1-controlled network of new synapse formation and cell cycle arrest in human neural stem cells. Mol. Cells. 43, 551-571.
  • Kramer, M.S., Cutler, N., Feighner, J., Shrivastava, R., Carman, J., Sramek, J.J., Reines, S.A., Liu, G., Snavely, D., Wyatt-Knowles, E. (1998). Distinct mechanism for antidepressant activity by blockade of central substance P receptors. Science. 281, 1640-1645.
  • Lamont, E.W., Legault-Coutu, D., Cermakian, N., Boivin, D.B. (2007). The role of circadian clock genes in mental disorders. Dialogues Clin. Neurosci.. 9, 333-342.
  • Lee, H.W., Yang, S.H., Kim, J.Y., Kim, H. (2019). The role of the medial habenula cholinergic system in addiction and emotion-associated behaviors. Front. Psychiatry. 10, 100.
  • Lee, M.R., Sheskier, M.B., Farokhnia, M., Feng, N., Marenco, S., Lipska, B.K., Leggio, L. (2018). Oxytocin receptor mRNA expression in dorsolateral prefrontal cortex in major psychiatric disorders: a human post-mortem study. Psychoneuroendocrinology. 96, 143-147.
  • Lee, S., Woo, J., Kim, Y.S., Im, H.I. (2015). Integrated miRNA-mRNA analysis in the habenula nuclei of mice intravenously self-administering nicotine. Sci. Rep.. 5, 12909.
  • Li, J.Z., Bunney, B.G., Meng, F., Hagenauer, M.H., Walsh, D.M., Vawter, M.P., Evans, S.J., Choudary, P.V., Cartagena, P., Barchas, J.D. (2013a). Circadian patterns of gene expression in the human brain and disruption in major depressive disorder. Proc. Natl. Acad. Sci. U. S. A.. 110, 9950-9955.
  • Li, K., Zhou, T., Liao, L., Yang, Z., Wong, C., Henn, F., Malinow, R., Yates, J.R., Hu, H. (2013b). βCaMKII in lateral habenula mediates core symptoms of depression. Science. 341, 1016-1020.
  • Li, Y., Li, G., Li, J., Cai, X., Sun, Y., Zhang, B., Zhao, H. (2021). Depression-like behavior is associated with lower Per2 mRNA expression in the lateral habenula of rats. Genes Brain Behav.. 20, e12702.
  • Liao, W., Liu, Y., Huang, H., Xie, H., Gong, W., Liu, D., Tian, F., Huang, R., Yi, F., Zhou, J. (2021). Intersectional analysis of chronic mild stress-induced lncRNA-mRNA interaction networks in rat hippocampus reveals potential anti-depression/anxiety drug targets. Neurobiol. Stress. 15, 100347.
  • Lin, L.C., Sibille, E. (2013). Reduced brain somatostatin in mood disorders: a common pathophysiological substrate and drug target?. Front. Pharmacol.. 4, 110.
  • McLaughlin, I., Dani, J.A., De Biasi, M. (2017). The medial habenula and interpeduncular nucleus circuitry is critical in addiction, anxiety, and mood regulation. J. Neurochem.. 142, 130-143.
  • Meynen, G., Unmehopa, U.A., Hofman, M.A., Swaab, D.F., Hoogendijk, W.J. (2007). Hypothalamic oxytocin mRNA expression and melancholic depression. Mol. Psychiatry. 12, 118-119.
  • Morris, J.H., Kuchinsky, A., Ferrin, T.E., Pico, A.R. (2014). enhancedGraphics: a Cytoscape app for enhanced node graphics. F1000Res. 3, 147.
  • Navarria, A., Wohleb, E.S., Voleti, B., Ota, K.T., Dutheil, S., Lepack, A.E., Dwyer, J.M., Fuchikami, M., Becker, A., Drago, F. (2015). Rapid antidepressant actions of scopolamine: role of medial prefrontal cortex and M1-subtype muscarinic acetylcholine receptors. Neurobiol. Dis.. 82, 254-261.
  • Nwokafor, C., Serova, L.I., Nahvi, R.J., McCloskey, J., Sabban, E.L. (2020). Activation of NPY receptor subtype 1 by [D-His(26)]NPY is sufficient to prevent development of anxiety and depressive like effects in the single prolonged stress rodent model of PTSD. Neuropeptides. 80, 102001.
  • Ozsoy, S., Esel, E., Kula, M. (2009). Serum oxytocin levels in patients with depression and the effects of gender and antidepressant treatment. Psychiatry Res.. 169, 249-252.
  • Purba, J.S., Hoogendijk, W.J., Hofman, M.A., Swaab, D.F. (1996). Increased number of vasopressin- and oxytocin-expressing neurons in the paraventricular nucleus of the hypothalamus in depression. Arch. Gen. Psychiatry. 53, 137-143.
  • Quina, L.A., Wang, S., Ng, L., Turner, E.E. (2009). Brn3a and Nurr1 mediate a gene regulatory pathway for habenula development. J. Neurosci.. 29, 14309-14322.
  • Raudvere, U., Kolberg, L., Kuzmin, I., Arak, T., Adler, P., Peterson, H., Vilo, J. (2019). g:Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res.. 47, W191-W198.
  • Redrobe, J.P., Dumont, Y., Fournier, A., Quirion, R. (2002). The neuropeptide Y (NPY) Y1 receptor subtype mediates NPY-induced antidepressant-like activity in the mouse forced swimming test. Neuropsychopharmacology. 26, 615-624.
  • Saunders, A., Macosko, E.Z., Wysoker, A., Goldman, M., Krienen, F.M., de Rivera, H., Bien, E., Baum, M., Bortolin, L., Wang, S. (2018). Molecular diversity and specializations among the cells of the adult mouse brain. Cell. 174, 1015-1030.e16.
  • Scantamburlo, G., Hansenne, M., Fuchs, S., Pitchot, W., Maréchal, P., Pequeux, C., Ansseau, M., Legros, J.J. (2007). Plasma oxytocin levels and anxiety in patients with major depression. Psychoneuroendocrinology. 32, 407-410.
  • Seo, J.S., Zhong, P., Liu, A., Yan, Z., Greengard, P. (2018). Elevation of p11 in lateral habenula mediates depression-like behavior. Mol. Psychiatry. 23, 1113-1119.
  • Szklarczyk, D., Gable, A.L., Lyon, D., Junge, A., Wyder, S., Huerta-Cepas, J., Simonovic, M., Doncheva, N.T., Morris, J.H., Bork, P. (2019). STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res.. 47, D607-D613.
  • Wallace, M.L., Huang, K.W., Hochbaum, D., Hyun, M., Radeljic, G., Sabatini, B.L. (2020). Anatomical and single-cell transcriptional profiling of the murine habenular complex. Elife. 9, .
  • Wishart, D.S., Knox, C., Guo, A.C., Shrivastava, S., Hassanali, M., Stothard, P., Chang, Z., Woolsey, J. (2006). DrugBank: a comprehensive resource for in silico drug discovery and exploration. Nucleic Acids Res.. 34, D668-D672.
  • Xu, P., Wang, K., Lu, C., Dong, L., Chen, Y., Wang, Q., Shi, Z., Yang, Y., Chen, S., Liu, X. (2017). Effects of the chronic restraint stress induced depression on reward-related learning in rats. Behav. Brain Res.. 321, 185-192.
  • Yang, Y., Wang, H., Hu, J., Hu, H. (2018). Lateral habenula in the pathophysiology of depression. Curr. Opin. Neurobiol.. 48, 90-96.
  • Yi, J.H., Jeon, J., Kwon, H., Cho, E., Yun, J., Lee, Y.C., Ryu, J.H., Park, S.J., Cho, J.H., Kim, D.H. (2020). Rubrofusarin attenuates chronic restraint stress-induced depressive symptoms. Int. J. Mol. Sci.. 21, 3454.
  • Yoo, H., Yang, S.H., Kim, J.Y., Yang, E., Park, H.S., Lee, S.J., Rhyu, I.J., Turecki, G., Lee, H.W., Kim, H. (2021). Down-regulation of habenular calcium-dependent secretion activator 2 induces despair-like behavior. Sci. Rep.. 11, 3700.
  • Yu, G., Wang, L.G., Han, Y., He, Q.Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 16, 284-287.
  • Zeisel, A., Hochgerner, H., Lönnerberg, P., Johnsson, A., Memic, F., van der Zwan, J., Häring, M., Braun, E., Borm, L.E., La Manno, G. (2018). Molecular architecture of the mouse nervous system. Cell. 174, 999-1014.e22.
  • Zelikowsky, M., Hui, M., Karigo, T., Choe, A., Yang, B., Blanco, M.R., Beadle, K., Gradinaru, V., Deverman, B.E., Anderson, D.J. (2018). The neuropeptide Tac2 controls a distributed brain state induced by chronic social isolation stress. Cell. 173, 1265-1279.e19.

Figure 1

Genome-wide gene expression analysis showed 379 DEGs, selected by P < 0.05 and |fold change| > 1.2 criteria. (A) A volcano plot revealed upregulated (red dots) and downregulated (blue dots) DEGs in the microarray analysis. The top 10 hits of the significantly upregulated and downregulated genes are labeled. (B) Heatmap of all DEGs among the control and CRS groups (control, n = 4; stress, n = 4).

Figure 2

A total of 379 DEGs were used to analyze the functionally enriched pathways. Only cases with an adjusted P < 0.05 were selected. (A) KEGG pathway and (B-D) GO analysis results. The bubble size indicates the number of genes related to the pathway, and the color indicates the significance of enrichment. The rich factor denotes the ratio of DEGs to the total number of genes annotated in the pathway.

Figure 3

The network indicates the genes affiliated with the top five significantly enriched pathways in order of their P value in (A) KEGG and (B) GO enrichment analyses. The size of the pathway dot indicates the number of related DEGs, and the color of each gene dot denotes the relative expression level of the genes in the CRS group compared to the control group.

Figure 4

(A) Enriched pathway and clustered modules are circled in the PPI network. (B) Results of the clustering analysis related to the enriched pathways. Neuroactive ligand-receptor, circadian oscillatory gene, and cholinergic synapse-related modules are identified. The thickness of edges indicates STRING DB score.

Figure 5

Violin plot indicating the expression levels and number of cells expressing genes in habenular cell-type clusters from GSE 137478 that are related to (A) neuroactive ligand-receptor pathway and (B) trans-synaptic signaling. The colors of the graph indicate log2 fold change (FC) of each DEGs.

Figure 6

(A) Changes in the expression levels from microarray results validated by qRT-PCR for DEGs. (B) Genes related to neuroactive ligand-receptor pathway. (C) Genes related to circadian entrainment. (D and E) Genes related to synaptic transmission and cholinergic synapses (*P < 0.05, **P < 0.01, Student’s t-test).