Abstract
Corticosteroid (CS) pharmacogenomics was studied using gene microarrays in rat liver. Methylprednisolone (MPL) was administered intravenously at 50 mg/kg. Rats were sacrificed and liver excised at 17 time points over 72 h. RNAs from individual livers were used to query Affymetrix GeneChips that contain sequences for 8000 genes. Cluster analysis revealed six temporal patterns consisting of 197 CS-responsive probes representing 143 genes. Based on our fifth-generation model of CS pharmacokinetics/pharmacodynamics (PK/PD), mechanistic models were developed to describe the time pattern for each CS-responsive gene. Two clusters showed increased expression with different effect duration. PK/PD models assuming CS stimulation of mRNA synthesis were applied. Another two clusters showed an initial decline followed by delayed increase, suggesting two mechanisms might be involved jointly. The initial suppression was captured by CS inhibition of mRNA synthesis or stimulation of degradation. CS may also stimulate the production of a biosignal (transcription factors or other hormones), which can cause secondary induction of the target mRNA. One cluster showed a very abrupt increase in message followed by rapid decrease. These genes were lymphocytic in origin and were modeled combining the fast gene induction effect of CS in lymphoid cells and its direct lymphocyte trafficking effect. Another cluster showed reduction persisting for 18 h, which was described by CS inhibition of mRNA synthesis. Our results reveal the marked diversity of genes regulated by CS via a limited array of mechanisms. These PK/PD models provide quantitation of CS pharmacogenomics and new hypotheses regarding understanding of diverse mechanisms of CS receptor-gene mediated action.
The most potent anti-inflammatory/immunosuppressive drugs are synthetic glucocorticoids collectively referred to as corticosteroids. They are widely used for diseases such as lupus erythematosus, bronchial asthma, leukemia, and organ transplantation. Beneficial effects derived from inhibition of immune system are accompanied by numerous metabolic side effects, including osteoporosis, muscle wasting, steroid diabetes, and others (Schimmer and Parker, 1995).
Both the therapeutic and adverse effects of CS are produced by binding to the glucocorticoid receptor (GR) (Perez et al., 2001). Based on their time courses, responses can be classified as rapid effects (cortisol suppression and cell trafficking) and delayed effects (gene-mediated enzyme/cytokine induction or repression) (Jusko, 1994). Virtually all tissues contain genes that are at least in part subject to transcriptional regulation by drug-activated GR. CS also influence the activity of transcription factors, such as nuclear factor-κB, activator protein-1, and CAAT enhancer-binding protein (C/EBP), via “transcriptional cross-talk” and/or change of their gene expressions (Gotoh et al., 1997; Gottlicher et al., 1998; Miesfeld, 2001). Alteration of a critical transcription factor can affect expression of a battery of downstream genes. Genes affected by CS include both immunosuppressive genes for therapeutic effects and metabolic genes for adverse effects. Better understanding of CS pharmacogenomics will not only provide more insights on their mechanism of action but also aid in optimization of clinical therapy.
In spite of the rich literature concerning CS effects on gene expression, most studies were done in vitro. Little information is available on the temporal pattern of mRNA changes after in vivo treatment with CS. In the past decade, our laboratory has been involved in exploring the time profiles of several individual genes after CS treatment in rats using quantitative Northern hybridization (Sun et al., 1998, 1999; Almon et al., 2002; Ramakrishnan et al., 2002b). Two distinct time profiles were observed after a single-dose of MPL: delayed induction and delayed repression.
Development of mechanistic PK/PD models for CS effects is essential for quantitative understanding of its molecular and cellular mechanisms. A series of PK/PD models have been proposed to describe the major physiological system and the rate-limiting steps for both rapid and gene-mediated effects of CS in animals and humans (Jusko, 1994; Sun et al., 1998; Chow et al., 1999; Ferron et al., 1999; Ramakrishnan et al., 2002a). Our most current “fifth-generation model” describes the pharmacokinetics, GR mRNA repression, receptor dynamics, tyrosine aminotransferase mRNA induction, and its enzyme induction in liver upon single and double doses of MPL in adrenalectomized (ADX) rats (Ramakrishnan et al., 2002a).
However, traditional message quantification methods such as Northern blot and reverse transcriptase-polymerase chain reaction only allow measurement of single genes. This greatly limits the number of genes that can be well explored. Because of the broad nature of CS effects, indications of efficacy and toxicity cannot be evaluated comprehensively based on a limited number of genomic changes. Lately, microarray-based assay technology enables measurements of thousands of genes in a single experiment. The two commonly used systems are cDNA and oligonucleotide microarrays.
Recently, we used cDNA microarrays to study the change of 5200 genes at 5.5 h in rat liver after 50 mg/kg MPL i.v. injection (Almon et al., 2002). The CS pharmacogenomics observed in this microarray dataset was limited by investigating only one time point.
In the current study, we describe the use of oligonucleotide microarrays (Affymetrix GeneChips; Affymetrix, Inc., Santa Clara, CA) to examine the entire profile of temporal changes of 8000 genes in rat liver after a single dose of 50 mg/kg MPL. The liver was selected because of its major role in both the efficacious and adverse effects of CS. In liver, these agents influence the expression of serum proteins that regulate immune/inflammatory responses (Morand and Leech, 1999). Expression of liver enzymes involved in metabolic effects, such as gluconeogenesis and lipid metabolism, is also altered (Berdanier, 1989; Andrews and Walker, 1999). In addition, CS modify the metabolism of many endogenous and exogenous substrates (Prough et al., 1996). This is possibly the first extensive microarray dataset describing the time dependence of in vivo responses of a tissue to any drug. We found 197 CS-responsive probes representing 143 different genes that were sorted into six clusters based on their temporal behaviors. Our major goal is to develop mechanism-based PK/PD models that can quantitatively describe CS pharmacogenomics, provide insights on their mechanisms of regulation, and understand the signaling network in this biological system. In this report, we expand our fifth-generation model that described simple gene induction to several pharmacogenomic models that may explain the unique expression time patterns of all the CS-responsive genes observed in our study.
Materials and Methods
Experimental Procedures
Liver samples were obtained from a previously performed animal study in our laboratory (Sun et al., 1998). In brief, 47 male ADX Wistar rats weighing 225 to 250 g (Harlan, Indianapolis, IN) underwent right jugular vein cannulation under light ether anesthesia 1 day before the study. Forty-three animals received a single intravenous bolus dose of 50 mg/kg MPL (Solu-Medrol; The Pharmacia-Upjohn Co., Kalamazoo, MI) via cannula. Rats were sacrificed by exsanguinations under anesthesia at 0.25, 0.5, 0.75, 1, 2, 4, 5, 5.5, 6, 7, 8, 12, 18, 30, 48, or 72 h after dosing. The sampling time points were selected based on preliminary studies describing GR dynamics and enzyme induction in liver. Four untreated rats were sacrificed at 0 h as controls. Liver was rapidly excised, flash frozen in liquid nitrogen, and stored at –80°C. Frozen liver tissues were ground into powder using a liquid nitrogen chilled mortar and pestle.
Microarrays. Liver powder (100 mg) from each individual animal was added to 1 ml of prechilled TRIzol Reagent (Invitrogen, Carlsbad, CA). Total RNA extractions were carried out according to manufacturer's directions. Extracted RNAs were further purified by passage through RNAeasy mini-columns (QIAGEN, Valencia, CA) according to manufacturer's protocols for RNA cleanup. Final RNA preparations were resuspended in diethyl pyrocarbonate-treated water and stored at –80°C. RNAs were quantified spectrophotometrically and purity assessed by gel electrophoresis.
Isolated RNA from each individual liver was used to prepare target according to manufacturer's protocols. The biotinylated cRNAs were hybridized to 47 individual Affymetrix GeneChips Rat Genome U34A (Affymetrix, Inc.), which contained 7000 full-length sequences and approximately 1000 expressed sequence tag clusters. Unlike cDNA arrays used in the previous study, the high reproducibility of in situ synthesis of oligonucleotide chips allows accurate comparison of signals generated by samples hybridized to separate arrays, which is applicable to study temporal patterns.
Cluster Analysis. Affymetrix Microarray Suite 4.0 (Affymetrix, Inc.) was used for initial data acquisition and basic analysis. In this first step, a “call” of Present, Absent, or Marginal was determined for each probe set on each chip. The results were also normalized for each chip using a distribution of all genes around the 50th percentile. The dataset was scrubbed to retain only those genes that indicated at least four Present calls in the Affymetrix analysis. The scrubbed dataset was then used for cluster analysis.
The results from the first step were inputted to another program, GeneSpring 4.1 (Silicon Genetics, Redwood City, CA). Unlike the Affymetrix software, GeneSpring was designed to accommodate temporal sequence analysis. Each individual gene was normalized as a ratio to the mean of the four controls, which had a distribution around 1. Three approaches were used to define clusters of genes with similar temporal patterns. The first approach was a supervised method that allowed the dataset to be queried using a particular pattern. Based on our previous work we inputted the temporal patterns exhibited by two enzymes, tyrosine aminotransferase in liver and glutamine synthetase in skeletal muscle (Sun et al., 1998, 1999). The other two approaches were nonsupervised methods: self-organizing maps and k-means clustering. Both tools divided genes into groups based on expression patterns. Self-organizing maps had the additional benefit of providing an indication of the relationship between clusters. A correlation coefficient of at least 0.95 among genes in each cluster was used as criteria for the cluster analysis.
Pharmacokinetic/Pharmacodynamic/Pharmacogenomic Models
Pharmacokinetics. The PK of MPL was described by a biexponential equation with the parameters reported previously (Sun et al., 1998) (Table 1): where Ci and λi are the coefficients for the intercepts and slopes and CMPL is the plasma concentration of MPL in nanograms per milliliter.
Receptor Dynamics. The cellular processes for CS pharmacogenomics are depicted in Fig. 1. Unbound CS in blood is moderately lipophilic and freely diffuses into the cytoplasm of liver cells. These steroids quickly bind to the cytosolic GR and cause activation of the receptor. This is accompanied by release of heat-shock proteins 90, 70, and 56, hyperphosphorylation, and conformational changes of the receptor. The activated steroid-receptor complex rapidly translocates into the nucleus where it can bind as a dimer to the glucocorticoid responsive element (GRE) in the target DNA and alter rates of transcription of target genes. CS are also known to cause homologous down-regulation of their own receptor (Oakley and Cidlowski, 1993). Binding of the activated steroid-receptor complex to the GRE results in decreased transcription and reduced levels of GR mRNA. The mRNA translocates to the cytoplasm and is translated into protein. This further leads to the decrease in the free GR density in the cytosol. After the transcriptional control of target genes, the steroid-receptor complex in nucleus may dissociate from the GRE and return to the cytosol. Part of the receptors may be degraded during the process, whereas the rest may be reassembled with heat-shock proteins and recycled.
Based on the cellular mechanisms, a PD model as depicted in Fig. 2 was used to describe the receptor dynamics in rat liver after MPL treatment. This was partially adapted from the most current model developed by our laboratory based on acute dosing of MPL in ADX rats (Ramakrishnan et al., 2002a). The differential equations for receptor dynamics include the following: where symbols represent the plasma molar concentration of MPL (D), the receptor mRNA (mRNAR), the free cytosolic GR density (R), cytosolic drug-receptor complex (DR), and drug-receptor complex in nucleus [DR(N)]. The rate constants in the equations include zero-order rate of GR mRNA synthesis (ks_Rm); the first-order rates of GR mRNA degradation (kd_Rm), receptor synthesis (ks_R), and degradation (kd_R), translocation of the drug-receptor complex into the nucleus (kT), and the overall turnover of DR(N) to cytosol (kre); as well as the second-order rate constant of drug-receptor association (kon). In addition, IC50_Rm is the concentration of DR(N) at which the synthesis rate of receptor mRNA drops to 50% of its baseline value, and Rf is the fraction of free receptor being recycled.
The baselines were defined using the following equations: where and R0 are the baseline values of receptor mRNA and free cytosolic GR density. These baseline values were fixed based on means of the control animals. Parameters reported in the fifth-generation model (Ramakrishnan et al., 2002a) were used for simulations of receptor dynamics in the present study (Table 1).
Pharmacogenomics. As depicted in Fig. 1, binding of the activated steroid-receptor complex to the GRE or negative GRE in target DNA leads to the induction or repression of several genes. Binding of the complex to GRE may directly interact with proteins associated with the proximal promoters and either enhance or repress DNA transcription. This binding may also cause rapid and local influences on chromatin structure, followed by conformational changes in the transcriptosome complex. Therefore, RNA polymerase action will be affected, causing either activated or blocked initiation of transcription. These manifest as an increased or decreased rate of transcription, resulting in changed levels of mRNA in the cell. (Miesfeld, 2001)
Furthermore, upon binding to GRE, CS may affect DNA transcription of functional biosignals/transcription factors (represented by BS), such as C/EBP, thereby causing changes in mRNA and protein concentrations of these regulators (Gotoh et al., 1997). Gene expression controlled by these factors at transcriptional and/or post-transcriptional levels will thus be secondarily affected by CS.
The CS may also affect gene expression without binding to GRE in DNA. Steroid-receptor complex can directly or indirectly affect mRNA stability in cytosol. The complex may replace some binding proteins on the mRNA, which are essential for protection of the message from ribonucleases in cytosol (Vedeckis et al., 1989). This leads to an increased degradation rate of the mRNA. In addition, it has been shown that CS can decrease ribonuclease activity in liver (Onyezili and Goodlad, 1981). This may cause repressed degradation of certain mRNAs that are sensitive to ribonuclease activity.
Based on the above-mentioned mechanisms, the following mathematical models were proposed to describe different gene expression patterns after MPL treatment in rat liver. In all these models except model C, target mRNA (normalized as ratio to control) was assumed to be synthesized at zero-order rate ks_m, and degraded at first-order rate kd_m without drug administration:
The message level was assumed to be at steady-state at time 0 (control animals), yielding the following baseline equation:
Baseline message level mRNA0 was fixed as 1 for most of the genes, unless estimation of this parameter showed significant improvement of model fitting.
Model A. Induced transcription. As depicted in Fig. 3A, mRNA with induced production was described as follows: where the enhancement of transcription rate ks_m is dependent on DR(N) concentration with a linear efficiency constant (S).
Model B. Repressed transcription. As depicted in Fig. 3B, mRNA with repressed production was described as follows: where the inhibition of transcription rate ks_m is dependent on DR(N) concentration, and IC50 represents the concentration of DR(N) at which mRNA synthesis rate drops to 50% of its baseline value.
Model C. Induced mRNA degradation in cytosol plus secondarily induced transcription by BS. As depicted in Fig. 3C, mRNA with induced degradation and secondarily induced production was described as follows: where symbols represent the message (mRNABS) and protein (BSr) levels of the intermediate regulator BS (both normalized as ratio to control). The extra rate constants in the equations include the zero-order rate of BS mRNA synthesis (ks_BSm); the first-order rates of BS mRNA degradation (kd_BSm), translation to BS protein (ks_BS), and protein degradation (kd_BS). The stimulation of BS transcription process ks_BSm is dependent on DR(N) concentration with a linear efficiency constant (SBSm). The stimulation of mRNA synthesis, ks_m, is dependent on the relative changes of regulator BS with a linear efficiency constant (Sm_s). This stimulation is present even at baseline conditions. The mRNA degradation kd_m in cytosol is regulated by DR concentration with a linear stimulation factor (Sm_d).
At time 0, eqs. 12, 13, and 14 yield the following baseline equations: where and are the baseline values of normalized BS mRNA and protein levels. These baseline values were fixed as 1.
Model D. Repressed transcription plus secondarily induced transcription by BS. As depicted in Fig. 3D, mRNA with repressed then secondarily induced production was described as follows: where the intermediate regulator BS is described in a simplified manner using a linear transduction model (Sun and Jusko, 1998). In this model, BSa represents the absolute change of regulator level from the control, and this change is produced by DR(N) via first-order rate (k1). Drug [DR(N)] and the regulator (BSa) both act on mRNA synthesis independently, characterized by a sigmoidal inhibition and a linear stimulation, respectively. The IC50_s represents the concentration of DR(N) at which mRNA synthesis rate drops to 50% of its baseline value, and Ss represents the efficiency of BS stimulation on transcription. The initial condition of eq. 18 () was fixed as 0.
Model E. Repressed transcription plus secondarily repressed mRNA degradation by BS. As depicted in Fig. 3E, mRNA with repressed production and secondarily repressed degradation was described as follows: where BSa represents the absolute change of the regulator level from the control, and is characterized as in model D by a first-order rate constant (k1). The inhibition of transcription rate, ks_m, is dependent on DR(N) concentration, and IC50_s represents the concentration of DR(N) at which mRNA synthesis rate drops to 50% of its baseline value. The inhibition of degradation rate kd_m is dependent on absolute changes of BS, and IC50_d represents the changes of BSa at which mRNA degradation rate drops to 50% of its baseline value. The initial condition of eq. 20 () was fixed as 0.
Model F. Induced transcription in immune cells. Because we are using the whole liver to study CS pharmacogenomics, gene expression detected in our current gene array study may not only come from liver tissues but also from blood contained in the sample. As an immunosuppressant, CS is well known to regulate immune-related genes in lymphatic cells. Recent data suggest that at least in lymphocytes there may be membrane receptors for CS that may mediate some gene regulation effects (Watson and Gametchu, 1999; Falkenstein et al., 2000). Another aspect of CS effects is altered lymphocyte trafficking, a rapid effect causing redistribution of lymphocytes in the body (Chow et al., 1999; Ferron et al., 1999). Specifically, CS cause a decrease of lymphocyte counts in blood and concentrates them in lymphoid tissues. Therefore, changes of immune-related message in the liver sample can result from two sources: change of the gene expression in lymphocytes and change of blood lymphocyte number contained in the sample. As depicted in Fig. 3F, immune-related mRNA with induced production was described as follows:
As reflected in eq. 22, a cell trafficking model was used to describe the change in blood lymphocyte counts. According to this model, lymphocytes from the lymphatic tissues enter into blood at a zero-order rate, kin_C, assuming the tissue compartment is very large to produce a constant input rate. The return of lymphocytes from the blood is controlled by a first-order rate constant, kout_C. It is hypothesized that CS instantaneously cause a change in the affinity of lymphocytes for lymphatic tissues thus accumulating the cells in tissues by inhibiting the movement of lymphocytes from tissues to blood. Therefore, an inhibition function is applied on the zero-order entry rate kin_C where IC50_C represents the drug concentration that inhibits kin_C by 50%. The lymphocyte counts in blood (CellB) were expressed as a ratio to the predose value. At time 0, the system is at steady-state, yielding the baseline equation as where is 1. The lymphocyte trafficking was fixed (kout_C = 0.643 h–1 and IC50_C = 6.15 ng/ml) based on literature estimates of MPL effect in ADX rats (Ferron et al., 1999; Ramakrishnan et al., 2002b).
The immediately enhanced message level in single lymphocyte (mRNAC) is described as induced production (ks_m) by CS upon binding to its receptor via a linear efficiency constant (S). The mRNAC is also normalized as ratio to control ( fixed as 1). The total message level in liver sample (mRNA) is represented by the combination of mRNA in a single lymphocyte and total lymphocyte counts from blood source in the sample (multiplication).
Data Analysis. Data in this study were generated from a so-called “giant rat” study in our laboratory. Animals were sacrificed to obtain serial blood and tissue samples. Each point represents the measurement from one individual rat and data from all these different rats were pooled to obtain a time profile as though it came from one giant rat. Mean data from two to three animals at each time point were used to fit the pharmacogenomic models using ADAPT II software (D'Argenio and Schumitzky, 1997). The maximum likelihood method was used with variance model specified as , where V(σ,θ,ti) is the variance for the ith point, Y(θ,ti) is the ith predicted value from pharmacogenomic model, θ represents the estimated structural parameters, and σ1, σ2 are the variance parameters that were estimated. Proposed models were fitted and compared for each individual gene that showed regulation by CS. Models were selected based on visual inspection of curve fitting, estimator criterion value, sum of squared residuals, Akaike information criterion, Schwartz criterion, and confidence of parameter estimations.
Results
Cluster Analysis
From the original 8000 genes on the Affymetrix GeneChips, the scrubbed dataset contained 5300 genes, which showed at least four Present calls in the Affymetrix analysis. Cluster analysis using GeneSpring yielded six different temporal patterns containing 197 CS-responsive probe sets representing 143 different genes (Fig. 4). These genes relate to a variety of processes, including immunosuppression, blood clotting, carbohydrate metabolism, amino acid metabolism, lipid metabolism, and xenobiotic metabolism. Preliminary description and discussion of these clusters are presented in another article from our laboratory (Almon et al., 2003). In brief, cluster 1 was found by the supervised approach as described under Materials and Methods. This cluster contained 23 probe sets representing 21 genes. As depicted in Fig. 4a, the typical pattern for this cluster showed a simple induction with maximum at around 5 h followed by a fast decline to baseline after 18 h. Clusters 2 to 6 were found by self-organizing maps and k-means clustering. Both methods defined fundamentally the same clusters. Cluster 2 contained 9 probes coding for six genes. As depicted in Fig. 4b, this cluster also showed a simple induction with a maximum at around 10 h followed by a slower return to baseline as late as 72 h, compared with the first cluster. Cluster 3 contained 21 probes representing 18 genes. As depicted in Fig. 4c, this cluster showed a more complex pattern with an initial decline for the first 30 min followed by a fast increase to peak at around 12 h. This delayed induction quickly returned to baseline at around 30 h. Cluster 4 contained 66 probes coding for 49 genes. As depicted in Fig. 4d, this cluster showed a longer down-regulation in the first 8 h followed by even more sustained message induction almost throughout the 72-h period. Cluster 5 contained 10 probe sets representing five genes. As depicted in Fig. 4e, this cluster was a unique cluster showing a very fast induction with peak at 30 min followed by an abrupt decrease to baseline by 4 h. Cluster 6 contained 68 probes coding for 44 identifiable genes. Genes in this cluster mostly had low abundance of message in liver, as illustrated by colors of blue and yellow in Fig. 4f. This was the only cluster showing simple down-regulation, which persisted for about 18 h. The decline phase of this cluster was slower compared with the down-regulation phase for cluster 4. This whole dataset is available online at http://microarray.cnmcresearch.org/(link Programs in Genomic Applications).
Pharmacokinetics and Receptor Dynamics
Simulations of various components from the MPL pharmacokinetic and receptor dynamic model are shown in Fig. 5. The MPL plasma concentrations after an i.v. dose exhibit a biexponential decline with a terminal half-life of 35 min. Also shown are the simulated profiles of drug-receptor complexes in cytosol and in the nucleus, which can act as the driving force for genomic effects of CS independently or jointly. Figure 5b illustrates the free GR density in hepatic cytosol and GR mRNA expressions, showing that CS can down-regulate its own receptor at the transcriptional level. The observed mRNA level was nicely captured by the PD model (Sun et al., 1998).
Pharmacogenomics
Cluster 1. Cluster 1 genes were all well captured by model A assuming induced transcription by DR(N), suggesting that all these genes are regulated by CS via similar mechanisms. Figure 6 shows three representative fittings in this cluster. The induction peak of these genes (5 h) exhibited a marked time shift compared with the profile of the driving force DR(N) as shown in Fig. 5a (peak at 2.2 h). Table 2 lists the estimated parameters for all 23 probes using model A. The mRNA degradation rate constant kd_m represents the drug-independent property of the physiological system. The limited range of kd_m estimates from 0.2 to 0.7 h–1 indicates that these cluster 1 genes have similar stability with half-lives ranging from 3.5 to 1 h. The linear stimulation factor S represents the drug specific property of the message. The limited range of S estimates from 0.008 to 0.053 (fmol/mg protein)–1 implies that the transcriptional machineries of cluster 1 genes have similar sensitivity to CS action. All parameters were estimated with reasonable precision (CV% mostly less than 30%). Three separate probe sets of ornithine decarboxylase were all sorted into this first cluster. Interestingly, they showed different magnitudes of induction as partly illustrated in Fig. 6, nos. 1 and 10. This difference in peak values was associated with different S estimates. These three probes showed similar rates of return to baseline yielding similar kd_m estimates.
Cluster 2. Cluster 2 genes were also all well captured by model A assuming induced transcription by DR(N), suggesting that all these genes are regulated by CS via similar mechanisms to cluster 1 genes. Figure 7 shows three representative fittings in this cluster. The induction peak of these genes (9 h) exhibited a further time shift compared to the profile of the driving force DR(N). Table 3 lists the estimated parameters for all nine probes using model A. The slower induction and much slower return to baseline of cluster 2 genes were explained by the smaller kd_m estimates ranging from 0.027 to 0.1 h–1 compared with cluster 1 genes. This suggests that these genes may have greater stability with half-lives ranging from 26 to 7 h. This cluster yielded a large range of S estimates from 0.015 to 2.17 (fmol/mg protein)–1, implying that the transcriptional controls of these genes by CS have a wide range of efficiencies. All parameters were estimated with reasonable precision (CV% mostly less than 30%). Four separate probe sets of α2-macroglobulin were all sorted into this second cluster. They also showed very different magnitude of induction with similar peak time (data not shown). This was described by the different S estimates and similar kd_m estimates from model fitting as well.
Cluster 3. Cluster 3 showed a transient decline in mRNA followed by a delayed induction. This pattern suggests that two mechanisms might be involved in CS action. Various models including models C, D, and E, and others that contained dual mechanisms of CS action were tested for genes in this cluster. All cluster 3 genes were best captured by model C assuming a rapid stimulation of mRNA degradation in cytosol by DR followed by a delayed enhanced transcription that is secondary to the steroid-enhanced expression of a presumed transcription factor BS. The mRNA for BS was assumed to be a cluster 1 gene with linearly induced transcription by DR(N). Thus kd_BSm and SBSm were fixed as 0.3 h–1 and 0.03 (fmol/mg of protein)–1 based on the typical property of cluster 1 genes. Figure 8a shows three representative fittings in this cluster. Figure 8b exhibits the simulated components of model C using arginase (9) as an example. The initial decline in mRNA resulted from a rapid effect of DR. The DR(N)-induced BS mRNA reached a peak at 5 h and returned to baseline after 24 h. The BSr showed a delayed induction with peak at 11 h, which represented the translation process of BS protein from its message. The further delayed induction peak of cluster 3 genes at 12 h resulted from the transcriptional control by BSr. Table 4 lists the estimated parameters for all 21 probes using model C. This cluster yielded a large range of Sm_d estimates from 0.0004 to 0.011 (fmol/mg protein)–1, which correlated with the difference in the nadir value of the initial inhibition of cluster 3 genes ranging from 0.9 to 0.2. The Sm_s estimates also showed a relatively large range from 0.2 to 4.3, which correlate with the difference in the peak value of the delayed enhancement ranging from 1.5 to 3. The turnover rate of message level was described by the combination effects of kd_BS and kd_m. Some parameter estimates were associated with relatively high CV%, implying that the model might be somewhat overparameterized for some genes.
Cluster 4. Cluster 4 showed a fast and prolonged decline in mRNA followed by a further delayed sustained induction. This pattern suggests that, similar to cluster 3 genes, two mechanisms might be involved in CS action. Various models, including models D, E, and others that contained dual mechanisms of CS action were also tested for genes in this cluster. About two-thirds of cluster 4 genes was best captured by model D, assuming a repressed transcription by DR(N) followed by an enhanced transcription that was initiated by the steroid-enhanced regulator BS. These genes typically showed a sustained induction with very slow return to baseline even until 72 h. Figure 9a shows three representative fittings of cluster 4 genes by model D. Fourteen probe sets in this cluster were best described by model E assuming the mRNA induction was caused by inhibition of mRNA degradation by the steroid-enhanced regulator BS. Compared with the genes characterized by model D, the induction profiles of these genes showed a faster return phase with message returning to baseline as early as 24 h. Although the short-lived induction peaks mimicked cluster 3 gene behavior, the long-lasting initial repression of these genes made them sort into cluster 4. Figure 9b shows three representative fittings of cluster 4 genes by model E. The residual 11 probes in this cluster exhibited simple down-regulation for the first 12 h with no induction. These genes were well captured by model B assuming repressed transcription. The very fast decline phase of these cluster 4 genes separated them from the cluster 6 genes, which had a slower decrease. Figure 9c shows three representative fittings of cluster 4 genes by model B. Table 5 lists the estimated parameters for the three subgroups of cluster 4 genes using model D, E, and B. This cluster yielded relatively high kd_m estimates ranging from 0.3 to 14 h–1, suggesting that messages in this cluster may have relatively low stability with half-lives ranging from 2.3 to 0.05 h. These high kd_m values produced the very steep initial slopes of mRNA down-regulation in this cluster. The low IC50_s values ranging from 5 to 91 fmol/mg protein for genes fitted by models D and E implied that the transcriptional machineries of these genes are sensitive to CS repression. With DR(N) as high as 248 fmol/mg protein, maximum inhibition was approached at the current dose of MPL. The rates of onset and return of the induction were partially controlled by the rate of transduction process, k1. In general, genes fitted by model D had smaller k1 than genes fitted by model E, which was consistent with their more sustained inductions. The magnitude of induction was controlled by Ss in model D, and IC50_d in model E. Some parameter estimates were associated with relatively high CV%, implying that the model might be somewhat overparameterized. Ten separate probe sets for hydroxysteroid sulfotransferase were all sorted into this fourth cluster and were all captured by model D. They showed very similar time patterns with a little difference in magnitude of induction (partially shown in Fig. 9a, nos. 1 and 7). This was explained by the similar k1, kd_m, and IC50_s estimates and different Ss estimates from model fitting.
Cluster 5. Cluster 5 was unusual in that the genes were immune-related and were unlikely to be expressed in hepatocytes. The fact that they were all sorted to the same cluster suggested that they had a common cellular origin and were regulated via the same mechanism. Because whole liver was used in our study, it is likely that these RNAs for these genes were from blood lymphocytes in the sample. This cluster showed a very abrupt increase in message followed by a very abrupt decrease. Neither characteristic suggested that the nuclear component DR(N) was involved. The almost immediate increase in message implied a rapid cytosolic effect. Because CS also cause an acute effect on lymphocyte trafficking, it is possible that lymphocyte efflux from the liver blood stream could account for the abrupt decrease of these immune-related message levels. Various models using different driving forces of drug action on mRNA via different mechanisms were tested for genes in this cluster with or without the complication of cell trafficking effects. All cluster 5 genes were best captured by model F assuming a rapid stimulation of mRNA production in lymphocyte by DR in combination with a decrease in blood lymphocyte counts by steroid. Figure 10a shows three representative fittings in this cluster, and Fig. 10b shows the simulation of the change of blood lymphocyte counts. The very steep incline was described by the fast action of DR. After MPL administration, lymphocytes leave the blood stream, yielded a nadir count at 4.5 h, and returned back to normal after 12 h. This drop of lymphocyte number abruptly reduced the total mRNA expression to even lower than baseline value for some of the genes around 5 h (Fig. 10a, no. 3). Once the cell counts were back to normal but the gene induction effect was still not fully diminished, a small second peak was observed in some of the genes around 10 h (Fig. 10a, no. 7). Table 6 lists the estimated parameters for all 10 probes using model F. Genes in this cluster yielded relative large kd_m values, which reflected the fast turnover rate of mRNA. Most parameters were estimated with reasonable precision.
Cluster 6. About 80% of cluster 6 genes exhibited simple down-regulation and were well captured by model B assuming repressed transcription by DR(N), suggesting that all these genes were regulated by CS via similar mechanisms. Figure 11a shows three representative fittings by model B. The residual 15 probe sets exhibited down-regulation for the first 12 h followed by a delayed induction persisting almost until 72 h. Similar to cluster 4 genes, various models, including models D and E that contained dual mechanisms of CS action, were tested for these 15 probes. Seven of these probes were best described by model D, and the other eight were best described by model E. Figure 11b shows three representative fittings of cluster 6 genes by model D, and Fig. 11c shows three representative fittings of cluster 6 genes by model E. Compared with the cluster 4 genes that also had three subgroups fitted by models B, D, and E, cluster 6 genes had a slower decline phase and the repression lasted longer. Table 7 lists the estimated parameters for the three subgroups of cluster 6 genes using models B, D, and E. This cluster yielded relatively low kd_m estimates ranging from 0.2 to 2 h–1, suggesting that messages in this cluster may have relatively high stability with half-lives ranging from 3.2 to 0.35 h. These low kd_m values produced the shallower initial slopes of mRNA down-regulation in this cluster. The IC50_s values ranging from 40 to 220 fmol/mg protein for genes fitted by model B implies that the transcriptional machineries of these genes may have intermediate sensitivity to CS repression. The low k1 estimates from genes fitted by models D and E produced the delayed and sustained induction profile. Most parameters were estimated with reasonable precision associated with CV% less than 30%.
Discussion
The CS can regulate gene expression at various steps from DNA to mRNA, alone or jointly with other transcription factors and/or hormones. CS may induce or inhibit the same gene depending on other regulatory proteins already present on DNA control regions (Alberts et al., 2001). Furthermore, drug may initiate cascades of events where changes in one gene product may become the controlling factor for changes of other genes. Our microarray dataset contains substantial information about the mechanism and extent of CS effects on various genes. The small number of distinct temporal patterns indicates that a limited number of mechanisms may mediate CS pharmacogenomics.
Cluster 1 had patterns similar to tyrosine aminotransferase and glutamine synthetase. It is widely accepted that both enzymes were induced at the transcriptional level by CS-receptor complex binding to GRE (Max et al., 1988) and could be well captured by our fifth-generation model. Cluster 1 genes were described using the same model. Hirvonen et al. (1988) reported that hepatic ornithine decarboxylase mRNA showed 60-fold induction at 4 h after 5 mg/kg dexamethasone i.p. injection in rats. This is comparable with our findings. The higher magnitude of induction was possibly due to the higher potency of dexamethasone than MPL.
Cluster 2 had strong similarity to cluster 1, and the same model was applied. Kurtz et al. (1978) reported that hepatic α2-macroglobulin mRNA was induced after 20 mg/kg hydrocortisone i.p. injection in rats. As demonstrated by Chen and Feigelson (1978), this induction by CS could be prevented by the transcriptional inhibitor α-amanitin in rat hepatocytes. Additionally, α2-macroglobulin mRNA is very stable with a half-life longer than most hepatic protein mRNAs, which is consistent with our small kd_m estimates. The α2-macroglobulin is a protease inhibitor that binds to proteolytic enzymes. Induction of its message may contribute to protein regulation by CS.
In cluster 3, the initial decline was short-lived, and the predominant induction was assumed to be secondary to a CS-enhanced transcription factor. There was strong evidence for at least two genes, arginase and carbamyl phosphate synthetase, that the delayed induction was secondary to the primary CS-enhanced transcription factor C/EBP. Gotoh et al. (1994, 1997) reported that the enhancer region of rat arginase gene contained at least two binding sites for C/EBP. After incubation with dexamethasone, rat hepatoma cells showed induced C/EBP mRNA at 0.5 h, whereas arginase mRNA had no induction until after 6 h. The C/EBP mRNA was primarily induced by CS at transcription without requiring de novo protein synthesis (Matsuno et al., 1996). However, CS effect on arginase mRNA disappeared after adding cycloheximide, a protein synthesis inhibitor (Takiguchi and Mori, 1995). In summary, CS first stimulate a certain primary-response-gene (C/EBP mRNA), whose protein product will subsequently induce transcription of secondary-responsegenes (arginase mRNA). Because of the sensitivity limitation in microarrays, C/EBP message was not detected in the current experiment. To decrease the number of estimated parameters, the primary-response-gene (mRNABS) was assumed to be a cluster 1 gene in model C. This hypothesis is testable by measuring C/EBP mRNA using real-time reverse transcription-polymerase chain reaction.
Most cluster 4 genes also showed dual regulation by CS with longer suppression. This suggests that instead of rapid influence on mRNA stability in cytosol, these genes may be down-regulated by CS-receptor binding to negative GRE on DNA. The secondary induction may be caused by a CS-enhanced biosignal similar to cluster 3, as described in model D. This biosignal may exhibit a delayed and prolonged message induction, similar to cluster 2 genes. Xiao et al. (1995) reported that aldehyde dehydrogenase expression was inhibited by 1 μM dexamethasone in rat hepatocytes. Moreover, arylhydrocarbon receptor acts as a transcription factor by binding to aldehyde dehydrogenase DNA and facilitates its transactivation. In hepatoma cells, arylhydrocarbon receptor could be induced 2.5-fold by dexamethasone (Wiebel and Cikryt, 1990). This evidence supports our model assumptions. Some cluster 4 genes were best fitted by model E assuming secondary inhibition of mRNA degradation. Hepatic ribonuclease activity was suppressed in rats receiving 10 mg/kg prednisolone i.p. injection (Onyezili and Goodlad, 1981). This may serve as general mechanism for stabilizing genes that are sensitive to ribonuclease degradation. The hypothetical compartment (BSa) in models D and E was the simplified form of the primary-response-product (BSr) in model C. This simplification allowed illustration of the underlying transduction process without assuming BSr mRNA behavior.
Cluster 5 contained immune-related genes that may originate from blood lymphocytes. The almost immediate message increase was unlikely to be GRE-mediated. One possibility was inhibition of cytosolic mRNA degradation. However, the very abrupt increase would yield unrealistically high estimation of mRNA synthesis (data not shown). Some recent literature focused on rapid membrane-initiated steroid effects (Watson and Gametchu, 1999; Falkenstein et al., 2000). This pathway involves specific steroid membrane receptors, intracellular second messengers, and effectors on plasma membrane. The intracellular signaling may contribute to genomic modulation of steroid. This information supported model F. The dramatic mRNA decrease might result from removal of blood lymphocytes, as discussed before. This hypothesis is testable by counting lymphocytes in whole liver after in vivo treatment. Normally, 45% of blood lymphocytes migrate to spleen, and 42% move into lymph nodes. Only about 10% distribute to tertiary extralymphoid tissues such as liver (Goldsby et al., 2000). Although hepatic interstitial fluid might contain few, if any, lymphocytes, their number was considered to be ignorable.
Cluster 6 was similar to cluster 4, but with a slower initial decrease. Most of the genes showed simple down-regulation presumably at transcription. For instance, cAMP phosphodiesterase showed message suppression in rat hepatocytes upon adding 0.1 μM dexamethasone (Hermsdorf et al., 1999). Inhibition of phosphodiesterase, the cAMP degradation enzyme, may lead to enhanced cAMP, which is essential for controlling numerous genes/proteins. In one of our studies, single injection of MPL increased hepatic cAMP 1.5-fold for 10 h (data not shown). Ornithine carbamoyltransferase was a cluster 6 gene with dual regulation by CS and captured by model D. Rats receiving 50 μg of hydrocortisone i.p. injection showed ornithine carbamoyltransferase suppression, and 50-μg dibutyryl-cAMP administration could induce this enzyme (Gautier et al., 1985). Therefore, CS-induced cAMP may act as BSa causing the delayed induction. In addition, CS can affect several endocrine factors. Insulin will be induced secondary to the gluconeogenic effect of CS. Genomic changes may involve regulation by several hormones.
Although the literature confirms that CS alter expression of many genes that we observed, it is sometimes difficult to compare our results with previous work quantitatively or even qualitatively, especially for cluster 3, 4, and 6 genes that are found to be regulated differently at different times. Questions of drug, dose, time, in vitro/in vivo differences, and message quantification methods arise. Of special note, several genes such as hydroxysteroid sulfotransferase (Table 5) and ornithine decarboxylase (Table 2) were monitored using different probes. Reasonable concordance in profiles and dynamic parameters indicated a degree of reproducibility in use of microarray probes.
Proper analysis and understanding of microarray results have been challenging since the technology evolved. Cluster analysis will group genes with high correlation coefficients assuming they share similar regulation mechanism. By performing mechanism-based PD modeling for each individual gene, we found some interesting results. Most genes in the same cluster were successfully described by the same mathematical model. However, cluster 4 and 6 each had three subgroups that were modeled differently, implying that the statistical approach was not sensitive enough to identify dissimilarity. Overall, six models were used to capture all CS-responsive genes. Although clusters 1 and 2 had statistically different patterns, they could be captured by the same model, implying that the difference was not caused by different mechanism but by different sensitivity of mRNA to endogenous/exogenous regulators. Clusters 4 and 6 also shared the same general models. At present, we analyzed genes individually. Further integrated models incorporating multiple gene inter-regulations will provide additional insights into signaling networks at molecular, cellular, and systemic levels. The present models serve to provide hypotheses on how mRNA turnover is controlled by direct and secondary factors. These models sometimes confirm known mechanisms and sometimes are only possibilities that will need further exploration with specific studies (Hargrove, 1993).
This study was limited by technical factors such as genes available on the commercial chip, limited assay sensitivity especially for low abundance messages, and clustering sensitivity especially for suppressed genes. Our pharmacogenomic dataset was generated from whole liver that contains various cell types. This may mimic the real biological situation, but extra complexities such as with cluster 5 were introduced. Only one dose level was studied, which obliged the use of the linear stimulation coefficient S that limited prediction capability. Chronic dosing conditions are expected to cause additional complications (Ramakrishnan et al., 2002b). The rats were adrenalectomized to eliminate endogenous corticosterone effects, but this might have distorted the natural CS response owing to the disturbed endocrine system.
In summary, we used microarray technology to identify time patterns of diverse CS-responsive genes to evaluate pharmacogenomic processes and develop generalized mathematical models for receptor-gene-transduction dynamics. This facilitates the understanding of the global picture of CS actions and provides new insights for microarray data analysis.
Acknowledgments
We thank Ms. Sukyung Woo for technical assistance in data analysis.
Footnotes
-
This study was supported by Grant GM24211 from the National Institute of Health.
-
DOI: 10.1124/jpet.103.053256.
-
ABBREVIATIONS: GR, glucocorticoid receptor; CS, corticosteroid(s); C/EBP, CAAT enhancer-binding protein; MPL, methylprednisolone; PK/PD, pharmacokinetics/pharmacodynamics; ADX, adrenalectomized; GRE, glucocorticoid responsive element; BS, intermediate biosignal; DR, cytosolic drug-receptor complex; DR(N), drug-receptor complex in nucleus.
- Received April 18, 2003.
- Accepted June 3, 2003.
- The American Society for Pharmacology and Experimental Therapeutics