Author response: Identification of human glucocorticoid response markers using integrated multi-omic analysis from a randomized crossover trial (2021)

Article Figures and data Abstract eLife digest Introduction Results Discussion Materials and methods Appendix 1 Appendix 2 Appendix 3 Data availability References Decision letter Author response Article and author information Metrics Abstract Background: Glucocorticoids are among the most commonly prescribed drugs, but there is no biomarker that can quantify their action. The aim of the study was to identify and validate circulating biomarkers of glucocorticoid action. Methods: In a randomized, crossover, single-blind, discovery study, 10 subjects with primary adrenal insufficiency (and no other endocrinopathies) were admitted at the in-patient clinic and studied during physiological glucocorticoid exposure and withdrawal. A randomization plan before the first intervention was used. Besides mild physical and/or mental fatigue and salt craving, no serious adverse events were observed. The transcriptome in peripheral blood mononuclear cells and adipose tissue, plasma miRNAomic, and serum metabolomics were compared between the interventions using integrated multi-omic analysis. Results: We identified a transcriptomic profile derived from two tissues and a multi-omic cluster, both predictive of glucocorticoid exposure. A microRNA (miR-122-5p) that was correlated with genes and metabolites regulated by glucocorticoid exposure was identified (p=0.009) and replicated in independent studies with varying glucocorticoid exposure (0.01 ≤ p≤0.05). Conclusions: We have generated results that construct the basis for successful discovery of biomarker(s) to measure effects of glucocorticoids, allowing strategies to individualize and optimize glucocorticoid therapy, and shedding light on disease etiology related to unphysiological glucocorticoid exposure, such as in cardiovascular disease and obesity. Funding: The Swedish Research Council (Grant 2015-02561 and 2019-01112); The Swedish federal government under the LUA/ALF agreement (Grant ALFGBG-719531); The Swedish Endocrinology Association; The Gothenburg Medical Society; Wellcome Trust; The Medical Research Council, UK; The Chief Scientist Office, UK; The Eva Madura's Foundation; The Research Foundation of Copenhagen University Hospital; and The Danish Rheumatism Association. Clinical trial number: NCT02152553. eLife digest Several diseases, including asthma, arthritis, some skin conditions, and cancer, are treated with medications called glucocorticoids, which are synthetic versions of human hormones. These drugs are also used to treat people with a condition call adrenal insufficiency who do not produce enough of an important hormone called cortisol. Use of glucocorticoids is very common, the proportion of people in a given country taking them can range from 0.5% to 21% of the population depending on the duration of the treatment. But, like any medication, glucocorticoids have both benefits and risks: people who take glucocorticoids for a long time have an increased risk of diabetes, obesity, cardiovascular disease, and death. Because of the risks associated with taking glucocorticoids, it is very important for physicians to tailor the dose to each patient's needs. Doing this can be tricky, because the levels of glucocorticoids in a patient's blood are not a good indicator of the medication's activity in the body. A test that can accurately measure the glucocorticoid activity could help physicians personalize treatment and reduce harmful side effects. As a first step towards developing such a test, Chantzichristos et al. identified a potential way to measure glucocorticoid activity in patient's blood. In the experiments, blood samples were collected from ten patients with adrenal insufficiency both when they were on no medication, and when they were taking a glucocorticoid to replace their missing hormones. Next, the blood samples were analyzed to determine which genes were turned on and off in each patient with and without the medication. They also compared small molecules in the blood called metabolites and tiny pieces of genetic material called microRNAs that turn genes on and off. The experiments revealed networks of genes, metabolites, and microRNAs that are associated with glucocorticoid activity, and one microRNA called miR-122-5p stood out as a potential way to measure glucocorticoid activity. To verify this microRNA's usefulness, Chantzichristos et al. looked at levels of miR-122-5p in people participating in three other studies and confirmed that it was a good indicator of the glucocorticoid activity. More research is needed to confirm Chantzichristos et al.'s findings and to develop a test that can be used by physicians to measure glucocorticoid activity. The microRNA identified, miR-122-5p, has been previously linked to diabetes, so studying it further may also help scientists understand how taking glucocorticoids may increase the risk of developing diabetes and related diseases. Introduction Glucocorticoids (GCs) have a key role in the metabolic, vascular, and immunological response to stress (Cain and Cidlowski, 2017; Oster et al., 2017). GC secretion from the adrenal gland is under tight dynamic control by the hypothalamic–pituitary–adrenal axis and is regulated in a classic circadian pattern (Cain and Cidlowski, 2017; Oster et al., 2017). Most actions of GCs are mediated by the ubiquitously expressed GC receptor (Cain and Cidlowski, 2017; Oster et al., 2017). The tissue-specific effects of GCs are regulated by many local factors, including pre-receptor metabolism of GCs and the interaction of the GC receptor with tissue-specific transcription factors, or through non-genomic mechanisms (Cain and Cidlowski, 2017; Oster et al., 2017). As a result of this complexity, circulating levels of cortisol relate poorly to tissue action of cortisol, and serum cortisol therefore has limited value as a biomarker for GC action (Karssen et al., 2001). GCs are among the most commonly prescribed drugs, and GC treatment remains a cornerstone in the management of many rheumatic and inflammatory diseases despite the introduction of modern disease-modifying antirheumatic drugs and biological immunomodulatory treatment (Smolen et al., 2017). GC replacement is essential for survival in patients with various forms of adrenal insufficiency (Johannsson et al., 2015). However, metabolic and other side effects of GC treatment or replacement are common (Björnsdottir et al., 2011; Fardet et al., 2012), indicating that current methods to monitor their action and tailor their treatment are inadequate. Unphysiological GC exposure has been implicated in the etiology of several common diseases such as type 2 diabetes mellitus, hypertension, abdominal obesity, and cardiovascular disease (Ragnarsson et al., 2019). Against this background, it is highly desirable to be able to measure and quantify GC action as this might be useful to refine current GC therapy. Biomarkers of GC action will also provide potential mechanistic understanding for the role of GCs in the etiology of many common diseases. Previous attempts to identify biomarkers using metabolomics have identified circulating metabolites associated with GC exposure (Alwashih et al., 2017a; Alwashih et al., 2017b). Integrated multi-omic analysis provides increased robustness over analysis of individual 'omic data sets (Ideker et al., 2011). In particular, the identification of groups within one 'omic 'layer' with shared co-regulation within another 'omic layer implies a functional relationship that can be used both to assess the mechanistical relevance and to support the identification of biomarkers (Karczewski and Snyder, 2018; Misra et al., 2018). The aim of this exploratory study was to define multi-omic patterns derived from independent tissues related to GC action and to use these patterns to search for clinically applicable circulating biomarkers of GC action. Subjects with primary adrenal insufficiency, Addison's disease, lack GC production from the adrenal cortex and can therefore be considered a human GC 'knock-down' model (Figure 1A). An experimental study design including subjects with Addison's disease, standardizing for diurnal variation and food intake, allowed a within-individual comparison between physiological GC exposure and GC withdrawal (Figure 1B). A multi-omic analysis strategy combining data from gene expression in circulation (peripheral blood mononuclear cells [PBMCs]) and an important metabolic tissue, adipose tissue, integrated with circulating microRNAs (miRNAs) and metabolites was used to identify putative biomarkers. The strongest putative biomarkers were then replicated in independent study groups with different GC exposure. Figure 1 Download asset Open asset Clinical and analytical part of the exploratory study and the replication step. (A) Subjects with Addison's disease (primary adrenal insufficiency, step 1) were studied in a random order during both physiological glucocorticoid (GC) exposure and GC withdrawal (step 2). Transcriptomics (whole-genome expression) in peripheral blood mononuclear cells (PBMCs) and adipose tissue (n = 28,869,869 genes), plasma miRNAomics (n = 252), and serum metabolomics in morning samples were analyzed (n = 164) (step 3). Integration of the multi-omic data derived a network including gene expression (derived from two independent tissues), microRNAs (miRNAs), and metabolites that were statistically differentiated between the two interventions (step 4). The miRNA findings, because of their centrality in the network, were replicated in subjects with different GC exposures (within the physiological range) from three independent studies (step 5). (B) Subjects with Addison's disease (primary adrenal insufficiency) received in a random order intravenous (i.v.) hydrocortisone (HC) infusion mixed in 0.9% saline in a circadian pattern (physiological GC exposure) or the same volume of 0.9% saline alone (GC withdrawal) during 22 hr starting at 9 AM more than 2 weeks apart. During the GC exposure, HC(Solu-Cortef) was administered at a dose of 0.024 mg/kg/hr between 9 AM and 12 PM (first day), 0.012 mg/kg/hr between 12 PM and 8 PM (first day), 0.008 mg/kg/hr between 8 PM and 12 AM (first day), and 0.030 mg/kg/hr between 12 AM and 7 AM (second day). Samples for the 'omics analyses were collected at 7 AM on day 2 of the intervention (morning samples). p.o.: oral. Results Clinical experimental study Patient characteristics Eleven subjects with well-defined Addison's disease and no other endocrinopathies were recruited and included in the study between September 2013 and September 2015. One subject discontinued the study after randomization and before the first intervention because of persistent orthostatic hypotension. Ten subjects (four women with three of them post-menopausal) with a median age of 50 years (range, 25–57) and a median disease duration of 23.5 years (range, 1–33) completed all aspects of the study between May 2014 and October 2015. The median daily replacement dose of hydrocortisone (HC) prior to the study was 30 mg (range, 20–30), and 9 out of 10 subjects had treatment with fludrocortisone (mineralocorticoid) at a median daily dose of 0.1 mg (range, 0.1–0.2). Clinical and biochemical outcomes The main time points for sample collection in each intervention were at 9 AM on the first intervention day ('before start') and at 7 AM on the second intervention day ('morning') (Figure 1B). The subjects' last ordinary oral HC dose was administered the day before admission to the study unit. Infusion of HC mixed with isotonic saline ('GC exposure') had no effect on systolic and diastolic blood pressure, body weight, serum sodium and potassium, or plasma glucose concentrations compared to the same amount of isotonic saline infusion alone ('GC withdrawal') (Table 1). HC and saline infusion achieved the intended differences in GC exposure. Both median morning serum cortisol and cortisone during the HC infusion were within the physiological range (298 and 81.2 nmol/L, respectively) and markedly lower during the saline infusion (44.4 and 42 nmol/L, respectively, both p<0.001) (Figure 2). Serum cortisol and cortisone were detected in all subjects' morning samples during the saline infusion, but both overnight (between 12 AM and 7 AM) urinary cortisol and cortisone excretion were below the limit of detection. Both HC and saline infusions were well-tolerated, and no serious adverse events were observed. Three subjects reported mild physical and/or mental fatigue, and one subject reported mild salt craving during the GC withdrawal period. Table 1 Clinical and biochemical tests assessed or collected immediately before each intervention and at 7 AM on the second intervention day during both interventions (n = 10). 9 AM before start of intervention (first intervention day)7 AM (second intervention day)Median (IQR)p-ValueMedian (IQR)p-ValueS-cortisol (nmol/L)HC43.2 (38.0–55.1)0.36298 (228–359)<0.001Saline46.7 (43.1–61.7)44.4 (36.8–52.5)Overnight U-free cortisol (µg/7 hr)*,†HC–678 (459–814)<0.001Saline–<0.01‡S-cortisone (nmol/L)HC40.7 (29.5–49.2)0.7681.2 (60.5–94.9)<0.001Saline42.1 (29.7–50.4)42.0 (28.9–47.3)Overnight U-free cortisone (µg/7 hr)*,†HC–136 (106–151)<0.001Saline–<0.01‡SBP (mmHg)HC122 (111–134)0.88123 (107–139)0.65Saline127 (112–131)124 (101–137)DBP (mmHg)HC76 (63–83)0.8869 (59–76)0.55Saline72 (67–80)65 (60–70)S-sodium (mmol/L)HC141 (139–142)0.85140 (138–141)0.97Saline141 (138–142)140 (138–141)S-potassium (mmol/L)HC4.4 (4.3–4.7)0.274.4 (4.3–4.8)0.26Saline4.2 (4.0–4.5)4.4 (4.3–4.5)P-glucose (mmol/L)HC5.0 (4.5–5.3)0.795.6 (5.2–5.8)0.08Saline5.0 (4.6–5.1)5.2 (5.2–5.4)Body weight (kg)HC72.2 (69.6–77.3)0.5372.8 (69.3–77.7)0.63Saline73.5 (70.5–79.5)74.0 (70.5–79.2) *Overnight U-free cortisol and cortisone were collected from midnight to morning (12 AM to 7 AM) during physiological GC exposure (HC infusion) and GC withdrawal (saline infusion). †One of the 10 subjects was not included in the analysis because of a problem during sample collection. ‡Below the limit of detection. DBP: diastolic blood pressure; GC: glucocorticoid; HC: hydrocortisone; IQR: interquartile range; P: plasma; S: serum; SBP: systolic blood pressure; U: urinary. Figure 2 Download asset Open asset Violin plots of serum and urinary glucocorticoids (GCs) during GC exposure and withdrawal. Individual data and changes for morning (7 AM) serum cortisol and cortisone (A, B), and overnight (12 AM to 7 AM) urinary-free (U-free) cortisol and cortisone (C, D) from both interventions. Boxes represent interquartile range (IQR); whiskers mark spread of points within 1.5 times IQR; violins demonstrate distribution of results with the maximum width representing the highest density within each violin. Difference of median values between interventions is presented with p-values. HC: hydrocortisone. Differentially regulated 'omic elements associated with response to GCs Similarity network fusion (SNF) was used to demonstrate overall similarity between subjects across and between 'omic layers, prior to analysis (Appendix 1 and Appendix 1—figure 1). Differential gene expression was associated with GC response in both PBMC and adipose tissue (Appendix 1). Differential expression of metabolites and miRNA was identified in blood in relation to GC response (Appendix 1). Differentially expressed 'omic elements (DEOEs) are presented in Table 2 and Supplementary file 1a–d. All DEOEs were used for integrated analysis, and false discovery rate (FDR)-corrected DEOEs were used for all other analyses (Table 2). DEOEs from the PBMC and adipose tissue transcriptomes were shown to have limited overlap in response to GC but were enriched for shared pathways, revealing an overlap that indicated shared mechanism in relation to GC exposure (Appendix 2 and Appendix 2—figure 3). Table 2 Summary of differentially expressed 'omic elements in association with response to glucocorticoids. 'Omic data setTotal number of 'omic elementsNumber of significant elements (p<0.05)Number of significant elements (FDR < 0.05)PBMC transcriptome28,86944263997Adipose tissue transcriptome28,86935203115Plasma miRNAome25299Serum metabolome1643814 FDR: false discovery rate; miRNA: microRNA; PBMC: peripheral blood mononuclear cell. We assessed the impact of differential expression on the entire interactome to aid in the identification of similar GC-related function. Interactome network models were generated using differentially expressed genes (DEGs) from both the PBMC transcriptome and the adipose tissue transcriptome. These were shown to be consistent with one another (Appendix 2 and Appendix 2—figures 1 and 2) despite the limited overlap of DEGs. GC-responsive genes were shown to have higher connectivity in the human interactome than expected by chance, demonstrated using 10,000 permutations of this network model (Appendix 2). Integration of PBMC and adipose tissue transcriptomes with plasma miRNAomic and serum metabolomic data Hypernetworks are network structures where edges are not restricted to defining a relationship between two nodes but may be shared between many nodes. As such, these structures can be used to describe complex relationships that link multiple elements. Hypernetworks also allow for the same pair of nodes to be connected by multiple edges. This means that relationships between nodes can be ranked by the number of edges shared between them. Hypernetworks allow for the summary of correlation matrices, compressing the high-dimensional relationships between data points (transcripts/miRNA/metabolites) into a single metric of similarity. Hypernetworks facilitate integration of 'omic data and can be used to define strongly associated elements. Elements with large numbers of shared edges are more similar and likely to be of functional relevance; clustering allows refinement of large 'omic data sets to highly associated elements (Figure 3A, B). Hypernetworks are robust to random error and act to filter out false-positive correlations as these will not have a uniform pattern of correlation across all 'omic elements. Figure 3 with 4 supplements see all Download asset Open asset Hypernetwork analyses of integrated 'omic data. (A) Hypernetworks differ from traditional networks in that edges can connect more than two nodes. Nodes are represented by black circles, edges by colored lines and surfaces. This demonstration shows how one edge can connect (i) two nodes as a one-dimensional line, (ii) three nodes as a two-dimensional surface, and (iii) four nodes as a three-dimensional structure. Hypernetworks of 'omic data can have edges shared between hundreds of nodes. (B) Hypernetwork diagram illustrating how a pair of nodes (a–h) can be connected by more than one edge. In this example, nodes e and d share two edges, as do b and d. (C) A hypernetwork plotted as a heat map can be used to investigate clustering of blood peripheral blood mononuclear cell (PBMC) transcripts, based on correlation to, for example, adipose tissue transcriptome. A central cluster, defined using hierarchical clustering, groups PBMC transcripts based on high numbers of shared edges (red square, n = 965). This approach was applied to define groups of PBMC transcripts with similar profiles when correlated against each other 'omic layer. (D) Gene probe level overlaps between PBMC transcriptome clusters identified by hypernetwork shared with the other 'omic data sets. PBMC transcriptomic changes are correlated with changes in miRNAome, adipose tissue transcriptome, and metabolome (gas chromatography-mass spectrometry and liquid chromatography-mass spectrometry overlaps combined); overlaps are common PBMC transcripts with correlation to the 'omic data sets. Values in brackets represent the size of PBMC transcriptomic clusters drawn from all differentially expressed PBMC transcripts (n = 4426, p<0.05). Data demonstrates a fundamental relationship in glucocorticoid response between PBMC and adipose tissue (965 genes) and reinforces the presence of common pathways in these two independent tissues. To assess similarity, we defined the correlation coefficient between each differentially expressed 'omic measurement and assessed as 'present' in the network model those correlations with an r-value >| 1.5 | were defined as PBMC transcripts with shared for example, two PBMC transcripts that are both correlated with the same three metabolites are connected by three edges. We the shared correlations as a measure of similarity between each pair of GC-responsive PBMC transcripts, correlations across the other 'omic data sets (Figure 1). The number of correlations shared was between PBMC and adipose tissue transcriptome genes, Figure the the overlap of differential expression was common pathways are in both tissues related to GC which similar networks of The order of the number of correlations shared with the GC-responsive PBMC transcriptome was adipose tissue transcriptome plasma serum and this was confirmed both by comparison of the heat (Figure 1) and by a diagram (Figure The diagram also a between the serum metabolome and both PBMC and adipose tissue and of a shared transcriptomic profile in both and adipose tissue GC response was in which were generated to model based on the of correlations with PBMC genes that were by these were from the predictive analysis. this we defined of PBMC transcripts with maximum predictive This of genes the and groups using analysis (Figure We identified of using and the experimental using analysis (Figure This identified a of genes as of with changes in the same in both transcriptomic data sets that HC from saline treatment file of these genes were differentially expressed in both PBMC and adipose tissue transcriptomes (Figure of these genes, were associated with GC response gene and (Figure The gene (n = that we identified, which both PBMC and adipose tissue transcriptomes in relation to GC exposure, was under the by further in other studies of GC action by other research groups in models (Table 3). robustness of the random was by that the at which the of in was small (Figure 1). Figure 4 with 1 see all Download asset Open asset An gene in peripheral blood mononuclear cell (PBMC) and adipose tissue transcriptome can be used to glucocorticoid (GC) These analyses were to the common predictive genes in and adipose (A) analysis of hydrocortisone (HC) infusion (GC exposure, from saline infusion (GC using of PBMC transcripts confirmed as robust in the hypernetwork by analysis of 1 and (B) of in using a approach to model experimental confirmed experimental transcripts were identified as important or in GC exposure from GC as as the same in both PBMC and adipose tissue transcriptomic data sets. (C) genes that are differentially expressed between GC exposure and GC withdrawal in both PBMC and adipose tissue transcriptomes and in the same in both tissues (n = (D) of predictive genes out of with GC response through gene Data demonstrates the presence of a robust transcriptomic profile GC response in two independent tissues. Table 3 of the predictive genes from the current exploratory study against studies GC response in error effect on in et al., effect on and cell et al., GC effect on peripheral blood et al., cell response to of GC receptor et al., et al., effect on et al., the approach of that is using of the that can be used for a to demonstrates a high of a sample from each error rate = error using the data. The gene that both PBMC and adipose tissue transcriptomes in relation to GC exposure with in the same Figure genes) was by further in other studies of GC action in under the of the GC: glucocorticoid; Gene study number with study number PBMC: peripheral blood mononuclear number of the the data. Integration of circulating 'omic data sets to miRNA and of GC action We further between the circulating 'omics data associated with GC exposure (Figure All of the circulating 'omics data was to a correlation and hierarchical clustering used to identify 'omic data points with similar correlation (Figure 1). Eleven clusters including miRNAomic, and metabolomic data were identified, and these clusters were shown to have within the interactome network model file and Appendix 2). We then the number of correlations between all the circulating 'omic data associated with GC exposure (n = using a This approach was used to define a of highly connected multi-omic elements with a relationship to GC exposure (Figure Figure with 1 see all Download asset Open asset Integration of all circulating 'omic data sets associated with glucocorticoid (GC) These analyses were to to putative biomarkers of GC action. (A) Hypernetwork summary heat map of shared correlations between all circulating 'omic elements (peripheral blood mononuclear cell transcriptome, plasma microRNA serum n = with differential expression between GC exposure and GC withdrawal. (B) of central (red in A = of which map to genes with differential differentially expressed differentially expressed one of the genes identified as a of GC response Figure (C) of the representing multiple highly connected genes) by the most central gene in each miR-122-5p is in the of two of centrality most central in the central in the A hypernetwork model of the of highly connected elements was generated (Figure was the gene shared with the adipose tissue transcriptome that also had predictive value with a in Figure of the gene has been shown to be associated with increased to GCs et al., an in with the in expression we in both PBMC and adipose tissue transcriptomes in association with GC exposure (Figure The hypernetwork model (Figure also a range of related and A hierarchical model of within the network was assessed using the measure of network centrality (Figure These revealed multi-omic relationships and demonstrated that miR-122-5p was the miRNA