Identification of human glucocorticoid response markers using integrated multi-omic analysis from a randomized crossover trial

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.

Keywords: biomarkers; computational biology; gene expression; glucocorticoids; human; medicine; metabolites; microRNA; multi-omics; systems biology.

Conflict of interest statement

DC DC has received lecture fees from Otsuka, Sanofi, and Shire. PS, TG, CG, RB, RS, SB, UF, PJ, SS No competing interests declared, BW BW is a consultant with Actinogen Medical, inventor on patents owned by the University of Edinburgh relating to HSD1 inhibitors and to the discovery of glucocorticoid-sensitive biomarkers (the HSD1 inhibitor patents have been licensed to Actinogen Medical, patent numbers: WO 2011/135276, WO 2011/033255, and WO 2009/112845); BW is Non-Executive Director of Roslin Technologies Ltd. OR OR has received lecture fees from Novo Nordisk, Ipsen, Sandoz, and Pfizer, an unrestricted research grant from HRA Pharma, and consultancy fees from Novartis and HRA Pharma. PT PT has received lecture fees from Novartis and Novo-Nordisk. AS AS has received lecture fees from Merck. GJ GJ has served as a consultant for Novo Nordisk, Shire, and Astra Zeneca and has received lecture fees from Eli Lilly, Ipsen, Novartis, Novo Nordisk, Merck Serono, Otsuka, and Pfizer.

© 2021, Chantzichristos et al.

Figures

Figure 1.. 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.

Figure 2.. 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.

Figure 3.. 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.

Figure 3—figure supplement 1.. Hypernetwork heat maps of integrated ‘omic data.

The analyses were performed to quantify groups of highly correlated ‘omic data sets. Hypernetwork heat maps show clusters of correlated peripheral blood mononuclear cell (PBMC) transcripts (total 4426, 3.9 × 10−8 < p<0.05) with differential expression (glucocorticoid [GC] exposure vs. GC withdrawal) shared with (A) adipose tissue transcriptome (red square, n = 965), (B) plasma miRNA (red square, n = 768), (C) serum metabolome measured gas chromatography-mass spectrometry (red square, n = 915), and (D) serum metabolome measured by liquid chromatography-mass spectrometry (red square, n = 312). Increasing numbers of shared correlations appear as a scale from blue to white as shown in the keys.

Figure 3—figure supplement 2.. Transcriptome quality control measures.

(A) Peripheral blood mononuclear cell (PBMC) transcriptome probe-level log-intensity distribution. (B) Adipose tissue transcriptome probe-level log-intensity distribution. Density distributions represent intensity per probe, normalized to proportion of overall measurements at each intensity. (C) PBMC transcriptome per-sample log-intensity distribution (D) Adipose tissue transcriptome per-sample log-intensity distribution. Boxplots represent distribution of intensity per sample; boxes represent interquartile range, with median values and whiskers represent the range. All figures were generated in R from measurements generated by the R package ‘oligo’.

Figure 3—figure supplement 3.. Metabolome quality control measures.

Distribution before and following normalization of (A) gas chromatography-mass spectrometry (GC-MS) samples, (B) GC-MS metabolites, (C) liquid chromatography-mass spectrometry (LC-MS) samples, and (D) LC-MS metabolites. For all figures, distribution of concentration values is provided as boxplots and density distributions for individual samples and metabolites. Boxes represent interquartile range with median values, and whiskers represent the range; outliers are denoted as individual points. Density distributions represent overall distribution of concentrations combined across samples or metabolites. All figures were produced using R package ‘MetaboanalystR’.

Figure 3—figure supplement 4.. Chi-squared distance distribution from a range of dichotomization thresholds.

Hypernetworks were generated, dichotomizing the data using a range of cut-off r-values, from 0.1 to 0.9 by steps of 0.05, in order to generate the hypernetwork incidence matrix M. Overall chi-squared distance was calculated for each hypernetwork as the sum of chi-squared distance across the hypernetwork matrix sM. A general additive model was used to fit a curve through these data using ggplot2. This approach determines how edges are assigned between pairs of ‘omic elements and, therefore, the number of shared edges represented by the hypernetwork matrix sM. Chi-squared distance represents the signal window present in a hypernetwork by measuring the spread of values (numbers of shared edges). The greatest spread in values here was found when dichotomizing around an r-value of 0.5. Above and below this threshold, the range of values in the hypernetwork was reduced, limiting the ability to refine a set of interacting ‘omic elements supported by most edges.

Figure 4.. An overlapping gene set in peripheral blood mononuclear cell (PBMC) and adipose tissue transcriptome can be used to classify glucocorticoid (GC) response.

These analyses were performed to depict the common predictive genes in PBMCs and adipose tissue. (A) Partial least squares discriminant analysis (PLSDA) showing complete separation of hydrocortisone (HC) infusion (GC exposure, blue points) from saline infusion (GC withdrawal, orange points) using 271 of 965 PBMC transcripts confirmed as robust in the hypernetwork by analysis of dissimilarity. X-variates 1 and 2: PLSDA components; expl. var: explained variance. (B) BORUTA feature selection identifies variables (genes) of importance in classification using a Random Forest approach to model experimental background noise (green: confirmed classification; yellow: tentative classification: red: rejected classification; blue: ‘shadow’ variable modeling experimental noise). Of 271 transcripts initially used, 59 were identified as important (confirmed [green] or tentative [yellow]) in separating GC exposure from GC withdrawal, as well as having the same direction fold change in both PBMC and adipose tissue transcriptomic data sets. (C) Predictive genes that are significantly differentially expressed between GC exposure and GC withdrawal in both PBMC and adipose tissue transcriptomes and display fold change in the same direction in both tissues (n = 9). (D) Association of predictive genes (six out of nine) with GC response through gene ontology. Data demonstrates the presence of a robust transcriptomic profile predicting GC response in two independent tissues.

Figure 4—figure supplement 1.. Random forest distribution of minimal depth.

Distribution of minimal depth among trees from random forest generated in Figure 4B. Minimal depth equals the depth of the node closest to the root, which is capable of splitting participants into groups. Minimal depth for variables in this tree is very low as all variables are significantly different between groups at an uncorrected p-value

Figure 5.. Integration of all circulating ‘omic data sets associated with glucocorticoid (GC) response.

These analyses were performed to lead 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 [miRNA], serum metabolome; n = 336) with differential expression between GC exposure and GC withdrawal. (B) Network representation of central cluster (red square in A [n = 139], of which 120 map to genes/miRNA/metabolites). Blue circles: genes with differential expression; red circle: differentially expressed metabolites; green: differentially expressed miRNA. Red box highlights DCK, one of the nine genes identified as a classifier of GC response (see Figure 4C). (C) Module decomposition of the hypernetwork. Genes modules (hexagons representing multiple highly connected genes) named by the most central gene in each module. miR-122-5p is present in the core of two modules (shown); color of modules represents centrality hierarchy: red: most central in the network; green: least central in the network.

Figure 5—figure supplement 1.. Heat map with clusters of circulating ‘omic data associated with glucocorticoid exposure identified using a correlation matrix.

Along the diagonal 11 clusters including peripheral blood mononuclear cell transcriptomic, miRNAomic and metabolomic data were identified (depicted in squares with their numbers). Color scale represents correlation coefficient. .

Figure 6.. Replication of miR-122-5p as a putative biomarker of glucocorticoid (GC) action in the current biomarker discovery study.

Targeted replication of plasma miR-122-5p fold change from the current study population between subjects with Addison’s disease during GC exposure and GC withdrawal showed a significant down-regulation of miR-122-5p expression with increased GC exposure (p=0.009) conducted using an independent RNA separation procedure in the same samples.

Figure 7.. Replication of miR-122-5p as a putative biomarker of glucocorticoid (GC) action in independent patient groups with different GC exposure.

(A) The expression of miR-122-5p was higher in subjects with rheumatoid arthritis and reduced GC exposure due to tertiary adrenal insufficiency after a short-term stop of the GC treatment (AI) than in those without tertiary adrenal insufficiency (Normal). (B) Subjects with Addison’s disease (AD) had higher expression of miR-122-5p than healthy matched controls (Control). In an experimental study in healthy subjects, the expression of miR-122-5p was increased both after low and excessive high GC exposure (LowGC and ExcessGC, respectively) compared to medium GC exposure (MedGC) at both (C) high and (D) low serum insulin levels. Diamond: mean. Box = median ± interquartile range. Whiskers = upper and lower quartiles. miR-122-5p axis is presented as normalized expression.

Appendix 1—figure 1.. Integration of unsupervised peripheral blood mononuclear cell (PBMC) transcriptomic, metabolomic, and miRNAomic data in subjects between glucocorticoid (GC) exposure and GC withdrawal at 7 AM using similarity network fusion.

Integration of all circulating ‘omic data (A) simultaneously and (B) stepwise. Each matrix represents a subject-by-subject comparison using similarity network fusion of change in signal of ‘omic data between GC exposure and GC withdrawal. The Euclidean distances between subjects based on their relative positions as determined by isomap dimensional scaling were used as a basis to determine similarity. Colors denote normalized continuous similarity score; red: dissimilar; yellow: similar. The complete data set for each specific ‘omic layer of data was used, transcriptome n = 53,617, GC-MS n = 82, LC-MS n = 82, miRNA n = 251. GC-MS, gas chromatography-mass spectrometry; LC-MS: liquid chromatography-mass spectrometry; miRNA, microRNA; S, serum.

Appendix 2—figure 1.. Gene ontology and interactome network model of differential gene expression in peripheral blood mononuclear cells (PBMCs) between glucocorticoid (GC) exposure and GC withdrawal.

Differentially expressed transcriptomic data in PBMCs was used to determine (A) the canonical pathways involved and (B) infer an interactome network model using BioGRID database, consisting of 2467 nodes and 64 modules. (C) Modules of genes (octagons) within the interactome network model were determined using the Moduland algorithm and ranked by a hierarchy based on network centrality; central modules are functionally relevant as they can influence the activity of a wider proportion of the network. The modules are arranged (left to right) in columns of the top 10 genes in each ranked by network centrality (top to bottom). (D) Gene expression marked within the network model (red: up; green: down in hydrocortisone treatment). Causal elements (known literature associations implying regulatory function) are marked in blue.

Appendix 2—figure 2.. Correlated networks related to glucocorticoid (GC) action in the peripheral blood mononuclear cell (PBMC) and adipose tissue transcriptome.

The analyses were performed to identify similarities and differences in transcriptome (gene expression) associated to GC exposure in two different tissues. Network associations of GC response to genes with differential expression between GC exposure and GC withdrawal in (A) PBMC transcriptome and (B) adipose tissue transcriptome. Gray circles represent genes involved in response to GCs that are not differentially expressed. Sixteen genes up-regulated and four down-regulated in PBMCs. Thirteen genes up-regulated and eight down-regulated in adipose tissue.

Appendix 2—figure 3.. Gene Ontology related to glucocorticoid (GC) action in the peripheral blood mononuclear cell (PBMC) and adipose tissue transcriptome.

(A) Numeric summary of biological pathways (Gene Ontology) associated with differential gene expression between GC exposure and GC withdrawal in the PBMC and adipose tissue transcriptomes. (B) Biological pathways (Gene Ontology) associated with differential gene expression in both PBMCs and adipose tissue. Data demonstrate that while different genes may be involved in the transcriptomic response to GCs in the two independent tissues, there are common pathways through which GC action is manifested.