Mol. Cells 2022; 45(5): 306-316
Published online April 25, 2022
https://doi.org/10.14348/molcells.2022.2257
© The Korean Society for Molecular and Cellular Biology
Correspondence to : kimhyun@korea.ac.kr (HK); biocais@korea.ac.kr (HWL)
This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License. To view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/3.0/.
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 (
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
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
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).
The habenula was isolated from the brain immediately after decapitation and stored at –80°C in RNA
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).
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.
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
Habenular single-cell RNA sequencing data (Gene Expression Omnibus accession No. GSE137478) were obtained from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/). 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 (
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
The datasets generated for this study can be found in the NCBI GEO database (GSE183386, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE183386) and the single-cell RNA sequencing data datasets analyzed for this study is GSE137478.
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
To investigate the biological implications of DEGs, we performed KEGG and GO pathway enrichment analyses. The KEGG pathway analysis identified 10 significant pathways (adjusted
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 (
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 (
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,
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
Several other neuropeptides among the selected DEGs have also been closely related to neuropsychiatric symptoms. For example, the activation of neuropeptide Y receptor 1 (
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,
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
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,
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.
This research was supported by the Ministry of Science and ICT through the National Research Foundation of Korea (NRF-2017M3C7A1079692 to H.K., NRF-2017R1D1A1B06032730 to H.W.L., and NRF-2019M3C7A1032764 to G.H.S.).
H.Y. and H.J.K. designed experiments and analyzed data. H.Y. and S.H.Y. conducted experiments and analysis. H.Y. and H.W.L. wrote the first draft of the manuscript. H.J.K., G.H.S., J.A.G., and H.K. reviewed and edited the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.
The authors have no potential conflicts of interest to disclose.
Mol. Cells 2022; 45(5): 306-316
Published online May 31, 2022 https://doi.org/10.14348/molcells.2022.2257
Copyright © The Korean Society for Molecular and Cellular Biology.
Hyeijung Yoo1,5,6 , Hyun Jung Kim1,4,6
, Soo Hyun Yang1
, Gi Hoon Son2
, Jeong-An Gim3
, Hyun Woo Lee1,5,*
, and Hyun Kim1,5,*
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, 6These authors contributed equally to this work.
Correspondence to:kimhyun@korea.ac.kr (HK); biocais@korea.ac.kr (HWL)
This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License. To view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/3.0/.
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 (
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
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
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).
The habenula was isolated from the brain immediately after decapitation and stored at –80°C in RNA
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).
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.
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
Habenular single-cell RNA sequencing data (Gene Expression Omnibus accession No. GSE137478) were obtained from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/). 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 (
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
The datasets generated for this study can be found in the NCBI GEO database (GSE183386, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE183386) and the single-cell RNA sequencing data datasets analyzed for this study is GSE137478.
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
To investigate the biological implications of DEGs, we performed KEGG and GO pathway enrichment analyses. The KEGG pathway analysis identified 10 significant pathways (adjusted
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 (
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 (
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,
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
Several other neuropeptides among the selected DEGs have also been closely related to neuropsychiatric symptoms. For example, the activation of neuropeptide Y receptor 1 (
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,
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
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,
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.
This research was supported by the Ministry of Science and ICT through the National Research Foundation of Korea (NRF-2017M3C7A1079692 to H.K., NRF-2017R1D1A1B06032730 to H.W.L., and NRF-2019M3C7A1032764 to G.H.S.).
H.Y. and H.J.K. designed experiments and analyzed data. H.Y. and S.H.Y. conducted experiments and analysis. H.Y. and H.W.L. wrote the first draft of the manuscript. H.J.K., G.H.S., J.A.G., and H.K. reviewed and edited the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.
The authors have no potential conflicts of interest to disclose.
Soo Min Kim, Soo Young Cho, Min Woong Kim, Seung Ryul Roh, Hee Sun Shin, Young Ho Suh, Dongho Geum, and Myung Ae Lee
Mol. Cells 2020; 43(6): 551-571 https://doi.org/10.14348/molcells.2020.0071