Immune regulation by glucocorticoids can be linked to cell type–dependent transcriptional responses

Nearly seven decades after their introduction to clinical practice (Hench and Kendall et al., 1949), glucocorticoids remain the most widely used class of anti-inflammatory and immunosuppressive agents. Despite their extensive clinical use, there are still substantial gaps in our understanding of glucocorticoid-mediated immunoregulation, particularly regarding their effects in specific cell types, their key cellular targets in particular disease states, and the actions that are broadly shared among cell types and tissues versus those that are unique to the immune system (Cain and Cidlowski, 2017).

Glucocorticoids act primarily by binding in the cytosol to the glucocorticoid receptor (GR; UniProt no. P04150), a nuclear receptor of the steroid/thyroid hormone receptor superfamily (Stahn and Buttgereit, 2008). The ligand-bound GR translocates into the nucleus and can dimerize and directly bind DNA at specific recognition sequences known as glucocorticoid response elements, increasing transcription rates. Monomeric GR can also bind DNA at a distinct set of recognition sequences known as negative glucocorticoid response elements (Surjit et al., 2011; Hudson et al., 2013), decreasing transcription rates. In addition, ligand-bound GR can be recruited to specific genomic sites without directly binding DNA, via protein–protein interactions with other DNA-bound transcription factors (Sacta et al., 2016). Genomic sites of direct GR binding represent glucocorticoid-induced enhancers, and genomic sites of indirect (tethered) GR binding appear to cluster around and amplify the activity of direct binding sites (Vockley et al., 2016). Composite sites of direct and tethered interactions with DNA have also been described (Sacta et al., 2016). Beyond the direct or tethered recruitment of ligand-bound GR to specific genomic sites, a key component of the mechanism of action of glucocorticoids involves interference with the activity of other transcription factors and signaling molecules, most notably NF-κB. This form of interference can be mediated by direct protein–protein interactions between the ligand-bound GR and other transcription factors (Ratman et al., 2013) but also by indirect effects via inhibitory long noncoding RNAs (Rapicavoli et al., 2013), proteins that dissociate from the GR–chaperone complex upon glucocorticoid binding (Croxtall et al., 2000), or competition for nuclear coactivators. Finally, some of the most rapid effects of glucocorticoids may occur independently of the cytosolic GR. These include alterations in ion transport across membranes (Buttgereit et al., 1993; Schmid et al., 2000), which have been hypothesized to result from intercalation of glucocorticoid molecules into the membrane (Buttgereit and Scheffold, 2002). They also include interactions with membrane-bound forms of GR (Gametchu, 1987; Gametchu et al., 1993; Bartholome et al., 2004).

While the mechanisms are diverse, a consistent outcome of glucocorticoid exposure is a significant reprogramming of a cell’s transcriptional state (Galon et al., 2002; Olnes et al., 2016). The genomic locations of GR binding have been shown to vary widely across cell types (Rao et al., 2011; Grøntved et al., 2013; Love et al., 2017), a phenomenon that is explained at least in part by differences in chromatin accessibility and expression differences of GR cofactors (John et al., 2011; Reddy et al., 2012; Grøntved et al., 2013). This, in turn, suggests that the transcriptional response to glucocorticoids could vary significantly across cell types. In this context, studies of specific cell subpopulations, in the species of interest, are necessary to gain a realistic view of the genomic effects of glucocorticoids in any system. Immortalized and tumor-derived cell lines have been valuable tools for the study of the molecular biology of GR signaling. However, their genomic composition and chromatin landscape are known to differ substantially from those of human primary cells. Similarly, complex cell mixtures, such as whole blood and peripheral blood mononuclear cells (PBMCs), have offered an initial glimpse of the genes and pathways affected by a glucocorticoid stimulus in primary human cells (Galon et al., 2002; Olnes et al., 2016), but they are limited in their ability to discern cell-specific effects.

To develop a greater understanding of how pharmacologic doses of glucocorticoids regulate immunity and the extent to which they differentially affect distinct cell subsets, we studied the genome-wide transcriptional response to glucocorticoids in nine primary human hematopoietic and nonhematopoietic cell types. We then employed this rich set of transcriptome data to generate a pathway-level map of glucocorticoid effects across cell types and tested the ability of this approach to uncover specific actions of glucocorticoids on individual cell types, focusing on a highly relevant population, namely B cells, that contribute to autoantibody formation.

We began by comparing the short-term transcriptional response to glucocorticoids across nine primary human hematopoietic and nonhematopoietic cell types obtained from healthy donors: B cells, CD4⁺ T cells, monocytes, neutrophils, endothelial cells, osteoblasts, myoblasts, fibroblasts, and preadipocytes. The hematopoietic cell types were selected because of their relevance to diseases often treated with glucocorticoids, whereas the nonhematopoietic cells were selected as representative of tissues showing undesirable off-target effects of glucocorticoids. For each cell type, we studied cells from four unrelated donors. Cells from each donor were cultured, treated, and sampled independently. We performed total RNA sequencing (RNA-seq) before, 2 h after, and 6 h after in vitro treatment with methylprednisolone or vehicle. We found evidence of differential expression (adjusted P ≤ 0.05) in at least one cell type, and at one or both time points following treatment, for 9,457 unique genes or ∼17% of the human transcriptome as annotated in GENCODE 28 (Harrow et al., 2012). This is consistent with prior estimates, based on protein-coding genes, which had suggested that glucocorticoids regulate ∼20% of the human genome (Galon et al., 2002). However, we found the transcriptional responses of each cell type to be quite distinct. First, the number of genes that responded to the glucocorticoid stimulus differed substantially among cell types (Fig. 1 a). Overall, the number of glucocorticoid-responsive genes was higher in hematopoietic than in nonhematopoietic cells, and neutrophils had the strongest transcriptional response of the nine cell types studied. Second, the specific genes that responded to the glucocorticoid stimulus were very different in each cell type. Of the 9,457 glucocorticoid-responsive genes identified, 4,806 (50.8%) were significantly differentially expressed after glucocorticoid treatment in only one of the nine cell types, and only 25 (0.3%) were significantly differentially expressed after such treatment in all cell types (Fig. 1 b). As expected, this short list of genes that were differentially expressed in all cell types includes those from classic glucocorticoid-responsive genes such as TSC22D3, DUSP1, and IRF1, which have been studied extensively in many human and animal cells. Most of the glucocorticoid-responsive genes (56.4%) were responsive only in hematopoietic cells, whereas only 21.5% of the glucocorticoid-responsive genes were unique to nonhematopoietic cells (Fig. 1 c). Importantly, however, the strongly cell type–dependent nature of the transcriptional response to the glucocorticoid is evident even when comparing hematopoietic or nonhematopoietic cells separately and when comparing cells of similar ontogenetic origin (Fig. 1 d). Varying the threshold for differential expression has the expected effect on the number of genes that are classified as differentially expressed, but it does not affect the observed cell type dependence of the transcriptional response to the glucocorticoid (Fig. S1).

We then assessed whether the observed differences in the transcriptional response to glucocorticoids among cell types were due to simple differences in the baseline transcriptome of the cells. If, at any given locus, glucocorticoids have a consistent transcriptional effect across cell types, then if two or more cell types express similar transcript levels from a gene at baseline, one would expect the glucocorticoid response of that gene to be similar in those cell types. In that case, the baseline transcriptional repertoire of each cell would determine its set of potential glucocorticoid-responsive genes. Alternatively, glucocorticoids could exert different effects at the same locus, even if the baseline expression levels are similar across cell types, in which case the baseline transcriptional repertoire of each cell type would not be the key determinant of the glucocorticoid response. To examine this issue, we first identified genes that were highly differentially expressed in one group of cells (hematopoietic or nonhematopoietic) but not in the other (Fig. 2 a, left) and assessed their baseline expression levels in each of the two groups (Fig. 2 a, right). Interestingly, we found no correlation between the average response to glucocorticoid and the average baseline level of expression of a gene. Many genes that are similarly expressed at baseline in the two groups of cells (i.e., along the diagonal of Fig. 2 a, right) are only glucocorticoid responsive in one group or the other. One example is the interferon-inducible gene TRIM22, which is similarly expressed at baseline in hematopoietic and nonhematopoietic cells (Fig. 2 a, right) but is only differentially expressed in hematopoietic cells upon glucocorticoid exposure (Fig. 2 a, left). A closer look at the expression of this gene across cell types (Fig. 2 b) reveals that it is similarly expressed at baseline in the nine cell types studied, but the magnitude of the response ranges from strong and highly significant up-regulation in neutrophils to the absence of any detectable change in any of the nonhematopoietic cells studied. An even more striking example of the cell type–specific nature of transcriptional regulation by glucocorticoids is a subset of genes in which not only the magnitude but also the direction of the response differs across cell types. One example is the gene encoding the integrin α subunit 5 (ITGA5), in which transcript levels increase significantly in response to the glucocorticoid in some cell types (most notably osteoblasts, myoblasts, and preadipocytes), whereas they decrease significantly in neutrophils, despite similarly high levels of expression at baseline in each of those cell types (Fig. 2 c).

To obtain a global view of the molecular pathways and biological processes affected by the glucocorticoid stimulus in human primary cells and of the similarities and differences in pathway-level effects across cell types, we performed a gene set enrichment analysis of glucocorticoid-responsive genes with gene sets reflective of molecular pathways and biological processes from the Molecular Signatures Database (MSigDB; Subramanian et al., 2005). Hierarchical clustering of cell types based on their pathway-level responses to glucocorticoid largely follows the known ontogenetic relationships among cell types (Fig. 3 a and Fig. S2). However, at a finer-grained level, each cell type has a distinct pattern of pathway enrichment, and differences are evident even among developmentally and phenotypically similar cells. We used k-means clustering to identify 12 modules of pathways defined by their pattern of response to the glucocorticoid stimulus across the nine cell types (Fig. 3 a and Fig. S2), so that the pathways within each module share a common pattern. Some pathway modules (M1 and M2) were strongly and similarly affected by the glucocorticoid stimulus across the nine cell types. These mostly correspond to pathways in which protein kinase signaling cascades play a central role, such as MAPK signaling and cytokine–cytokine receptor interactions. A second set of modules (M3–M5) involves pathways that are more strongly affected in hematopoietic cells. These include cell adhesion molecules, chemokine signaling, extracellular membrane-receptor interaction, TLR, and B cell receptor (BCR) signaling pathways. A third set of modules (M6 and M7) involves pathways that are more strongly affected in nonhematopoietic cells. These include calcium signaling and glucose, fatty acid, and steroid hormone metabolic pathways. A fourth set of modules (M8–M10; Fig. S2) involves pathways that are similarly affected in hematopoietic and nonhematopoietic cells with a moderate magnitude of enrichment, including the ERBB, WNT, and TGF-β signaling pathways. The two modules with the highest number of pathways (M11 and M12; Fig. S2) include a large number of metabolic pathways in which the magnitude of enrichment is low but similar across all cell types examined.

The cell type specificity of the pathway-level response to the glucocorticoid is even more evident if one considers the predominant direction of change in expression for the genes that are driving the pathway enrichment. In modules with strong pathway enrichment in hematopoietic but not in nonhematopoietic cells (M3–M5), down-regulation of expression was the predominant effect of the glucocorticoid (Fig. S3 a). In contrast, in the larger of the two modules in which pathway enrichment was stronger in nonhematopoietic than in hematopoietic cells (M6), the predominant effect of the glucocorticoid was up-regulation of expression (Fig. S3 b).

We then evaluated whether this pathway-level map could be used to uncover mechanisms of glucocorticoid action in specific cell types. We focused on B cells, which play a key role in antibody-mediated autoimmune diseases, for which glucocorticoids are used extensively. An initial view of the top glucocorticoid-responsive genes in our data set revealed significant changes in multiple genes known to be important in B cell biology. One example is up-regulation of PRDM1 (Fig. 4), which encodes BLIMP-1, a zinc finger protein with a PRDI-BF1 and RIZ homology domain, known to be involved in terminal differentiation and reduced proliferation of B cells. Another example is B cell–intrinsic up-regulation of IL10 (Fig. 5), which has been implicated in immunoregulation by B cells (Mauri and Menon, 2017) and is a known target of glucocorticoids in human monocytes and macrophages (Mozo et al., 2004; Frankenberger et al., 2005).

B cells express both Toll-like receptors and antigen-specific BCRs, which together participate in modulating B cell function and autoantibody production (Rawlings et al., 2012; Suthers and Sarantopoulos, 2017). In our pathway-level map of the effects of glucocorticoids across cell types, these two pathways were in one of the modules with evidence of differential enrichment between hematopoietic and nonhematopoietic cells (M4). Therefore, we selected them for further analysis.

There was a strong suppressive effect of glucocorticoids on transcriptional output from the three human immunoglobulin loci (Fig. 6) and on the expression of several genes that encode key proximal BCR signaling proteins (Fig. 3 b). Some of the observed effects were shared with other hematopoietic cells, while others were unique to B cells. Notably, glucocorticoid exposure led to a significant decrease in the expression of CD79B, which encodes Igβ, a protein required for BCR assembly and for signal initiation after antigen stimulation (Venkitaraman et al., 1991; Matsuuchi et al., 1992). In contrast, expression of CD79A, which encodes Igα, was not significantly affected. The genes CR2 and CD19, which encode the two components of the B cell co-receptor complex that serves as an enhancer of BCR-mediated signaling (Fearon and Carroll, 2000), were significantly reduced in expression in response to the glucocorticoid in vitro. We also observed decreased expression of SYK and BTK, which encode key tyrosine kinases immediately downstream of the BCR complex, and of BLNK, which encodes a B cell adaptor protein immediately downstream of SYK. The gene encoding the critical B cell effector protein phospholipase C-γ2 (PLCG2) and those encoding the phosphoinositide-3-kinase (PI3K) catalytic subunit delta (PIK3CD) and the B cell adaptor of PI3K (PIK3AP1) were also reduced in response to the glucocorticoid. Finally, we observed increased expression of the genes encoding the PI3K regulatory subunit 1 (PIK3R1) and the negative regulators SHIP1 and CD72. Other known negative regulators of BCR signaling showed either reduced expression (FCGR2B and CD22) or no transcriptional change (SHIP2) in response to the glucocorticoid.

To determine whether these in vitro findings were consistent with transcriptional changes in B cells exposed to glucocorticoids in vivo, we measured the expression of the above genes over time in a cohort of 20 unrelated healthy volunteers who received a single intravenous dose of methylprednisolone. B cells were isolated before and 2 and 4 h after the methylprednisolone infusion. RNA was purified and gene expression was measured by real-time PCR. We observed a significant drop in the expression of BLNK, BTK, CD79B, and CR2 and no change in the expression of CD79A (Fig. 7). The genes CD19 and SYK had lower mean expression values after glucocorticoid treatment in vivo, but the difference did not reach statistical significance, whereas the gene that encodes the key upstream kinase LYN was significantly reduced in expression after glucocorticoid administration in vivo, although it had not reached statistical significance in vitro.

For the key B cell genes identified as glucocorticoid responsive in the in vitro RNA-seq experiments, we then investigated the kinetics of the transcriptional response to glucocorticoid exposure over a longer period. Circulating B cells were isolated from a new cohort of five unrelated healthy donors and incubated in the presence of methylprednisolone or vehicle. RNA was purified after 4, 24, and 48 h, and expression was measured by real-time PCR. We found that the kinetics of the transcriptional response to the glucocorticoid varied by locus (Fig. 8). Most of the down-regulated genes had a nadir at the 4-h time point, with return to expression values similar to those of vehicle-treated cells by 48 h. In contrast, the strong up-regulation of IL10 and PRDM1 expression was sustained for the 48-h period of incubation.

Based on the above observations, we hypothesized that glucocorticoids would lead to a functional defect in BCR-dependent activation of human B cells following antigen stimulation. We examined this using an anti-IgM antibody to trigger B cells expressing IgM-BCRs. Given the observation that glucocorticoid exposure leads to decreased expression of immunoglobulin genes and of CD79B, which is required for expression of the BCR on the cell surface, we first assessed surface IgM-BCR levels over time in primary human B cells exposed in vitro to methylprednisolone or vehicle. This was done to ensure that if glucocorticoid exposure led to the expected reduction in the expression of surface IgM-BCR, the anti-IgM stimulation was performed at a time when sufficient levels of IgM-BCR were still present. In parallel with the 48-h incubation experiment of gene expression in B cells from five unrelated healthy donors, we measured surface expression of IgM and IgD proteins by flow cytometry after 4, 24, and 48 h of in vitro exposure to glucocorticoid or vehicle (Fig. 9, a and b). At 4 h, the mean surface IgM signal in methylprednisolone-treated cells was 84.3% (SD = 3.9%) that of vehicle-treated cells. This percentage had dropped to 42.4% (SD = 7%) at 24 h and 32.1% (SD = 4.1%) at 48 h (Fig. 9 b). Therefore, to test the effect of glucocorticoids on signaling downstream of the IgM-BCR, we stimulated purified circulating B cells from a new cohort of five unrelated healthy volunteers with an anti-IgM antibody after 4 h of in vitro glucocorticoid exposure, when sufficient levels of surface IgM-BCR were still present for adequate stimulation. To differentiate changes in signal-induced phosphorylation from changes in overall protein abundance induced by the glucocorticoid, we selected phospho-CD79A as the primary readout of IgM-BCR signaling (Fig. 9 c), given that the gene encoding this protein was not subject to glucocorticoid-mediated transcriptional regulation (Fig. 3 b and Fig. 7). We also tested the downstream signaling molecule phospho-PLCγ2, which was down-regulated at the transcript level in vitro (Fig. 3 b). We combined phosphoprotein measurements with surface staining (Fig. 9 d), to identify possible differences in the effect of glucocorticoids on IgM-BCR signaling in the two main subsets of circulating human B cells that express this BCR: naive B cells (CD19⁺IgD⁺CD27⁻) and unswitched memory B cells (CD19⁺IgD⁺CD27⁺). IgM-BCR signaling after anti-IgM stimulation was reduced in methylprednisolone-treated cells when compared with vehicle-treated cells, as measured by both phospho-CD79A (Fig. 9 e) and phospho-PLCγ2 (Fig. 9 f). For both readouts, the observed reduction in BCR signaling in total B cells appeared to be primarily driven by a stronger and significant effect on unswitched memory B cells (Fig. 9, e and f). Of note, the glucocorticoid-induced decrease in surface IgM seen in Fig. 9 b was greater on naive B cells than on unswitched memory B cells, although the difference was not statistically significant (P = 0.069).

TLR signaling is involved in critical crosstalk with BCR signaling in the activation of these lymphocytes (Suthers and Sarantopoulos, 2017), and variation in expression of intracellular TLRs is linked to development of autoimmunity (Pisitkun et al., 2006; Fairhurst et al., 2008). We therefore examined our data set specifically for glucocorticoid-induced changes in TLRs and associated signaling molecules. As predicted from the pathway-level map, the transcriptional effect of the glucocorticoid on genes whose products are related to TLR signaling is stronger in hematopoietic cells (Fig. 3 c and Fig. S5). The known glucocorticoid-induced increases in transcript abundance of key negative regulators downstream of TLRs, such as DUSP1 (Imasato et al., 2002), IL1RL1 (Díaz-Jiménez et al., 2017), IRAK3 (Miyata et al., 2015), NFKBIA (Scheinman et al., 1995), and TNFAIP3 (Altonsy et al., 2014), were evident in our data. In B cells, we observed a significant decrease in transcript abundance for multiple TLRs, which was strongest for TLR1, TLR6, and TLR7. We focused on TLR7 given the magnitude of the glucocorticoid effect and the strong association between TLR7 copy number and the development of autoimmunity (Pisitkun et al., 2006; Deane et al., 2007; Fairhurst et al., 2008; García-Ortiz et al., 2010; Suthers and Sarantopoulos, 2017; Souyris et al., 2018). We observed significant down-regulation of B cell TLR7 expression in vivo at 2 and 4 h after methylprednisolone administration (Fig. 10 a), suggesting that this would result in a functional defect in TLR7 signaling. To examine this possibility, we separately stimulated purified circulating total B cells (Fig. 10 b) from each of six unrelated healthy volunteers with the TLR7 ligand imiquimod, in the presence of methylprednisolone or vehicle. To assess whether any observed effect of the glucocorticoid was selective for TLR7, and thus likely related to the transcriptional effect on the TLR7 gene, we stimulated the cells in parallel with motolimod, a ligand for TLR8, which has similar downstream signaling properties to TLR7 but was not transcriptionally down-regulated by the glucocorticoid (Fig. 3 c). Using phospho-p38 MAPK as a readout, stimulation of total B cells with either TLR ligand led to a bimodal distribution of the signal (Fig. 10 c), suggesting that at any given concentration of ligand, a certain proportion of the total B cells has a measurable response to TLR stimulation. We generated a dose–response curve and found that the proportion of cells that respond to TLR7 or TLR8 stimulation increased as the ligand concentration increased, until a saturating concentration was reached (Fig. 10 d). At saturating or higher concentrations of the ligand, the proportion of cells that respond to the TLR7 stimulus was significantly lower in glucocorticoid-treated than in vehicle-treated cells. In contrast, the proportion of cells that responded to the TLR8 stimulus was not significantly different between glucocorticoid-treated and vehicle-treated cells (Fig. 10, d and e). The dose–response curve for phospho-p38 MAPK signal showed a similar result: at saturating or higher concentrations of the ligand, phospho-p38 MAPK signal in response to TLR7 stimulation was significantly lower in glucocorticoid-treated than in vehicle-treated cells, whereas no significant difference was observed in response to TLR8 stimulation (Fig. 10, f and g). To extend these findings, we measured expression of CCL3 and CCL4, two genes that are responsive to TLR7 and TLR8 stimulation but whose expression is unaffected by glucocorticoid treatment of B cells. Consistent with our prior observations, glucocorticoids significantly decreased the transcriptional response of CCL3 (Fig. 10 h) and CCL4 (Fig. 10 i) after TLR7 stimulation, whereas no significant difference was observed after TLR8 stimulation.

Our results illustrate the ability of a functional genomics approach to uncover previously undescribed cell type–dependent transcriptional responses that can be linked to specific immunoregulatory actions of glucocorticoids.

The strong cell type dependence of the transcriptional response has important implications for our understanding of glucocorticoid immunoregulation. For decades, nearly all the work with human material in this field has been performed with immortalized or cancer cell lines or with mixed primary cell populations such as PBMCs. While this has provided valuable insight into the molecular biology of GR signaling and GR–DNA interactions, our results suggest that additional information at the level of individual cell types is needed to develop a more accurate understanding of the effects of glucocorticoids on the immune system, as observations made in one cell type are unlikely to apply to others. Our work broadens the range of human cells in which glucocorticoid responses have been studied and focuses on highly purified populations of human primary cells, which are more likely to yield medically relevant information.

In principle, the observed cell type dependence of the transcriptional response to glucocorticoids could be explained by differences in the baseline: cells have different transcriptomes, and therefore, they have different transcriptional responses to a stimulus. However, we found that the magnitude and direction of transcriptional regulation by glucocorticoids can differ across cell types, even for genes that are expressed at similar levels at baseline. This suggests that classifying genes as glucocorticoid-induced or glucocorticoid-repressed, as is common in the literature, should only be done in the context of a specific cell type. Such findings also suggest that baseline expression, chromatin accessibility, and GR binding can only partially explain how a gene responds to a glucocorticoid stimulus. The mechanisms that allow genes with similar levels of baseline expression to respond differently to glucocorticoids across cell types require additional study, and they are likely specific to particular cell types and loci.

Many of the medically relevant effects of glucocorticoids are known to exhibit phenotypic differences among ontogenetically related cell types or tissues. These include the preferential induction of apoptosis in certain cell types (Weinstein et al., 1998; Kim et al., 2001; Herold et al., 2006), opposing effects on circulating human neutrophil and eosinophil counts (Hills et al., 1948; Dale et al., 1975), and variable effects on adipose tissue depending on its type and anatomical location (John et al., 2016). Some of these differences may be the result of variation in the levels of the type II isozyme of 11β-hydroxysteroid dehydrogenase (11βHSD2) or of the β isoform of the GR (GRβ) across cell types, both of which are known to decrease the overall sensitivity of cells to glucocorticoids (Funder et al., 1988; Lu et al., 2007). Our results suggest that these phenotypic differences could also be explained by cell type–dependent transcriptional responses, some of which are likely unrelated to the differences in overall sensitivity to glucocorticoids across different cell types. Thus, integrating information obtained from our functional genomics approach with phenotypic differences from relevant primary cell types could lead to a better understanding of the connection between transcriptional and phenotypic outputs in response to glucocorticoids.

Integrating the transcriptome data sets, we have generated a pathway-level map of the effects of glucocorticoids across nine primary human cell types. We observed substantial differences in the way that glucocorticoids affect specific pathways in individual cell types and in hematopoietic versus nonhematopoietic cells. In pathway modules with stronger enrichment in hematopoietic cells, the predominant effect of the glucocorticoid was down-regulation of genes, while in nonhematopoietic cells the predominant effect was up-regulation. It is clear from previous work that glucocorticoid-mediated immunosuppression goes beyond a simple model of down-regulation of immune effector genes, as increased expression of negative regulators of immune activation is also common and is evident in the results presented here. That said, our pathway-level analysis suggests that the characteristics of GR signaling and activation could be broadly different between hematopoietic and nonhematopoietic cells and that this could contribute to the contrasting directionality trends in gene expression.

We demonstrate the utility of this approach in uncovering relevant functional effects of glucocorticoids by focusing on B cells, which play a key role in antibody-mediated autoimmune diseases. Although glucocorticoids are essential components of the treatment of this group of diseases, it has been pointed out recently that surprisingly little is known about the effect of glucocorticoids on mature human B cells (Cain and Cidlowski, 2017). Crosstalk between the BCR and TLR signaling pathways plays a central role in the integration of functional B cell responses (Rawlings et al., 2012) and in the development of autoreactive B cells (Suthers and Sarantopoulos, 2017). Our findings indicate that pharmacologic doses of glucocorticoids functionally impair both signaling pathways, through rapid transcriptional effects on key genes.

With respect to antigen receptor signaling, in vivo or in vitro exposure of B cells to glucocorticoids led to significant down-regulation of the genes encoding two of the three components of the BCR complex (surface immunoglobulin and Igβ); the B cell co-receptor CR2 (CD21); and the upstream kinases BLNK, BTK, and LYN. There was reduced IgM-BCR signaling following stimulation with anti-IgM as early as 4 h after glucocorticoid exposure, when levels of surface IgM were still sufficient to allow response to anti-IgM stimulation of glucocorticoid-treated cells, as evidenced by the stronger effect of the treatment on naive than unswitched memory B cells for surface expression of IgM and the reverse effect for BCR signaling. The kinetics of the transcriptional effect varied by locus. For most of the BCR signaling genes in which transcript abundance decreased in response to glucocorticoids, this reduction occurred rapidly, being evident after 2 h of in vivo or in vitro glucocorticoid administration. In this group of genes, serial sampling after in vitro treatment suggests that the nadir of expression occurs between 4 and 24 h, with return to expression values similar to those of vehicle-treated cells by 48 h. The magnitude of the transcriptional response was also locus dependent but remarkably consistent across subjects and when comparing the two methods used in the different cohorts (RNA-seq and real-time PCR). For example, the nadir of expression for CD79B had a mean fold change value of 0.38 (62% lower in glucocorticoid-treated than in vehicle-treated cells at that time point), whereas the nadir of expression for BTK had a mean fold change value of 0.68 (32% lower in glucocorticoid-treated cells). It seems likely that the early blunting of BCR signaling that we have documented results from multiple rapid effects on the abundance of key upstream components of the signaling system, which individually may not have been sufficiently strong to have a measurable effect. In clinical practice, glucocorticoids are often given as a daily dose, which may over time potentiate the transcriptional and functional effects on B cells that we have observed in this study of the early response.

At 24 and 48 h after in vitro exposure, we found a striking difference in surface IgM-BCR and IgD-BCR between glucocorticoid-treated and vehicle-treated cells. By 24 h, B cells had ∼60% less surface IgM and 50% less surface IgD than vehicle-treated cells, and the difference was even more pronounced among the cells that were still viable at 48 h. In the context of this study, the purpose of looking at surface IgM-BCR over time was to identify a kinetic window for investigating signaling downstream of the receptor. Because we stimulated the cells with an anti-IgM antibody, it was important to ensure that the stimulation experiments were performed at a time when sufficient IgM-BCR was still present on the cell surface to allow adequate stimulation. However, the striking difference in receptor abundance between vehicle-treated and glucocorticoid-treated cells could have functional consequences that go beyond the blunting of BCR signaling that we documented after short-term glucocorticoid exposure. It seems likely that this difference results from a combination of the transcript-level effects on IGHM and CD79B that we identified, and some level of glucocorticoid-induced cell death. The latter may affect distinct B cell subsets in a differential manner both in vivo and in vitro, and further study of this possibility is certainly warranted, though it is a technically fraught subject due to effects of apoptosis on marker expression and to in vivo glucocorticoid-induced changes in cell migration and, hence, representation in the blood.

Glucocorticoids also induced a transcriptional program that can be predicted to impair TLR signaling by decreased expression of receptor genes and increased expression of known negative regulators. These changes include early down-regulation of TLR7 expression, which results in selective impairment of signaling through that receptor, as revealed by follow up experimental assessment. TLR7 signaling is involved in the recognition of microbial-derived nucleic acids but also in the activation of B cells and plasmacytoid dendritic cells (PDCs) by endogenous immune complexes containing nucleic acids. This results in the production of anti-nuclear antibodies and type I interferon, both important to the pathogenesis of systemic lupus erythematosus. High doses of glucocorticoids have been shown to normalize an interferon-α signature in cells from patients with systemic lupus erythematosus, which is correlated with a reduction in PDCs (Guiducci et al., 2010). Our results suggest glucocorticoid-induced down-regulation of TLR7 in B cells as a complementary mechanism. Taken together, our data indicate that glucocorticoids can functionally constrain TLR-driven, BCR-activated human B cells and that this effect may play an important role in their therapeutic effects in autoimmune diseases.

Glucocorticoid-mediated down-regulation of TLR7 expression and functional impairment of TLR7 signaling could also have implications on gender differences in the incidence of autoimmune diseases. Most human autoimmune diseases have a strong gender bias, with higher incidence in females. Human TLR7 resides in the X chromosome and has been proposed as a candidate gene for the observed gender bias. This idea has received additional support from the recent demonstration that TLR7 escapes X-inactivation in human B cells, monocytes, and PDCs, so that a proportion of each cell type in females has two effective copies of the gene (Souyris et al., 2018). While we only studied pharmacologic concentrations of glucocorticoids, it is intriguing to speculate that physiological concentrations may also play a role in restricting TLR7 expression in human immune cells, thus modulating the threshold for autoimmunity.

An important technical aspect of our assessment of intracellular TLR signaling in primary human B cells was the selection of ligand concentrations. The EC50 values that are often cited for the TLR7 ligand imiquimod and the TLR8 ligand motolimod, and the estimates of the relationship between ligand dose and receptor specificity, are based on experiments performed on HEK293 cells transfected with the respective human TLR gene and with a plasmid containing a reporter gene under the control of an NFκB promoter (Shukla et al., 2010; Lu et al., 2012). We chose not to make assumptions about the direct applicability of the dose–response characteristics and saturation dynamics of such a system to our upstream signaling assay in primary B cells, performing instead a dose–response curve for each ligand and displaying the results of the full range of concentrations tested. To ensure that the assay could be sensitive to modest changes in receptor expression, it was important to use saturating concentrations of each ligand, selected from the respective dose–response curve.

In addition to BCR signaling and the levels of surface immunoglobulins, the observed drop in transcript levels of the three human immunoglobulin loci induced by glucocorticoids could also affect the production of secreted antibodies. Clinical administration of glucocorticoids has been shown to result in lower levels of plasma immunoglobulins, except for IgE and IgG4 (Posey et al., 1978; Klaustermeyer et al., 1992; Akdis et al., 1997). Whether the clinical observations of previous studies are indeed related to the transcriptional effects we have observed here remains to be investigated. Specifically, it will be important to establish whether the differences in immunoglobulins reflect selective effects of glucocorticoids on B cell subpopulations, including long-lived plasma cells in the bone marrow or B cells in secondary lymphoid tissues. Such selective effects would be consistent with our observation of a stronger effect of methylprednisolone on unswitched memory than on naive B cells in response to BCR stimulation. In this regard, our findings of a selective effect on unswitched memory B cells are also consistent with the increased susceptibility to nontuberculous mycobacterial infections that have been independently linked to selective IgM deficiencies and use of glucocorticoids (Hojo et al., 2012; Gharib et al., 2015). Finally, given the evidence that unswitched memory B cells are critical for longevity of protective immunity and can generate new germinal center reactions (Pape et al., 2011), the latter being important for responses to mutating pathogens such as influenza, a closer assessment of glucocorticoid use in these contexts is also warranted.

Our results also indicate that B cell–intrinsic transcriptional up-regulation of IL10 could represent another mechanism of glucocorticoid action in autoimmune diseases. IL10 is a known target of glucocorticoids in macrophages, but we also observed increased expression in response to glucocorticoid treatment in B cells, both in vivo and in vitro. IL-10 secretion is an important mechanism by which B cells exert immunoregulatory functions, and circulating IL-10–producing regulatory B cells have been found to be significantly lower in patients with rheumatoid arthritis than in healthy controls (Bankó et al., 2017).

Another finding that deserves additional study is the observation that glucocorticoids lead to up-regulation of PRDM1, which encodes BLIMP-1 and was induced nearly fivefold within 2 h of methylprednisolone administration in our in vivo study. BLIMP-1 plays a central role in the terminal differentiation of B cells into plasma cells (Yu and Lin, 2016). The functional consequences of increased PRDM1 expression are unclear, but it is possible that glucocorticoids induce a transcriptional program of terminal differentiation and reduced proliferation via increased PRDM1 expression and related transcriptional programs, while dampening the responsiveness of the cells and reducing their antibody-producing capacity via the mechanisms described above.

Beyond immunoregulation, our findings could offer insights into the mechanisms behind glucocorticoid actions in human B cell malignancies. For decades, glucocorticoids have been a constant in the treatment regimens for multiple myeloma (Burwick and Sharma, 2018) and B cell lymphomas (National Guideline Alliance, 2016; Chaganti et al., 2016; Burwick and Sharma, 2018). However, their specific modes of action in this group of disorders remain as unclear as those behind their immunoregulatory effects. It is interesting to note, in the context of our findings, that some B cell lymphomas are known to be strongly dependent for survival on nonantigen-mediated BCR signaling (Dunleavy et al., 2018) and that recent whole-genome sequencing studies in multiple myeloma have suggested the possibility of PRDM1 having a tumor-suppressor role in that condition (Bolli et al., 2018). Studies in primary tumor cells will be necessary to establish the extent to which our observations of the transcriptional effects of glucocorticoids in circulating B cells from healthy volunteers are also seen in specific B cell malignancies.

Taken together, these findings provide new mechanistic understanding of glucocorticoid action and emphasize the multifactorial, cell type–dependent effects of this class of drugs. Glucocorticoids remain a mainstay of therapy in a broad range of conditions, despite the well-recognized adverse effects that high doses and prolonged use engender and despite a large and rapidly growing number of other drugs with immunotherapeutic potential. Identifying the drugs or drug combinations that best mimic the clinically beneficial effects of glucocorticoids on relevant cell types while limiting off-target action on nonhematopoietic cells will be facilitated by greater insight into the mechanisms underlying desirable versus untoward actions of glucocorticoids in distinct human cell types. The pathway-level map of glucocorticoid effects across nine primary human cell types reported here and the rich transcriptome data on which it is built provide an initial step toward achieving this knowledge.

Human peripheral blood hematopoietic cells for all in vitro experiments were obtained from the Department of Transfusion Medicine at the National Institutes of Health (NIH) Clinical Center, under NIH study 99-CC-0168, Collection and Distribution of Blood Components from Healthy Donors for In Vitro Research Use, which was approved by the Clinical Center’s Institutional Review Board.

Mononuclear cell subsets were obtained by isolation of PBMCs, followed by immunomagnetic enrichment for the specific cell subset with EasySep Human cell enrichment kits (STEMCELL Technologies). PBMCs were isolated from a leukapheresis sample for biological replicate 1 of the RNA-seq experiments and from peripheral blood collected in Vacutainer EDTA tubes (Becton Dickinson; cat. no. 366643) for all other subjects. In both cases, the PBMCs were isolated by gradient centrifugation in SepMate tubes (STEMCELL Technologies; cat. no. 15460), with Ficoll-Paque PLUS (GE Healthcare Life Sciences; cat. no. 17-1440-03).

B lymphocytes and CD4⁺ T lymphocytes were isolated from PBMCs by negative selection (STEMCELL Technologies; cat. nos. 19054 and 19052, respectively). Monocytes were isolated from PBMCs by positive selection (STEMCELL Technologies; cat. nos. 18058 or 17858), to ensure inclusion of the CD14⁺/CD16⁺ fraction, which would be excluded with the use of a negative-selection kit. Neutrophils were isolated directly from whole blood by negative selection (STEMCELL Technologies; cat. no. 19666).

Immediately after isolation, and before treatment, mononuclear cells were incubated overnight in 12-well polystyrene plates (Corning; cat. no. 3512), in Advanced RPMI 1640 medium (Thermo Fisher Scientific; cat. no. 12633-012) supplemented with 1 x L-glutamine, 10 mM Hepes, and 10% autologous serum, at 37°C and 5% CO₂. Immediately after isolation, and before treatment, neutrophils were incubated for 4 h in 12-well polystyrene plates (Corning; cat. no. 3512) pre-coated with 20% autologous plasma, in RPMI 1640 medium without phenol red (Thermo Fisher Scientific; cat. no. 11835030) supplemented with 10 mM Hepes, at 37°C and 5% CO₂.

Human primary skeletal muscle myoblasts from adult, isolated from four unrelated healthy human donors (Lonza; cat. no. CC-2580; lot nos. 0000419228, 0000650386, 0000657512, and 0000583849), were cultured in SkGM-2 medium (Lonza; cat. no. CC-3245) without antibiotics or glucocorticoids. Cells were incubated at 37°C, 5% CO₂, and a growth curve was generated to ensure glucocorticoid treatment was performed in the early plateau phase of growth. At the time of treatment, skeletal myoblasts were on day 7 (lot no. 0000419228) or day 9 (lot nos. 0000650386, 0000657512, and 0000583849) of incubation and on passage 3 or 4.

Human primary subcutaneous preadipocytes from adult, isolated from four unrelated healthy human donors (Lonza; cat. no. PT-5020, lot nos. 0000399826, 0000629514, and 0000645827; Cell Applications, Inc.; cat. no. 802s-05a, lot no. 1687), were cultured in PGM-2 medium (Lonza; cat. no. PT-8002) without antibiotics or glucocorticoids. Cells were incubated at 37°C, 5% CO₂, and a growth curve was generated to ensure glucocorticoid treatment was performed in the early plateau phase of growth. At the time of treatment, preadipocytes were on day 10 (Lonza; lot no. 0000399826) or day 8 (Lonza; lot nos. 0000629514 and 0000645827; Cell Applications; lot no. 1687) of incubation and on passage 3.

Human dermal microvascular endothelial cells from adult, isolated from four unrelated healthy human donors (Lonza; cat. no. CC-2543; lot nos. 0000442486, 0000577921, 0000550175, and 0000664503), were cultured in EGM-2MV medium (Lonza; cat. no. CC-3202) without antibiotics or glucocorticoids. Cells were incubated at 37°C, 5% CO₂, and a growth curve was generated to ensure glucocorticoid treatment was performed in the early plateau phase of growth. At the time of treatment, endothelial cells were on day 5 (lot nos. 0000442486 and 0000664503), day 9 (lot no. 0000550175), or day 14 (lot no. 0000577921) of incubation and on passage 4 or 5.

Human primary dermal fibroblasts from adult, isolated from four unrelated healthy human donors (Lonza; cat. no. CC-2511; lot nos. 0000409270, 0000509796, 0000540991, and 0000545147), were cultured in FGM-CD medium (Lonza; cat. no. 00199041) or HyClone Minimum Essential Medium (Thermo Fisher Scientific; cat. no. SH30265.FS). The medium was supplemented with 10% FBS, without antibiotics or glucocorticoids. Cells were incubated at 37°C, 5% CO₂, and a growth curve was generated to ensure glucocorticoid treatment was performed in the early plateau phase of growth. At the time of treatment, fibroblasts were on day 7 of incubation and on passage 2 or 3.

Human primary osteoblasts from adult (Lonza; cat. no. CC-2538, lot no. 0000435102) or child (Lonza; cat. no. CC-2538, lot nos. 0000336963 and 0000426160; iXCells Biotechnologies; cat. no. 10HU-179, lot no. 200211), isolated from four unrelated healthy human donors, were cultured in OGM medium (Lonza; cat. no. CC-3207) without antibiotics or glucocorticoids. Cells were incubated at 37°C, 5% CO₂, and a growth curve was generated to ensure glucocorticoid treatment was performed in the early plateau phase of growth. At the time of treatment, osteoblasts were on day 7 of incubation and on passage 3 or 4.

Purity and viability of the enriched immune cell populations was assessed by flow cytometry. Cells were stained in PBS containing 1% BSA. Surface staining was performed with a panel containing the following monoclonal antibodies: CD45 clone HI30 (BD Biosciences; cat. no. 564357), CD66b clone G10F5 (BioLegend; cat. no. 305104), CD16 clone 3G8 (Beckman Coulter; cat. no. A33098), CD14 clone M5E2 (BD Biosciences; cat. no. 561390), CD3 clone HIT3a (BioLegend; cat. no. 300312), CD4 clone OKT4 (BioLegend; cat. no. 317417), CD8 clone SK1 (BioLegend; cat. no. 344712), and CD19 clone HIB19 (BioLegend; cat. no. 302210). Cell viability was assessed with the LIVE/DEAD Fixable Blue Dead Cell Stain Kit (Thermo Fisher Scientific; cat. no. L-23105). Evidence of early apoptosis was assessed by fluorochrome-labeled Annexin V staining (BioLegend; cat. no. 640908). Spectral compensation was performed with UltraComp eBeads (Thermo Fisher Scientific; cat. no. 01-2222-42). Data acquisition was performed with a BD Biosciences LSR Fortessa flow cytometer. Flow cytometry data were analyzed with FlowJo X software. Purity was defined as the proportion of cells with the specified cell-lineage markers, out of the total CD45⁺ cells in the purified cell preparation. Cell lineages were classified as follows: CD66b⁺/CD16⁺ (neutrophils), CD14⁺/CD66b^– (monocytes), CD3⁺/CD4⁺/CD8^– (CD4⁺ T cells), and CD19⁺/CD3^– (B cells).

The choice of methylprednisolone was based on two facts. First, this glucocorticoid is commonly used in clinical practice when rapid immunosuppression is required. Second, it has dose-linear pharmacokinetics, which makes it easier to estimate plasma concentrations from the doses that are commonly used in the clinic (Möllmann et al., 1989; Derendorf et al., 1991). For the initial in vitro studies, in which cells were treated with methylprednisolone and the transcriptional response was studied by RNA-seq, we chose to treat the cells in vitro with a concentration of 22.7 µM, which was estimated to be equivalent to the peak plasma concentration after an intravenous dose of 1 g (a dose commonly used in acute presentations of autoimmune diseases). When we performed the in vivo studies of glucocorticoid response in circulating human B cells, the study volunteers were given a dose of 250 mg. We measured methylprednisolone levels in the plasma of each volunteer, and the mean value of the methylprednisolone concentration at 4 h was 5.34 µM, which is very close to what we had predicted based on previous data (Möllmann et al., 1989; Derendorf et al., 1991) and on the dose-linear pharmacokinetics of the drug. To bring the in vitro conditions as close as possible to those of our in vitro study, all subsequent experiments used a methylprednisolone concentration of 5.34 µM.

For in vitro glucocorticoid treatment and RNA-seq, nine primary human cell types were purified and cultured as described under the Cell purification and cell culture conditions for in vitro treatment section. For this set of experiments, each cell type was obtained from four unrelated healthy human donors. The experiments corresponding to each of the four biological replicates were performed independently. For each cell type and biological replicate, the cells were cultured in four wells of a 12-well plate. On the day of treatment, methylprednisolone 22.7 µM (Sigma; cat. no. M0639) was added to two of the wells and vehicle (ethanol, 0.08%) to the other two. At 2 and 6 h after the stimulus, cells were harvested from one of the wells that received vehicle and one of the wells that received methylprednisolone. Immediately after harvesting, the cells originating from each well were independently separated from the supernatant by centrifugation, resuspended in 500 μl of TRIzol Reagent (Thermo Fisher Scientific; cat. no. 15596018), and stored at −80°C until the time of RNA purification. All downstream processing steps (RNA purification, RNA-seq library preparation, and sequencing) were performed separately for each well.

Total RNA was isolated by extraction with TRIzol Reagent (Thermo Fisher Scientific; cat. no. 15596018), followed by column-based purification with the RNA Clean & Concentrator-5 kit (Zymo Research; cat. no. R1016). RNA quantity was measured on a Qubit 2.0 fluorometer (Thermo Fisher Scientific; cat. no. Q32866), with RNA BR quantitation assays (Thermo Fisher Scientific; cat. no. Q10211). RNA quality was assessed by microfluidic electrophoresis on an Agilent 2100 Bioanalyzer system (Agilent; cat. no. G2939A), with RNA 6000 Nano chips (Agilent; cat. no. 5067-1511).

Sequencing libraries were prepared with the TruSeq Stranded Total RNA with Ribo-Zero Gold kit (Illumina; cat. no. RS-122-2303). Indexed libraries were normalized and pooled. For biological replicate 1, the cBot system (Illumina; cat. no. SY-301-2002) was used for paired-end cluster generation with a TruSeq PE Cluster Kit v3-cBot-HS (Illumina; cat. no. PE-401-3001). Paired-end sequencing (2 × 94 bp) was performed on an Illumina HiSeq 2000 sequencer (Illumina; cat. no. SY-401-1001), with the TruSeq SBS v3-HS kit (Illumina; cat. no. FC-401-3001). For biological replicates 2–4, paired-end cluster generation was performed with the HiSeq 3000 PE Cluster Kit (Illumina; cat. no. PE-410-1001). Paired-end sequencing (2 × 75 bp) was performed on an Illumina HiSeq 3000 sequencer (Illumina; cat. no. SY-401-3001), with the HiSeq 3000 SBS kit (Illumina; cat. no. FC-410-1001).

Illumina base call (.bcl) files were converted to FASTQ format with bcl2fastq2 v.2.20 (Illumina, Inc.). Adapter sequences were trimmed with Cutadapt v.1.10 (Martin, 2011) in Python v.2.7.9, using the following adapter sequences as input: Read 1: 5′-AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC-3′; Read 2: 5′-AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT-3′. Adapter-trimmed reads under 20 bp were discarded. The adapter-trimmed FASTQ files were aligned to the reference human genome assembly (GRCh38) with STAR v.2.6 (Dobin et al., 2013). The transcript annotation (GTF) file was obtained from GENCODE, release 28 (Harrow et al., 2012). The binary alignment files (.bam) were then used for generation of a matrix of read counts with the featureCounts program of the package Subread v.1.5.3 (Liao et al., 2014). Paired-end exonic fragments were grouped at the level of genes, based on the GENCODE 28 annotation file.

Normalization and differential expression analysis for RNA-seq data were performed with the DESeq2 package v.1.18.1 (Love et al., 2014) in R v.3.5.0 (R Development Core Team, 2009). Our choice of this widely used statistical package was based on the fact that it provides appropriate tools to effectively model differential expression in RNA-seq data by accounting for features like non-normal distribution of data and overdispersion (increasing variance as the mean count increases). In the specific case of our study design, in which we compared data from methylprednisolone-treated versus vehicle-treated cells in each of four subjects, it was also important to choose a method that allows a paired analysis.

To effectively account for discreteness and the dependence of the variability on the mean, DESeq2 employs a generalized linear model, where the read count in gene i from sample j is assumed to follow a negative binomial distribution with mean S_ijq_ij and variance as a function of the mean and a dispersion parameter. The normalizing factor S_ij accounts for differences in sequencing depth between samples. The quantity q_ij depends on the gene and the biological condition of the sample, and it is related to these through the loglinear model

In the case of our analysis comparing methylprednisolone-treated versus vehicle-treated cells at each time point, with vehicle-treated cells as the reference, r takes the values 0 or 1: X_j0 = 1, X_j1 = 0 for the vehicle-treated cells and X_j0 = 0, X_j1 = 1 for the methylprednisolone-treated cells. The log fold change in gene i is reflected by β_i1 − β_i1. To improve efficiency in assessing the effect of methylprednisolone treatment, where gene expression is measured in cells from the same subjects treated with methylprednisolone or vehicle, we performed a paired analysis by adding a subject effect into the above negative binomial model, following the method of Love et al. (2014).

Another important consideration in the analysis of RNA-seq data is that it relies on count data. In that context, ratios (such as the log fold change that is often used to compare gene expression in a sample against that in the reference) are inherently noisier when the read counts are low. This can result in differential expression estimates that are much stronger for weakly expressed than for highly expressed genes. DESeq2 overcomes this limitation by using an empirical Bayes procedure, whereby the log fold change estimates are shrunk toward zero in a manner such that shrinkage is stronger when the available information for a gene is low as a result of low read count, high dispersion, or a small sample size (Love et al., 2014). The empirical Bayes estimates of log fold change are used for differential expression analysis with the Wald test: for each gene, a z statistic is calculated as the empirical Bayes estimate of log fold change in that gene divided by its standard error, which is compared with a standard normal distribution. This is followed by multiple-testing adjustment with the method of Benjamini and Hochberg (1995).

To functionally interpret the results of differential expression analysis, we tested for enrichment of differentially expressed genes across functional gene sets corresponding to annotated molecular pathways and biological processes. From the Broad Institute’s MSigDB (Subramanian et al., 2005), version 6.2, we downloaded gene set definitions from three groups: H, hallmark gene sets; C2/CP, canonical pathways, which include gene sets from KEGG, Biocarta, Pathway Interaction Database, Reactome, Signaling Gateway, Signal Transduction KE, SuperArray, and Sigma-Aldrich; and C5/BP group, Gene Ontology Biological Processes.

We then applied a gene set enrichment test, implemented in the R package limma (Smyth, 2005) as function geneSetTest. This test first ranks all genes by a differential expression statistic, in this case, log2 fold change in glucocorticoid-treated versus vehicle-treated cells. For each gene set, it then tests whether the genes in the set tend to be located near the top of the ranked list of differentially expressed genes. We used the ranks of the differentially expressed genes [ranks.only = TRUE], in which case geneSetTest calculates the enrichment P values by a Wilcoxon test. Multiple-testing correction was performed separately for each category of gene sets (GO, KEGG, Reactome, etc.) with the method of Benjamini and Hochberg (1995). We used geneSetTest instead of the classical hypergeometric tests to avoid potential noise associated with setting a cutoff for calling differential expression.

For the in vivo glucocorticoid administration experiments, healthy volunteers were enrolled in NIH study 16-I-0126, Genomic Responses of Human Immune and Non-Immune Cells to Glucocorticoids (ClinicalTrials.gov Identifier NCT02798523), which was approved by the Institutional Review Board of the National Institute of Allergy and Infectious Diseases (NIAID). Informed consent was obtained from each subject before enrollment. On the day of the screening visit (day −30 to day −1), each subject underwent a medical history, physical examination, and baseline blood collection. Subjects who passed the clinical and laboratory criteria for enrollment were scheduled for a glucocorticoid infusion visit. 20 healthy volunteers were deemed eligible for the study and subsequently enrolled. 10 of the 20 volunteers were female. The mean age was 37 yr. On the day of the glucocorticoid infusion visit (day 0), volunteers underwent a targeted history and physical examination. Each volunteer then received a single intravenous dose of 250 mg of methylprednisolone sodium succinate (SOLU-MEDROL; Pfizer, Inc.) over 30 min. Blood was collected 2 and 4 h after the start of the infusion. Blood was collected in Vacutainer EDTA tubes (Becton Dickinson; cat. no. 366643) and immediately processed for B cell immunomagnetic isolation, which was performed as described in the Cell purification and cell culture conditions for in vitro treatment section. Immediately after isolation, purified B cells were lysed in TRIzol reagent (Thermo Fisher Scientific; cat. no. 15596018) and stored at −80°C. RNA purification was performed as described in the RNA purification and quality control section. RNA-level validation of the RNA-seq findings was then performed on these samples, as described in the Real-time PCR section.

After the glucocorticoid infusion visit, each study participant received two follow up phone calls, on day 1 and day 5. Total participation time was 1–5 wk.

Human peripheral blood from healthy donors was obtained in Vacutainer EDTA tubes (Becton Dickinson; cat. no. 366643) and immediately processed for B cell immunomagnetic isolation, which was performed as described in the Cell purification and cell culture conditions for in vitro treatment section. Isolated B cells were incubated overnight at 37°C, 5% CO₂, in RPMI 1640 with 10% FBS, then treated once with methylprednisolone (5.34 µM) or vehicle (0.08% ethanol). Immunophenotyping was performed after 4, 24, and 48 h of in vitro treatment. Phosphorylation assays were performed after 4 h of in vitro treatment.

For B cell immunophenotyping, multicolor flow cytometry analyses of B cell subpopulations were performed using the following monoclonal antibodies (mAbs): CD3 clone OKT3 (BioLegend; cat. no. 317324), CD10 clone HI10a (BD Biosciences; cat. no. 340923), CD19 clone SJ25C1 (Thermo Fisher; cat. no. 45-0198-42), CD27 clone L128 (BD Biosciences; cat. no. 563815), CD27 clone O323 (Thermo Fisher; cat. no. 25-0279-42), IgD clone IA6-2 (BioLegend; cat. no. 348232), IgD clone IA6-2 (BD Biosciences; cat. no. 555778), and IgM clone MHM-88 (BioLegend; cat. no. 314517). Cell viability was assessed with the Zombie Aqua Fixable Viability Kit (BioLegend; cat. no. 423102). Flow cytometry data acquisition was performed on a LSR Fortessa flow cytometer (BD Biosciences), and data were analyzed using FlowJo software.

Five unrelated healthy adult volunteers were studied, each on a different day. Peripheral blood B cells were isolated as described in the Cell purification and cell culture conditions for in vitro treatment section. For phosphorylation assays, glucocorticoid- or vehicle-treated B cells were stained with mAbs against CD19, CD27, CD10, and IgD, as detailed in the Assessment of B cell surface markers section, then stimulated with 10 µg/ml goat F(ab′)₂ anti-human IgM (Jackson ImmunoResearch Laboratories) at 37°C for 2 min. For the detection of phosphorylated signaling intermediates, cells were fixed and permeabilized using BD Cytofix and Phosflow Perm/Wash buffers (BD Biosciences) and stained separately with a PE-conjugated mAb against phosphorylated PLC-γ2 (pY759) clone K86-689.37 (BD Biosciences), or an unconjugated polyclonal antibody against phosphorylated CD79A(Tyr182) (Cell Signaling Technology; cat. no. 5173) followed by staining with secondary goat anti-rabbit IgG-PE (Thermo Fisher; cat. no. P-2771MP). Flow cytometric analyses were performed as described above.

Six unrelated healthy adult volunteers were studied, each on a different day. Peripheral blood B cells were isolated as described in the Cell purification and cell culture conditions for in vitro treatment section. For phospho-p38 MAPK staining, B cells (2 × 10⁵ per well) were cultured in RPMI 1640, 10% FBS at 37°C, 5% CO₂, in 96-well U-bottom plates. Cells were treated overnight with methylprednisolone (5.34 µM) or vehicle (0.08% ethanol), then stimulated for 15 min at 37°C with either the TLR7 ligand imiquimod-HCl, molecular weight 276.8 g/mol (InvivoGen; cat. no. tlrl-imqs), or the TLR8 ligand motolimod, molecular weight 458.60 g/mol (MedChem Express; cat. no. HY-13773). After TLR stimulation, B cells were fixed by addition of paraformaldehyde (Electron Microscopy Sciences; cat. no. 15710) to a final concentration of 1.6%, followed by a 15-min incubation at room temperature (RT). Dose–response curves for imiquimod were generated by stimulating cells at five ligand concentrations: 0, 5, 10, 20, and 40 µg/ml (0, 18, 36.1, 72.2, and 144.4 µM). Dose–response curves for motolimod were generated by stimulating cells at six ligand concentrations: 0, 1.25, 2.5, 5, 10, and 20 µg/ml (0, 2.7, 5.5, 10.9, 21.8, and 43.6 µM). After one wash with PBS containing 1% FBS and 2 mM EDTA, samples were permeabilized using ice-cold methanol (Fisher Chemical; cat. no. A452-1) for 30 min at 4°C, blocked using 5% goat serum and Fc receptor specific antibody Trustain FcX (BioLegend; cat. no. 422302) for 15 min at RT, and stained for 1 h at RT with Alexa Fluor 647 mouse anti-p38 MAPK (pT180/pY182) clone 36 (BD Phosflow; cat. no. 612595). Data were acquired on a BD Fortessa flow cytometer and analyzed using FlowJo software.

For CCL3 and CCL4 expression, cells were treated overnight with methylprednisolone (5.34 µM) or vehicle (0.08% ethanol) in RPMI 1640, 10% FBS at 37°C, 5% CO₂, then stimulated with either the TLR7 ligand imiquimod-HCl (40 µg/ml; 144.4 µM) or the TLR8 ligand motolimod (20 µg/ml; 43.6 µM) for 1 h at 37°C. Cells were pelleted, then resuspended in 500 μl of TRIzol Reagent and stored at −80°C until the time of RNA purification.

Two types of variable were used in the statistical analysis of flow cytometry data. For the quantification of surface proteins and phosphorylated intracellular proteins, the variable was the mean fluorescence intensity of the protein or phosphoprotein being assayed. In the specific case of the analysis of the phosphoprotein signals in B cells treated with TLR7 or TLR8 agonists, where the response was bimodal, a second variable was the percentage of cells that responded to the stimulus (those in the upper mode). For both types of variable, statistical inference was performed with a paired t test. The choice of a paired t test was based on the fact that the two groups we were comparing involved cells purified from each individual and treated with either vehicle or glucocorticoid. Therefore, the two groups are not independent, and a paired design is appropriate.

cDNA was synthesized with the SuperScript IV VILO Master Mix with ezDNase Enzyme kit (Thermo Fisher Scientific; cat. no. 11766050). Enrichment for the targets of interest was performed with the Fluidigm Preamp Master Mix (Fluidigm; cat. no. 100-5581) and the same set of primers to be used for real-time PCR. TaqMan quantitative PCR assays were then performed on a BioMark HD system (Fluidigm; cat. no. BMKHD-BMKHD), with the TaqMan Fast Advanced Master Mix (Thermo Fisher Scientific; cat. no. 4444556) and the FLEXsix IFC Gene Expression Kit (Fluidigm; cat. no. 100-6309) for BLNK, BTK, CCL3, CCL4, CD19, CD79A, CD79B, CR2, IL10, LYN, PRDM1 (BLIMP-1), SYK, and TLR7.

We calculated the ΔCt values (ΔCt = Ct of the target gene − Ct of the reference gene TBP) for each subject, at each time point. For ease of visualization, the real-time PCR data in Figs. 4 c, 5 c, 7, and 8 are displayed as fold changes in gene expression: time point x over baseline for the in vivo data and glucocorticoid-treated over vehicle-treated for the in vitro data. For this, we employed the 2^–ΔΔCt method (Schmittgen and Livak, 2008), where the first Δ is the ΔCt at time point x and the second is the ΔCt for the sample that was not treated with glucocorticoid (baseline or vehicle for in vivo or in vitro data, respectively). This method allows the display of the fold difference or fold change between two conditions, which is an intuitive way to visualize and understand real-time PCR gene expression data. The error bars in the above figures correspond to the geometric mean of the fold change values and its standard error. The standard error of the geometric mean was calculated by the Delta method, as

where FC is fold change, $log 2 FC ¯$ is the sample mean of fold changes in log₂ scale, and n is number of biological replicates.

To assess the statistical significance of gene expression changes measured by real-time PCR after in vivo glucocorticoid administration, we performed paired (signed-rank) Wilcoxon tests at each time point (time point x versus baseline), with the 2^−ΔCt values as input. The 2^−ΔCt values represent the relative expression of a gene relative to the reference gene. This value can be used as input for statistical inference (Schmittgen and Livak, 2008), but normality cannot be assumed. The choice of a paired test was based on the fact that this experiment had a paired design: the two groups we were comparing involved cells from each individual, sampled before or after glucocorticoid exposure, and they are therefore not independent. The choice of a nonparametric test was based on the fact that the input values cannot be assumed to be normally distributed. The 2^−ΔCt values that were used as input for statistical testing from the in vivo data are displayed in Fig. S4.

The RNA-seq data set for this study has been uploaded to the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus, under accession no. GSE112101. It includes links to the NCBI Sequence Read Archive database, from which the raw data will be accessible in FASTQ format, under accession no. SRP136108.

Fig. S1 shows that varying the threshold for differential expression does not affect the observed cell type dependence of the transcriptional response to glucocorticoid. Fig. S2 shows a pathway-level map of the transcriptional response to glucocorticoids in nine human primary cell types. Fig. S3 shows enrichment for down-regulated or up-regulated genes among pathways affected by the glucocorticoid stimulus. Fig. S4 shows gene expression over time in key BCR signaling genes before and after in vivo treatment with methylprednisolone. Fig. S5 shows transcriptional effects of glucocorticoids on TLR signaling genes.

This study used the High-Performance Computing cluster of the NIAID Office of Cyber Infrastructure and Computational Biology. We thank Gustavo Gutierrez-Cruz at the National Institute of Arthritis, Musculoskeletal and Skin Disorders and Yuesheng Li at the DNA Sequencing and Genomics Core of the National Heart, Lung, and Blood Institute for sequence data generation. We thank the Flow Cytometry Section at the NIAID Research Technologies Branch, specifically Kevin Holmes and David Stephany, for their technical expertise and assistance. We thank Angélique Biancotto and Yuri Kotliarov at the NIH Center for Human Immunology for their help in obtaining data from previous studies of glucocorticoid action. We thank Andrew Oler and Manikandan Narayanan at NIAID for bioinformatics support. We thank Silvia Bolland at NIAID for sharing her expertise in intracellular TLR biology. We thank the volunteers who participated in the clinical study and the nursing staff at the NIH Clinical Center Day Hospital for their dedicated clinical support. We thank Derek Einhaus at the Laboratory of Immune System Biology, NIAID, for organizational support.

This research was supported by the Division of Intramural Research, NIAID.

The authors declare no competing financial interests.

Author contributions: L.M. Franco, J.S. Tsang, and R.N. Germain conceived and designed the project. L.M. Franco, K. Howe, and M. Gadkari carried out the clinical study. M. Gadkari and S. Kumari performed the in vitro experiments with primary human cells. L.M. Franco performed the data analysis with guidance from J.S. Tsang and Z. Hu. M. Gadkari, L. Kardava, S. Moir, and L.M. Franco performed the BCR signaling experiments. P. Kumar provided pharmacokinetics expertise in the clinical and in vitro studies. M. Gadkari, J. Sun, I.D.C. Fraser, and L.M. Franco performed the TLR signaling experiments. S. Kumari performed the experiments while employed at the NIH; American Type Culture Collection was not involved in the funding or performance of any part of this work. All authors participated in drafting or revising the manuscript and gave final approval of the version submitted.