Abstract
Introduction
Despite their potent anti-inflammatory effects,1–3 corticosteroids (CSs) are associated with numerous side-effects, particularly related to long-term treatment, eg hyperglycemia, dyslipidemia, arteriosclerosis, muscle wasting, and osteoporosis.4–7 These complex side-effects, which manifest as systematic transcriptional changes in multiple tissues,8,9 necessitate approaches to better decipher the pharmacogenomic effects of CS to properly assess the balance between their advantageous and harmful effects.
There has been an increasing number of recent studies that explored the tissue-specificity of gene expression and regulation,10,11 especially its interplay with pathology.12–14 Several databases were developed to establish a knowledgebase of large-scale data regarding tissue-specific gene expression, regulation, and disease association in a variety of human tissues.15,16 Furthermore, tissue-specificity has also been shown to link with many significant outcomes including expression-quantitative trait loci, 17 evolution,14,18 and disease-association,13,19 eg tissue-specific effects of insulin signaling in diabetes. 20 Consequently, in this study we placed emphasis on exploring the tissue-specificity of chronic CS administration in liver and muscle to obtain a better understanding of CS-responsive functions and its side-effects. We first identified tissue-specific transcriptional dynamics using a novel clustering approach proposed in our previous study. 21 Subsequently, with the hypothesis that common functions activated across multiple tissues will be more likely to be important CS-responsive functions, the KEGG database was utilized to identify such common functions.
However, it has been noted that genes affected by CS include both immunosuppressive genes, mostly associated with therapeutic effects, and metabolic genes often associated with adverse effects whose regulation is mainly controlled by glucocorticoid receptor (GR) through gene-mediated pathways. 6 However, chronic infusion of CS causes a sustained down-regulation of the receptor (mRNA and thus protein).22,23 Although several alternative mechanisms have been proposed,24–26 it is still not understood why drug effects remain strong even after the down-regulation of GR mRNA to the point of almost being eliminated. A plausible explanation is that, in addition to direct regulation through glucocorticoid response element (GRE) binding sites in the 5′ regulatory regions of genes, there are changes in expression that are also the indirect results of CS effects because of changes in expression of transcription factors (TFs) that act as secondary bio-signals directly or indirectly modulating the transcription of genes.23,27,28 Therefore, in addition to identifying critical transcriptional responses and CS-responsive functions, finding relevant TFs and understanding their relationship with CS-responsive functions is also an important aspect in this analysis.
Materials and Methods
Datasets
Adrenalectomized (ADX) male Wistar rats with body weights of 339 ± 28 g (SD) were housed and maintained under constant temperature (22 °C) and humidity with a controlled 12-hour light/dark cycle for a period of at least two weeks before surgery.22,29 Rats had free access to rat chow and 0.9% NaCl drinking water. Forty rats were given 0.3 mg/kg/hour infusions of methylprednisolone (MPL) sodium succinate for more than 168 hours via an Azlet pump. The drug solutions were prepared for each rat based on its predose body weight. Pumps were subcutaneously implanted between the shoulder blades. Animals were sacrificed at various times up to seven days, specifically at time-points 6, 10, 13, 18, 24, 36, 48, 72, 96, and 168 hours. A control group of four animals was implanted with a saline-filled pump and killed at various times throughout the seven-day study period. Liver and gastrocnemius muscle samples from each animal were ground into a fine powder in a mortar cooled by liquid nitrogen, and RNA was extracted with TRIzol reagent (Invitrogen, Carlsbad, CA, USA). The gene expression was obtained via the Affymetrix RAE230A array, which consists of 15,923 probesets. The data are publicly available through the GEO database under accession numbers GDS972 and GDS2688 corresponding to the liver and muscle datasets during chronic MPL infusion.
Clustering
Utilizing the concept of the agreement matrix (AM) in consensus clustering, we recently proposed a novel method to identify the core group of probesets showing the most agreement indicating they belong to the same patterns of gene expression.
21
To produce the AM, a number of different clustering methods (hierarchical clustering, divisive analysis clustering, partitioning around medoids, fuzzy analysis clustering,
CS-responsive function identification
We first identified the pathways from the KEGG database that were enriched in each of the gene expression patterns from each tissue (
Promoter extraction and processing
The promoters of genes including all transcript-relevant alternative promoters were extracted from a rich database of promoter information with a default length (500 bp upstream and 100 bp downstream of the transcription start site) if there is no experimentally defined length suggested by Genomatix. 32 To accelerate the process of identifying putative transcriptional regulators, promoters were pre-processed as in Nguyen et al. 33 Specifically, MatInspector 34 was applied to scan position weight matrix (PWM) matches on those promoter sequences using optimal parameters from MatBase, which ensures that the minimum number of matches found in non-regulatory sequences, ie, the false-positive matches, was minimized. 32 Each gene promoter was then re-modeled to become a list of TF binding sites (TFBSs) ordered by their local positions on the promoter sequences and represented by the corresponding TF names along with their binding orientations. The conversion supports fast search for the presence of a TFBS or a m-regulatory module (CRM) on promoter sequences.
Putative transcriptional regulator prediction
To predict putative transcriptional regulators, we utilized the context-specific CRM search technique to identify over-represented CRMs in the promoter set of a gene battery. 33 Each gene battery contains a certain set of genes that are hypothetically co-regulated, ie co-expressed and co-functional in this study. With the hypothesis that common functions activated across multiple tissues may play important roles in response to CS dosing, we applied our previous tool 33 to identify TFs relevant to CS transcriptional responses. In brief, we computationally defined a CRM as a list of non-overlapping TFBSs ordered by their positions on the promoter sequence and characterized with their corresponding binding strand orientation. The procedure first identified all potential TFBSs that are common in the corresponding promoter set and then searched for all possible combinations of all commonly found TFBSs above using the breadth-first search technique. Owing to the fact that a CRM can be present on promoters of many genes in the background set, we estimated the statistical significance of commonly identified CRMs for each gene battery vs. those for the background set to select those that are significantly over-represented. Subsequently, the selected CRMs are decomposed to obtain a list of TFs that are associated with the corresponding TFBSs present on CRMs.
Regulatory closeness (RC) definition
Using the hypothesis that with the more relevant transcriptional regulators, two gene batteries will share higher closeness of regulatory mechanisms, the RC was estimated to quantify the relationship between gene expression profiles, function, and transcriptional regulation. In this study, we defined the RC as the common ratio (CR) of relevant TFs shared between two gene batteries CR =
Results
This study involved infusion of MPL for more than 168 hours in ADX rats. Three animals were sacrificed at each of 12 time-points, and diverse tissues harvested and processed to assess gene expression using Affymetrix arrays. The liver and muscle tissues were compared in this analysis (see Materials and Methods). The pharmacokinetics of MPL for this study was reported previously 35 showing steady-state concentrations in plasma at six hours and thereafter.
Critical transcriptional responses
To analyze gene expression data during chronic MPL infusion, a pre-processing step was performed by filtering differentially expressed probesets using an ANOVA test (
For responses in liver, we obtained 12 significant patterns with 2,285 genes in total. The expression patterns as

Critical dynamic transcriptional responses within individual tissues under chronic corticosteroid administration. Each pattern is characterized by the average gene expression profile of the corresponding cluster of genes in the liver dataset. The error bar shows the standard deviation of all probeset transcript levels at each time-point in each corresponding pattern.
In muscle responses, there are seven significant expression patterns with 2,291 genes in total. However, almost all genes follow two main transcriptional responses as characterized by patterns 1 and 2 in Figure 2. Pattern 1 including 718 genes exhibited an early down-regulation during the first 6 hours followed by a robust up-regulation to reach the peak at about 36 hours with expression levels maintained until 96 hours before returning to the baseline. In contrast, 1,012 genes in cluster 2 expressed an opposite pattern that shows a transient up-regulation in expression at six hours during MLP infusion followed by a strong decline reaching the nadir at about 36 hours. These genes show sustained repression until about 96 hours and then gradually return to the baseline. Genes in cluster 3 (261 genes) exhibited a similar pattern, but there is no early transient up-regulation in the first 6 hours, and the nadir is reached at about 18 hours. Pattern 4 (76 genes) showed a robust, fast down-regulation during the first 10 hours reaching a new steady state without returning to the baseline. Pattern 5 (63 genes) characterizes an initial up-regulation to the peak at about 10 hours followed by a down-regulation to reach the nadir at 96 hours. Patterns 6 (65 genes) and 7 (96 genes) exhibited a simple up-regulation followed by a down-regulation to the baseline. However, genes in pattern 6 reach the peak faster at about 18 hours, which was maintained until 96 hours, whereas genes in pattern 7 reach their peak at about 96 hours and then simply return to the baseline (see Supplemental Material 1).

Critical dynamic transcriptional responses in muscle under chronic corticosteroid administration. Each pattern is characterized by the average gene expression profile of the corresponding cluster of genes in the muscle dataset. The error bar shows the standard deviation of all probeset transcript levels at each time-point in each corresponding pattern.
CS-responsive genes and functions
To identify specific genes and functions caused by drug dosing, we assume that co-expressed genes involved in pathways activated by MPL are more likely to be important genes and functions. Using the KEGG database and ArrayTrack,
31
we searched for enriched pathways in these co-expressed clusters (
Tissue-specific regulation by functional characterization.
Detailed information of selected gene batteries in liver and muscle.
A number of pathways are significantly enriched in both tissues. Because these co-expressed genes reflect critical transcriptional responses during MPL infusion, we hypothesize that common functions activated across multiple tissues are important CS-responsive functions. Furthermore, along with the concept of “gene battery,”38,39 we define gene sets associated with each function in each corresponding cluster as gene batteries. Table 2 shows 18 common functions activated in both tissues during MPL infusion, which include 23 gene batteries in liver and 19 in muscle. These functions mainly belong to four pathway categories including signaling pathways (insulin signaling, MAPK signaling, neurotrophin signaling, PPAR signaling, and Wnt signaling), metabolic pathways (oxidative phosphorylation, pentose phosphate pathway, and purine metabolism), and pathways relevant to cellular process (adherens junction, focal adhesion, tight junction, lysosome, and cell cycle relevant processes) and information processing (aminoacyl-tRNA biosynthesis, proteasome, spliceosome, ubiquitin-mediated proteolysis, and ECM–-receptor interaction). Additionally, it has been noted that expression levels of many CS-affected genes are mediated through the binding motifs located on their control regions, called GREs. We thus examined the presence of this binding site on the promoter of genes in each of the enriched pathways to assess the potential effect of GRE on common functions. Following the approach in our previous study, 30 genes shown in bold are those whose promoters contain GRE motifs located on conserved regions identified from the corresponding sets of orthologous promoters (details in Supplemental Material 2). The average ratio between genes containing GRE motifs vs. genes without GRE motifs in liver gene batteries is 0.35 where 7 out of 23 gene batteries have ratios ≥0.5. The average ratio in muscle is 0.47 where there are 10 out of 19 gene batteries with ratios ≥0.5. Also, we find that many CS-responsive genes in signaling cascades and information processing pathways (eg, M APK signaling, Wnt signaling, aminoacyl-tRNA biosynthesis, and proteasome) are more likely to be directly regulated by the CSs and GR complex in both tissues. Some other functions, eg insulin signaling, PPAR signaling, focal adhesion, and cell cycle relevant processes, contain a greater number of genes directly regulated by the CS complex in muscle (Table 2).
Potential transcriptional regulators of CS-responsive common functions
During chronic infusion, MPL concentrations reach and remain at a stable steady state after six hours. 35 The drug binds to cytosolic GRs and then rapidly translocates into the nucleus to alter the expression of target-genes. As GRs are greatly diminished in response to CSs,22,23,27,28 mRNA levels of CS target-genes should quickly return to the baseline. However, as observed in Table 2 many gene data sets exhibit a long-term response before returning to baseline or even reach a new steady state without returning. In addition to a number of possibilities (eg, multiple GR isoforms, multiple GREs with different affinities to the drug–-receptor complex, multiple receptors that can mediate the effect of CSs),35,40 we hypothesize that these effects can be caused by the regulation of secondary bio-signals where TFs are the most likely candidates. After activation by CSs, they in turn further modulate the expression of their target-genes as a continuous cascade of events that were initiated by the drug.
Consequently, with the hypothesis that common functions activated across multiple tissues are more likely to be important as CS-responsive, the gene batteries selected in Table 2 were further analyzed to predict putative transcriptional regulators that are potential secondary bio-signals relevant to the regulation of transcriptional responses. For 18 common functions expressed across multiple patterns, we found 23 critical gene batteries in liver responses and 19 in muscle responses. Using the context-specific CRM search technique, a set of TFs associated with TFBSs in CRMs that are statistically significantly present on the corresponding promoter set was extracted for each gene battery (see Materials and Methods). Frequent TFs present in more than 20% of all gene batteries in each tissue are shown in bold (Tables 3 and 4). Specifically, we have 22 frequent TFs relevant to the regulation of transcriptional responses in liver and 20 in muscle. Among them, 17 relevant TFs are common including CLOX, CREB, CTCF, E2FF, EGRF, ETSF, FKHD, HOMF, HOXF, NKXH, NR2F, OCT1, RXRF, SORY, SP1F, STAT, and ZBPF. Almost all TF families consist of TF members that are recognized as differentially expressed genes in one or both tissues (see Supplemental Material 2). This finding highlights the possibility that secondary bio-signals are involved in the regulatory complexities of expression changes for CS-affected genes.
Functions and transcriptional regulators of gene batteries in liver.
Functions and transcriptional regulators of gene batteries in muscle.
Putative functional regulatory networks
As biological processes do not happen in isolation, we defined a hypothetical quantity, called the “regulatory closeness,” which putatively reflects the similarity of transcriptional regulatory mechanisms between two gene batteries to characterize the relationship among expression, function, and transcriptional regulation (see Materials and Methods). We graphed these relationships in a so-called functional regulatory network with nodes representing gene batteries and edges high transcriptional regulatory similarities (RC ≥0.8 in this case). Figure 3 exhibits corresponding functional regulatory networks in liver (top) and in muscle (bottom). Gene batteries having no relationship of high RC with any other gene battery are not displayed in the graph. The edges present on both graphs are displayed in “red” color. As a result, we have 19 gene batteries with 25 edges for the functional regulatory network in liver and 19 gene batteries with 38 edges in muscle. The networks show that gene batteries with similar expression patterns are better connected to each other than batteries with different expression patterns. This is the main reason why muscle that has the same number of gene batteries as liver contains more connections. Moreover, under the same condition and for the same function, involved genes can exhibit different expression patterns in liver than muscle. Also, even in the same tissue, different sets of genes with the same function can have different expression patterns (eg, gene battery “tight junction” in liver). Additionally, gene expression patterns in liver are more diverse than those in muscle, which may result in the potential capability of sharing common regulatory mechanisms among gene batteries, which is higher in muscle. These results provide an overview landscape of tissue-specific expression and regulation under chronic MPL infusion through critically identified transcriptional responses and CS-affected functions as well as their transcriptional regulatory relationships.

Tissue-specific regulation represented by putative functional regulatory networks (top: liver; bottom: muscle). Each node represents a gene battery which is a set of coexpressed genes sharing a common function (pathway). Edges characterize the regulatory closeness between two gene batteries if the ratio of sharing common transcriptional regulators is greater than 0.8; ‘red’ edges express high regulatory closeness that occurs in both liver and muscle. Although these gene batteries are involved in similar functions under chronic corticosteroid administration, their expression patterns and transcriptional regulatory relationships are specific in liver vs. in muscle.
Discussion
Although CSs affect gene expression in multiple tissues, the array of genes that are regulated by these steroids is diverse, highly tissue specific, and depends on their functions in the tissues. Liver is an organ that plays a critical role in performing and regulating important physiological processes such as energy metabolism, detoxification, inflammation, and immune responses against pathogens. Muscle plays a key role in maintaining systemic energy homeostasis and accounts for 80% of insulin-directed glucose disposal. In addition, both liver and muscle are important sites for immune reactions caused by different immune cell types. Therefore, we explored tissue-specific CS-responsive gene expression profiles, their functional significance, and the effect of secondary bio-signals in the form of TFs to gain better insight into CS pharmacogenomics.
In this study, 2,285 differentially expressed genes in liver during chronic MPL infusion were parsed into 12 distinct temporal clusters. Functional analysis of the pathways in which these clustered genes are involved suggests that each cluster is enriched with specific biological functions. For example, most of the genes from cluster 6 are present in pathways involved in performing or regulating energy metabolism (oxidative phosphorylation, pentose phosphate pathway, and starch and sucrose metabolism). Similarly, genes from cluster 7 are primarily involved in pathways regulating immune response and inflammation (complement and coagulation cascades, leukocyte transendothelial migration, antigen processing and presentation) or associated processes (cell adhesion, ECM–-receptor interaction, focal adhesion). It is interesting to see that the expression of genes in cluster 6 is up-regulated over time with continuous MPL infusion reaching a plateau. In contrast, expression of genes in cluster 7 is continuously down-regulated, ultimately reaching a plateau. These observations go well with the fact that CSs repress inflammation and immune response cascades that lead to their therapeutic applications while steroid effects on metabolic processes lead to their adverse effects. The CSs, being potent anti-inflammatory agents, not only repress inflammatory signaling but also affect immune cell trafficking. This is confirmed by the continuous down-regulation of genes involved in leukocyte migration and associated processes including cell adhesion. Another interesting observation is the development of tolerance, which is a commonly observed phenomenon in chronic CS treatment because of the down-regulation in the expression of GRs that mediate their genomic effects. For example, in liver gene expression, clusters 1, 2, 3, 4, and 5 show clear tolerance, which may occur because of the down-regulation in receptor expression following CS dosing.
In muscle responses during CS infusion, 2,291 differentially expressed genes were parsed into seven distinct temporal clusters. Similar to liver, genes present in pathways regulating immune response (B-cell receptor signaling pathway, Fc gamma R-mediated phagocytosis, chemokine signaling pathway) and related processes (ECM–-receptor interaction, focal adhesion) are down-regulated in muscle during MPL infusion (clusters 2 and 3), although the specific immune regulating pathways are different between these two tissues. In the case of energy metabolism in muscle, genes involved in oxidative phosphorylation show a robust up-regulation followed by the development of tolerance, genes involved in fatty acid metabolism show a transient down-regulation followed by a robust up-regulation (cluster 1), and genes involved in glycolysis, pyruvate metabolism, and pentose phosphate pathway show exactly the opposite expression profile with initial transient up-regulation followed by robust down-regulation (cluster 2). This observation is consistent with the fact that glucocorticoids direct the muscle to use free fatty acids as the primary energy source, producing NADPH which is used for the production of ATPs through oxidative phosphorylation. Furthermore, the observation that genes involved in glycolysis and pyruvate metabolism are down-regulated confirms that GC decrease the dependence of muscle on glucose, which helps in the maintenance of proper plasma glucose concentrations. It is interesting to see that, even though the genes involved in oxidative phosphorylation are up-regulated both in liver and muscle at some point in time after MPL dosing, the expression profiles are very different between the two tissues with the liver showing three different sets of expression patterns and muscle showing the development of tolerance by 10 hours after the start of MPL infusion. Furthermore, expression of genes involved in ubiquitin-mediated proteolysis falls under cluster 1, where the genes show a robust up-regulation in muscle suggesting that the proteasomal machinery is activated by MPL. This is one of the important causes of clinical side-effects of CS, which is muscle wasting caused by the breakdown of muscle proteins to provide the amino acid carbon for gluconeogenesis in liver.
Apart from the above-mentioned processes that are affected by CSs in both liver and muscle, CSs have tissue-specific effects on certain processes and pathways. For example, they cause up-regulation in the expression of ribosomal genes in liver, but have no effects on muscle ribosomal gene expression. In liver, clusters 1 and 5 consist of ribosomal genes that show an up-regulation of expression with MPL infusion with the peak expression occurring at 18 and 10 hours after the start of drug infusion and the development of tolerance thereafter. CSs are known to cause hypertrophy in liver, and up-regulation of expression ribosomal genes is a part of the process of increasing the cell mass, which is not observed in muscle. Similarly, genes involved in the metabolism of drugs and retinol show initial down-regulation followed by a rebound in expression in liver after the drug infusion, while the expression of these genes is unaffected in muscle. Xeno-biotic and retinol metabolism primarily occur in liver, and the muscle has no role in these processes; hence, the expression of these genes is not affected in muscle. However, the expression of genes present in pathways involved in adipocytokine signaling, calcium signaling, ErbB signaling, GnRH signaling, and phosphatidylinositol signaling are all affected by CS dosing in muscle, but the drug has no effects on these genes in liver, at least in part because of the higher relevance for these signaling pathways in muscle compared to liver.
Although the plasma concentrations of MPL attained steady state and remained almost constant from six hours after the start of the infusion, the GR mRNA expression and the receptor density are diminished because of the negative feedback regulation. This will cause a drastic reduction in the activated drug receptor complex, which is the primary driving force for the changes in gene expression. Hence, if the gene expression changes in these animals are caused only by the activated drug-receptor complex in the nucleus, then the expression profile for all the genes would go to a steady-state expression close to the baseline because of the marked tolerance caused by the drastic reduction in receptor density. However, from the expression profiles in both tissues, there are many genes that show long-term responses. It is well known that CS can regulate the expression of many TFs, and hence, we hypothesized that many of the long-term responses in the gene expression profiles are caused by other TFs that are affected and differentially expressed by CS. As shown in Tables 3 and 4, we were able to computationally identify the important TFs that could be involved in CS-responsive transcriptional profiles. Furthermore, from this analysis we could see that CSs regulate many of these TFs even at genomic levels (see Supplemental Material 2). For example, important regulators of inflammation and immune response (eg, nuclear factor kappa-B (NF-κB), interferon regulatory factor (IRF)) are found to be involved in pathogenesis of inflammatory disorders when CSs are used as therapeutic drugs. In this analysis, both are identified as potential TFs that are relevant to the CS-responsive transcription profiles, and their members are also down-regulated at the genomic level by CS in both liver and muscle. Similarly, many other TFs including STAT, CREB, and RXRF, which play important roles in energy metabolism and many signaling pathways, are also identified as potential TFs with distinct expression profiles during CS dosing. These results suggest that many TFs act as secondary bio-signals in controlling CS-responsive transcriptional behaviors, which in turn regulate the processes that are involved in both the therapeutic effects (caused by anti-inflammatory effects and immune suppression) and adverse effects (impact on metabolic processes).
Although similar pathways are affected by CS in both tissues, the genes and their expression patterns can be very different depending on significance of that pathway in each tissue. For example, as shown in Figure 4, the insulin signaling pathway is affected by MPL in both liver and muscle, but the individual genes that are CS-responsive and the expression patterns are very different between the two tissues. Among these, there is a tissue-specific coherent set of genes that are significantly coexpressed and regulate a certain set of different/similar downstream functions dependent on the tissue. Insulin signaling plays an important role in regulating many aspects of energy metabolism (protein synthesis and lipid and carbohydrate metabolism) and cell cycle processes (proliferation, differentiation, and apoptosis). Components of insulin signaling pathway that are involved in protein synthesis are responsive to CS only in liver, whereas genes involved in anti-lipolysis, proliferation, and differentiation are responsive to CS only in muscle. However, genes involved in glycogenesis, glycolysis, and lipogenesis are CS-responsive in both liver and muscle, although the expression patterns are very different between the two tissues. The mTOR and EIF-4EBP1 genes are continuously up-regulated only in liver with MPL infusion. They play an important role in activating liver protein synthesis, which is required to support the increase in liver mass caused by CS. As this phenomenon is not observed in muscle, these proteins may not be relevant and hence are not CS-responsive. Similarly, muscle is the major consumer of lipids as the energy source with high CS concentrations, and hence PRKAR1A, the gene involved in anti-lipolysis, is down-regulated only in muscle and is not affected in liver. It is also important to remember that TFs in the form of secondary bio-signals play an important role in regulating tissue-specific expression. As shown in Figure 4, all TFs involved in insulin signaling except for EBOX are different between the two tissues and hence result in regulating different genes involved in the same pathway and also cause differences in the expression profiles. This observation reinforces the tissue-specificity of transcriptional responses during MPL infusion (see more in Supplemental Material 3).

Tissue-specific expression within individual functions–-a case of insulin signaling pathway. The left subfigure shows an abstract of the insulin signaling pathway where expressed genes in liver and muscle and their corresponding downstream affected functions are included. The right panel displays the tissue-specific regulation of genes within the same functional category; putative TFs are those significantly overrepresented on the promoters of corresponding genes.
