A549 human lung carcinoma cells were obtained from the American Type Culture Collection (ATCC: Manassas, VA, United States) and cultured in complete medium (DMEM with 14.5 g/L glucose, supplemented with 10% fetal bovine serum and 1% penicillin and streptomycin, all from Invitrogen). Cells were cultured at 37°C in a humidified 5% CO₂ atmosphere. For cell passaging or plating, cells were first washed out with 1X phosphate buffered saline (Invitrogen) and then trypsinized using 1X trypsin-EDTA (10X 0.5% Trypsin, Invitrogen). For all ligand treatments, A549 cells (200,000 cells/well) were seeded into six-well plates (Corning International, NY, United States) and cultured for 24 h in complete medium and maintained in steroid-free medium prepared as complete medium but using charcoal-stripped serum (Adams et al., 2003) for 24 h before ligand or vehicle (ethanol) treatments. DEX, RU486, and CYP (Sigma) stock solutions (1 mM) were dissolved in 100% ethanol and kept at −20°C and further diluted in steroid-free medium before usage.

Total cellular RNA was extracted using the TRIzol^® reagent (Invitrogen) following the supplier’s manual (Invitrogen). Total RNA was dissolved in RNAase free water, denatured for 5 min at 65°C and RNA concentration was quantified by spectrophotometric OD260 measurement using the Bioanalyzer (Agilent Technologies, Palo Alto, CA, United States). RNA samples were stored at −80°C until further use. 1 μg of total RNA was used for cDNA synthesis. In order to remove genomic DNA carry-over, RNA samples were treated with 1.5 u of DNAase I (Invitrogen) for 15 min at 25°C. DNAase I treated samples were then incubated at 65°C for 10 min following addition of 25 nM of EDTA (Invitrogen). Finally, they were reverse transcribed using the iSCRIPT^TM cDNA Synthesis Kit according to the manufacturer’s instructions (Bio-Rad). From each DNAase I treated RNA sample, a non-reverse transcribed (-RT) sample was similarly generated (reverse transcriptase was replaced with water). cDNA as well as –RT samples were kept at −20°C.

Forward and reverse primer pairs against reference (b-actin) and GR were generated using the primer3 Input on line software¹ and designs were based on publicly available human mRNA sequences. Primers were designed to have approximately 50% G/C content and to generate 150–250 bp amplicons. Primer pair specificity against target sequence was checked in the NCBI GenBank database using BLAST². The sequences of the primers used to detect glucocorticoid-induced leucine zipper (GILZ), Solute Carrier 19A (SLC19A), and thrombomodulin (THBD) were provided by Dr. J. C. Wang and have been used before to detect gene expression in A549 cells (Wang et al., 2004). The sequences of the primers used to detect secretory leukocyte protease inhibitor (SLPI) and nuclear receptor subfamily 0 group B member 1 (NR0B1) were as follows: SLPI forward 5′-TCAAATGCCTGGATCCTGTTGA-3′; SLPI reverse 5′-GCATCAAACATTGGCCATAAGTC-3′; NR0B1 forward 5′-TGCTCTTTAACCCGGACGTG-3′; NR0B1 reverse 5′-GCGTCATCCTGGTGTGTTCA-3′. In all cases, primers were supplied by Isogen Life Sciences (Netherlands) and dissolved in water according to the supplier’s instructions and kept at −20°C until use.

qPCR monitoring and analysis was performed using the LightCycler^® Carousel-Based Detection System 2.0 (Roche). PCR reactions were performed in a total volume of 10 μl containing 2 μl of LightCycler^® FastStart DNA Master^PLUS SYBR Green I master mix (Roche), 2 μl undiluted cDNA and 1 μl of each forward (5 pmol/μl) and reverse primer (5 pmol/μl). Every PCR reaction mix was filled in the LightCycler glass capillaries which were subsequently closed and centrifuged using the LC Carousel Centrifuge 2.0 (Roche). Cycling conditions were a single pre-incubation step at 95°C for 10 min followed by 45 cycles of 10 s at 95°C, 10 s at 60°C and 10 s at 72°C. To verify that the primer pairs used yielded single PCR products, a dissociation protocol was added after thermocycling, determining dissociation of the PCR products from 65 to 95°C for 15 s. Finally, a cooling step was set for 20 s at 40°C.

To estimate the efficiency of the amplification reaction, serial half logarithm unit dilutions of cDNA from the A549 cells were used and standard curves were generated. The linear slope of the standard curve for each primer pair was estimated using GraphPad Prism 4 software and the efficiency was calculated based on the following equation (1).

Additionally, the -RT samples and a water-template were included in the analysis to confirm the absence of any residual DNA or contamination. All cDNA samples were analyzed in triplicates. Finally, the following equation (2) was used to calculate the fold induction of gene expression.

Transfections with siRNAs against the GR were performed as described before (Fitzsimons et al., 2013). Briefly, a total of 100 pmol of a GR targeting (Hs_NR3C1_6_HP validated siRNA, Qiagen and siGENOME NR0B1 siRNA, Dharmacon) or a non-targeting control siRNA (AllStrars Neg. siRNA AF 546, Qiagen or siGENOME Non-Targeting Control siRNAs #1, Dharmacon) were transfected into A549 cells using the Nucleofector I (Lonza) and the Cell Line Nucleofector^® Kit T for the A549 cell line, according to the supplier’s instructions (Lonza) using the U-29 Nucleofector program. The control siRNA was tagged with a red fluorophore, to monitor transfection efficacy, which was always higher than 90%. The cell medium was refreshed 1 day after transfection and all ligand treatments were performed 3 days after transfection.

Total cellular proteins were extracted from A549 cells with fresh RIPA lysis buffer (200 μl/well) on ice, transferred into 1.5 ml eppendorfs, mixed and kept at −20°C until use. Total protein concentration was measured using the BCA^TM Protein Assay Kit according to the supplier’s guidelines (Pierce). 3X of sample buffer was added to 10 μg of total protein, denatured for 5 min at 95°C, spinned shortly and then loaded onto 10% SDS-polyacrylamide gels. Samples were run through stacking gel at 100 V for 10 min and separating gel at 200 V for approximately 1 h. Protein transfer onto methanol-activated Immobilon^TM – P^SQ Transfer membranes (Millipore) was performed overnight at 4°C at 125 mA.

For the detection of the GR protein levels, the blots were incubated in blocking buffer, consisted of 5% low fat milk powder in TBST solution for 1 h at RT (10 ml/membrane) and subsequently with a primary antibody against the GR (H-300 rabbit polyclonal IgG, Santa Cruz), or against NR0B1 [Anti-NR0B1/Dax1 antibody (EP13786) – N-terminal (ab196649), abcam], or against alpha-tubulin antibody (clone DM1A, Sigma), or against GAPDH antibody (Santa Cruz). Primary antibodies were added in blocking buffer (1:2000 and 1:1000 dilution, respectively) for 1 h at RT (5 ml/membrane). Following 3X washing with TBST, the blots were probed with species-specific horseradish peroxidase-conjugated secondary antibodies (Santa Cruz) in blocking buffer (1:5000 dilution) for another 1 h at RT (10 ml/membrane) and finally washed 5X with TBST. For all incubations and washes rolling shakers were used. For luminescent signal detection, membranes were incubated with 10 ml of luminol solution, supplemented with 100 μl of enhancer solution and 3.1 μl of 30% H₂O₂ for approximately 1 min at RT in the dark. Following film exposure, development and fixation, GR protein levels among samples were quantified relatively to a-tubulin signal using the ImageJ software³ (Rasband, W.S., ImageJ, United States National Institutes of Health, Bethesda, MD, United States, 1997–2012).

Interactions between the GR-LBD and coregulator NR-box peptides were determined using a MARCoNI assay (PamChip no. 88011; PamGene International) as described before (Koppen et al., 2009; Zalachoras et al., 2013). Each array was incubated with a reaction mixture of 1 nM Purified Glucocorticoid Receptor Recombinant Human Protein, Ligand Binding Domain, (Thermo Fisher Scientific, cat # A15668), ALEXA488-conjugated anti-GST antibody and buffer F (PV4689, A-11131, and PV4547; Invitrogen). For ligand induced peptide interaction profiling experiments 1 μM DEX, RU486, CYP or solvent (2% DMSO in water) were added. Incubation was performed at 20°C in a PamStation96 (PamGene International). GR binding to each peptide on the array, reflected by fluorescent signal, was quantified by image analysis using BioNavigator software (PamGene International).

GR translocation to the nucleus was studied using a YFP-GR construct kindly supplied by Dr. Cidlowski (National Institute of Environmental Health Sciences, National Institutes of Health) as previously described, with some modifications (Fitzsimons et al., 2008). Briefly, the previously described protocol was scaled down to a 96-wells plate format and semi-automated. 6000 A549 cells/well were plated 24 h prior to transfection. The cells were transfected with YFP-GR plasmid using a Nuclefector I (Lonza) as described before. Complete medium was refreshed 24 h after transfection and 48-h after transfection cells were incubated for 6 h in steroid free medium. After this procedure, the GR localized to the cytosol in all cells as described before (Fitzsimons et al., 2008). The test compounds or vehicle (ethanol) were manually dispensed into the corresponding wells and cells were incubated with the compounds for 30 min. Subsequently, cells were fixed with 80% acetone in water and stained with Hoechst 3342 (1:10.000) for nuclear staining. All experiments were performed in triplicates. Three non-overlapping images were taken from each well by using a Zeiss Axiovert 200/200M inverted microscope and 10× magnification. DAPI and FITC filters with excitation wavelength 409 and 487 nm were used to excite Hoechst (blue emission) and YFP (green emission), respectively. All images were collected with the same settings in Microsoft Window’s BMP format. For post-acquisition image analysis Images were opened in ImageJ software and image backgrounds were subtracted by using the built in subtract background command with Rolling ball radius parameter set at 50. Subsequently, the images were opened with the CellProfiler software⁴ to automatically identify the nuclear and cytoplasm compartments of cells and the green fluorescence intensity in each compartment (Carpenter et al., 2006).

For comparisons between groups a two-tailed Student’s t-test was applied using GraphPad’s Prism 5 Software. For multigroup comparisons, a one-way ANOVA test with a Tukey’s post test was performed using the same software package. Dose response curves were fit to a sigmoidal (four-parameter logistic) curve using GraphPad’s Prism 5.

To model ligand-specific effects on GR-mediated gene transcription, we used three GR ligands with different pharmacological characteristics: DEX is a well-characterized GR agonist, while RU486 is usually described as a partial agonist (often used as antagonist), and CYP as a passive antagonist (Rousseau and Baxter, 1979; Honer et al., 2003; Meijer et al., 2005; Matthews et al., 2009). We studied the effects of increasing concentrations of DEX on the expression of three GR responsive genes in A549 cells. Glucocorticoid Induced Leucine Zipper (GILZ) is an important mediator of the anti-inflammatory effects of GCs (Ronchetti et al., 2015), THBD is an endothelial cell surface glycoprotein that controls thrombosis by downgrading thrombin-mediated fibrin generation and promoting protein C activation (Loghmani and Conway, 2018) and SLC19A2 is a thiamine transporter associated with the thiamine-responsive megaloblastic anemia syndrome (TRMA) (Aoyagi and Archer, 2011). In this cellular system, DEX dose-dependently induced the expression of the three genes tested, albeit with different pharmacological parameters. We found significant differences in calculated effective concentration 50 (EC₅₀) and maximal response (R_max) values. Specifically, DEX was significantly more potent in inducing GILZ than THBD or SLC19A2, and significantly less efficacious in inducing SLC19A2 than THBD or GILZ (Figure 1A and Table 1).

These differences are difficult to explain from the standpoint of theoretical models of ligand-induced gene expression because the ligand, the amount of receptors and all other components of the cellular environment were the same in all experiments. One possible explanation could be in the structure of the gene promoter and its GRE composition (Wang et al., 2004; So et al., 2008), suggesting that response magnitude depends on the strength of GR binding to GREs. However, numerous previous studies argue to the contrary (So et al., 2008; Pfaff and Fletterick, 2010; Dougherty et al., 2012; Chen et al., 2013), and previously published chromosome immunoprecipitation (ChIP) data analyzing these promoters do not predict all the differences in sensitivity to the GR agonist DEX (Wang et al., 2004). The promoter sequence data predicted GILZ to be strongly bound by the occupied GR with two high GR binding score (GBS) sites in its promoter, followed in order of predicted responsiveness by SLC19A, with one high GBS sites, and THBD with one low GBS site (Figure 1B), while DEX dose response curves showed a different order in both response efficacy (THBD∼GILZ>>SLC; Figure 1A) and potency (GILZ>>THBD∼SLC; Figure 1A). Similarly, RU486 induced significant changes in GILZ, SLC19A2 and THBD, expression in A549 cells, albeit only at high doses (Figure 1C and Table 2).

Importantly, the observed order of gene-transcription potency (SLC19A∼THBD>GILZ; Table 2) and efficacy (GILZ>SLC19A∼THBD; Table 2) induced by RU486 were different than those induced by DEX (Table 1). Noteworthy, although CYP alone did not have any detectable effect on the expression of the three genes analyzed (Figure 1D), DEX, RU486 and CYP induced GR translocation from the cytosol to the nucleus, indicating that the three ligands promote active changes on GR behavior and therefore cannot be ascribed as passive (or inactive) antagonists (Table 3).

To understand the differential, ligand- and gene-specific, biological behavior observed in A549 cells in more detail, we aimed to develop an MTC model that could integrate ligand–receptor–DNA–coregulator interactions.

An alternative hypothesis to explain the gene-specific response to DEX observed in A549 cells involves the differential recruitment of specific coregulators to the MTC active at each gene promoter, which could induce gene-specific expression. To model this situation, we postulated a theoretical equilibrium model of receptor action that explicitly includes the receptor (R), the ligand (L), the coregulator(s) (C), and the DNA (D) and four parameters that govern receptor species equilibria: α represents the effect of ligand binding on the binding of coregulator, β describes the effect of coregulator binding on the binding of receptor to DNA and γ the effect of ligand binding on the binding of receptor to DNA. In turn, δ represents the extent to which the joint effect of any two of ligand binding, coregulator binding or DNA binding varies conditional on the level of the third. A detailed description of all these parameters is presented as Supplementary Material and schematically in Figure 2. Our model assumes that the GR can spontaneously couple to the coregulators or the DNA even in the absence of the ligand (Power et al., 1991; Mani et al., 1994). Thus, the model is fully described by three basic equilibrium constants that account for ligand binding, coregulator coupling and DNA binding and four parameters that illustrate the interaction effect between them (Figure 2A and Supplementary Information S1).

Assuming that the MTC formed by the receptor, the ligand, the coregulator(s) and the DNA is responsible for the final induction of gene transcription, dose–response curves can be simulated with equations describing how the relative composition of the different components of the MTC affects ligand-dependent gene induction (Figure 2B–E). The model predicted a series of system characteristics that could be experimentally validated. According to the model, the concentration–response curve to a ligand can be modified adjusting the values assumed for each component of the MTC (Supplementary Information S1). Indeed, a first validation of the model using siRNAs to reduce GR expression in A459 cells resulted, as indicated by our model, in significant changes in R_max but not EC₅₀ (Figure 3). Transfection of A549 cells with increasing siRNA concentrations resulted in concomitant decreases in GR expression (Figure 3A,B) and in GILZ’s R_max to DEX, without affecting its EC₅₀ (Figure 3C and Table 4).

Importantly, when changes in the α, β, or γ parameters are simulated, the results predict changes in either R_max or EC₅₀, but never both simultaneously (Figure 2B–D and Supplementary Information S1). This prediction leaves out the possibility that differences in coregulator recruitment or in DNA binding alone could explain the simultaneous change in EC₅₀ and R_max observed in DEX concentration–response curves in A549 cells (Figure 1). However, when the δ parameter is taken into consideration, simultaneous variation in both EC₅₀ and R_max can be explained theoretically using our model (Figure 2E). δ describes how two binding events impact on a third within the MTC proposed by the model, reflecting interaction between all the components of the complex. Therefore, our model suggests that non-independent binding events between the different components of the multifactor complex could be responsible for the simultaneous differences in EC₅₀ and R_max to DEX observed in A549 cells (Figure 1A and Table 1).

The sequence of binding events leading to transcriptional activation by the GR has been intensively studied but not fully established yet (Rogatsky et al., 2003; John et al., 2009; Dougherty et al., 2012; Chen et al., 2013). It is considered to be a dynamic process and the most parsimonious hypothesis is that transcriptional activation by the GR, and NRs in general, involves multiple factors that act in both a sequential and combinatorial manner to reorganize chromatin templates (Pollard and Peterson, 1998; Glass and Rosenfeld, 2000; Voss et al., 2006; Stavreva et al., 2012). Moreover, the temporal order of the events leading to the formation and composition of the MTC that leads to transcription activation can take place in a gene- and cell-specific manner (Gronemeyer et al., 2004). Potential dissimilarities in the structure of promoters for GILZ, SLC19A, and THBD genes do not explain the differential expression patterns obtained with DEX and RU486 (Figure 1A vs. Figure 1C). This divergence can be explained by our MTC model, considering that each ligand induces a specific pattern of coregulator binding to GR forming the ternary complex LRC, which in turn may display differential affinity for DNA. Therefore, using a DNA-free system represented by the MARCoNI peptide array, we focused on understanding how the interactions between ligand and GR could cooperate to modulate coregulator binding.

Our model includes the existence of ligand-specific parameters governing multifactor complex formation (Zalachoras et al., 2013; Atucha et al., 2015) that can be tested experimentally in DNA-free conditions (Supplementary Information S1). Ligand affinity constant Ka (specific for each ligand), coregulator affinity constant Kc (specific for each individual coregulator NR-box), and the parameter α (characteristic of each ligand/coregulator pair) indicate how ligand and coregulator affects each other’s binding to the receptor and are key components of the model when only receptor, ligand, and coregulators are present. According to the model’s prediction, GR ligands should induce a characteristic coregulator binding profile as previously demonstrated (Zalachoras et al., 2013; Atucha et al., 2015).

Ligand-independent and -dependent coregulator binding profiles were measured using a MARCoNI peptide array. This array contained 53 coregulator-derived peptides representing a wide range of coregulator NR-boxes known to interact with the GR and other nuclear receptors, which sequence details, Gene Name and UniProt Knowledge Base accession numbers are shown in Supplementary Table S2. First, NR box peptides were titrated against the GR LBD in the absence of ligand using a customized MARCoNI array (Figure 4). A concentration series of each NR box peptide was immobilized on the array and incubated with 3 nM of the GR LBD and binding isotherms were calculated from these data. We observed a significant basal binding of the GR LBD to several NR-box peptides in the absence of any ligand (Figure 4A). This basal coregulator recruitment in the absence of ligand is predicted by the model (RC species in Figure 2) and is in agreement with similar observations done using fluorescence polarization assays (Pfaff and Fletterick, 2010). Twelve NR box peptides showed positive ligand-independent binding to the GR LBD, although with variable binding profiles (Supplementary Table S3). Ten of these NR box peptides showed lower affinity binding profiles, while two NR box peptides showed higher affinity binding. These higher affinity binding peptides corresponded to the coactivator PRGC1 and the corepressor NRIP1 (Figure 4B and Supplementary Table S3). Secondly, we used a MARCoNI array in which 1 mM of each NR box peptide was immobilized and incubated with 3 nM of the apo GR LBD in the presence of a receptor-saturating concentration of three selected GR ligands or a vehicle control, to study the effect of GR ligands on the basal binding of GR LBD to NR box peptides. We analyzed the coregulator binding profile induced by DEX, RU486, and CYP. These three ligands induced characteristic coregulator binding profiles, in some cases favoring and, in some others, disfavoring the basal binding of the apo GR LBD to NR box peptides observed in the absence of ligand (Figure 5). Notably, the NR box peptide binding profile induced by RU486 resembled for some peptides the effect induced by DEX, although RU486’s effects were substantially weaker in many cases, reflecting its partial agonist activity (Figure 5A,B). Interestingly, GR LBD binding to NR box peptides from the coactivators NCOA2, NCOA3 and NR0B1 was strongly favored by DEX, while it was disfavored by RU486 (Figure 5A,B), suggesting that this differential effect on coregulator binding is involved in the differential pharmacological effects of the two ligands. Supporting this preliminary conclusion, we did not find any example of the opposite pattern in our dataset, this is, binding to NR box peptides that were disfavored by DEX but promoted by RU486.

In contrast, the NR box peptide binding profile induced by CYP was more divergent from that induced by DEX (Figure 5B,C). We did not find any example in our dataset of NR box peptides for which binding to GR was favored by CYP but disfavored by DEX (Figure 5B). Interestingly, from the 10 NR box peptides whose binding was disfavored by CYP, 7 were favored by DEX (Figure 5B), suggesting that these differences in coregulator binding may be crucial to understand the pharmacological differences between DEX and CYP.

In summary, the pharmacological behaviors observed for the three GR ligands used in experiments in A549 cells (Figure 1) can be interpreted using our MTC model considering the interactions between receptor–ligand–coregulator–DNA binding events. Our model suggests that there exist an allosteric phenomenon involving the joint effect of the three GR partners within the MTC (L; C and D), reflected by the δ factor, i.e., a specific ligand induces the recruitment of a specific set of coregulators, that differentially affects the expression of a particular gene. In fact, our experimental observations in vitro with the MARCoNI, support the induction of ligand-specific binding profiles between the GR LBD and NR box peptides (Zalachoras et al., 2013; Atucha et al., 2015) that can explain the differential gene-specific behaviors of GC ligands. The model indicates that the distinctive ligand-specific binding of GR to coregulators would differentially impact on the expression of specific genes when cells are co-incubated with DEX and RU486 or CYP. Indeed, Figure 6 shows that RU486 can block DEX-induced GILZ expression without affecting SLC and THBD expression significantly, while CYP specifically blocked DEX-induced SLC and THBD expression without affecting GILZ maximal expression levels (Figure 6).

According to the MTC, each ligand may induce a specific response by recruiting a defined array of cofactors. Among the coregulators that were differentially bound by GR ligands, we chose NR0B1 as a proof-of-concept experimental validation. NR0B1 displayed a maximum difference between DEX-favored and RU486-disfavored binding to the GR LBD in vitro (Figure 5). SLPI, a protease inhibitor expressed by cells at mucosal surfaces, is stimulated by IL-1β. IL-1 β-induced SLPI expression is increased by DEX and inhibited by RU486 in A549 cells (Ito et al., 2001). We found that SLPI expression is strongly induced by DEX but not by RU486 in non-transfected A549 cells (Figure 7A). While NR0B1 downregulation using specific siRNAs (Figure 7B–D) significantly diminished the maximal DEX-induced expression of SLPI (Figure 7E) and simultaneously increased DEX’s EC₅₀ for SLPI induction, it did not affect RU486-induced SLPI expression (Figure 7F), demonstrating its specific role in DEX-induced SLPI expression in A549 cells.