A cross-tissue transcriptomic approach decodes glucocorticoid receptor-dependent links to human metabolic phenotypes

Abstract

Glucocorticoids, acting through the glucocorticoid receptor (GR), control metabolism, maintain homeostasis, and enable adaptive responses to environmental challenges. Their function has been comprehensively studied, leading to identification of numerous tissue-specific GR-dependent mechanisms. Abundant evidence shows that GR-triggered responses differ across tissues, however, the extent of this specificity was not comprehensively explored. It is also unknown how particular GR-induced molecular patterns are translated into profile of higher-level human traits. Here, we examine cross-tissue effects of GR activation on gene expression. We assessed changes induced by stimulation with GR agonist, dexamethasone in nine tissues (adrenal cortex, perigonadal adipose tissue, hypothalamus, liver, kidney, anterior thigh muscle, pituitary gland, spleen, and lungs) in adult male C57BL/6 mice, using whole-genome microarrays. Dexamethasone induced balanced transcriptional responses across all examined tissues with 585 identified dexamethasone-regulated transcripts, including 446 with significant treatment-tissue interaction effects. Clustering analysis revealed sixteen GR-dependent patterns, including those universal across tissues and tissue-specific. We leveraged existing gene annotations and created new annotation sets based on chromatin immunoprecipitation sequencing, recent large-scale genome-wide association studies, and human transcriptome collections. As expected, GR-dependent transcripts were associated with essential metabolic processes (glycolysis/gluconeogenesis, lipid-metabolism) and inflammation-related pathways. Beyond these, we found novel links between regulated gene patterns and human phenotypic traits, like reticulocyte count or blood triglyceride levels. Overall effects of GR stimulation are well coordinated and closely linked to biological roles of tissues and organs. Our findings provide novel insights into complex systemic and tissue-specific actions of glucocorticoids and their potential impacts on human physiology and pathology.

Background

Glucocorticoids are steroid hormones that regulate the biology of organs from the early development, through cell differentiation, to adult physiology [1,2,3]. Physiologically, glucocorticoids are secreted in a circadian mode as effectors of the hypothalamus – pituitary – adrenal (HPA) axis, acting to synchronize metabolism across the organism [4]. In response to increased glucocorticoid levels, the liver increases gluconeogenesis, while muscles and adipose tissue increase protein catabolism and lipolysis, respectively, to provide substrates for glucose production by the liver. This sequence leads to an increase in circulating glucose levels that sustains the functioning of the central nervous system (CNS) [5]. The activation of the HPA axis is particularly important for the coordination of stress response. A repetitive or robust increase in glucocorticoid signaling elicited by stress is considered a general risk factor contributing to the development of pathologies, notably metabolic and psychiatric disorders. Glucocorticoids regulate vascular tone, bone mineralization, and activation of the immune system. Due to their potent immunosuppressant ability, synthetic glucocorticoids are some of the most commonly prescribed drugs in the world today [6]. However, prolonged treatment with glucocorticoids is associated with serious adverse effects, which may include weight gain, impaired skin healing, mood disorders, euphoria, hypertension, increased risk of infection, and high intraocular pressure. Therefore, understanding of the molecular rules of systemic GR stimulation to tissue-specific actions is mandatory for better management of multiple disorders associated with distinct organs.

Glucocorticoids act through glucocorticoid receptors (GRs) and mineralocorticoid receptors (MRs). GR belongs to the nuclear receptor superfamily of transcription factors and canonical effects of glucocorticoids are mediated through transactivation or transrepression of gene expression [7, 8]. Transactivation is based on GR dimerization, translocation to the nucleus, and binding to proximal or distal response elements. Examples of genes activated in this mechanism include Fkbp5 [9] and Sgk1 [10]. An alternative mechanism, independent of direct GR binding to DNA, involves protein–protein interactions with other transcription factors [11], such as STAT5 as observed in the case of Igf1 [12]. Conversely, GR has also been reported to repress transcription through various mechanisms that include direct DNA binding, protein–protein interactions, employment of GR isoforms, and secondary, indirect molecular mechanisms [13,14,15].

While GR is ubiquitously expressed in almost every cell in an organism, the effects of its activation are highly tissue- and organ-specific [16, 17]. Consequently, the vast majority of published data focus on the effects of GR activation in discrete cell types and tissues. The information including the systemic context is scarce and may be found in meta-analyses and reviews on tissue-specific effects of GR elimination on general physiology and inflammation control [18], cross-tissue analysis of GR cofactors [19], cross-tissue glucocorticoid metabolism [20] or GR-induced chromatin accessibility [21]. Comprehensive studies of glucocorticoid-induced gene expression profiles permit the identification of tissue-type specific expression patterns, and open the possibility of developing synthetic glucocorticoids with targeted, organ-specific actions and reduced adverse effects [22]. However, these efforts have thus far been limited (a notable exception includes a study on the primary human blood cell types) [23].

To fill this gap, we investigated parallel transcriptional effects of selective GR stimulation in vivo across nine mammalian tissues, including key elements of the HPA axis. Using rigorous statistical analyses, we identified shared and tissue-specific gene clusters regulated by GR activation. Subsequent analyses revealed that the engagement of distinct co-factors shapes local alterations in the transcriptional landscape upon GR stimulation. Finally, we investigated the implications of these discoveries for human health and disease. To our knowledge, this is the first comprehensive gene expression profiling study across multiple tissues, linking local profiles of GR-dependent responses with traits of human disorders.

Materials and methods

Animals

Adult male (8 to 10 weeks old) C57BL/6 J mice were housed in groups of 5–7 per cage (31 × 16 × 14 cm), under a 12-h dark/light cycle (lights on 7:00 AM, lights off 7:00 PM), with free access to food and water. Animals aged 8–10 weeks and weighing 20 to 30 g were used for all the experiments. While a limitation, using only male mice controlled for sex-based differences and enabled direct comparisons with previous male-focused studies. Future experiments should be performed in female mice to assess the differences in transcriptional response to dexamethasone. Experiments were performed on 24 C57BL/6 J male mice bred at the Maj Institute of Pharmacology of the Polish Academy of Sciences animal facility. The animal protocols used in the study were approved by the II local ethics committee at the Maj Institute of Pharmacology PAS (1156/2015, Krakow, Poland). The experiments were planned and executed in accordance with the ARRIVE guidelines [24], European and Polish laws concerning the use and welfare of laboratory animals (Directive 2010/63/UE, European Convention for the Protection of Vertebrate Animals Used for Experimental and other Scientific Purposes ETS No.123, and Polish Law Dz.U. 2015 poz. 266).

Drug treatment and tissue collection

Mice were killed by decapitation 4 h after a single dexamethasone (DEX, 10 mg/kg) i.p. injection, while saline-treated mice (10 ml/kg) served as the control group. Injections were performed between ZT0 and ZT2 (early light phase, corresponding to the nadir of endogenous corticosterone) to minimize interference from the circadian rhythm, and all the mice were killed at the same time of the day, 4 h later (ZT4 to ZT6). The time point was selected based on our previous studies investigating acute transcriptional responses in vivo, as it typically captures both primary GR target gene induction and the initiation of secondary effects [25, 26]. All the mice were killed at the same time of the day. The dose of DEX was based on our [27] and others [28, 29] previous experiments, paying particular attention to avoiding systemic toxic effects and overcoming the blockade of the blood–brain barrier. Tissue samples (adrenal cortex, perigonadal adipose tissue, hypothalamus, liver, kidney, anterior thigh muscle, pituitary gland, spleen, and lungs) were fixed using the RNAlater reagent (Qiagen Inc., Valencia, CA, USA) and stored at −70 °C. Left–right symmetric organs were pooled within one sample. Pituitary tissue included both the anterior and posterior pituitary, while adrenal glands comprised both the cortex and medulla.

RNA preparation

RNA was isolated using the RNeasy Mini Kit (Qiagen Inc.) and further purified following the manufacturer’s protocol. Total RNA concentration was measured using a NanoDrop ND-1000 Spectrometer (NanoDrop Technologies Inc., Montchanin, DE, USA). RNA quality was determined by chip-based capillary electrophoresis using RNA 6000 Nano LabChip Kits and an Agilent Bioanalyzer 2100 device (Agilent, Palo Alto, CA, USA). RNA from two mice was randomly pooled to prepare samples for each microarray to reduce biological variability. For each tissue and treatment condition, 4–6 independent pooled samples were processed, resulting in a total of 83 microarrays passing quality control.

Gene expression profiling

A starting amount of 200 ng high-quality total RNA was used to generate cDNA and cRNA with the Illumina TotalPrep RNA Amplification Kit (Illumina Inc., San Diego, CA, USA). The procedure consisted of reverse transcription with an oligo(dT) primer that also included a T7 promoter sequence using Array-Script. The obtained cDNA was used as a template for in vitro transcription with T7 RNA polymerase and biotin-labeledUTP, which generated multiple copies of biotinylated cRNA. The purity and concentration of the cRNA were checked using an ND-1000 Spectrometer. Validated cRNA was then hybridized with Illumina’s direct hybridization array kit (Illumina). Each cRNA sample (1.5 μg) was hybridized overnight to the MouseWG-6 v2 BeadChip arrays (Illumina) in a multiple-step procedure according to the manufacturer’s instructions; the chips were washed, dried, and scanned on the BeadArray Reader (Illumina). Raw microarray data was generated using BeadStudio v3.0 (Illumina).

Microarray data analysis

Microarray quality control was performed using BeadArray R package v2.2.0. The following parameters were checked on all the arrays: number of outliers, number of beads, and percent of detected probes. After background subtraction, the data were normalized using quantile normalization and log₂-transformed. The obtained signal was taken as the measure of mRNA abundance derived from the gene expression level. Statistical analysis of the results was performed using two-way ANOVA (for tissue and treatment) followed by Tukey’s HSD post-hoc tests (where appropriate). The false discovery rate (FDR) was estimated using the Benjamini and Hochberg method [30]. All statistical analyses were performed using the R software (version 4.1.1). For clustering and visualization, the gene expression data underwent several preprocessing steps. Initially, the data was scaled using the R scale function to standardize the expression values across genes. The saline condition was then established as the baseline (control, CTR) by subtracting its median value from all samples, effectively highlighting the treatment-induced changes. To mitigate the impact of extreme outliers, expression values were capped at a threshold of ± 5 standard deviations from the mean. For clustering, 1 minus the correlation distance was used. Hierarchical clustering was then performed using the R hclust function with the complete linkage method. Finally, the resulting dendrogram was partitioned into 16 distinct clusters using the R cutree function (specifying k = 16). This number was chosen empirically to balance granularity and interpretability in the gene expression patterns observed in the heatmap (Fig. 2) and ensure most clusters were sufficiently large for meaningful enrichment analysis.

Evaluation of DEX impact on tissue-level transcriptional response

We employed a comprehensive statistical approach to evaluate the differential impact of DEX across various tissues. First, we calculated the log₂ ratio of gene expression levels between DEX-treated and CTR samples for each tissue. We selected the top 50 significantly altered transcripts (identified by probe IDs) with the highest absolute log₂ ratios (DEX vs. CTR within that tissue, among probes significant in the ANOVA) for each tissue to focus on the most responsive genes. We then computed the median log₂ ratio for these top 50 transcripts in each tissue to represent the overall transcriptional response magnitude. To assess whether the differences in response across tissues were statistically significant, we performed a one-way ANOVA with tissue as a factor and the log₂ ratios as the dependent variable. Following the ANOVA, we conducted post-hoc analyses using Tukey’s Honest Significant Difference (HSD) test to identify specific pairwise differences between tissues. This test allowed us to determine which tissues exhibited significantly different responses from others while controlling for multiple comparisons. The ANOVA was conducted using the R aov function, and Tukey’s HSD test was performed using the R TukeyHSD function. Statistical significance was set at nominal p < 0.05 for impact on tissue analyses. The results of these analyses were visualized using a box plot created with ggplot2 (version 3.3.5).

Gene enrichment analysis

Overrepresentation analyses were performed using Enrichr [31], a comprehensive gene set enrichment analysis tool that integrates data from multiple genomic resources. Lists of unique Entrez Gene IDs derived from the probes in each cluster were used as an input. Transcripts were assigned to clusters using the complete agglomeration method. Gene lists derived from each cluster were input into Enrichr to identify overrepresented biological processes (track GO Biological Process 2023) drug effects (DSigDB, Drug Perturbations from GEO 2014), and transcription factor binding sites based on ChIP-seq datasets (ChEA 2022). Enrichr utilizes statistical measures, including Fisher’s exact test, to determine the significance of enrichment. The resulting p-values were adjusted for multiple testing using the Benjamini–Hochberg procedure to control for false discovery rate (FDR). The level of significance was set at FDR < 0.01 and at least two DEX-regulated genes in a gene set. For tissue-selectivity confirmation of enriched TFs we analyzed expression levels of all TFs identified through Enrichr analysis of ChIP-seq datasets across relevant tissues using two complementary approaches. First, we utilized RNA-seq data from the GTEx project, where gene expression levels were averaged across all available samples for each tissue and normalized as counts per million (CPM). TFs were considered expressed in a tissue if their mean expression exceeded 50 CPM. Second, we examined microarray data generated in our study, where raw expression values were log₂-transformed and averaged for each tissue following saline treatment. Here, TFs were deemed expressed when their mean signal intensity surpassed a threshold of log₂ > 8.

Gene overlap analyses

Three datasets containing lists of (I) GR-dependent genes, (II) genes linked to metabolic traits, and (III) genes associated with various human phenotypes were prepared based on external resources. The lists of GR-regulated genes were assembled from a PubMed database search for publications on GR-dependent gene expression in different tissues of mice, rats, or humans (see Supplementary Table 2 for a full list of PMIDs). Only lists with five or more genes were included in the database. Gene lists linked to human metabolite levels were extracted from Metabolon, and Nightingale (https://www.omicspred.org/Scores) databases. The resulting set of lists was filtered to exclude all instances with fewer than 3 genes, resulting in a total of 520 and 137 lists from Metabolon and Nightingale, respectively. The third set of lists was derived from genes associated with 729 human phenotypes obtained from the Pan-UK Biobank database [32]. Lists of genes containing variants showing GWAS association with phenotypic traits were extracted at p-value thresholds of p < 10^–8. We performed separate overlap analyses for each of the datasets with the lists of genes included in the clusters of GR-regulated transcripts probes identified in the current study. Mouse Entrez Gene IDs corresponding to the probes in each cluster were mapped to their human orthologs using current annotations from Ensembl. Statistical analyses were performed using the chi-square test, followed by FDR correction for multiple tests. The total number of, 19,437 human genes (protein-coding genes downloaded from Biomart 110 [33] corresponding to the universe of potential orthologs) was used as a background test set for these human phenotype/metabolite overlap analyses. The level of significance was set at FDR < 0.01.

Results

Dexamethasone regulated transcripts

To determine tissue-specific profile of GR activation, we administered i.p. a synthetic GR agonist, dexamethasone (DEX) or vehicle (CTR) to adult mice. Four hours after the injection, 9 tissues were isolated, RNA was extracted, and processed for transcriptomic profiling. Two-way ANOVA was performed with tissue and treatment factors to identify transcripts with altered abundance after DEX treatment. We found: (1) 470 DEX-regulated transcripts (p_adj < 0.01 for treatment factor in ANOVA, with absolute log₂ ratio > 1 in at least one tissue), (2) 446 transcripts regulated by DEX in a tissue-specific manner (p_adj < 0.01 for tissue-treatment interaction, with absolute log₂ ratio > 1 in at least one tissue). There was a strong overlap between these two lists, with 331 transcripts meeting both criteria. We used a comprehensive set of 585 differentially expressed transcripts identified by probe IDs (mapped to 362 unique Entrez gene IDs) for further analyses (Supplementary Table 1), representing all genes regulated by DEX that were present in either or both lists.

Impact of dexamethasone on transcription in different tissues

To assess the specific impact of DEX on different tissues, we calculated the median log₂ ratio of the expression level (DEX vs. CTR) for the top 50 transcripts within each tissue (Fig. 1). The analysis revealed significant differences in the transcriptional response across tissues (ANOVA, F(8, 441) = 13.47, p < 1e⁻¹⁵). The strongest response was observed in the kidney (KID, 1.31 mean log₂ ratio), with a significantly higher median fold change compared to all other tissues (Tukey’s HSD, p < 0.05). The effects of DEX were noticeable in other tissues, with the following descending order: adipose tissue (FAT, 1.16), pituitary (PIT, 1.12), muscle (MUS, 1.08), lung (LUN, 1.07), liver (LIV, 1.01), spleen (SPL, 1.01), and the adrenal gland (ADR, 0.93). The hypothalamus (HTH, 0.67) exhibited the lowest transcriptional changes, which were significantly lower than those in all other tissues (Tukey’s HSD, p < 0.05).

Fig. 1

Dexamethasone treatment affects gene expression in different tissues. Barplots (mean ± SD) of log₂ ratio for the top 50 transcripts are shown for each tissue to represent overall tissue impact. ADR—adrenal cortex, FAT—perigonadal adipose tissue, HTH—hypothalamus, KID—kidneys, LIV—liver, LUN—lungs, MUS—anterior thigh muscle, PIT—pituitary gland, SPL—spleen. Significant differences in fold changes between tissues obtained by a Tukey HSD test are indicated by # (vs. Kidney; p < 0.05) and by * (vs. Hypothalamus; p < 0.05)

Dexamethasone-regulated gene expression patterns

To identify gene expression patterns, we used a list of 585 unique microarray probe values modified in response to DEX administration. Upon standardization for tissue-specific basal expression, hierarchical clustering of these probe values revealed 16 clusters (Fig. 2, clusters named from A to P). These clusters were selected for further analysis based on containing at least 6 unique corresponding Entrez Gene IDs to allow for meaningful enrichment analysis. The number of probes belonging to individual clusters ranged from 8 to 117, with a median of 23. The corresponding number of unique genes per cluster ranged from 6 to 93 (Supplementary Table 1). Six clusters (A to F) contained transcripts mostly downregulated by DEX, while 10 clusters (G-P) contained probes generally upregulated by DEX. Notably, only a few clusters (D, K, I, O, P) featured transcripts regulated similarly across all tissues. In these sets, we found numerous known GR-dependent transcripts, such as members of the clock machinery (including Bhlhe40, Dpd, Per1), regulators of inflammatory response (Cdkn1a, Cxcl10, Sphk1), MAPK signaling pathway (Dusp1, Dusp4, Map3k6) and many other bona fide GR targets (Tsd22d3, Ddit4, Sgk1, Fkbp5). The majority of clusters (A, B, C, E, F, G, H, J, K, M, N) contained transcripts preferentially regulated in selected tissues. This data point to a highly tissue-specific pattern of transcriptional regulation by DEX.

Fig. 2

Hierarchical clustering of dexamethasone-induced transcriptional alterations. Relative levels of all transcripts (585) with differential expression are shown as a heat map (Supplementary Table 1). Data were standardized for the expression level in each tissue. Colored rectangles represent transcript abundance 4 h after injection of DEX or CTR in the specific tissue, as indicated below the heatmap. The intensity of the color is proportional to the standardized values (between −4 and 4) from each microarray (according to the scale below the heatmap). Gene clusters are depicted as colors and letters (A-P) on the left. Clustering was performed using Pearson correlation as distance, the results are summarized on the dendrogram shown on the right. Clusters A-F are classified as down-regulated (DOWN). Clusters G-P are classified as up-regulated (UP) in response to DEX

In silico analysis of promoter regions of differentially-regulated genes

Next, we examined the transcriptional mechanisms influenced by the GR. To this end, we analyzed the overrepresentation of GR binding sites across identified clusters, using ChIP-seq data processed through Enrichr (ChEA 2022) (Supplementary Table 3). The analysis revealed a significant overrepresentation of GR (encoded by NR3C1 gene) binding sites in clusters P (15 genes, p_adj = 8.3*10^–6), K (9, p_adj = 3.0*10^–5), O (10, p_adj = 1.1*10^–3), I (7, p_adj = 2.9*10^–3), and D (14, p_adj = 0.037), that is, exclusively in clusters containing transcripts regulated by DEX throughout multiple tissues (compare with Fig. 2). To explore the transcriptional control beyond the GR, we identified the transcription factor binding sites (TFBS) overrepresented (p_adj < 1*10^–1) in each cluster (see Fig. 3). NR3C1 was found in the top 10 significantly enriched TFBS of four clusters (I, K, O, P). For remaining clusters, a variety of TFBS were identified, representing known GR transcriptional cofactors and effectors: LXR (clusters H, I, L, upregulated), RXR (clusters H, I, L, upregulated), CLOCK (clusters K, P), EZH2 (clusters F, O), PPARA (cluster I, upregulated), RNF2 (cluster O), SOX2 (clusters D, M), and TP53 (cluster O). This data indicate that tissue-specific transcriptional effects of GR stimulation engage discrete, tissue-specific cofactors. To further explore tissue specificity, we quantified baseline mRNA expression of genes encoding these enriched TFs across the nine examined tissues using control samples (Fig. 3, Supplementary Table 3). We discovered that 94.6% of overrepresented TFs can be found (mean log₂ expression value > 8) in tissues where they are enriched.

Fig. 3

The spectrum of putative transcriptional regulators among the DEX-induced gene expression patterns. The heatmap presents the statistical significance of enrichment of binding sites for transcription factors (TF) in the regulatory regions of genes from the clusters identified in this work. The information about TFs binding data were obtained from the ChEA 2022 database. The columns represent GR-responsive genes grouped in clusters labeled A to P. UP and DOWN are combined groups of all DEX-increased and DEX-decreased transcripts, respectively. For each cluster, the number of genes and the observed direction of the transcriptional regulation (‘ + ’ increase, ‘-’ decrease) are indicated. Statistically significant results are represented by colored squares: p_adj. < 0.2 (gray); p_adj. < 0.1 (light red); p_adj. < 0.05 (medium red); p_adj. < 0.01 (dark red); p_adj. < 0.001 (deep red). TFs were annotated as specified in the legend. The dots represent tissue co-localization between DEX-induced gene expression patterns and corresponding TF expression in at least one tissue associated with the gene cluster (as detailed in Supplementary Table 3). Up to 10 top statistically significant results per cluster with a minimum of two genes per term were included. Light green dots represent TF expression level in tissues based on GTEx data (mean CPM > 50) [34]. Dark green dots represent TF expression measured in our dataset (mean log₂ value for the control group > 8). Full results are listed in Supplementary Table 3

Functional classification of regulated transcripts

To gain insight into the biological processes associated with the observed gene expression patterns, we performed gene ontology (GO) enrichment analysis, using Enrichr. Restricting the analysis to terms that included at least two genes per term, we identified 336 GO term overrepresentations in all 16 clusters (280 unique terms, p < 0.05; Supplementary Table 4). This analysis revealed several key biological processes with broad, high-level terms, overrepresented in multiple clusters and low-level terms, associated with individual clusters. The criteria of being overrepresented in at least four clusters was met by terms: “Negative Regulation Of Cellular Process” (GO:0048523; clusters D, K, O, P) and “Positive Regulation Of Cell Population Proliferation” (GO:0008284; clusters B, D, K, O). Interestingly, in both cases, there are different transcripts contributing to the overlap in the particular clusters. Genes involved in “Negative Regulation of Cellular Processes” include Irf1, Il1b, Sox7, Ifit3 (cluster D), Bmp4, Txnip, Axp11, Tob2 (cluster K), Cdkn1a, Tex, Rgcc, Sox17, Dcun1d3 (cluster O), Dusp1, Tp53inp1, Rhob (cluster P). Genes involved in “Positive Regulation of Cell Population Proliferation” include Wnt3a, Sox11 (cluster B—tissue-selective), Fgf7, Esm1, Clec7a, Il1b, S1pr2 (cluster D), Bmp4, Tsc22d1, Aqp11 (cluster K), Cdkn1a, Sphk1, Slc25a33, Fgfr2 (cluster O).

Gene ontology terms overrepresented in at least three clusters include: Apoptotic Process (GO:0006915), Negative Regulation Of Multicellular Organismal Process (GO:0051241), Regulation Of p38 MAPK Cascade (GO:1,900,744), Positive Regulation Of Cold-Induced Thermogenesis (GO:0120162), Positive Regulation Of Multicellular Organismal Process (GO:0051240), Regulation Of Gene Expression (GO:0010468), Regulation Of Phagocytosis (GO:0050764). In line with expectations, terms associated with metabolism and inflammatory response were enriched both in broad or tissue-selective clusters, including: Glucose Homeostasis (GO:0042593) in clusters O (Pdk4, Lrrc8a) and L (Obp2a, Pck1), Regulation Of Glucose Metabolic Process in cluster K (GO:0010906; Slc45a3, Gnmt), Acute Inflammatory Response in cluster L (GO:0002526; Crp, Itih4, Hp), Cellular Response To Cytokine Stimulus in cluster D (GO:0071345; Irf1, Il1b, Ccl4), “Steroid Biosynthetic Process in cluster L (GO:0006694; Hsd3b2, Cyp8b1, Slc27a5). An example of tissue-specific functional cell reprogramming is the presence of genes connected to “Fat Cell Differentiation” (GO:0045444) and “Regulation of Cold-Induced Thermogenesis” (GO:0120151) in cluster C, specifically downregulated by DEX in the adipose tissue. The results indicate that GR through activation of tissue-specific molecular programs is involved in regulation of biological pathways in multiple organs and systems.

Overlap between observed and previously reported glucocorticoid-induced gene transcription patterns and published drug perturbation datasets

To further investigate the GR-dependent gene expression patterns, we performed analyses of functional overlap between the clustered genes and a set of external gene lists. The first set of gene lists was obtained from transcriptional studies involving GR stimulation. The dataset was compiled for this research and included 91 gene lists from 37 studies, which investigated GR-dependent transcriptional changes in various tissues and cell types. In total, 45 lists from 26 studies were found to have a statistically significant overlap with at least one of the gene clusters identified in this study. Out of these, 35 lists significantly overlapped with the full list of 585 DEX-induced transcripts. Genes included in clusters O and P overlapped respectively with 19 and 23 lists of genes obtained from the literature (Fig. 4, Supplementary Table 2), confirming that those clusters contain bona fide GR targets regulated ubiquitously across tissues.

Fig. 4

Heatmap of gene overlaps between the GR-dependent transcriptional patterns and literature-based GR-dependent gene lists for various tissues and cells. Lists of genes for tissues and cell types were extracted from transcriptomic studies of GR-dependent gene expression (PMIDs are provided in the rightmost column). The chi-square test was used to assess the overlap between the clusters of GR-regulated genes (x-axis) and literature (y-axis). χ² values were transformed in the following way: log₂(χ² + 1) and are displayed in each rectangle of the heatmap (white blocks—no overlap; red blocks—higher overlap). Values in brackets within each square indicate the number of overlapping genes. Statistically significant results are indicated as colored frames: p_adj < 0.01 (green frames); p_adj < 0.0001 (purple frames). The gene clusters are labeled at the top of the columns. For each cluster, the number of genes and the observed direction of the transcriptional regulation (‘ + ’ increase, ‘ − ’ decrease) are indicated. Cluster UP and cluster DOWN represent combined groups of all DEX-induced and DEX-decreased transcripts, respectively. The row description on the right of the heatmap provides details on the gene lists obtained from the literature ('+'represents upregulation and'−'represents downregulation). The gene lists presented on the heatmap are organized according to hierarchical clustering, indicated by both row and column dendrograms. An extended table summarizing gene list overlaps is included in Supplementary Table 2

Next, we exploited our resource of DEX-regulated clusters to inspect drug gene-expression signatures, based on published studies available through Enrichr—Drug Perturbations from GEO and DSigDB [35]. Using the Drug Perturbations dataset, we identified drugs and conditions overlapping with particular clusters. Similar to the manually curated DEX literature database, in the Enrichr performed analysis showed that DEX regulation was found (at adjusted p < 0.1) predominantly in broadly regulated clusters O (GSE484 mouse lung; GSE7683 mouse chondrocytes; GSE34313 human airway smooth muscle; GSE44208 mouse skeletal muscle; GSE37474 human trabecular meshwork), P (GSE484 mouse lung; GSE34313 human airway smooth muscle; GSE34313 mouse placenta; GSE37474 human trabecular meshwork; GSE2342 mouse B-cells; GSE7683 mouse chondrocytes) and K (GSE29912 rat liver; GSE484 mouse lung; GSE4165 mouse placenta). Additionally, DEX was previously found to regulate genes from clusters B (GSE7683 mouse chondrocytes), and D (GSE37474 human trabecular meshwork). The search also included other steroid drugs, which were previously found to regulate genes included in several of the clusters, namely hydrocortisone (GR/MR agonist, GDS3071, cluster P), methylprednisolone (GR agonist, GDS964, clusters D, I, K, L, O, P), estradiol (Estrogen Receptor agonist, GSE11567, clusters B, D, I, K, L, P), mifepristone (GR/Progesterone Receptor antagonist, GSE39270, cluster G). Other metabolic drugs were also found as regulators of particular clusters: troglitazone (clusters I, K, O), glipizide (clusters I, L), cholecalciferol (cluster I), and probucol (clusters B, K, L, P).

To further elucidate the potential environmental and pharmacological factors influencing the observed gene expression patterns, we analyzed the association between our gene clusters and known chemical perturbations using the Comparative Toxicogenomics Database (CTD) implemented in Enrichr. We identified the top overrepresented chemical perturbations across various gene clusters, focusing on those affecting multiple clusters (see Supplementary Table 5). Notably, aflatoxin B1 (CTD:00007128) and calcitriol (CTD:00005558) emerged as the most widely influential, each associated with 10 different clusters including both up- and downregulated genes. Other significant perturbations included anisomycin (associated with 9 clusters), benzo[a]pyrene (CTD:00005488, 9 clusters), estradiol (CTD:00005920, 9 clusters), and tetradioxin (CTD:00006848, 9 clusters).

This analysis revealed several important insights: First, substantial overlap between identified broad clusters (particularly O and P) with existing literature-derived gene lists provides robust validation that these clusters contain genuine GR target genes with cross-tissue relevance. Second, the fact that 45 lists from 26 studies showed significant overlap demonstrates the biological reproducibility of these findings across different experimental contexts. Third, beyond DEX, the analysis revealed significant overlap with other steroid drugs (hydrocortisone, methylprednisolone, estradiol, mifepristone), suggesting common regulatory pathways. Fourth, The overlap with patterns associated with metabolic drugs (troglitazone, glipizide, cholecalciferol, probucol) points to potential crosstalk between glucocorticoid signaling and metabolic pathways.

GR-dependent transcriptional signatures show specific associations with human health-related traits

We explored the associations between GR-dependent transcriptional patterns and gene lists associated with human phenotypic traits. Gene-phenotype associations were obtained from the Pan-UK Biobank study [36]. Among the 729 phenotypes examined, 11 were identified as showing a statistically significant overlap (p_adj < 0.01). Detailed examination revealed that cluster I particularly overlapped with phenotypic traits. This cluster, containing transcripts upregulated predominantly in the kidney and the pituitary gland, was identified as potentially connected to four phenotypes: three associated with reticulocyte count, and one linked to blood triglyceride levels. Additionally, cluster D turned out to contained transcripts related to varicose veins and alanine aminotransferase levels. Another finding was the identification of genes connected to systolic blood pressure among transcripts in cluster O. Notably, oppositely-DEX-regulated gene clusters showed associations with distinct sets of health-related phenotypes (Fig. 5A, Supplementary Table 6).

Fig. 5

GR-dependent transcriptional patterns overlap with groups of genes associated with (A) human phenotype- or (B) metabolism-related traits. A The heatmap on the left presents an overlap between the GR-regulated genes and human phenotype-associated gene lists. Lists of genes associated with phenotypes were extracted from the Pan-UK Biobank database (phenocodes are provided). B The right panel presents the overlap between the clusters of GR-dependent genes and gene lists associated with metabolite levels (based on Metabolon and Nightingale databases). The chi-square test was used to examine the overlap of the lists of the GR-regulated (columns) and phenotype- or metabolism-associated genes (rows). χ² values were transformed using the formula: log₂(χ² + 1) and are presented in each rectangle of the heatmap. The intensity of the red color is proportional to the level of overlap, as indicated in the legend (white, no overlap; red, high-level overlap). The statistical significance of the overlap is pinpointed using two thresholds: p_adj < 0.01 (green frames) and p_adj < 0.0001 (purple frames). Values in brackets in each box show the number of overlapping genes. The results with a minimum of three genes were included. The identified gene clusters are indicated on top column labels with the number of genes inside brackets. Cluster UP and cluster DOWN represent combined groups of DEX-induced and DEX-decreased transcripts, respectively. For each cluster, the number of genes and the observed direction of the transcriptional regulation (‘ + ’ increase, ‘ − ’ decrease) are indicated. The row description presented on the right of the heatmaps provides details on lists of human phenotypes and metabolomes from genome-wide association studies (GWAS). The gene lists presented on the heatmaps are organized according to hierarchical clustering, indicated by both row and column dendrograms. To reduce redundancy, for cases where the overlap contained the same genes for more than one phenotype, only the phenotype with the lowest p_adj is displayed. Full results are available in Supplementary Table 6

Subsequently, we investigated the relationship between GR-dependent transcripts and genes associated with human metabolite levels. This analysis was performed on 654 lists of genes extracted from human metabolome GWAS results [37]. We found statistically significant associations with lists of genes connected to six metabolic traits. Gene lists connected to the four traits overlapped with the group of genes induced by DEX (Fig. 5B, cluster UP). Transcripts induced by DEX in the kidney clusters I and L overlapped with four and two traits, respectively. Transcripts included in cluster L were associated with LDL and HDL diameter, whereas cluster I overlapped with gene lists related to triglyceride levels and HDL diameter (Fig. 5B, Supplementary Table 6). These data uncovered several examples of complex functional interactions between GR-regulated genes and physiological parameters of the organism.

Discussion

We have comprehensively analyzed DEX-induced gene expression patterns across multiple tissues within the context of the functional annotations and human phenotype associations. Globally, the magnitude of the effect on various tissues is similar, with one exception—the hypothalamus, where the impact of DEX is less prominent than in other tissues. This observation may be explained by limited DEX penetration through the blood–brain barrier [38].

Examination of the transcriptional effects of systemic GR activation revealed complex mechanisms behind tissue-specific and non-tissue-specific responses to DEX. Notably, four upregulated (P, K, O, and I) and one downregulated (D) cluster seemingly represent a set of core GR-responsive genes that may play foundational roles in coordinated glucocorticoid action across tissues. Accordingly, these clusters show significant enrichment for direct GR binding sites, unlike the clusters that define tissue-specific responses. The shared transcripts include well-established GR-responsive genes like Fkbp5, Tsc22d3, and Dusp1 [25, 39,40,41]. Most of these core genes were found previously to be regulated by glucocorticoids in specific contexts [25, 42,43,44,45].

A key strength of this study is the parallel profiling of GR responses across multiple tissues in vivo under identical experimental conditions. While numerous studies have investigated GR effects in specific tissues [25], or even multiple cell types [23], and meta-analyses attempt to synthesize these [17, 18], our approach minimizes confounding variables (such as differing treatments, time points, genetic backgrounds, or platforms) that complicate cross-study comparisons. This allows for a more direct and reliable assessment of both shared systemic responses and genuine tissue-specific transcriptional signatures orchestrated by GR activation, providing a unique systemic perspective.

Strikingly, more than two-thirds of DEX-induced changes were controlled in a tissue-selective manner (significant interaction between treatment and tissue factors), in agreement with the previously postulated tissue specificity of glucocorticoid response [46, 47]. This heterogeneity suggests that rather than pursuing a single, unified biological mechanism, glucocorticoids orchestrate a diverse array of tissue-specific responses that together contribute to systemic homeostasis.

An intriguing finding was the lack of overrepresentation of GR binding sites in clusters showing high regional selectivity (A, B, C, E, F, G, H, J, L, M, N). This pattern suggests that DEX-induced changes in these clusters may be mediated through indirect transcriptional mechanisms. Indeed, our analysis revealed enrichment of binding sites for other transcription factors in specific clusters. For instance, LXR and RXR, known GR interactors [48], were prominently enriched in clusters H, I, and L. CLOCK, which is a known regulator of GR activity [49], showed significant enrichment in clusters K and P. The enrichment of PPARA, which tends to occupy the same sites on chromatin as GR [50, 51], was observed in cluster I and EZH2 enrichment was observed in clusters F and O. This data point to an important crosstalk between glucocorticoid signaling and other regulatory pathways involved in metabolism and epigenetic regulation, respectively [52]. The SOX2 enrichment, a known inhibitor of GR binding to chromatin [53], in clusters D and M, and TP53 in clusters N and O, further highlights the diversity of regulatory networks potentially influencing selected aspects of GR activity. We observed strong concordance between tissue-specific enrichment patterns and mRNA expression levels for several TFs identified in specific clusters. For instance, Nr1h3 (encoding LXRα) showed high expression in the liver, adipose tissue, and spleen, aligning with the enrichment of LXR sites in clusters H, I, and L which feature upregulated genes in these tissues such as LIV, FAT, SPL, KID, PIT, ADR, HTH, and MUS. Similarly, Ppara expression was highest in the liver, kidney, and adipose tissue, consistent with its enrichment in cluster I (upregulated mainly in the kidney and pituitary). Conversely, Sox2 showed prominent expression in the pituitary and hypothalamus, correlating with its enrichment in clusters D and M where tissue-specific effects were noted. While mRNA levels do not perfectly equate to protein activity, these correlations suggest that the differential availability of specific transcription factors across tissues contributes to the tissue-specific transcriptional outcomes following GR stimulation. These findings show that the tissue-specific effects of glucocorticoids are achieved through a complex interplay of direct GR-mediated transcription and indirect regulation via other transcription factors and signaling pathways [53, 54]. Our data highlight the key role of recruitment of tissue-specific transactivators or transrepressors in shaping the GR-induced transcriptional reprogramming at the single-organ resolution.

There are several consequences of the observed tissue-selectivity of GR responsiveness for the biology of the system. First, the multidirectional nature of these effects, namely simultaneous upregulation of gluconeogenesis genes (cluster K), lipid metabolism (cluster I) and the modulation of inflammatory response genes (cluster D) may reflect an body's attempt to increase the availability of energy resources while simultaneously suppressing excessive inflammation. This pattern is highly relevant for increased periods of physiological activity, such as the beginning of the active phase (peak of plasma glucocorticoid concentration) and is reminiscent of the acute stress response. In the context of chronic stress or prolonged glucocorticoid exposure, however, these adaptive responses may become maladaptive. For instance, the persistent upregulation of genes involved in glucose metabolism could contribute to the development of insulin resistance and type 2 diabetes, while chronic suppression of immune-related genes might increase susceptibility to infections [55].

The data prompt us to also reconsider our current model of GR action, particularly in metabolic regulation. For example, while the GR-dependent regulation of glucose metabolism is thought to be uniform across tissues, our data demonstrate highly tissue-specific patterns of key metabolic gene regulation by DEX (Igfbp1 in cluster H, Slc45a3 and Gnmt in cluster K, Pck1 in cluster L, Pdk4 and Lrrc8a in cluster O) (see [56,57,58]). Regulation of Eno2 [59], central to the glycolytic pathway, and Pck1 [60], a key enzyme in gluconeogenesis, exemplifies the direct modulation of glucose production and utilization. Meanwhile, Igfbp1 [61], which influences insulin sensitivity, and Insig2 [62], involved in lipid metabolism, showcase the indirect routes through which DEX affects glucose homeostasis.

The consequences of GR regulation for human disease were further revealed through the overlap analysis with data from the Pan-UK Biobank [32]. Notably, the discovery that cluster (I) of genes upregulated predominantly in the kidney and pituitary gland, correlates with phenotypes related to blood triglyceride levels, suggests a potential role for these tissues in mediating the systemic effects of glucocorticoids on lipid metabolism. The genes involved in this action are Angptl4, Apoa1, Apoc1, Cdkal1, Ppp2r3a, Snx10, and Tsku [63]. The analysis exploiting drug perturbation datasets and metabolite-associated gene lists corroborates this notion. The association with metabolic drug signatures (such as troglitazone, glipizide) may indicate the molecular basis for glucocorticoids’ diabetogenic effects. DEX-induced gene clusters showing connections to human metabolic traits, particularly lipid metabolism (cluster I with triglyceride levels and HDL diameter and cluster L—LDL and HDL diameter) provide insights into how glucocorticoids influence systemic lipid profiles. The associations with physiological traits appear to be mediated through tissue-specific interactions with other transcription factors, such as PPARs, rather than through direct GR homodimer activation. Hence, our findings highlight the complex effects of glucocorticoids on metabolism, as evidenced by the diverse regulation profiles of key genes such as Pck1 in gluconeogenesis, Insig2 in lipid metabolism, and Apoc4 in triglyceride levels. These findings also suggest potential biomarkers and targets for managing metabolic side effects of glucocorticoid therapy (such as Apoc4, Asgr1).

The observed associations highlight potential mechanisms by which chronic stress might contribute to metabolic syndrome and cardiovascular disease. The link between systolic blood pressure and genes from cluster O suggests a potential contribution of glucocorticoids in the etiology of cardiovascular disorders. Significant overlaps with other steroid drugs (such as hydrocortisone, methylprednisolone, estradiol) indicate that shared molecular pathways of steroid hormones might be exploited to counteract these effects [23, 64]. Overlaps with genes responsive to various chemical perturbations (such as aflatoxin B1 and benzo[a]pyrene) in multiple tissues may also suggest activation of antitoxicological mechanisms by stress. Thus, our findings provide a molecular framework for understanding how the initially adaptive effects of glucocorticoids under acute stress can, may lead to the development of stress-related pathologies, emphasizing the delicate balance that glucocorticoid signaling must maintain for optimal physiological function.

Several limitations of this study should be acknowledged. The transcriptomic data were generated using whole-genome microarrays. While providing a broad overview, this technology has inherent limitations compared to RNA-sequencing, including a narrower dynamic range, reduced sensitivity for detecting low-abundance transcripts, and an inability to identify novel transcripts or splice variants not covered by the array probes. Consequently, some DEX-regulated transcripts, particularly those expressed at very low levels or representing novel isoforms, may not have been captured in our analysis. Furthermore, the study focused on a single time point post-DEX administration (4 h), capturing a snapshot of transcriptional changes but not the full dynamic response profile. Finally, only male mice were used to control for sex-based variability; future studies should include female mice to investigate potential sex differences in GR-mediated transcriptional responses across tissues.

In conclusion, our comprehensive analysis bridges the intricate molecular mechanisms of GR-mediated transcriptional regulation with clinically relevant phenotypes. Moreover, we demonstrate that the dichotomy between direct GR-dependent regulation of shared responses and indirect regulation of tissue-specific effects represents a fundamental principle of glucocorticoid action. This new mechanistic understanding of how a single nuclear receptor can orchestrate diverse tissue-specific responses through differential engagement with local transcriptional networks provides a framework for developing more targeted therapeutic approaches that could potentially harness tissue-specific regulatory mechanisms while minimizing systemic effects. Future research could build upon these findings to further elucidate the complex interplay between GR and tissue-specific transcriptional machinery, potentially leading to more selective glucocorticoid-based therapies.