Single-cell molecular profiling of all three components of the HPA axis reveals adrenal ABCB1 as a regulator of stress adaptation – Science Advances
♫ Monday, February 1st, 2021Abstract
Chronic activation and dysregulation of the neuroendocrine stress response have severe physiological and psychological consequences, including the development of metabolic and stress-related psychiatric disorders. We provide the first unbiased, cell typespecific, molecular characterization of all three components of the hypothalamic-pituitary-adrenal axis, under baseline and chronic stress conditions. Among others, we identified a previously unreported subpopulation of Abcb1b+ cells involved in stress adaptation in the adrenal gland. We validated our findings in a mouse stress model, adrenal tissues from patients with Cushings syndrome, adrenocortical cell lines, and peripheral cortisol and genotyping data from depressed patients. This extensive dataset provides a valuable resource for researchers and clinicians interested in the organisms nervous and endocrine responses to stress and the interplay between these tissues. Our findings raise the possibility that modulating ABCB1 function may be important in the development of treatment strategies for patients suffering from metabolic and stress-related psychiatric disorders.
The hypothalamic-pituitary-adrenal (HPA) axis is pivotal for the maintenance of homeostasis in the presence of real or perceived challenges (1, 2). This process requires numerous adaptive responses involving those of the neuroendocrine and central nervous systems (3). When a situation is perceived as stressful, the paraventricular nucleus (PVN) of the hypothalamus releases corticotropin-releasing factor (CRF) to the hypophyseal portal system, connecting the hypothalamus with the anterior pituitary gland, where it stimulates the secretion of adrenocorticotropic hormone (ACTH) into the peripheral bloodstream. In turn, upon binding to the melanocortin 2 receptor, ACTH stimulates the production and secretion of glucocorticoids (GCs) from the adrenal cortex that bind to corticosteroid receptors (4). These act as transcriptional regulators providing the necessary energy resources and behavioral (emotional and cognitive) adaptations to cope with the stressful challenge and also to exert the main negative feedback at different levels of the HPA axis. While necessary for immediate response, prolonged GC exposure can increase morbidity and mortality (5, 6). Dysregulation of the neuroendocrine stress response can have severe psychological and physiological consequences, and chronic activation of the HPA axis has been linked to stress-related disorders such as anxiety disorders, major depression, posttraumatic stress disorder, and metabolic syndrome (7). As exemplified in Cushings syndrome, endogenous overproduction of GCs has detrimental effects such as impaired glucose metabolism; infectious, musculoskeletal, and cardiovascular complications; and neuropsychiatric comorbidities (8). However, despite decades of research, the molecular underpinnings of HPA dysfunction after prolonged exposure to stress are still not fully understood. Furthermore, most of the work in the field has largely focused on investigating chronic stress effects in the brain, yet much less is known about how chronic stress exposure affects the peripheral components of the HPA axis at the molecular level (9). Recent advances in the field of genomics now allow us to obtain genome-wide data on an individual cell level. Single-cell transcriptomics thereby provide powerful insight into the complexity of different tissues by enabling the identification and characterization of molecular signatures at extraordinary resolution, which can ultimately reveal previously unidentified dimensions of cell identities and their relationships with disease (10).
In this study, using single-cell RNA sequencing (scRNA-seq), we comprehensively cataloged transcriptional changes associated with chronic stress exposure in all three components of the HPA axis. We analyzed 21,723 single cells from the PVN, the pituitary, and the adrenal gland from 10 mice across two conditions (controls, n = 5; stress, n = 5). We found cell typespecific transcriptional signatures of chronic stress adaptation across the HPA axis. We identified a novel subpopulation of stress-responsive adrenocortical cells, which play an important role in the plasticity and adaptation process associated with chronic stress exposure in the adrenal cortex. We validated our findings using mouse tissues, human adrenal samples from patients with ACTH-dependent Cushings syndrome, in vitro adrenal cell models, and peripheral cortisol and genotyping data from treatment-nave, depressed patients. Our study provides the first unbiased and systematic characterization of cell typespecific signatures of the HPA axis under baseline (unstressed) and chronic stress conditions. Furthermore, our results allow a deeper understanding of HPA axis activity and its association with stress-related and metabolic disorders. Ultimately, these findings could lead to more accurate, and more reliable, molecular signatures to monitor disease progression and efficacy of treatment.
To induce chronic activation of the HPA axis, we used the chronic social defeat stress (CSDS) model, a validated, commonly used paradigm to induce long-lasting, depression- and anxiety-like endophenotypes in mice (Fig. 1A) (11). Stress exposure resulted in hallmark features of chronically stressed mice, including reduced social interaction, as demonstrated by the social avoidance test (SAT), a significant increase in basal corticosterone (CORT) levels, enhanced adrenal weight, and reduced fur quality, which is a measure associated with decreased grooming behavior (Fig. 1, B to F) (12). Body weight was not significantly different across groups after CSDS (Fig. 1G). Notably, the natural variability shown by control (unstressed) mice in the SAT did not correlate or was indicative of any of the hallmark features of chronically stressed mice (fig. S1). Five mice from each group (controls versus stressed) were selected for molecular characterization. The PVN, pituitary, and adrenal gland from these mice were used for scRNA-seq experiments (Fig. 1H).
(A) Experimental timeline of chronic social defeat stress (CSDS) paradigm for control (n = 15) and stressed (n = 15) mice. (B) CSDS reduced interaction ratios in stressed mice during the social avoidance test. Bigger dots represent the five mice from each group selected for molecular characterization (P = 0.0084, unpaired t test, two-tailed). (C and D) Twenty-one days of social defeat exposure significantly increased (a.m.) basal corticosterone (CORT) levels (P < 0.0001, unpaired t test, two-tailed) and enhanced adrenal weight (P < 0.0001, unpaired t test, two-tailed). (E) Representative adrenal glands from control and stressed mice. Scale bars, 0.5 mm. (F) CSDS significantly reduced fur quality in stressed mice [two-way analysis of variance (ANOVA), P < 0.0001]. Coat state score: (0) no wounds, well-groomed and bright coat, and clean eyes; (1) no wounds, less groomed and shiny coat OR unclean eyes; (2) small wounds, AND/OR dull and dirty coat and not clear eyes; (3) extensive wounds, OR broad piloerection, alopecia, or crusted eyes. (G) Body weight was not significantly affected by chronic stress (two-way ANOVA, P > 0.05). (H) Experimental design for scRNA-seq experiment. Individual cell suspensions were prepared from the PVN, pituitary gland (PG), and adrenal gland (AG) from selected control (n = 5) and stressed mice (n = 5). **P < 0.01, ****P < 0.0001.
To characterize inter-and intratissue heterogeneity of the HPA axis, we sequenced the transcriptome of 21,723 single cells from the PVN, pituitary, and adrenal, obtained from both unstressed (n = 5) and chronically stressed (n = 5) mice. We systematically cataloged cell identities using Scanpy, a scalable toolkit for analyzing single-cell gene expression data (13) following best practices. Graph-based clustering was performed to group cells according to their unique gene expression profiles, and dimension reduction (UMAP, Uniform Manifold Approximation and Projection) plots were used for visualization (Fig. 2) (14). In the PVN, unsupervised cluster analysis revealed a total of 18 cell clusters with distinct gene expression signatures (Fig. 2A). We determined the identity of each cluster based on the expression of established cell typespecific markers from the literature (1520). Expression of these markers across all PVN clusters can be found in fig. S2 (A to D). The 18 clusters from the PVN were further subdivided into eight major cell types as neurons, oligodendrocytes, astrocytes, microglia, endothelial, ependymal, tanycytes, and vascular cells (Fig. 2B). In the pituitary, we identified 22 unique cell clusters across 12 populations, which were grouped into somatotropes, lactotropes, corticotropes, melanotropes, gonadotropes, thyrotropes, stem cells, Pou1f1-expressing mixed cells, macrophages, endothelial cells, vascular cells, and posterior pituitary cells (Fig. 2, C and D, and fig. S3, A and B). Last, in the adrenal gland, we identified 16 unique clusters grouped into eight major groups of cells from the zona glomerulosa, zona fasciculata, a transition zone of cortical cells, medullar cells, capsular and vascular cells, macrophages, endothelial cells, and a small cluster of unknown cells (Fig. 2, E and F, and fig. S4, A and B). Expression of the top 100 genes defining the individual clusters in each of the three tissues can be found in tables S1 to S3.
(A) Dimensionality reduction Uniform Manifold Approximation and Projection (UMAP) plot depicting 6966 single cells from the PVN of the hypothalamus. Colors represent each of the 18 Louvain groups of individual cell types labeled with an abbreviation as follows: glutamatergic neurons (nGLUT1 and nGLUT2), GABAergic neurons (nGABA1 and nGABA2), mixed neurons (nMixed), vasopressin neurons (nAVP), neuropeptides (nNeuP), oligodendrocytes (Oligo1 and Oligo2), committed oligodendrocyte progenitor cells (COPs), oligodendrocyte progenitor cells (OPCs), astrocytes, endothelial, microglia, macrophages, ependymal, tanycytes, and vascular cells. (B) Distribution and percentage of eight major cell types in the PVN (purple). (C) UMAP plot depicting 9879 single cells from the pituitary. Colors represent each of the 22 Louvain groups representing individual cell types labeled with an abbreviation as follows: somatotropes (Somato1 to Somato8), lactotropes (Lacto1 and Lacto2), corticotropes (Cortico1 and Cortico2), melanotropes, gonadotropes (Gonado1 and Gonado2), thyrotropes (Thyro), endothelial, macrophages, vascular cells, stem cells, Pou1f1 mixed cells (Pou1g1 MCs), and posterior pituitary cells (PPCs). (D) Distribution and percentage of 12 major cell types in the pituitary (green). (E) UMAP plot depicting 4878 single cells from the adrenal. Colors represent each of the 16 Louvain groups representing individual cell types labeled with an abbreviation as follows: zona fasciculata (zFasc1 to zFasc5), zona glomerulosa (zGlom1 and zGlom2), transition zone of mixed fasciculata and glomerulosa cells (tZone1 and tZone2), cycling adrenocortical cells (cACCs), macrophages 1 and 2, endothelial, medullar cells, capsular and vascular cells, and unknown cells. (F) Distribution and percentage of eight major cell types in the adrenal (blue).
Next, we performed inter- and intratissue analyses to characterize cell typespecific molecular signatures of chronic stress in all three tissues of the HPA axis. First, we assessed the distribution of cell numbers for each cluster by comparing the total number of cells from the stressed group to controls (Fig. 3, A to C, and fig. S5, A to C). In the PVN, we observed a significant decrease in the number of cells from the Glut2 (32%) and neuropeptide (25%) neuronal clusters (Fig. 3A). In the pituitary, we found a significant increase in two subclusters of somatotropes (Somato6 and Somato8, 67 and 69%, respectively) (Fig. 3B). Last, in the adrenal gland, we identified the largest and most significant changes in cell distribution between the two groups. Specifically, we observed a significant increase in the number of zona fasciculata 1 cells (82%) and macrophages 2 (70%), as well as a significant decrease in the number of zona glomerulosa 1 cells (40%) (Fig. 3C).
(A to C) Distribution of cell numbers by cluster in each condition (control versus stress). Bars represent the percentage of cells from the control and stressed group per cluster (0 to 100%). All controls (gray), PVN (purple), pituitary (green), and adrenal (blue). Fishers exact test *P < 0.05, **P < 0.01, ***P < 0.001. (D) Sixty-six DEGs in 10 clusters of the PVN. Dark purple represents neurons, purple represents glial cells, and light purple represents vascular cells. (E) Six hundred ninety-two DEGs in 17 clusters of the pituitary. Dark green represents endocrine cells, green represents support cells, and light green represents stem/progenitor cells. (F) Nine hundred twenty-two DEGs in 10 clusters of the adrenal gland. Dark blue represents endocrine cells and light blue represents support cells. Size of the circle represents the number of DEGs in each cluster for all three tissues. (G) DEGs across tissues (intertissue analysis). Sixteen DEGs in common (PVN, pituitary, and adrenal), 3 DEGs (PVN and pituitary), 6 DEGs (PVN and adrenal), and 97 DEGs (pituitary and adrenal). Fourteen DEGs exclusively in the PVN (purple), 162 DEGs exclusively in the pituitary (green), and 343 DEGs exclusively in the adrenal gland (blue). Size of the circle represents the total number of DEGs in each cluster for all three tissues. (H) Expression patterns of dysregulation across DEGs per tissue. Heatmaps represent the percentage of up- and down-regulated DEGs per cluster within the PVN (purple), pituitary (green), and the adrenal (blue). Heatmap scale, 0% (gray); 50% (white); 100% (dark purple/green/blue).
Subsequently, we performed differential expression analyses to evaluate cell typespecific molecular signatures of chronic stress. We compared differentially expressed genes (DEGs) within tissues (intratissue analysis) and found that no single gene was differentially expressed (DE) across all cell types for any of the three tissues (tables S4 to S9), suggesting that cell typespecific effects of chronic stress could be masked or diluted in alternative studies using bulk RNA-seq. In contrast, when gene expression was analyzed within cell types, interesting effects emerged. In the PVN, we identified a total of 66 DEGs in 10 of the 18 cell types (Fig. 3D). In the pituitary, our analysis revealed a total of 692 DEGs in 17 of the 22 pituitary clusters (Fig. 3E). Consistent with cell distribution changes by condition, in the adrenal gland, we also observed the largest number of DEGs. Specifically, we identified 922 DEGs in 10 of the 16 adrenal clusters, ranging from 21 to 171 DE transcripts per cell type (Fig. 3F). A full list of DEGs per cell type across all three tissues can be found in tables S4 to S6.
We further compared DEGs across tissues (intertissue analysis). First, we collapsed all DEGs per tissue and identified 39 unique DEGs in the PVN, 278 in the pituitary, and 462 in the adrenal. We then looked for common genes and found 16 DEGs across all tissues (Fig. 3G). There were also 6 DEGs in common between the PVN and the adrenal, 3 DEGs between the PVN and the pituitary, and 97 DEGs between the pituitary and the adrenal glands. In addition, there were 14 genes exclusively DE in the PVN, 162 in the pituitary, and 343 in the adrenal gland (Fig. 3G and table S10). Among the genes dysregulated across the three tissues, we found several genes coding for protein members of the GC receptor (GR) chaperone complex known to play key roles in the stress response (21). Among these were HSP90 (Hsp90aa1 and Hsp90ab1), which is responsible for the direct binding of GR to the chaperone complex; HSP70 (Hspa1a and Hspa8), which encodes the first protein that recognizes and binds newly synthesized GR; and HSP40 (Dnaja1 and Dnajb1), which mediates the interaction between GR and its chaperones (22). We also observed consistent differences between the PVN and the adrenal gland for the transcription factor Nfkbia (NFB), known to interact with GCs due to their strong anti-inflammatory properties (22, 23), and Fkbp4 (encoding for the FKBP52 protein), a major regulator of GR activity (table S10) (23, 24).
Moreover, we found that most cell populations in the PVN and the pituitary showed an up-regulation of DEGs after exposure to chronic stress, except for microglial cells (PVN), macrophages, and vascular cells (pituitary), where DEGs were down-regulated (Fig. 3H). In the adrenal gland, we noticed a different pattern of regulation with several cell types, including macrophages and adrenocortical cells, showing a down-regulation of DEGs after exposure to chronic stress (Fig. 3H), suggesting a larger range of transcriptional plasticity after chronic stress at the adrenal level. Overall, these results suggest that the most profound differences due to chronic stress in the HPA axis occur in the adrenal gland, where our intra- and intertissue analyses identified the largest number of DEGs and the most significant changes at the cell population level.
The adrenal gland is a highly dynamic organ, which can quickly adapt and regenerate in response to different types of stimuli (25). For example, the adrenal significantly increases its weight in response to chronic stress, a phenomenon that has been documented in rodents, as well as human psychiatric patients (2628). In our study, we confirmed a significant increase of the adrenal weight of mice exposed to chronic social stress (Fig. 1, D and E), and single-cell transcriptomic analyses of the adrenal gland revealed a specific population of overrepresented zona fasciculata cells within the stressed group (zFasc1) (Fig. 3C). In an attempt to further investigate zFasc1 cells and identify what makes them unique, we compared their molecular profiles against all other cells in the adrenal. Because this population was so strongly driven by stress, we reasoned that the genes defining this cell type are also important responders to chronic stress. We found that the top three genes that defined the zFasc1 population were the adenosine 5-triphosphate (ATP)binding cassette subfamily B member 1B (Abcb1b) [qval: 3.27 10146; fold change (FC): 7.4], Suprabasin (Sbsn) (qval: 4.08 1084; FC: 6.6), and the 5-reductase (Srd5a2) (qval: 7.78 1082; FC: 5.8) (Fig. 4A, table S3, and fig. S5D). These genes have been previously associated with GC transport (29, 30), cell proliferation (31), and glucose metabolism (32). To validate our findings and to rule out any potential bias introduced by single-cell dissociation methods that can affect the proportions of cells in the original intact tissue, we performed mRNA in situ hybridization of Abcb1b, Sbsn, and Srd5a2 using RNAscope in adrenal glands obtained from nave or chronically stressed mice. Consistent with our single-cell results, the expression of these genes was restricted to adrenocortical cells from the zona fasciculata (Fig. 4B). Moreover, we observed a significant increase of Abcb1b and Sbsn, but not Srd5a2, mRNA expression in the zona fasciculata of stressed mice as compared to controls (Fig. 4, C to E).
(A) UMAP plot showing the expression of the top three genes that differentiate zFasc1 from other zFasc clusters: Abcb1b, Sbsn, and Srd5a2. Cyp11b1 is expressed in all zona fasciculata cells (zFasc1 to zFasc5). (B) Expression of Abcb1b, Sbsn, and Srd5a2 is restricted to adrenocortical cells from zona fasciculata. Representative adrenal glands from control and stressed mice, showing mRNA expression (brown) of Abcb1b, Sbsn, and Srd5a2 by RNAscope. Nuclei were stained with vector hematoxylin QS (purple). Scale bars, 500 m. (C to E) Chronic stress leads to a significant increase of Abcb1b (P < 0.0002) and Sbsn (P < 0.0005), but not Srd5a2 (P = 0.9715), mRNA expression in the zona fasciculata of stressed (n = 14) as compared to control (n = 14) mice. Representative images show the percentage of mRNA expression (brown) and nuclei (purple). Scale bars, 50 m. (F) CSDS leads to cellular hypertrophy in the adrenal cortex of chronically stressed mice (P < 0.0001). Bar graphs represent the average number of nuclei from the zona fasciculata. Average cell area was calculated by dividing the number of nuclei by the total area. Values are multiplied by 1000 for graphical representation. (G and H) Hypertrophy in the adrenal cortex is associated with higher levels of Abcb1b mRNA expression. Bar graphs represent the average number of nuclei present in areas of high Abcb1b (P < 0.0001) and Sbsn (P = 0.9628) mRNA expression as compared to low expressing regions in zona fasciculata of stressed mice (n = 11). All unpaired t tests, two-tailed. ***P < 0.001, ****P < 0.0001.
Subsequently, we tested whether the increase in adrenal weight after chronic stress exposure was due to an increase in the number of cells (hyperplasia) or an increase in the size of cells (hypertrophy) at the adrenal cortex. Our analysis revealed that, in stressed mice, the number of nuclei present in the zona fasciculata was significantly lower as compared to controls (Fig. 4F), suggesting cellular hypertrophy in the adrenal cortex of chronically stressed mice. Last, we evaluated whether growth characteristics of zona fasciculata cells with high expression levels of Abcb1b or Sbsn were different from low-expressing cells. Unexpectedly, our analysis revealed that the number of nuclei present in areas with high Abcb1b expression was significantly lower than in regions with low Abcb1b expression (Fig. 4G). We did not find any differences in nuclei density between regions of high or low Sbsn expression (Fig. 4H). These results suggest that hypertrophy in the adrenal cortex is associated with higher levels of Abcb1b mRNA expression.
Next, we investigated how the adrenal expression levels of Abcb1b, Sbsn, and Srd5a2 change over time, during 21 days of chronic stress exposure. Therefore, we exposed six groups of mice to a different number of social defeat sessions (0, 3, 5, 10, or 21 days). An additional group of mice received 21 days of social defeat, followed by 48 hours of recovery time to match the end point of our original CSDS paradigm (23 days) (Fig. 5A). We observed a significant and gradual increase in adrenal weight across time points (Fig. 5B). Three days of social defeat were sufficient to stimulate a significant increase in adrenal weight, which continued steadily and plateaued between days 10 and 21. We then quantified bulk mRNA expression levels of Abcb1b, Sbsn, and Srd5a2 in the adrenal cortex from these mice using quantitative real-time polymerase chain reaction (qRT-PCR). We found a significant increase of Abcb1b mRNA levels after 5 days of social defeat, while an increase for Sbsn was present only after 21 days (Fig. 5C). Consistent with our in situ nuclei quantification, we also identified a significant correlation between adrenal weight gain and the expression levels of Abcb1b (r = 0.73; P < 0.0001) and Sbsn (r = 0.51; P = 0.004) (Fig. 5D), suggesting that increases in the expression of these genes were proportional to increases in adrenal weight. In contrast, Srd5a2 did not yield any significant results in these experiments (Fig. 5, C and D). We did not find any significant differences in adrenal weight or mRNA expression levels of these genes between days 21 and 23, suggesting that the long-lasting effects of the CSDS paradigm are still present 48 hours after the last defeat session. Last, our results suggest that chronic stress exposure causes zona fasciculata cells to enlarge and increase their expression of Abcb1b, perhaps as a mechanism to cope with the increased production of GCs in the system.
(A) Experimental timeline. Six different groups of mice (n = 5) were exposed to a different number of social defeat sessions. (i) Controlno defeat, (ii) 3 days, (iii) 5 days, (iv) 10 days, (v) 21 days, and (vi) 21 days, followed by 48 hours of recovery time. (B) Three days of social defeat are sufficient to stimulate a significant increase in adrenal weight, which continued steadily and plateaued between days 10 and 21 (P < 0.0001). (C) Social defeat exposure leads to a significant increase of Abcb1b mRNA levels after 5 days of social defeat (P < 0.0001), an increase of Sbsn after 21 days (P < 0.001), while no significant changes in Srd5a2 expression (P = 0.12). Bar graphs represent mRNA levels of Abcb1b, Sbsn, and Srd5a2 normalized to Hprt. qRT-PCR, quantitative real-time polymerase chain reaction. (D) Social defeat leads to a significant correlation between adrenal weight gain and mRNA levels of Abcb1b (Pearson r = 0.73; P < 0.0001) and Sbsn (r = 0.51; P = 0.004), but no correlation with Srd5a2 expression (r = 0.29; N.S., no significance. P > 0.05). ***P < 0.001, ****P < 0.0001. CTRL, controls; SD, social defeat; D, day.
The Abcb1 gene, also known as multidrug resistance protein 1 (MDR1) or P-glycoprotein 1 (P-gp), is a well-characterized, ATP-dependent efflux pump, whose role is to transport xenobiotics and endogenous cellular metabolites across cellular membranes (33). The protein product of Abcb1 is encoded by two gene variants in mice (Abcb1a and Abcb1b) but only one gene in humans (ABCB1) (34). Moreover, it has been hypothesized that this gene modulates HPA axis activity and mediates antidepressant treatment response by regulating access of GCs and antidepressants into the brain (35). Most of the current literature in biological psychiatry has been primarily focused to understand the activity of Abcb1a in the brain, based on early observations that, in humans, the ABCB1 gene is highly expressed in endothelial cells of the blood-brain barrier (36). However, translational studies in rodents have not been successful in explaining how Abcb1 regulates the response to stress or antidepressant treatment (37). One of the reasons might be that most of these studies were carried out under the assumption that Abcb1a and Abcb1b have similar patterns of expression in the brain. Our single-cell analysis shows a very different picture with limited coexpression among the two variants. Abcb1a is specifically expressed in endothelial cells from the PVN and the pituitary (Fig. 6A), while Abcb1b is expressed in microglia and macrophages of all three tissues, in lactotropes and somatotropes of the pituitary, and in a subsection of zona fasciculata cells (zFasc1), where it shows its highest expression (Fig. 6A). Furthermore, we quantified the expression of Abcb1a and Abcb1b using publicly available bulk RNA-seq data from 35 different mouse tissues (38) and found that their expression also differs considerably in other peripheral organs. Abcb1a is lowly expressed in the periphery, while Abcb1b is the predominant variant showing high expression levels among multiple tissues, particularly in the adrenal gland where the expression of Abcb1b is several magnitudes higher as compared to every other tissue tested (Fig. 6B). These findings suggest that the adrenal gland is an important site for Abcb1 activity.
(A) UMAP plots representing cell typespecific mRNA expression of Abcb1a and Abcb1b in the PVN, pituitary, and adrenal gland of mice. (B) Bulk RNA sequencing data from 35 different mouse tissues showing mRNA expression levels of Abcb1a and Abcb1b. Heatmaps represent expression levels (0 to 12). Red, high expression; white, low expression. Expression values are displayed as Transcripts per Kilobase Million (TPM) and are log2-transformed.
Previous studies in rodents have shown that in vivo inhibition of Abcb1 by intraperitoneal injection of tariquidar, a highly specific and potent Abcb1a/b inhibitor (39), leads to a decrease in CORT levels after acute stress (40). Others have shown that mutant mice lacking both variants (Abcb1a/b) have lower baseline CORT levels as compared to wild-type controls (41). However, these studies could not attribute changes in CORT to a specific Abcb1 variant (Abcb1a or Abcb1b), nor could they conclude that the effects are modulated at the level of the brain or any of the peripheral tissues where Abcb1a and Abcb1b are expressed. To specifically explore the function of Abcb1b in the adrenal gland, we examined whether pharmacological inhibition by tariquidar modulates secretion of CORT in vitro, using an adrenocortical cell line. Mouse Y1 cells were stimulated for 24 hours with 10 nM forskolin alone, or in combination with different concentrations of tariquidar. Forskolin induces secretion of CORT by stimulation of adenylate cyclase (42). While 24 hours of forskolin treatment significantly increased CORT levels as compared to controls, we found a dose-dependent decrease of CORT with increasing concentrations of tariquidar (Fig. 7A), suggesting that GC secretion from adrenocortical cells might be dependent on Abcb1b function. In an attempt to translate our findings to humans, we assessed the modulatory role of ABCB1 on GCs, using human NCI-H295R adrenocortical cells, a validated in vitro model for steroid profiling based on their ability to produce and secrete the major steroidogenic enzymes of the adrenal cortex (43). In line with our previous results, treatment of NCI-H295R cells with 10 nM forskolin led to a significant increase of medium cortisol levels, as compared to vehicle-treated controls (Fig. 7B). Treatment with increasing concentrations of tariquidar led to a significant decrease of media cortisol levels (Fig. 7B). Together, our results show that in vitro pharmacological manipulation of Abcb1 in adrenocortical cell lines leads to a decrease in GC secretion, suggesting that modulation of Abcb1b in adrenocortical cells affects GC secretion in both mice and humans.
(A) Pharmacological inhibition of Abcb1 by tariquidar in mouse Y1 adrenocortical cells. CORT levels (ng/ml) after 24-hour treatment with forskolin (0 and 10 nM) or tariquidar (0, 10, 50, 125, 250, 500, and 1000 nM). (B) Pharmacological inhibition of ABCB1 by tariquidar in human NCI-H295R adrenocortical cells. Cortisol levels (ng/ml) after 24-hour treatment with forskolin (0 and 10 nM) or tariquidar (0, 10, 50, 125, 250, 500, and 1000 nM). One-way ANOVA, ****P < 0.0001.
In humans, chronic endogenous oversecretion of ACTH due to a pituitary or ectopic tumor results in excessive GC secretion and enlargement of the adrenal glands (Fig. 8A) (44). Thereby, this disease stage overlaps with the chronic activation of the HPA axis and hypersecretion of GCs in stress-related disorders. In cases of unsuccessful pituitary surgery or in those patients in whom the ectopic source of ACTH remains obscure, bilateral adrenalectomy is required to treat steroid excess. This opened the possibility to study adrenal glands that had been chronically stimulated and to compare those with controls in the absence of ACTH oversecretion. We quantified mRNA expression of ABCB1 and SBSN, using RNAscope in cases (n = 8) and controls (n = 6). SBSN mRNA was not detectable in human adrenocortical samples (Fig. 8B). Following this approach, we identified a significant up-regulation of ABCB1 mRNA in Cushings adrenocortical samples, as compared to controls (Fig. 8C). These results are consistent with our initial findings in chronically stressed mice and reinforce our evidence for a role of ABCB1 as a modulator of GC activity in the adrenal gland. In addition, our results highlight ABCB1 as a potential regulator of the detrimental effects of impaired glucose metabolism associated with patients with Cushings syndrome.
(A) Graphical representation of the effects of Cushings disease on the adrenal gland. (B and C) Expression of ABCB1 (P = 0.0056) and SBSN (P = 1.0) mRNA, using RNAscope in adrenal glands from patients with ACTH-dependent Cushings syndrome (n = 8) and controls (n = 6). Representative images show the percentage of mRNA expression (brown) and nuclei (purple) normalized by total area. Scale bars, 50 m, **P < 0.01.
Individuals who lack ABCB1, as it occurs in some breeds of dogs with the ABCB1-1 mutation (45), have severe adverse reactions to common medications that act as substrates of this transporter, such as immunosuppressants and steroid hormones (46). Previous studies have shown that dogs and rodents lacking a functional Abcb1/ABCB1 gene have a blunted HPA axis response compared to wild-type animals (41, 47). In humans, multiple single-nucleotide polymorphisms (SNPs) map to the ABCB1 gene locus, and some of these variants have been associated with reduced protein function and activity (48). One of the most studied ABCB1 polymorphisms is the rs2032582 (G2677T), which is a nonsynonymous variant on exon 21 that has been linked to major depressive disorder and treatment response (49). To investigate the relevance of our findings in depressed human patients, we examined whether the ABCB1 polymorphism rs2032582 is associated with an altered HPA axis response, using peripheral plasma samples from 154 treatment-nave, depressed patients. We measured plasma ACTH and cortisol concentrations in depressed patients at baseline, following CRF stimulation, and 15-min intervals for the following hour (Fig. 9A). The genotype and allele distributions of rs2032582 in patients are shown in Fig. 9A. At baseline, patients with the minor allele (TT) showed a decrease in cortisol levels as compared to the major (GG) and heterozygote (TG) alleles; however, this effect did not reach statistical significance after Bonferroni correction (Fig. 9B). After CRF stimulation, we found a significant genotype-by-time interaction in patients cortisol levels (qval: 0.033). More specifically, patients with the minor (TT) allele showed a dampened cortisol response after CRF stimulation (Fig. 9C). We did not find any statistical differences in ACTH levels after CRF stimulation (Fig. 9D), suggesting that the effects of rs2032582 on the ABCB1 gene might be taking place at the level of the adrenal gland. These results are consistent with our mouse and cell culture findings and support the idea that Abcb1/ABCB1 function may regulate HPA axis response.
(A) Experimental design. Predictors of remission in depression to individual and combined treatments (PReDICT) cohort (N = 154) to investigate effects of the ABCB1 variant rs2032582 on HPA axis function. CRF stim, CRF stimulation test; SNP, single-nucleotide polymorphism; HWE-P, Hardy-Weinberg equilibrium P value. (B) Baseline cortisol levels (g/ml) for treatment-nave, depressed patients carrying the rs2032582 SNP genotype. (C and D) Cortisol (g/ml) and ACTH (pg/ml) levels after CRF stimulation (log-transformed). Only completers were included in the analysis. There are no dropouts in sample sizes over time for cortisol or ACTH. GG = major homozygotes (n = 56), GT = heterozygotes (n = 74), TT = minor homozygotes (n = 24). Mixed effects models, Bonferroni-corrected *P < 0.05.
Despite decades of research, the molecular and cellular identity of the HPA axis components, their inter-relationships, and their function after chronic stress exposure are still only partially understood. Here, using scRNA-seq, we describe cell typespecific molecular signatures of chronic stress in all three components of the HPA axis, providing a level of resolution never before reached.
The PVN integrates and coordinates the neuroendocrine HPA axis response to stressful stimuli. However, aside from containing the neuroendocrine neurons that control the synthesis and release of CRF, the PVN also exhibits a significant degree of cellular and molecular complexity, with multiple types of neuronal and nonneuronal subtypes. In this study, we characterized and described the cellular heterogeneity and identity of all cell types in the PVN. We identified many DEGs that are involved in the intracellular trafficking of the GC and mineralocorticoid receptors and play key roles in the response to chronic stress, such as heat shock proteins and Fkbp4 across multiple cell types (24). We also found groups of genes that were DE in unique cell types, such as the cysteine-rich angiogenic inducer 61 gene (Cyr61), which was only found DE in ependymal cells. Cyr61 is a target gene of the hippo signaling pathway, which regulates tissue homeostasis, regeneration, proliferation, and growth and has recently been linked to the pathophysiology of stress-related psychiatric disorders (50). In the neuropeptide cluster, we found a down-regulation of corticotropin-releasing factor (Crf) and vasopressin (Avp), as well as an up-regulation of oxytocin (Oxt) and somatostatin (Sst); however, these changes did not survive correction for multiple testing. We did not find any significant dysregulation of GR (Nr3c1) mRNA in any of the clusters of the PVN. Nevertheless, we did find a significant difference in the total number of Nr3c1+ cells in some of the cell clusters of the PVN (fig. S2D), suggesting that the GR mRNA differences reported in the literature (4, 6, 9) could be due to a decrease in the total number of Nr3c1+ cells after chronic stress, rather than lower expression levels of the existing cells. A decrease in the total number of Nr3c1+ cells is not found across all cell types of the PVN but is rather limited to specific cell populations. These populations could represent the cell types where stress exerts its main effects in the PVN via GR. However, these findings would need to be further validated and replicated in other studies. Last, we found that most cell populations across the PVN showed an up-regulation of DEGs after exposure to chronic stress, except for microglial cells where most DEGs were down-regulated. These changes in microglial cells in combination with gene expression differences (across multiple cell types) of several genes involved in the intracellular trafficking of GCs are possibly the result of overexposure to GCs during a prolonged (chronic) stress paradigm. GCs are released during the stress response and are well known for their immunosuppressive and anti-inflammatory properties. In addition, growing evidence suggests that changes in neuroendocrine function and metabolism are significant triggers of inflammation, which has been linked to the development of neuropsychiatric disorders. Ultimately, while this is an important issue in the field, it is logistically challenging to address considering that the effects or stress, GCs, and inflammation are closely intertwined, likely powering each other in a bidirectional way.
The second component of the HPA axis, the pituitary gland, is a complex organ and an important regulator of major physiological processes, including the neuroendocrine stress response (51). It is composed of a heterogeneous mix of endocrine, general support, and stem cells (17, 18). Despite a significant body of research characterizing attributes of pituitary activity, the cell typespecific regulation of chronic stress at the pituitary level is still poorly understood. Here, we characterized cell typespecific molecular signatures of chronic stress in the pituitary gland. Among many, our DE analysis revealed several genes that were consistently dysregulated in multiple endocrine cells, such as somatotropes, gonadotropes, lactotropes, and corticotropes. Specifically, we found an up-regulation of Cd63, Hsp90aa1, and Hsp90ab1, as well as a down-regulation of several ribosomal genes in all four cell types, suggesting altered GC and ribosomal activity. Moreover, corticotropes are directly stimulated by CRF and are responsible for the release of ACTH into circulation. In our analysis, we found 32 DEGs in this population. However, we did not find any significant differences in the expression of the corticotropin-releasing hormone receptor 1 (Crhr1) or the GR (Nr3c1) gene. Furthermore, and consistent with our findings in the PVN, we found that most cell types across the pituitary showed an up-regulation of DEGs after exposure to chronic stress, except for macrophages and vascular cells, where most DEGs are down-regulated. In our single-cell data, we found a large number of DEGs across multiple types of endocrine cells, suggesting that the stress response in the pituitary gland is a dynamic and complex process that is not limited to the effect that CRF exerts on corticotropes. In our analyses, somatotropes were the population of pituitary cells that showed the biggest changes after chronic stress, both in terms of number of DEGs and changes in proportions of cells. Somatotropes produce and release growth hormone, and they play an important role in the regulation of GC synthesis and adrenal growth and have been shown to positively affect adrenal cell size and number of adrenocortical cells (52). However, the role that somatotropes play in chronic stress and the development of stress-related psychiatric disorders are still poorly understood. Our high-throughput, cell typespecific findings of the effects of chronic stress on the pituitary and somatotropes are both novel and a significant advancement to our understanding of the mechanisms of stress adaptation in the pituitary gland.
Last, the adrenal gland is a major effector of the HPA axis, where interplay between several types of specialized cells takes place to coordinate a complex endocrine, immune, and metabolic response to stress. It is composed of the adrenal medulla and the adrenal cortex, two embryonically different endocrine tissues (25). The adrenal cortex is further divided into three major zones: zona glomerulosa (zG), zona fasciculata (zF), and zona reticularis (zR), each responsible for the synthesis and release of mineralocorticoids, GCs, and androgens, respectively (27). Zona reticularis has been shown to be absent in mice (25). Until now, our understanding of the mechanisms responsible for chronic stress adaptation in the adrenal has been limited. Our study is the first to provide a cell typespecific, unbiased, molecular characterization of the adult adrenal gland (under baseline or chronic stress conditions). Across several cell types, we found a significant dysregulation of genes coding for steroidogenic enzymes responsible for the biosynthesis of corticosteroids, such as GCs and mineralocorticoids. More specifically, we found a dysregulation of Star, Fdx1, Cyp11b1, Cyp21a1, Cyp11a1, Hsd3b1, Nr4a1, and Agtr1a after exposure to chronic stress. In contrast to what we found in the PVN and pituitary, cell types in the adrenal showed both up-regulation and down-regulation of DEGs after exposure to chronic stress, suggesting a larger range of transcriptional plasticity after chronic stress at the adrenal level. Although the changes in the expression of genes coding for steroidogenic enzymes are consistent with the current literature (27), our results offer a new level of resolution by describing the specific cell types where these changes take place in the adrenal. Furthermore, our results highlight that changes after chronic stress in the adrenal are not limited to the endocrine cells of the adrenal cortex or adrenal medulla. In our data, we also find significant changes in the number of macrophages, as well as the number of DEGs in this cluster, after chronic stress. Macrophages are modulated by GCs to secrete cytokines and regulate inflammation and the immune system (53). To the best of our knowledge, this study is the first to show a significant effect of chronic stress in macrophages of the adrenal gland. In addition, our results show a global dysregulation of transcriptional activity in macrophages across all three components of the HPA axis (PVN, pituitary gland, and adrenal gland), suggesting that this cell population is part of a common, multilevel and multitissue signaling network that regulates adaptation to chronic stress.
One of the main findings from our study is the identification of a novel population of overrepresented Abcb1b+ cells within the zona fasciculata of the stressed group. The identification of this novel and specialized cell type in the adrenal gland could not have been possible using standard bulk RNA-seq methods. All previous transcriptomic studies examining the effects of chronic stress in the adrenal gland have been limited to adrenocortical, adreno-medullar, or whole tissue homogenates that average out the signature of thousands of cells, which can mask, dilute, or even distort signals of interest coming from specialized cell populations. Hence, one can expect that any cell typespecific signature of chronic stress (as is the case for zFasc1 cells) has been diluted or even lost in these studies. Here, through a series of complementary experiments, we validated this novel subpopulation of Abcb1b+ cells in the adrenal cortex, which play an important role in stress adaptation. Our experiments showed that increased mRNA expression of Abcb1b+ cells in the adrenal gland is associated with increased adrenal weight and cellular hypertrophy in the adrenal cortex of stressed mice, suggesting that chronic stress exposure causes zona fasciculata cells to enlarge and increase their expression of Abcb1b, perhaps as a mechanism to cope with the increased and sustained production of GCs in the system. The Abcb1 gene is a well-characterized efflux pump whose role is to transport substances, deemed as harmful, across membranes. However, most of the work to study this gene in psychiatry has been primarily focused on understanding the activity of the variant Abcb1a in the brain. Our single-cell analysis in combination with bulk RNA-seq data from 35 different mouse tissues showed that Abcb1b is the predominant variant in the periphery showing high expression levels among multiple tissues, particularly in the adrenal gland, suggesting that the adrenal is an important site for Abcb1 activity. Furthermore, to disentangle the effects of Abcb1a versus Abcb1b in the response to stress, we performed a series of in vitro experiments in mouse and human adrenocortical cells. Our results showed that pharmacological inhibition of Abcb1b in adrenocortical cell lines leads to a decrease in GC secretion, suggesting that modulation of Abcb1b in adrenocortical cells affects GC activity in both mice and humans. Moreover, in an attempt to translate our findings to humans, we investigated the expression of ABCB1 in adrenal cortical tissues from patients diagnosed with ACTH-dependent Cushings syndrome. These patients suffer from excessive GC secretion and adrenal hypertrophy due to a pituitary or ectopic tumor. Thus, this disease stage overlaps with the chronic activation of the HPA axis and hypersecretion of GCs in stress-related disorders. We found a significant up-regulation of ABCB1 in cases, as compared to controls. In addition to being consistent with our findings in chronically stressed mice, these results highlight the role of ABCB1 as a modulator of GC activity in the adrenal gland and postulate ABCB1 as a potential regulator of the impaired glucose metabolism associated with Cushings syndrome. Last, we investigated the relevance of our findings in depressed human patients by examining whether the ABCB1 polymorphism rs2032582 (G2677T) is associated with an altered HPA axis response in peripheral plasma samples from treatment-nave, depressed patients. In humans, the rs2032582 polymorphism has been associated with reduced protein function and activity and has been linked to major depressive disorder, suicidal ideation, and treatment response (49). Consistent with our findings in mice, adrenocortical cell lines, and adrenocortical samples from human Cushings patients, we found that, after CRF stimulation, patients with the minor (TT) allele showed a dampened cortisol but normal ACTH response, suggesting that the effects of rs2032582 on the ABCB1 gene might be taking place at the level of the adrenal gland. In addition, our results support the idea that Abcb1/ABCB1 function may regulate HPA axis response.
Together, our data offer new insights into how chronic stress regulates transcriptional activity in a multilevel, cell typespecific fashion. We identified hundreds of novel genes that are dysregulated across all tissues and levels of the HPA axis. On the basis of our intra- and intertissue analyses, we found the most profound differences due to chronic stress in the adrenal gland, which had the highest number of DEGs and the most significant changes at the cell population level. Through a series of complementary behavioral, molecular, cellular, and functional experiments, we identified a novel subpopulation of Abcb1b+ cells in the adrenal cortex, which play an important role in the adaptation process and plasticity associated with chronic stress exposure. The exact mechanism underlying the effect of ABCB1 on GC regulation and secretion in the adrenal cortex still needs to be further explored. However, previous studies have shown that transcriptional regulation of the Abcb1 genes can be mediated through a putative GC response element (GRE) identified in the promoter region of these genes in both rodents and humans (54, 55). At least in mice, this GRE binding site is only found in the promoter region of the Abcb1b variant, suggesting that Abcb1b is directly regulated by GCs in the periphery, predominantly in the adrenal glands. Therefore, we propose that the Abcb1b/ABCB1 gene and protein are involved in mediating chronic stress adaptation through regulation and control of GCs in the adrenal gland. Our findings raise the prospect that modulating ABCB1 function may be important in the treatment of patients suffering from neuropsychiatric and metabolic disorders, such as stress-related disorders and Cushings syndrome. They further suggest that adrenal ABCB1 activity could be used to stratify patients and tailor treatment strategies. Ultimately, our results provide a deeper understanding of the complex mechanisms of HPA axis regulation.
All experiments were performed in accordance with the European Communities Council Directive 2010/63/EU. All protocols were approved by the Ethics Committee for the Care and Use of Laboratory Animals of the government of Upper Bavaria, Germany. Male mice aged between 7 and 10 weeks old were used for all experiments. Mice were bred in the animal facility of the Max Planck Institute of Biochemistry (Martinsried, Germany) and group-housed (four to five mice per cage) until 1 week before the start of the experiments, when mice were single-housed. Mice were kept in individually ventilated cages (IVCs; 30 cm by 16 cm by 16 cm; 501 cm2), serviced by a central airflow system (Tecniplast, IVC Green LineGM500), according to institutional guidelines. IVCs had sufficient bedding and nesting material as well as a wooden tunnel for environmental enrichment. Animals were maintained under pathogen-free, temperature-controlled (23 1C), and constant humidity (55 10%) conditions on a 12-hour light/12-hour dark cycle (lights on at 7:00 a.m.) with food and water provided ad libitum, at the Max Planck Institute of Psychiatry (Munich, Germany).
C57BL/6N males (7 weeks old) were exposed to the CSDS paradigm for 21 consecutive days, as previously described (12). Briefly, experimental mice were introduced daily into the home cage of a dominant CD1 resident mouse, which rapidly recognized and attacked the intruders within 2 min. To avoid serious injuries, the subordinate mouse was separated immediately after being attacked by the CD1 aggressor. After the physical encounter, mice were separated by a perforated metal partition, allowing the mice to keep continuous sensory but not physical contact for the next 24 hours. Every day, for a total of 21 days, mice were defeated by another unfamiliar, CD1 resident mouse, to exclude a repeated encounter throughout the experiment. Defeat encounters were randomized, with variations in starting time (between 9:00 a.m. and 6:00 p.m.) to decrease the predictability to the stressor and minimize habituation effects. Control mice were single-housed, in the same room as the stressed mice, throughout the course of the experiment. All animals were handled daily and weighed every 4 days. Coat state was scored on a scale of 0 to 3 according to the following criteria: (0) No wounds, well-groomed and bright coat, and clean eyes; (1) no wounds, less groomed and shiny coat, OR unclean eyes; (2) small wounds, AND/OR dull and dirty coat, and not clear eyes; (3) extensive wounds, OR broad piloerection, alopecia, or crusted eyes. End point and tissue collection were performed in the morning (9:00 a.m.) and 48 hours after the last social defeat session (day 23). This was done to capture the cumulative effects of chronic stress, rather than the acute effects of the last defeat session. The SAT was conducted during the last week of the CSDS paradigm, and based on their performance, five mice from each group were selected for molecular characterization, thus avoiding potentially stress-resilient animals. For end point, all mice were deeply anesthetized with isoflurane and perfused with cold phosphate-buffered saline (PBS), and target tissues were quickly dissected for molecular experiments. Cardiac blood was collected for the assessment of basal CORT levels. Adrenal glands were dissected from fat and weighed. The brains, pituitary, and adrenal glands from selected mice were immediately processed for RNA single-cell analysis. Tissues from all remaining mice were collected for downstream validation experiments.
Social avoidance behavior was assessed with a novel CD1 mouse in a two-stage social interaction test. In the first 2.5-min test (nontarget), the experimental mouse was allowed to freely explore the open-field arena containing an empty wire mash cage against one wall of the arena (labeled as the interaction zone). In the second 2.5-min test (target), the experimental mouse was returned to the arena with an unfamiliar male CD1 mouse enclosed in the wire mash cage. The ratio between the time in the interaction zone of the nontarget trial and the time in the interaction zone of the target trial was calculated and deemed as the interaction time ratio.
Blood sampling was performed during end point (9:00 a.m.) by collecting blood from the heart of each mouse before perfusion with PBS. All blood samples were kept on ice and centrifuged at 4C, and 10 l of plasma was removed for measurement of CORT. All plasma samples were stored at 20C until CORT measurement. CORT concentrations were quantified by radioimmunoassay (RIA) using a CORT double antibody 125I RIA kit (sensitivity: 25 ng/ml; MP Biomedicals Inc.) following the manufacturers instructions. Radioactivity of the pellet was measured with a gamma counter (Wizard2 2470 Automatic Gamma Counter; Perkin Elmer). All samples were measured in duplicate, and the intra- and interassay coefficients of variation were both below 10%. Final CORT levels were derived from the standard curve.
Six different groups (each N = 5) of C57BL/6N males were exposed to a different number of social defeat sessions to assess the cumulative effects of stress and its correlation with changes in mRNA levels of Abcb1b, Sbsn, and Srd5a2. The groups were defined as follows: (i) controlno defeat, (ii) 3 days of social defeat, (iii) 5 days of social defeat, (iv) 10 days of social defeat, (v) 21 days of social defeat, and (vi) 21 days of social defeat, followed by 48 hours of recovery time. The last group was introduced to match the end point of our original chronic social defeat paradigm cohort (23 days). All mice were 7 weeks old at the beginning of the experiment. Individual social defeat encounters were carried out exactly as previously explained in the CSDS paradigm section of the methods. End point and tissue collection were performed in the morning (9:00 a.m.) of day 23. All mice were deeply anesthetized with isoflurane, and the adrenal glands were quickly dissected for molecular experiments. Adrenal glands were further dissected from fat and weighed. Trunk blood was collected for the assessment of basal CORT levels.
Tissue dissociation. Mice were anesthetized lethally using isoflurane and perfused with cold PBS to get rid of undesired blood cells in target tissues. Brains, pituitaries, and adrenal glands were quickly dissected and immediately transferred to ice-cold oxygenated artificial cerebral spinal fluid (aCSF) (brains), ice-cold Hanks balanced salt solution (HBSS) (pituitaries), or ice-cold PBS (adrenals) and kept in the same solutions during dissection and dissociation. The aCSF was oxygenated throughout the experiment with a mixture of 5% CO2 in O2. Sectioning of the brain was performed using a 0.5-mm stainless steel adult mouse brain matrix (Kent Scientific) and a Personna Double Edge Prep Razor Blade. A slide (approximately 0.58 mm Bregma to 1.22 mm Bregma) was obtained from each brain, and the extended PVN was manually dissected under the microscope. Two cell suspensions were prepared for each of the three tissues with one pool for controls and one pool for stressed mice. The PVN from five different mice was pooled and dissociated for 35 min using the Papain Dissociation System (Worthington) following the manufacturers instructions. Similarly, the pituitaries from five mice were pooled and dissociated for 15 min using papain. Last, the adrenal glands from five mice were pooled and dissociated for 55 min using a 0.2% collagenase II solution. All cell suspensions were incubated at 37C using a shaking water bath. After this, cell suspensions were filtered with 30-m filters (Partec) and kept in cold aCSF, HBSS, or PBS.
Cell capture, library preparation, and high-throughput sequencing. Cell suspensions with approximately 1,000,000 cells/l were used. Each pool was loaded onto individual lanes of a 10X Genomics Chromium chip, as per factory recommendations. For all three tissues, the control and stress cell suspensions were loaded and processed together in the same chip to avoid batch effects by condition. Reverse transcription and library preparation were performed using the 10X Genomics Single-Cell v2.0 kit following the 10X Genomics protocol. Molar concentration and fragment length of libraries were quantified by qPCR using KAPA Library Quant (Kapa Biosystems) and Bioanalyzer (Agilent High Sensitivity DNA kit), respectively. Each library was sequenced on a single lane of an Illumina HiSeq4000 System generating 100base pair paired-end reads at a depth of ~340 million reads per sample.
Preprocessing and quality control. Preprocessing of the data was done using the 10X Genomics Cell Ranger software version 2.1.1 in default mode. The 10X Genomics supplied reference data for the mm10 assembly and corresponding gene annotation was used for alignment and quantification. All further analyses were performed using Scanpy (version 1.4.4.post1) (13), following guidelines from an established best practices workflow (14). For quality control, we looked at the distribution of count depth, number of genes, and mitochondrial read fraction per sample. Because distributions were fairly homogeneous, we chose to pick the same thresholds for all samples (tissues and conditions). Specifically, we filtered out (i) cells with less than 1000 counts, (ii) less than 400 genes detected, and (iii) percentage of mitochondrial gene counts higher than 25%. In addition, genes expressed in less than 20 cells were removed as well. Quality control (QC) plots can be found in fig. S5 (E to G). This resulted in a dataset of 21,723 single cells, of which 6966 cells and 16,168 genes were from the PVN, 9879 cells and 15,437 genes were from the pituitary, and 4878 cells and 13,997 genes were from the adrenal gland. The size factors used for normalization were obtained using Scran (version 1.14.5) (56), and the data were log1 Ptransformed. Each dataset was batch-corrected using Combat (57), available in Scanpy. For each tissue, we selected the top 4000 highly variable genes using the highly_variable_genes function. Dimensionality reduction was performed using principal components analysis computed on highly variable genes and taking the first 50 PCs. Last, we computed a k-nearest neighbor graph (KNN)graph (k = 15) on the low-dimensional embedding, necessary for UMAP visualization.
Clustering, marker gene identification, and cluster annotation. Data were clustered using the Louvain (version 0.6.1) algorithm implemented in Scanpy (13). This is a graph-based clustering method that relies on the KNN-graph discussed above. We clustered at two different resolution levels (r = 0.5 and r = 1). After inspection of the cell clusters, we observed that those obtained using a finer resolution (r = 1) aligned better with our annotations and therefore used them for visualization and downstream analyses. Marker genes for each cluster were detected using a Welchs t test between cells in the cluster and all cells outside of it as reference. This was done using the rank_genes_groups function implemented in Scanpy and computed on log-normalized nonbatch-corrected data. Cell types were determined using a combination of marker genes identified from the literature and Gene Ontology for cell types using the web-based tool: mousebrain.org (58).
Differential expression analysis. Differential expression analyses were performed using MAST (59) implemented in R, which models scRNA-seq data using a generalized linear model (GLM). The computation was performed on log-normalized nonbatch-corrected data, and, for each cell cluster, we fit the following model: [Y ~ 1 + condition + n_genes], where Y is the log-normalized nonbatch-corrected data, 1 is the intercept term, condition is the covariate that accounts whether the mouse was stressed or not, and n_genes is used as a technical covariate as a proxy for technical and biological factor that might influence gene expression. The test produced a P value for each gene in each cell cluster and a q value, which is the P value after adjustment for multiple testing, using false discovery rate (FDR) correction. Furthermore, the mean expression of each gene for the two different conditions was computed.
Ambient RNA assessment. After QC analyses, we noticed the presence of highly expressed genes across all cells, despite being known marker genes of specific cell types. We noticed that most of these genes coded for neuropeptides or hormones and decided to assess whether we could explain this as ambient RNA contamination. To investigate which genes were expressed as ambient RNA, we analyzed the unfiltered datasets for the three tissues. We once again looked at the count depth distribution for what we conclude are empty droplets and selected cells with counts between 100 and 300 for the PVN, between 300 and 600 for the pituitary, and between 50 and 200 for the adrenal gland. We also removed genes that are expressed in less than 20 cells. After these steps, we obtained a dataset of 120,320 droplets and 14,129 genes for the PVN, 113,043 droplets and 13,777 genes for the pituitary, and 107,698 droplets and 11,684 genes for the adrenal gland. Because these are empty droplets and we do not expect any meaningful clustering of the data, we used a less sophisticated normalization technique, normalizing each cell by total counts over all genes, thus obtaining the same total count per cell after this step. The number of counts per cell to obtain was selected automatically as the median count per cell before normalization. For all tissues, we computed the mean expression of each gene across all droplets by condition (stress versus control). Furthermore, to exclude from our list of significantly DEGs those that are detected owing to differential ambient RNA expression, we performed differential expression testing using MAST across all droplets using the same GLM formulation defined above (note that, in our previous analysis, we tested within each cluster).
For paraffin embedding, adrenal glands were dissected and the surrounding fat was removed and fixed in 10% neutral buffered formalin (Sigma-Aldrich, HT501128) overnight at room temperature. Tissue was embedded manually over 3 days. All washes were carried out for 1 hour at room temperature unless indicated. Day 1: three times PBS, 25% EtOH, 50% EtOH, 70% EtOH, and 70% EtOH overnight at 4C. Day 2: 80% EtOH, 90% EtOH, 95% EtOH, and 100% EtOH overnight at 4C. Day 3: 100% EtOH, Neoclear (Sigma-Aldrich, 109843) I for 10 min at room temperature; Neoclear II for 10 min at 60C; Neoclear: paraffin 1:1 for 15 min at 60C, paraffin I for 1 hour at 60C, paraffin II for 1 hour at 60C, and paraffin III for 1 hour at 60C. Samples were sectioned at 5 m. RNAscope was carried out on paraffin-embedded sections with the RNAscope 2.5 HD Kit-BROWN (ACD bio 322300) assay following the manufacturers protocols, with Standard timings for retrieval and protease treatment. The following probes were used (all ACD bio): Mm-Abcb1b (422191), Mm-Sbsn (564441), Mm-Srd5a2 (431361), Hs-ABCB1 (401191), and Hs-SBSN (447411). Positive control Mm-Ppib (313911), Hs-UBC (310041), and negative control dapB (310043) were also used. Nuclei were stained with Vector Hematoxylin QS (Vector Laboratories, H-3404), and slides were mounted in VectaMount Permanent Mounting Medium (Vector Laboratories, H-5000).
Paraffin sections were deparaffinized and rehydrated as per immunohistochemistry. Slides were incubated for 30 s with Hematoxylin QS, washed with running water, incubated for 30 s with eosin, washed with running water, and mounted in VectaMount Permanent Mounting Medium.
Hematoxylin and eosin and RNAscope slides were scanned with a NanoZoomer-XR digital slide scanner (Hamamatsu). Images were processed with NanoZoomer digital pathology view (Hamamatsu), and quantification was done with Fiji.
Four areas of the same dimensions (252 252 pixels) were selected from the zona fasciculata of the cortex. Nuclei were counted, and the average cell area was calculated by dividing the number of nuclei by the total area. Values were multiplied by 1000 for graphical representation.
Quantification of messenger RNA levels of Abcb1b, Sbsn, and Srd5a2 in the adrenal glands was carried out using qRT-PCR. Total RNA was reverse-transcribed using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems). RT-PCR reactions were run in triplicate using the ABI QuantStudio 6 Flex RT-PCR System and data were collected using the QuantStudio RT-PCR software (Applied Biosystems). Expression levels were calculated using the standard curve, absolute quantification method. The endogenous expressed gene Hprt was used to normalize the data. The following Taqman probes were used: Abcb1b: Mm00440736_m1, Sbsn: Mm00552057_m1, Srd5a2: Mm00446421_m1, and Hprt: Mm00446968_m1.
Mouse Y1 cells and human NCI H295R adrenocortical cells were seeded into 12-well plates and incubated overnight using Dulbeccos modified Eagles medium high glucose (4.5 g/liter) (Gibco) with 7.5% horse serum (Gibco), 2.5% fetal bovine serum (FBS) (Gibco), and 1% penicillin-streptomycin (Gibco) and RPMI 16/40 + GlutaMax (Gibco) with 10% FBS (Gibco), 1% Insulin-Transferrin-Selenium-Ethanolamine (ITS) (Thermo Fisher Scientific), and 1% penicillin-streptomycin (Gibco), respectively. In this experiment, 100,000 Y1 and NCI H295R adrenocortical cells per well were used. Cells were then stimulated for 24 hours with 10 nM forskolin and subsequently treated with different concentrations of tariquidar (0, 10, 50, 125, 250, 500, and 1000 nM) and incubated for 24 hours. Last, supernatants and cell pellets were collected and harvested for further analyses and measurement of CORT (ng/ml) and cortisol (g/liter) levels. Media CORT levels in Y1 cells were quantified by RIA using a CORT double antibody 125I RIA kit, as previously described in the animal experiments. Media cortisol levels in NCI H295R adrenocortical cells were determined using an enzyme-linked immunosorbent assay (ELISA) kit (RE52061, TECAN, IBL Hamburg, Germany). The standard range was 20 to 800 ng/ml.
The study was approved by the Ethics Committee of the University of Wuerzburg (Germany) (#88/11), and written informed consent was obtained from all subjects. Eight patients with biochemically confirmed persistent ACTH-dependent Cushings syndrome were studied. Cushings syndrome was established according to current guidelines (60). Half of the patients (n = 4) had pituitary-dependent Cushings syndrome, while, in the other patients (n = 4), ectopic Cushings syndrome had been diagnosed. The patients underwent bilateral adrenalectomy as ultima ratio to control life-threatening hypercortisolism after other therapies had failed. Formalin-fixed paraffin-embedded sections were stained as described above. The normal adrenal tissue was derived from adrenal glands removed as part of tumor nephrectomy (n = 6). They were histologically proven normal adrenal glands without neoplastic tissue.
Data of the Emory Predictors of Remission in Depression to Individual and Combined Treatments (PReDICT) (61, 62) study was used to investigate effects of the ABCB1 variant rs2032582 on HPA axis function in 154 unmedicated patients with a current Diagnostic and Statistical Manual of Mental Disorders (DSM)IV diagnosis of major depressive disorder. The PReDICT study was designed and conducted in accord with the latest version of the Declaration of Helsinki. The Emory Institutional Review Board (IRB) and the Grady Hospital Research Oversight Committee gave ethical approval for the study design, procedures, and recruitment strategies (Emory IRB numbers 00024975 and 00004719). The PReDICT study is registered at ClinicalTrials.gov Identifier: NCT03226912 and NCT00360399. DNA was extracted from whole blood, and genome-wide SNP genotyping was performed using HumanOmniExpress BeadChips. Quality control was performed in PLINK. Samples with low genotyping rate (<98%) were removed. SNPs with a high rate of missing data (>2%), significant deviation from the Hardy-Weinberg equilibrium (HWE, P < 105), or a low minor allele frequency (<5%) were excluded from further analyses. SNP genotypes were coded as 0 for major homozygotes (GG, n = 56), 1 for heterozygotes (TG, n = 74), and 2 for minor homozygotes (TT, n = 24) and did not deviate from HWE (2 = 0.27, P = 0.60). HPA axis function was assessed using the dexamethasone/corticotropin-releasing hormone (Dex/CRF) test, consisting of an oral administration of 1.5 mg of Dex at 11:00 p.m. and an infusion of ovine CRF (1 g/kg) at 3:00 p.m. on the next day. Cortisol and ACTH levels were measured from plasma samples taken immediately before CRF administration (pre-CRF) (i.e., at 3:00 p.m.) and again at 3:30 p.m. (30 min), 3:45 p.m. (45 min), 4:00 p.m. (60 min), and 4:15 p.m. (75 min). Baseline cortisol levels were available for all 154 patients with genotype data for the SNP rs2032582. Only completers were included in the analysis, so there are no dropouts in sample sizes over time. Linear regression models were used to test for effects of the SNP genotype on baseline cortisol levels using R version 3.6.2. To assess differences in cortisol and ACTH levels after the Dex/CRF test over time, linear mixed effects models with a random intercept for each patient were applied. All models included gender, age, and baseline depression severity sum scores on the 17-item Hamilton Depression Rating Scale (63) and the first five genetic ancestry (multidimensional scaling) components as covariates.
See the original post here:
Single-cell molecular profiling of all three components of the HPA axis reveals adrenal ABCB1 as a regulator of stress adaptation - Science Advances