Social status alters chromatin accessibility and the gene regulatory response to glucocorticoid stimulation in rhesus macaques.

Low social status is an important predictor of disease susceptibility and mortality risk in humans and other social mammals. These effects are thought to stem in part from dysregulation of the glucocorticoid (GC)-mediated stress response. However, the molecular mechanisms that connect low social status and GC dysregulation to downstream health outcomes remain elusive. Here, we used an in vitro GC challenge to investigate the consequences of experimentally manipulated social status (i.e., dominance rank) for immune cell gene regulation in female rhesus macaques, using paired control and GC-treated peripheral blood mononuclear cell samples. We show that social status not only influences immune cell gene expression but also chromatin accessibility at hundreds of regions in the genome. Social status effects on gene expression were less pronounced following GC treatment than under control conditions. In contrast, social status effects on chromatin accessibility were stable across conditions, resulting in an attenuated relationship between social status, chromatin accessibility, and gene expression after GC exposure. Regions that were more accessible in high-status animals and regions that become more accessible following GC treatment were enriched for a highly concordant set of transcription factor binding motifs, including motifs for the GC receptor cofactor AP-1. Together, our findings support the hypothesis that social status alters the dynamics of GC-mediated gene regulation and identify chromatin accessibility as a mechanism involved in social stress-driven GC resistance. More broadly, they emphasize the context-dependent nature of social status effects on gene regulation and implicate epigenetic remodeling of chromatin accessibility as a contributing factor.

Social status is an important predictor of health and survival in social animals, including humans (1–4). In the United States, for instance, men in the lowest income percentile live 15-y-shorter lifespans, on average, than men in the highest income percentile (5). Similarly, British civil servants in the lowest employment grades are twice as likely to develop cardiovascular disease compared with those in the highest employment grades, despite universal access to the British National Health Service (6). Social status gradients in health can be explained in part by correlated variation in health risk behaviors, resource access, and the long-term effects of social adversity earlier in life (7, 8). However, increasing evidence indicates that low social status also leads to direct physiological effects (9–18). For example, experimentally induced low social status in cynomolgus macaques and rhesus macaques leads to increased rates of atherosclerosis and hyperinsulinemia and an exaggerated gene regulatory response to bacterial LPS (9, 10, 19–26).

The direct effects of low status may be mediated, at least in part, by the chronic stress of social subordination (27, 28). In humans and other hierarchically organized primates, low status is associated with greater social unpredictability and reduced social control, fewer affiliative interactions, and reduced social integration (29). These experiences are thought to increase exposure to stressors (e.g., via increased rates of harassment) and reduce the capacity to buffer against stressful experiences (e.g., via lack of social support) (29). Such changes in turn lead to physiological alterations in the regulation of the hypothalamic–pituitary–adrenal (HPA) axis, which controls the production of glucocorticoids (GCs), a class of steroid hormones involved in metabolic and immune regulation (30). Chronic social stress has been associated with both elevated baseline GC levels (e.g., refs. 31 and 32, but see ref. 33) and impaired GC negative feedback (i.e., GC resistance; refs. 28 and 34), suggesting a physiological basis for social gradients in metabolic and inflammation-related disorders (35). In support of this hypothesis, the psychological stress of caring for a chronically ill child is associated with both GC resistance and impaired GC suppression of proinflammatory signaling (34). Similarly, major social stressors such as job loss or bereavement predict both GC resistance and susceptibility to the common cold (28). Causal effects of social stress on GC function have also been demonstrated in experimental animal models. For instance, red squirrel mothers exposed to cues for increased resource competition (but not increased resource competition itself) exhibit elevated GC levels, a change that in turn alters offspring growth (36).

At the molecular level, GCs operate as signaling molecules that affect transcriptional regulation, thus linking multiorgan HPA axis regulation to changes in metabolic and immune function within cells. Circulating GCs canonically influence gene regulation by diffusing across the cell membrane and binding to the cytosolic form of the nuclear receptor GC receptor (GR, encoded by the gene NR3C1). Once bound, the GC–GR complex translocates into the nucleus, where it binds directly or indirectly [via tethering to other transcription factors (TFs); ref. 37] to target DNA sequences [e.g., GC response elements (GREs); refs. 30 and 38]. There are multiple potential mechanisms through which GCs exert antiinflammatory effects, and GCs have been associated with changes in gene regulation at several levels (39–42), including (i) GR induction of IκBα, which binds to the proinflammatory TF complex NF-κΒ and prevents it from entering the nucleus, (ii) inhibition of histone acetylation at the promoters of proinflammatory genes, (iii) posttranscriptional degradation of inflammation-related mRNAs, (iv) recruitment of negative elongation factor at promoters of inflammatory genes, which inhibits activation, and (v) GR tethering to the NF-κB component p65 (RelA) at NF-κB binding sites to inhibit proinflammatory gene transcription. GC resistance, which impairs termination of proinflammatory signaling, is thus also linked to multiple regulatory mechanisms, including defective histone acetylation, increased expression of proinflammatory TFs, and disinhibition of proinflammatory signaling by the alternative GR receptor, GRβ (reviewed in ref. 43).

HPA axis signaling via GC regulation therefore provides a plausible mechanism for linking the central nervous system, which senses and processes social experiences, to gene regulation in the periphery, where social stress is translated into poor health outcomes. Indeed, observational studies of human subjects (44–50) and experimental animal models for social subordination-induced stress (9, 10) have identified a strong relationship between the chronic social stress associated with low social status and proinflammatory gene expression in peripheral immune cells. In support of a role for GC signaling, predicted GR and NF-κB binding sites are enriched near genes that are differentially expressed in association with social isolation (51), low social status (9, 49, 52), and early-life social adversity (53, 54). In addition, experimental manipulation of social status in rhesus macaques results in a more pronounced proinflammatory gene expression response to LPS exposure in low-ranking animals. Notably, TF binding sites for NF-κB are enriched in accessible chromatin near genes affected by social status, suggesting that changes in GC-related NF-κB transcriptional activity may explain the low status-driven proinflammatory phenotype (9). However, although low-status macaques exhibit some evidence of GC resistance (22–24), both the degree to which GC signaling is responsible for social status-driven gene expression patterns and the gene regulatory mechanisms involved remain unclear.

To address these gaps, here we combined a powerful nonhuman primate model for investigating the consequences of chronic social stress (9, 10) with a well-established model for investigating the gene regulatory response to GC signaling: in vitro treatment with the synthetic GC dexamethasone (Dex) (41, 55, 56). Specifically, we sequentially introduced 45 female rhesus macaques into nine new social groups of five females each, a manipulation that produces a strong relationship between introduction order and dominance rank (earlier introduction predicts higher social status/dominance rank; refs. 10 and 57). This approach allows the causal effects of social subordination-induced stress to be disambiguated from alternative explanations, while taking advantage of the natural tendency of female rhesus macaques to form stable status hierarchies (58). To test how social status affects the response to GC signaling, we then cultured peripheral blood mononuclear cells (PBMCs) from each study subject using a paired control (untreated) and Dex-treated study design (following ref. 55).

Because the majority of GR–DNA binding events are thought to occur in regions of chromatin accessible to TF binding before GR nuclear translocation (59), we measured both genome-wide gene expression and chromatin accessibility in our samples. We hypothesized that social status would not only predict variation in PBMC gene expression [as previously shown in rhesus macaques and humans (10, 49)] but also patterns of open chromatin and TF binding, which have not previously been investigated. In support of this idea, social experiences have been shown to modulate patterns of epigenetic marks and TF binding in other vertebrates [e.g., three-spined sticklebacks (60), mice (61), and prairie voles (62)]. Further, as shown for exposure to LPS (9), we hypothesized that the effects of social status on gene regulation would be conditional on GC exposure. Finally, we hypothesized that low-ranking animals, who exhibit signatures of social stress-mediated GC resistance (28, 34), would also exhibit gene regulatory patterns consistent with reduced responsiveness to GC signaling, particularly at GR- and NF-κB–regulated genes.

Results

Social Status and GC Stimulation Alter Chromatin Accessibility and Gene Expression.

We manipulated the dominance rank positions of 45 adult female rhesus macaques by sequentially introducing them into newly constructed social groups of five females each (as reported in ref. 9) (Fig. 1, , and Dataset S1). This approach produces a strong correlation between order of introduction and resulting dominance rank, in which females introduced earlier become higher-ranking [Pearson’s r = −0.67, P = 8.8 × 10−7 at time of sample collection; rank was measured using a continuous Elo rating score (63); Fig. 1] (see also refs. 9 and 57). Consistent with dominance rank dynamics typical for female rhesus macaques, dominance rank values were highly stable from group construction until sample collection.

Fig. 1.

Dominance rank and Dex treatment both induce widespread changes in gene expression and chromatin accessibility. (A) Social group formation for one study group. Each line represents a different female and starts at the date of her introduction. Females introduced earlier were higher-ranking at the time of sampling (Pearson’s r = −0.67, P = 8.8 × 10−7; n = 9 groups and 43 females). (B) Schematic of the in vitro GC challenge. Paired PBMC samples from each female were treated with either 1 μM Dex (Dex+), a synthetic GC, or 0.02% ethanol (the vehicle for Dex, Dex−) and subjected to a 90-min incubation period before RNA and DNA extraction. (C) Example response to Dex challenge for FKBP5, a known GC-induced gene, for chromatin accessibility in intron 4 of isoform FKBP5-201 (P = 7.36 × 10−18; this region also falls in intron 5 of isoform FKBP5-202) and for gene expression (P < ×10−100). (D) Dex treatment is strongly associated with the first PC of gene expression data (paired T42 = 80.61, P = 1.12 × 10−47) and (E) with the first three PCs of chromatin accessibility (PC1 T42 = −3.85, P = 3.95 × 10−4; PC2 T42 = 4.49, P = 5.44 × 10−4; PC3 T42 = −4.93, P = 1.35 × 10−5). Social status is significantly associated with PC3 for both (F) gene expression (Pearson’s r = 0.77, P = 3.33 × 10−18) and (G) chromatin accessibility (Pearson’s r = 0.48, P = 3.41 × 10−6).

We isolated PBMCs from 43 females (all members of seven social groups and four of the five members of the remaining two groups; ). We divided each sample into two aliquots of 1 million cells each and simultaneously treated one aliquot with 1 μM Dex (the “Dex+” sample) and the second aliquot with 0.02% EtOH, the vehicle for Dex (the “Dex−” control sample). Following a 90-min incubation at 37° C (an incubation period that results in a strong gene expression response to Dex; ref. 55), we extracted RNA and intact nuclei from each sample for RNA sequencing (RNA-seq) and assay for transposase-accessible chromatin using sequencing (ATAC-seq) analysis, respectively (64, 65). After sequencing and quality control, we retained all 43 control and 43 Dex+ samples for both gene expression and chromatin accessibility analyses (Fig. 1 ).

As expected (55), Dex treatment resulted in substantial changes to PBMC gene expression levels. Principal component (PC) analysis of all expressed genes (n = 10,385) clearly separated Dex− from Dex+ samples on PC1 (paired T42 = 80.61, P = 1.12 × 10−47; Fig. 1). At the individual gene level, we identified 7,359 differentially expressed genes (70.9% of those tested) between Dex+ and Dex− samples at a false discovery rate (FDR) of 10% (n = 5,032 differentially expressed genes at FDR < 1%) controlling for dominance rank, PBMC composition, and relatedness among study subjects. This set of 7,359 genes ( and Dataset S1) included major targets of Dex induction (e.g., NFKBIA, FKBP5, and PER1; Fig. 1) and repression (CSF1, IER2, and NR3C1, the gene that encodes GR) (55). It is also enriched for genes that are both near GREs and identified as Dex-responsive in an independent dataset (hypergeometric test, P = 9.6 × 10−3; ref. 55) and exhibits effect sizes for the Dex response that are well correlated to those reported in several other studies [Pearson’s r = 0.56, P = 3.81 × 10−11, n = 119 genes from the A549 lung epithelial cell line (55); Pearson’s r = 0.27, P = 6.62 × 10−25, n = 1450 genes from PBMCs (66); and Pearson’s r = 0.40, P = 7.19 × 10−125, n = 3,232 genes from PBMCs (67)].

Dex treatment also remodeled the PBMC chromatin landscape, although with less pronounced effects than on gene expression. Control and Dex+ samples separated on PCs 1, 2, and 3 of the complete ATAC-seq dataset (PC1 T42 = −3.85, P = 3.95 × 10−4; PC2 T42 = 4.49, P = 5.44 × 10−4; PC3 T42 = −4.93, P = 1.35 × 10−5; Fig. 1). Of the 19,859 300-bp windows we analyzed (), we identified 1,480 (7.45%; full list in and Dataset S1) that became significantly more accessible in the Dex+ condition and 1,602 (8.07%) that become significantly less accessible after Dex treatment, at an FDR threshold of 10% (215 and 190 at an FDR of 1%). Further, although GR is thought to primarily bind DNA in regions that are already accessible before nuclear translocation (59), regions with increased accessibility after Dex administration were significantly enriched for two predicted GR binding motifs (log2odds = 1.75 and 2.04, P = 1.86 × 10−16 and 9.56 × 10−30) and for binding sites for the GR cofactor activator protein 1 (Jun-AP1 complex; log2odds =1.45, P = 3.45 × 10−23; and Dataset S1). Conversely, Dex-repressed regions were enriched for TF binding sites for a component of the proinflammatory TF complex, NF-κB (the NF-κB–p50 binding site: log2enrichment = 1.20, P = 1.57 × 10−15; and Dataset S1).

Social status also exhibited a detectable relationship with PBMC gene regulation, although these effects were modest compared with the effects of Dex treatment. Dominance rank was strongly associated with PC3 of gene expression, which explains 7.6% of the total variance in gene expression (Pearson’s r = 0.77, P = 3.33 × 10−18; Fig. 1), and with PC3 of chromatin accessibility, which explains 5.7% of the total variance in chromatin accessibility (Pearson’s r = 0.48, P = 3.41 × 10−6; Fig. 1). For the gene expression data, the effects of social status for individual genes in the Dex− condition were broadly correlated with those described in previous work on rhesus macaque PBMCs [Pearson’s r = 0.39, P = 4.08 × 10−24 (10)]. However, the effects of rank on gene expression levels were substantially weaker than those described for either purified individual cell types or for all pooled white blood cells (including granulocytes) in the same study subjects (9) (), potentially because here we performed more sample processing after a longer time delay following blood sample collection than for previous mixed cell populations (9, 10). Thus, at an FDR of 10%, we identified only 69 social rank-associated genes in the Dex− samples (n = 1,162 genes at an FDR of 20%). However, consistent with previous findings, genes more highly expressed in low-status animals were enriched for proinflammatory signaling (e.g., inflammatory response, P = 4.2 × 10−4, TNF signaling, P = 1.6 × 10−4; complete enrichment results presented in and Dataset S1). Weighted correlation network analysis (68) also revealed seven coregulated gene modules that carried a signature of rank ( and Dataset S1), including one of the largest modules (“module 3,” composed of n = 1,011 genes), for which the module eigengene was significantly negatively associated with rank (β = −0.050, P = 0.042). This module was enriched for genes involved in the inflammatory response, including genes related to IL-6 production (4.4-fold enrichment, P = 1.8 × 10−6) and MyD88 TLR signaling (4.2-fold enrichment, P = 6.8 × 10−4; and Dataset S1). Finally, we identified 159 rank-associated chromatin accessibility windows in the Dex− samples (0.8% of those tested at an FDR < 10%; 1,084 at an FDR < 20%), controlling for kinship and tissue composition—the first evidence to date for an effect of social status on locus-specific chromatin accessibility in humans or other primates.

GC Treatment-Dependent Effects of Social Status.

Previous work indicates that LPS stimulation magnifies the effects of social status on gene expression levels. In contrast, our results show that Dex treatment dampens the social status–gene expression relationship. Among genes that were nominally associated with status in either the Dex+ or Dex− conditions (P < 0.05; n = 2,279 genes), social status effects on gene expression were, on average, four times stronger in the Dex− condition compared with Dex+ (paired T2278 = 17.33, P = 2.63 × 10−63; Fig. 2). Additionally, with increasing effect size, the number of social status-associated genes detected in the Dex− versus the Dex+ condition became increasingly more biased toward the Dex− condition (Fig. 2). This observation is not explained by increased noise in the Dex+ condition (). However, while social status effect sizes were smaller after Dex treatment, they were generally well correlated between conditions (Pearson’s r = 0.75; P < 10−300; ). Thus, although the change in effect sizes between Dex− and Dex+ conditions implies an interaction between social status and GC exposure, we did not identify individual genes for which the log-fold change gene expression response to Dex was significantly associated with social status (10% FDR).

Fig. 2.

Differences in the stability of social status effects on gene expression and chromatin accessibility after Dex treatment. (A) The effects of social status on gene expression are stronger before Dex treatment (paired T2278 = 17.33, P = 2.63 × 10−63), but status effects on chromatin accessibility are similar before and after Dex treatment (paired T4082 = 2.69, P = 0.007). “Scaled difference” is the difference between the absolute value of the effect of social status in the Dex− and Dex+ conditions scaled by the sum of the effect sizes: . Solid vertical lines represent the median difference in scaled effect sizes. (B) We observed an excess of genes with large social status effects on gene expression in the Dex− versus the Dex+ condition. In contrast, we observed no consistent bias in the number of social status-associated genes in the Dex− versus Dex+ condition as effect size increased. The x axis depicts mutually exclusive bins of effect sizes for the standardized effects of status on gene expression (blue) or chromatin accessibility (pink). The y axis represents the log2-transformed ratio (±SE) of the number of genes falling in each bin for the Dex− compared with the Dex+ condition.

Compared with social status effects on gene expression, social status effects on chromatin accessibility were stable between conditions. In regions nominally associated with rank in either condition (P < 0.05; n = 4,083 regions), the effect of social status was significantly different across Dex− and Dex+ conditions, but much more similar in direction and magnitude than for gene expression (paired T4082 = 2.69, P = 0.007; Fig. 2). Additionally, unlike for gene expression, we observed little evidence for stronger rank-chromatin accessibility associations in the Dex− relative to the Dex+ condition with increasing effect size (Fig. 2), despite higher variance in the Dex+ condition (). Thus, while GC treatment rapidly altered social status effects on gene expression levels it did not induce systematic changes in social status-associated chromatin accessibility on the same timescale.

Dex Treatment Decouples Social Status-Associated Chromatin Regions from Social Status-Associated Gene Expression.

To test whether social status effects on gene expression are explained by social status-driven differences in chromatin accessibility, we next examined the relationship between the RNA-seq and ATAC-seq datasets. In the Dex− condition, we found that regions that were more accessible in high-ranking females were associated with genes that were also more highly expressed in high-ranking females (and vice versa) (Pearson’s r = 0.23, P = 1.51 × 10−199; n = 16,677 region-nearest gene pairs; Fig. 3). This relationship becomes stronger as either the threshold for classifying social status-associated genes (Fig. 3) or the threshold for classifying social status-associated chromatin accessibility regions becomes more stringent (Fig. 3). Thus, while social status effects on chromatin accessibility explained only 5.3% (95% CI: 4.6–6.0%) of variance in social status effects on gene expression for all region–gene pairs, they explained 12.7% (95% CI: 10.7–14.9%) of variance for genes where rank was nominally associated with gene expression (P < 0.05). This estimate is almost certainly conservative, as the regulatory elements that influence a gene’s expression are not well predicted by proximity alone (69). Thus, under non-GC-stimulated conditions, social environment-driven changes in the chromatin landscape may mechanistically contribute to social environment-driven variation in gene expression levels.

Fig. 3.

Dex treatment affects the correlation between rank-associated chromatin accessibility and rank-associated gene expression levels. Standardized rank effects on chromatin accessibility and standardized rank effects on the expression of the closest gene are correlated in the Dex− condition (Pearson’s r = 0.23, P = 1.51 × 10−199; n = 16,677 region–gene pairs) (A) but negligibly so after Dex treatment (r = 0.10, P = 2.76 × 10−9; n = 16,677) (B). (C and D) Pearson’s correlation (±SE) between the standardized effect of social status on chromatin accessibility and the standardized effect of social status on gene expression for the closest gene (y axis) is stronger in the Dex− condition compared with the Dex+ condition, whether conditioning on the gene expression association with dominance rank (C, x axis) or the chromatin accessibility association with dominance rank (D, x axis). CA, chromatin accessibility; GE, gene expression.

Social status effects on chromatin accessibility and social status effects on the expression levels of nearby genes were also significantly correlated after Dex stimulation (r = 0.10, P = 2.76 × 10−9; n = 16,677 region–gene pairs; Fig. 3). However, the magnitude of this correlation was reduced more than twofold relative to before treatment (Dex−). Social status effects on chromatin accessibility explained only 1.07% of the variance in social status effects on gene expression in the Dex+ data, and this relationship was not appreciably strengthened by altering the statistical threshold for identifying status-associated genes or regions (Fig. 3 , dark colors). Thus, social status effects on different levels of gene regulation are partially decoupled following GC treatment, at least on a short-term timescale.

Shared Chromatin Accessibility Signatures of Social Status and GC Treatment.

Previous studies in female rhesus macaques have identified evidence for higher GC levels and increased GC resistance in low-status animals (reviewed in ref. 70), and chronic social stress effects on immune gene expression have been hypothesized to originate in part from GC dysregulation (28). If social status effects operate through a GC-mediated pathway, then gene regulatory phenotypes affected by social status should overlap with those affected by Dex treatment.

This prediction was better supported by the chromatin accessibility data than the gene expression data. The effects of rank on gene expression in unchallenged cells and the effects of Dex on gene expression were not correlated (Pearson’s r = 0.007, P = 0.47; Fig. 4). Further, rank effects were actually stronger in Dex-unresponsive genes (Dex effect FDR > 20%) than in Dex-responsive genes (Dex effect FDR < 1%; T2014 = 5.10, P = 3.66 × 10−7). In contrast, regions that were more accessible after Dex treatment were also more accessible in high-ranking females, and vice versa (Pearson’s r = 0.18, P = 5.78 × 10−147; Fig. 4). Further, TF binding site (TFBS) motifs enriched in regions that were more accessible in high-ranking individuals were similarly enriched in regions that opened after Dex treatment (Pearson’s r = 0.64, P = 1.95 × 10−41; Fig. 4, , and Dataset S1). Compared with PBMCs from low-status individuals, cells from high-status individuals thus appear epigenetically primed to respond to GCs. As was the case with gene expression, we found little evidence for rank-by-treatment interactions on chromatin accessibility: Rank predicted the log2-foldchange chromatin accessibility response to Dex at three loci (FDR < 10%; n = 31 loci at an FDR < 20%).

Fig. 4.

The relationship between rank effects and Dex effects on gene expression, chromatin accessibility, and TF accessibility. (A) The effects of rank and Dex treatment on gene expression were uncorrelated (Pearson’s r = 0.007, P = 0.47), unlike (B) the effects of rank and Dex treatment on chromatin accessibility (Pearson’s r = 0.18, P = 5.78 × 10−147). (C) TF binding sites enriched in Dex-induced regions (n = 2,874 regions, FDR < 20%) also were enriched in regions that were more accessible in high-ranking females (n = 544 regions, FDR < 20%; Pearson’s r = 0.64, P = 1.95 × 10−41). Conversely, TF binding sites enriched in Dex-repressed regions (negative values on the x axis; n = 2,905 regions, FDR < 20%) were also enriched in regions that were more accessible in low-ranking females (negative values on the y axis; n = 540 regions, FDR < 20%). Each point (n = 369 TFBS) represents the log2-fold enrichment of one TF motif in Dex-induced compared with Dex-repressed regions (x axis) and rank-induced versus rank-repressed regions (y axis). TF motifs closely associated with GC regulation (GRE, AP-1, and NF-κB) are shown in colored circles.

Discussion

Together, our results dovetail with those from eusocial insects [honey bees and ants: (71–74)], sticklebacks and cichlid fish (18, 60, 75, 76), and rodents (61, 77, 78) to emphasize the importance of social modulation of gene regulation, including both gene expression levels and epigenetic states. Our findings specifically reinforce the role of social interactions in immune cell gene regulation, building substantially on our prior work (9, 10) by revealing a key role for chromatin dynamics and accessibility to TF binding. Here, we provide support for systematic, status-induced changes in chromatin accessibility in primary immune cells, adding to the growing literature on social status effects on gene expression (9, 10, 44, 47, 49, 52, 79, 80), DNA methylation (10, 81–84), and H3K27 acetylation (60). Rank-driven variation in chromatin accessibility may also contribute to previously reported rank-driven variation in the gene expression response to LPS (9). Our results show that low-ranking animals present chromatin landscapes more accessible at some TFBSs for NF-κB, a key TF in the inflammatory response. In contrast, high-ranking females exhibit more accessible binding sites for the GR cofactor AP-1, which is involved in antiinflammatory activity and NF-κB repression (Fig. 4). These effects are likely downstream of changes in the central nervous system, where studies in mice and voles have shown gene regulatory sensitivity to social status (61, 62); indeed, GC release is partially controlled by activity of the hypothalamus. However, recent evidence suggests that social environmental effects in the peripheral immune system also reciprocally influence the brain (78), raising the possibility that effects reported here subsequently feed back to influence social perception and behavior (85).

While social status effects are exaggerated after exposure to LPS (9), here we found that social status effects are attenuated after GC treatment (Fig. 2). Our results are broadly consistent with the evidence that chronic social stress leads to poor control of the inflammatory response (28). Specifically, by administering a high dosage of antiinflammatory GCs we appear to have attenuated, and in some cases eliminated, the signature of low social status. However, we note that Dex treatment did not entirely eliminate the gene regulatory signature of social status. High-ranking and low-ranking females presented consistently different patterns of chromatin accessibility regardless of treatment condition, resulting in a detectable association between social status effects on gene expression and chromatin accessibility at baseline, but not following Dex administration. These observations suggest that the imprint of chronic stress exposure may persist in the epigenome even when gene expression levels shift in response to short-term environmental cues. Such a model is appealing because the effects of social adversity are sometimes observable over the span of years, suggesting that they remain “biologically embedded” even as individuals experience and respond to other environmental stimuli (86). Some degree of decoupling between the gene regulatory substrates for biological embedding and gene expression levels themselves could help satisfy this requirement for simultaneous stability and plasticity.

A nonmutually exclusive alternative explanation is that we stopped the experiment too early to observe treatment effects on social status-associated chromatin accessibility comparable to those observed for gene expression. We chose our sample collection time because it yields a strong gene expression response to Dex (55). However, recent evidence suggests that gene expression changes occur rapidly in response to a stimulus and often precede detectable changes in the epigenome (e.g., CpG methylation; refs. 87 and 88). After stimulation by GCs, many GR–DNA binding events occur at sites that are already accessible before GC stimulation (37). Thus, changes in social status effects on chromatin accessibility might not manifest until a later time point (e.g., as a consequence of the secondary effects of GC signaling and GR activation). Such a possibility would be consistent with our observation that the gene expression data are more clearly separated by treatment condition than the chromatin accessibility data after 90 min (Fig. 1 ). Future studies with more time points would enable finer dissection of social status effects on the gene regulatory response to GCs, including whether social status affects how quickly cells return to a baseline state.

Our findings combine with those of others (67, 87, 89) to emphasize that the detectability and magnitude of predictors of gene expression variation are often context-dependent. Here we demonstrate that this context dependence depends, in turn, on the stimulus and cell type assayed: In our system, social status effects increase after LPS exposure (9), which is detected by TLR4 receptors on monocytes and macrophages, but decrease after GC exposure, which is detected by GR, a much more broadly expressed gene. Context-dependent effects on gene regulation have attracted substantial interest in research on the genetic contribution to gene regulation, where an increasing number of studies have identified genetic variants that influence the response to environmental stimuli (i.e., “response-QTL”) as opposed to, or in addition to, gene expression levels at baseline (e.g., refs. 67 and 90–95). Recent work suggests that response-QTL can emerge as a consequence of genetic variation that influences chromatin accessibility in enhancer regions before stimulation (96). Under this “enhancer priming” model, changes in the cellular environment poststimulation alter the complement of active TFs, thus translating previously silent effects of genotype on chromatin architecture to detectable effects on gene expression. Our results suggest that social environmental effects on gene expression can be masked or unmasked by environmental stimuli in a similar fashion. Here, status-driven differentially accessible chromatin predicts status-driven differential gene expression under some conditions (Dex−) but not others (Dex+), likely in relationship to the activity of GC-associated TFs.

Taken as a whole, our results indicate a subtle but potentially important role for GC physiology in mediating the effects of social status on PBMC function. Key aspects of PBMC chromatin accessibility in high-status females—particularly increased accessibility to the GR cofactor AP-1 and decreased accessibility to NF-κB—overlap with Dex-induced changes. Given that GR binding is largely predicted by the chromatin accessibility landscape before GC stimulation (59, 97, 98), we speculate that increased accessibility at GR cofactor binding sites in high-ranking animals at baseline could drive more efficient GC negative feedback (99), which in turn would result in tighter control of the inflammatory response. If so, increasing chromatin accessibility at key GC-sensitive targets (or decreasing accessibility at key NF-κB targets) in cells from low-ranking animals [e.g., by using epigenome editing approaches (100)] would be expected to reconcile some of the gene expression differences between low- and high-status individuals. Notably, our results also suggest that, while most research on the gene expression response to GC stimulation has focused on its genetic determinants (56, 101, 102), social environmental factors may also be an important factor to consider in human population samples.

Finally, we note that the current study is limited to female rhesus macaques, where social hierarchies are stable, strictly enforced, and determined largely by order of introduction (or, in more naturalistic settings, by kin-based nepotism) rather than physical ability. Given that social environmental conditions can also affect gonadal function (103, 104), the effects of social status reported here may differ in males. Indeed, strong sex-dependent relationships between dominance rank and immune gene expression have recently been reported in wild baboons, with much more pronounced signatures of rank in males compared with females (80). Further, relatively weak effects of social status were observed in female baboons, potentially due to the stress-buffering effects of living with close kin. Indeed, comparative studies across primates indicate that dominance rank is most closely linked to GC levels when hierarchical relationships are strictly enforced and subordinates have few avenues for social support (105). If GC-mediated gene regulation plays an important role in social status-driven changes in immune gene regulation, as our findings here suggest, similar heterogeneity across sexes and social contexts should be expected.

Materials and Methods

Study Population.

Samples were collected from 43 captive female rhesus macaques (Macaca mulatta) at the Yerkes National Primate Research Center in March 2015. These study subjects were part of a multiyear study on the effect of dominance rank (the primary measure of social status in rhesus macaques) on immune function. Demographic details for each study subject are provided in and Dataset S1. The study commenced in January–June 2013 using 45 females (five females assigned to each of nine social groups), as described in ref. 9. In March–June 2014, social group membership was rearranged to form nine new groups in which females that occupied the same or adjacent ordinal ranks in the initial social groups were cohoused together. This approach maximized all possible changes in dominance rank between the initial and rearranged groups, a strategy that allows us to establish the causal relationship between dominance rank and downstream outcomes. Social groups were constructed by sequentially introducing females into indoor–outdoor run housing (25 m by 25 m for each area) over the course of 2–15 wk, and behavioral data collection started after the fifth and final female was introduced into each group.

For the present analysis, we used only samples collected in the second phase of the study (study groups formed in March–June 2014). One of the original study subjects died before we could sample her, and a second female was ill at the time of sampling so we excluded her from the current analysis, leaving a total of 43 study subjects. Samples were collected 38–50 wk after the fifth and final female was introduced to each social group, after rank hierarchies had stabilized following the group rearrangement. Behavioral data in this study included a mean of 19.5 h of observational data per group (range = 14.5–23.0 h), collected using focal sampling (106). We used the EloRatings package in R to assign social status using the Elo method, a continuous measure of dominance rank that is based on an animal’s history of dominance interactions with its groupmates (63, 107, 108). The Elo stability index ranged between 0.995–1.00 (n = 9 groups), where 1 represents a hierarchy in which higher-ranking females always win competitive encounters with lower-ranking females, with no rank reversals, and 0 represents one in which rank does not predict the outcomes of competitive interactions (109). Data included in this study were obtained in accordance with Institutional Animal Care and Use Committee protocols approved by Emory University (YER-2001677-00715) and Duke University (A079-12-03).

In Vitro Dex Treatment.

For each study subject, we drew 4 mL of blood into BD Vacutainer CPT tubes and shipped the blood sample overnight on ice from the Yerkes National Primate Research Center to Duke University, where PBMCs were purified via gradient centrifugation. To control for variation in cell type, an aliquot of 50,000 PBMCs was analyzed by flow cytometry. We quantified the proportion of 10 different cell types in the sample: classical monocytes (CD14+/CD16−), CD14+ intermediate monocytes (CD14+/CD16+), CD14− nonclassical monocytes (CD14−/CD16+), helper T cells (CD3+/CD4+), cytotoxic T cells (CD3+/CD8+), double-positive T cells (CD3+/CD4+/CD8+), CD8− B cells (CD3−/CD20+/CD8−), CD8+ B cells (CD3−/CD20+/CD8+), natural killer T lymphocytes (CD3+/CD16+), and natural killer cells (CD3−/CD16+). Cell type composition was summarized using PC analysis, where the first three PCs of cell composition together explained 97.7% of the variance in whole blood composition. We did not obtain flow cytometry data for two animals, so imputed the values of PCs one to three for these two samples using elastic net regression fit to the Dex− gene expression data (110).

For each study subject, we seeded 1 million PBMCs each into two wells of a tissue culture plate containing 2 mL of media (RPMI; 15% FBS; 1% penicillin–streptomycin). We treated one of the wells with 1 μM Dex (Dex+) and the other with 0.02% EtOH (Dex−), the vehicle for Dex. Cells were then incubated for 90 min (37 °C and 5% CO2), washed with 1× PBS, and then harvested for gene expression and chromatin accessibility analyses.

Chromatin Accessibility.

We measured chromatin accessibility using the ATAC-seq method (64, 65) and 50,000 cells from each treatment–sample combination (n = 43 Dex− and 43 matched Dex+). We followed the protocol from ref. 65 but skipped the cell lysis step, a modification that reduces the proportion of reads mapping to the mitochondrial genome. ATAC-seq libraries were amplified for a total of 12 PCR cycles, barcoded, multiplexed, and sequenced on an Illumina NextSeq 500 using paired-end, 38-bp reads. Reads were mapped to the rhesus macaque genome (Mmul_8.0.1; Ensembl gene build 87) using the bwa-mem algorithm with default settings (111). The resulting dataset exhibited the typical periodic pattern of insert sizes associated with nucleosome positioning (). For downstream analyses, we only considered reads that mapped uniquely to the nuclear genome with a mapping quality score ≥10. For differential accessibility analysis, we then counted the number of Tn5 insertion sites (5′ of the reads aligned to the positive strand offset by +4 bp and to the negative strands offset by −5 bp) that mapped to nonoverlapping 300-bp windows of the genome using the bedtools “coverage” function. Before normalization, we filtered all regions that did not pass our coverage threshold (median cpm < 0.5 in the 86 samples) or had poor mappability [<0.90 as measured using GEM (112)].

Gene Expression.

Total RNA was extracted from each PBMC sample using Qiagen’s RNeasy Mini kit. RNA-seq libraries were prepared from 200 ng of total RNA using the NEBNext Ultra RNA Library Prep Kit and sequenced on an Illumina NextSeq 500 using paired-end, 38-bp reads. We quantified gene expression levels from the resulting reads by aligning them to the rhesus macaque transcriptome (Mmul_8.0.1; Ensembl gene build 87) using kallisto (113) and analyzed gene-wise estimates for genes with reads per kilobase per million mapped reads ≥2 in at least half the samples from each condition (n = 10,385 genes). We obtained 40.0 million (±1.2 SE) and 53.1 million (±2.1 SE) paired-end reads per RNA-seq and ATAC-seq library, respectively (see also and Dataset S1).

Associating Genomic Regions with the Nearest Gene.

We paired each of the 19,859 nonoverlapping windows from the chromatin accessibility analysis to the transcription start site of the nearest gene. For a subset of these region–gene pairs (n = 4,129), the nearest gene was not included in the gene expression dataset; these pairs were excluded from downstream analysis. We thus analyzed a total set of 16,677 region–gene pairs, which consisted of 16,102 regions and 7,547 genes.

Read Count Normalization and Correction for Batch Effects.

We normalized the filtered RNA-seq and ATAC-seq read count matrices using the function voom from the R package limma (114). Next, we modeled these normalized values as a function of the social group membership of the sample donor (nine total social groups). This allowed us to remove biological variation related to differences in group dynamics and, because we sampled all females in each group on the same day, most technical batch effects related to sample collection and processing. We then used the residuals from these models as our primary outcome measure in all subsequent analyses.

PC Analyses.

To examine global patterns of gene expression and chromatin accessibility, we carried out a PC analysis on the residual gene expression and chromatin accessibility data, after removing effects of study subject and cell composition. For both datasets, we then performed PC analysis on the samplewise correlation matrices.

Modeling Rank and Treatment Effects on Gene Regulation.

We used a linear mixed-effects model that controls for relatedness to identify genes or regions that were significantly associated with Dex treatment or dominance rank (115–117). We analyzed the normalized chromatin accessibility and normalized gene expression datasets separately using the R package EMMREML (118). We estimated the effects of rank and treatment on chromatin accessibility in each region and expression of each gene using the following model:

where y is the 86-sample vector of normalized gene expression levels or chromatin accessibility levels for a given gene/region, μ is the intercept, t is the treatment (Dex− = 0, Dex+ = 1), r is the Elo rating (i.e., social status), i is an indicator variable for the treatment (Dex− = 0, Dex+ = 1), β is the effect size of status in the Dex− condition, β is the status effect in the Dex+ condition, and is the effect of Dex treatment. The m by 1 vector u is a random effects term to control for kinship and other sources of genetic structure. Here, m is the number of unique females in the analysis (m = 43), the m by m matrix K contains estimates of pairwise relatedness derived from a 54,165-variant genotype dataset described in ref. 9, is the genetic variance component (0 for a nonheritable trait; note that most gene expression levels are heritable; refs. 119 and 120), and Z is an incidence matrix of 1s and 0s that maps measurements in the Dex+ and Dex− conditions to individuals in the random effects term, thus accounting for repeated measurements for the same individual. Residual errors are represented by ε, an n by 1 vector, where represents the environmental variance component (unstructured by genetic relatedness), I is the identity matrix, and MVN denotes the multivariate normal distribution. Other fixed effects covariates were the first three PCs of the cell type proportion data, which together explain 97.7% of the variance in whole-blood composition (matrix X, with effect sizes V). For each dataset and gene, we tested the null hypothesis that the effect of interest (β or δ) was equal to 0 versus the alternative hypothesis that the effect was not equal to 0.

For all models, we produced empirical null distributions for rank effects by permuting the Elo rating values associated with pairs of samples (Dex− and Dex-treated) as a block, 1,000 times, and rerunning our analyses. To produce a corresponding empirical null for the effect of treatment, we randomized treatment (Dex− versus Dex+) within the pair of samples for each female.

For all of our analyses, we controlled for multiple testing using a generalization of the FDR method of Storey and Tibshirani (121) to empirical null P value distributions generated via permutation. Gene Ontology enrichment analyses were conducted using the R package, topGO (weight01 algorithm), using either a Kolmogorov–Smirnov test (for gene-by-gene rank effects) or Fisher’s exact test (for gene coexpression modules) (122). Gene coexpression analyses were carried out using the WGCNA package (68) on the untreated (Dex−) samples after removal of technical batch effects as described above, following the WGCNA author recommendations (https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials).

TFBS Enrichment.

To test if binding sites for specific TFs were consistently found in social status- or Dex-associated chromatin regions, we used the motif analysis program HOMER (123) and the function “findMotifsGenome.pl.” We conducted two comparisons to test for the enrichment of 309 known vertebrate binding motifs: (i) regions that were more open in high-status animals compared with regions that were more open in low-status animals (FDR < 20%; n = 544 regions more open with higher rank, n = 540 regions more open with lower rank) and (ii) regions that were more open after Dex treatment compared with regions that were less open after Dex treatment (FDR < 20%; n = 2,854 Dex-induced regions, n = 2,902 Dex-repressed regions).