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.

Animals

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

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.

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.

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 (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 (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 C_t method (ΔΔC_t).

Statistics

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, 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.

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).

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.

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).

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.

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.

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.

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.