Identification of glucocorticoid-related molecular signature by whole blood methylome analysis

Objective Cushing’s syndrome represents a state of excessive glucocorticoids related to glucocorticoid treatments or to endogenous hypercortisolism. Cushing’s syndrome is associated with high morbidity, with significant inter-individual variability. Likewise, adrenal insufficiency is a life-threatening condition of cortisol deprivation. Currently, hormone assays contribute to identify Cushing’s syndrome or adrenal insufficiency. However, no biomarker directly quantifies the biological glucocorticoid action. The aim of this study was to identify such markers. Design We evaluated whole blood DNA methylome in 94 samples obtained from patients with different glucocorticoid states (Cushing’s syndrome, eucortisolism, adrenal insufficiency). We used an independent cohort of 91 samples for validation. Methods Leukocyte DNA was obtained from whole blood samples. Methylome was determined using the Illumina methylation chip array (~850 000 CpG sites). Both unsupervised (principal component analysis) and supervised (Limma) methods were used to explore methylome profiles. A Lasso-penalized regression was used to select optimal discriminating features. Results Whole blood methylation profile was able to discriminate samples by their glucocorticoid status: glucocorticoid excess was associated with DNA hypomethylation, recovering within months after Cushing’s syndrome correction. In Cushing’s syndrome, an enrichment in hypomethylated CpG sites was observed in the region of FKBP5 gene locus. A methylation predictor of glucocorticoid excess was built on a training cohort and validated on two independent cohorts. Potential CpG sites associated with the risk for specific complications, such as glucocorticoid-related hypertension or osteoporosis, were identified, needing now to be confirmed on independent cohorts. Conclusions Whole blood DNA methylome is dynamically impacted by glucocorticoids. This biomarker could contribute to better assessment of glucocorticoid action beyond hormone assays.


Introduction
Cushing's syndrome is a state of glucocorticoid excess related either to glucocorticoid treatment (exogenous Cushing's syndrome) or to excessive secretion of adrenocortical glucocorticoids. While effective in suppressing inflammatory states (1), the prolonged use of administered glucocorticoids is associated with potentially serious adverse effects, restricting their widespread and chronic usage (2). Clinical consequences of systemic glucocorticoid treatment are mirrored by states of endogenous cortisol excess. While overt manifestations of Cushing's syndrome are rare and most often related to pituitary adenomas (3), mild autonomous cortisol hypersecretion is more common and mostly caused by adrenal adenomas (4).
Cushing's syndrome is associated with high morbidity and mortality, and impaired quality of life (5) through numerous systemic manifestations, including diabetes mellitus, hypertension, osteoporosis, cutaneous bruising, muscular atrophy, neuropsychiatric disorders, and immune deficiency (6). Duration and level of glucocorticoid excess are undisputedly the main determinants of Cushing's syndrome severity. However, individual susceptibility highly affects the likelihood to develop each type of complication and modulates their severity (7).
On the other side of the clinical spectrum, adrenal insufficiency is a state of cortisol deprivation, inducing fatigue and acute decompensations of metabolism or electrolyte balance with potentially lethal outcome (8). Adrenal insufficiency can be caused by structural or functional damage of adrenal glands, pituitary, or hypothalamus, with autoimmunity-, tumour-, or treatment-related reasons being the most common causes.
Quantification of glucocorticoid action on peripheral tissues is difficult to assess by clinical means or on the basis of hormonal evaluations. Indeed, for exogenous glucocorticoid administration, pharmacokinetics may importantly influence the level of glucocorticoid excess, especially for low-dose systemic treatments or in case of local administration. For endogenous Cushing's syndrome, increased morbidity and mortality related to mild autonomous cortisol excess are well demonstrated on the population level (9, 10) but cannot be estimated properly on an individual level using classical hormone assays. For patients with adrenal insufficiency, titration of glucocorticoid supplementation relies mainly on clinical assessment. This shortcoming underlines the need for specific biomarkers quantifying glucocorticoid action, with potential impact on diagnosis, treatment decision, and prediction of the individual risk for specific complications.
DNA methylation is a chemically stable yet dynamic biological hallmark, playing a key role in epigenetic regulation of gene expression in both health and disease (11). Several studies have suggested an association between hypothalamic-pituitary-adrenal axis dysregulation and specific blood DNA methylation profiles, particularly in post-traumatic stress disorders (12,13,14,15,16). In addition, the association between stress and DNA methylation has been explored for some targeted genes (17). Among them, methylation of FKBP5 -encoding a co-chaperone of HSP90 protein involved in the regulation of glucocorticoid receptor activity (18) -and NR3C1encoding the glucocorticoid receptor -are impacted by stress. Furthermore, a recent study showed a correlation between FKBP5 expression and cortisol levels in patients with Cushing's syndrome (19).
The present study explores the impact of glucocorticoids on leukocyte DNA methylation. Specifically, we analysed whole blood methylome in patients with endogenous Cushing's syndrome, eucortisolism, or adrenal insufficiency, and we identified a methylome signature reflecting glucocorticoid excess.

Patients and samples
Ninety-four blood samples were collected from 47 patients with a confirmed diagnosis of endogenous Cushing's syndrome. Patients were followed in two expert centres, Cochin hospital (APHP, Paris, France) and LMU hospital (Ludwig-Maximilians-University, Munich, Germany). Diagnostic criteria of Cushing's syndrome included: increased 24-h urine-free cortisol, abnormal cortisol after 1 mg dexamethasone suppression, and altered circadian cortisol rhythm, following consensus guidelines (20).
Blood samples were collected either before correction of Cushing's syndrome or at least 3 months after (3-41 months; median: 12 months). At the time of blood sampling, patients were classified as overt Cushing's syndrome, mild Cushing's syndrome, and eucortisolism or adrenal insufficiency, depending on clinical evaluation and hormone assays. Briefly, overt Cushing's syndrome patients presented clinical signs and increased 24-h urine-free cortisol (>240 nmol/24 h), increased salivary midnight cortisol (>6 nmol/L), and insufficient cortisol suppression after 1 mg dexamethasone (>50 nmol/L). Mild Cushing's syndrome patients lacked signs of clinically overt Cushing's syndrome but had mild alterations of cortisol secretion, including either a slightly increased 24-h urine-free cortisol, or increased midnight salivary cortisol, or insufficient cortisol suppression after 1 mg dexamethasone. Adrenal insufficiency was based on low plasma morning cortisol (<160 nmol/L) and on insufficient response to corticotropin stimulation (<500 nmol/L). For two patients in eucortisolism, exact values were not available. Detailed hormone values for each sample are provided in Supplementary Table 1 (see section on supplementary materials given at the end of this article).
Signed informed consent for molecular analysis of blood samples and for access to clinical data was obtained from all patients, and the study was approved by a local Ethic Committee (for Cochin hospital: Comité de Protection de Personnes Ile de France 1, project 13495; for Munich: project 152-10).
Ninety-one additional samples were available and suitable for methylome analysis from patients enrolled in five specialized centres of the ENSAT-HT consortium (http://www.ensat-ht.eu). They included 26 patients with endogenous Cushing's syndrome (20) and 65 healthy volunteers (eucortisolism) (Supplementary Table 2). Signed informed consent was obtained from all patients, and the study was approved by the Ethic Committee of each participating centre.

Whole-genome DNA methylation measurement
Leukocyte DNA was extracted from EDTA blood samples, using the DNA Isolation kit for Mammalian Blood (Roche). DNA quality was assessed on a Genomic DNA ScreenTape system (Agilent) and quantified using a Qubit 3.0 Fluorometer (Thermofisher). DNA was treated by bisulphite and then hybridized to the Infinium MethylationEPIC BeadChip (Illumina; ~850 000 sites), starting from 500 ng of DNA. All experiments were performed following the manufacturer's instructions at the P3S Post-Genomic Platform of Sorbonne University (Paris, France).

Bioinformatics and statistics
All samples passed the quality controls provided by the Genome Studio software (v. 2011.1; Illumina). Data were exported as Intensity Data (IDAT) format and then processed using the minfi package (v. 1.32.0) (21) in R software environment (v. 3.6.3) (https://cran.r-project.org/). Data were normalized using the stratified quantile normalization procedure implemented in the preprocessQuantile minfi function (22) and the methylation score for each CpG probe was extracted as a β-value. The ChAMP package (v. 2.16.1) was used to filter the probes (23). A total of 731 635 probes passed the following criteria: detection P -value < 0.01, presence of the targeted CpG, absence of frequent SNPs in the probe, single hybridization hit, and autosomal target.
The significant components of variation in the data set were assessed using the singular value decomposition method (SVD) for methylation data (24) and a detected batch effect (Slide) was corrected using the ComBat method (25), as implemented in the ChAMP package.
White blood cell count of subpopulations (neutrophils, lymphocytes B, lymphocytes T4, lymphocytes T8, lymphocytes NK, and monocytes) was estimated by the reference-based RefbaseEWAS method (26) implemented in the ChAMP package. To confirm the reliability of the inferred white blood cell counts, we compared the estimated and available measured proportions of both neutrophils and lymphocytes, obtaining high correlation (Pearson's r = 0.81 and r = 0.87, respectively) (Supplementary Fig.  1 and Supplementary Table 3). Since neutrophils were the most represented cell type in all samples, and since the proportions of neutrophils and lymphocytes were negatively correlated (Pearson's r = −0.97), the estimated proportion of neutrophils was chosen as the unique proxy reflecting variations in white blood cell count.
M-values, used for statistical analyses, were calculated from β-values (log 2 ratio of the intensities of methylated vs unmethylated probes) using the lumi package (v. 2.36.0) (27).
Global data structure was assessed on β-values by principal component analysis (PCA), using all CpG probes. Probe variability was calculated on M-values as s.d. for each CpG site among samples. The most variable CpG probes (n = 52 727 with s.d. > 0.4) were selected for subsequent analyses.
Differentially methylated CpG sites were identified starting from the whole data set using the Limma package (v. 3.40.6) (28), including the estimated neutrophils count as covariate, and considering a Benjamin-Hochberg adjusted P -value < 0.05. Gene set enrichment analysis of genes associated with differentially methylated CpG sites was performed using the gometh method implemented in the missmethyl package (v. 1.18.0) (29), adjusting for the number of CpG sites associated to each gene (30). Differentially methylated regions were identified using the DMRcate package (v. 1.20.0) (31), comparing overt Cushing's syndrome vs eucortisolism or adrenal insufficiency, and using estimated neutrophil's count as covariate. Default parameters were applied to smooth the differential methylation signal and to define regions.
For predicting the glucocorticoid status from methylation, a training cohort of 60 samples was selected, randomly including 30 samples corresponding to overt Cushing's syndrome with no anti-cortisolic treatment at the time of sampling and 30 samples corresponding to either eucortisolism or adrenal insufficiency after Cushing's syndrome treatment. Remaining samples (n = 34) were used as a first validation cohort. CpG site selection was performed on the training cohort using a penalized Lasso regression on the most variable CpG probes (M-value s.d. > 0.4), with a ten-fold cross-validation, using the glmnet package (v. 4.0-2) (32). The predictive model, including 29 discriminating CpG sites, was assessed on the validation cohort, graphically using a PCA projection of samples based on the 29 CpGs methylation level and statistically using an ordinal logistic regression model including a 29-CpGs predictor, calculated by adding the 29 CpGs M-values weighted by their Lasso coefficients. Similarly, the performance of the 29-CpGs methylation predictor was tested on the ENSAT-HT cohort, a second independent validation cohort.
In order to identify potential CpG sites specifically associated to glucocorticoid-related complications (hypertension, diabetes, and osteoporosis), CpG site selection was performed on 47 Cushing's samples, starting from the most variable CpG probes (M-value s.d. > 0.4), using a penalized Lasso regression with a ten-fold crossvalidation. A methylation predictor for each model was calculated as described for the 29-CpGs predictor.
Quantitative variable comparisons between groups were performed using two-tailed t-test or Wilcoxon's test, depending on variable distribution. Quantitative variable correlations were performed using Pearson's test. Multivariate analysis was performed using a logistic regression model including the 29-CpGs methylation predictor and the estimated proportion of neutrophils as covariates. All tests were computed in R software environment.
An additional independent cohort of 91 samples, part of the European ENSAT-HT consortium, was collected, including 26 patients with endogenous Cushing's syndrome and 65 healthy volunteers (eucortisolism; Supplementary Table 2).

Glucocorticoid levels impact blood methylome
Whole blood DNA methylome was determined for the 94 samples, with 731 635 informative CpG sites in all samples. Unsupervized PCA showed a discrimination of samples according to their glucocorticoid status, with a specific profile of overt Cushing's syndrome (Fig. 1A). This discrimination was mainly related to the global methylation level. Indeed, overt Cushing's syndrome status was associated with overall decreased methylation among the most variable CpG sites (t-test P -value < 0.05 for 52 727 CpG sites; Fig. 1B). Another significant determinant was the white blood cell count variation ( Supplementary Fig.  2), related to the well-established effect of glucocorticoids on white blood cell composition, inducing granulocytosis and lymphopenia (33,34).

Exploration of glucocorticoid-related blood methylome profile
The specific effect of glucocorticoids on blood methylome was evaluated by comparing the methylation level in overt Cushing's syndrome samples (n = 42) vs each of the other three groups individually -mild Cushing's syndrome (n = 13), eucortisolism (n = 14), and adrenal insufficiency (n = 25) samples (Supplementary Table 4).
The most significant difference was observed in the comparisons of overt Cushing's syndrome vs eucortisolism (n = 1290 differentially methylated CpG sites) and overt Cushing's syndrome vs adrenal insufficiency (n = 7120 differentially methylated CpG sites). Both comparisons showed a prevalence of hypomethylated CpG sites in overt Cushing's syndrome (80 and 73%, respectively). Hypomethylated CpG sites were observed both in 'Open Sea' and 'Island' regions, showing the independence of glucocorticoid-related hypomethylation from CpG enrichment in the genome ( Fig. 2A and Supplementary  Fig. 3A). Glucocorticoid-related hypomethylation was not related to any specific gene locus structure either ( Fig. 2B and Supplementary Fig. 3B).
GSEA of genes associated with the differentially methylated CpG sites in the two comparisons revealed an enrichment in immunity-related signalling pathways (Gene Ontology gene sets; FDR < 0.05), particularly those relating to neutrophils degranulation (Supplementary  Tables 5 and 6).
Differentially methylated CpG sites were distributed all along the genome. One gene locus was strongly enriched in differential CpG sites, on chromosome 6, corresponding to the FKBP5 gene locus (Fig. 2C, Supplementary Fig.  3C and Supplementary Table 4). Beyond the analysis of individual CpG sites, a specific analysis of differentially methylated regions identified 99 and 418 differentially methylated regions in overt Cushing's syndrome vs eucortisolism and vs adrenal insufficiency, respectively (Supplementary Tables 7 and 8). Again, one of the most significant differentially methylated regions associated with the FKBP5 gene promoter (Stouffer's Z-score < 0.001) (Supplementary Table 7). This FKBP5 promoter region included 5 CpG sites, whose methylation level properly discriminated overt Cushing's syndrome samples from all the others (t-test P -value < 0.05; Fig. 3A). This effect was tested on the ENSAT-HT cohort, confirming the lower methylation level of these FKBP5-associated CpG sites in Cushing's syndrome (Fig. 3B). Since promoter methylation usually negatively correlates with gene expression (11), we measured FKBP5 gene expression in a subset of 37 samples (14 overt Cushing's syndrome, 10 mild Cushing's syndrome, 7 eucortisolism, 6 adrenal insufficiency), for which whole blood RNA was available (Supplementary Methods). FKBP5 expression was negatively correlated with the five FKBP5 promoter-associated CpG sites (r = −0.55, P -value < 0.001), and positively with the 24-h urine-free cortisol (Pearson's r = 0.62, P -value < 0.001), demonstrating the potential interest of using FKBP5 expression as a biomarker of glucocorticoid excess.
We next explored the kinetics of glucocorticoidrelated methylome modification by comparing the methylation profile of four different samples available for one of the patients (patient P30), collected before, 4 days, 7 months, and 35 months after Cushing's syndrome correction. Hierarchical clustering of methylome profiles well discriminated the overt Cushing's syndrome sample. During the months following Cushing's syndrome correction, the three samples were properly ordered, showing a progressive overall re-increase of methylation (t-test P -value < 0.05; Fig. 4).

Predicting glucocorticoid status by blood DNA methylation
To select a limited set of CpG sites predicting the glucocorticoid status, we performed a Lasso-penalized linear regression on the training cohort, starting from the 52 727 most variable CpG sites. Twenty-nine CpG sites were selected (Supplementary Table 9), properly discriminating overt Cushing's syndrome in the training cohort, with confirmation in the validation cohort (Fig.  5A). A 29-CpGs methylation predictor was generated by combining the M-values of the 29 selected CpG sites weighted by their Lasso coefficients. This predictor was significantly associated with the glucocorticoid status in the validation cohort (odds ratio: 1.58; 95% CI: 1.25-2.08; P -value < 0.001).
The 29-CpGs methylation predictor was then tested on the ENSAT-HT cohort, a second independent validation cohort (26 overt Cushing's syndrome and 65 eucortisolism samples). Samples were properly classified as Cushing's syndrome and eucortisolism respectively (Fig. 5B) and the prediction value could be confirmed (odds ratio: 1.10; 95% CI: 1.07-1.12; P -value < 0.001), corresponding to an accuracy of 0.84.  One of the 29 selected CpG sites was located in the FKBP5 gene locus. Among the CpG sites in this locus, the one with the highest correlation to the 29-CpGs methylation predictor (Pearson's r = −0.89) (Supplementary Table  10) properly discriminated on its own overt Cushing's syndrome from other samples (Supplementary Fig. 4).
We finally tested to which extent the glucocorticoid effect on blood methylome was related to glucocorticoidinduced white blood cell count variations. In a multivariate model combining the 29-CpGs methylation predictor and the neutrophils proportion, the 29-CpGs methylation predictor remained significant (logistic regression P -value < 0.001; Table 2).

Candidate CpG sites predicting Cushing's syndrome-related complications
In order to address whether blood DNA methylation is associated with specific glucocorticoid-related complications, we performed an exploratory Lasso regression analysis on Cushing's syndrome samples from the 47 patients. A combination of four CpG sites was able to discriminate Cushing's syndrome patients with and without hypertension. Similarly, 14 CpG sites discriminated patients with and without osteoporosis ( Fig.  6 and Supplementary Tables 11, 12). The combination of the methylation level of selected CpG sites for hypertension and osteoporosis was not correlated with 24-h urine-free cortisol (Pearson's r = 0.06 and r = −0.05 for hypertension and osteoporosis, respectively). No combination of CpG sites was able to discriminate Cushing's syndrome patients with and without diabetes.

Discussion
In this study, we demonstrated that whole blood methylome quantified biological glucocorticoid action. This biomarker was able to discriminate glucocorticoid excess from eucortisolism and adrenal insufficiency, independently from hormone assays. This new insight may contribute to overcome common pitfalls in Cushing's syndrome diagnosis and management (35,36). In clinical practice, such a tool would be of limited benefit in case of overt Cushing's syndrome, when clinical signs and hormone assays straightforwardly establish the diagnosis. However, a nonhormonal biomarker, directly measuring glucocorticoid action, could particularly help in three conditions: (i) in patients with mild autonomous cortisol secretion, to decide between surveillance and surgical correction of glucocorticoid excess; (ii) in patients under local or lowdose systemic glucocorticoid treatments, to assess the global glucocorticoid level. Indeed, morning plasma cortisol is often low in these patients, and cannot properly assess the glucocorticoid level; this low glucocorticoid level may either reflect mild glucocorticoid excess with negative feedback on endogenous cortisol production, or adrenal insufficiency resulting from prolonged adrenal blockade (37); and (iii) in patients with adrenal insufficiency, to determine the optimal glucocorticoid supplementation (38). At this stage, though, the performance of our biomarker in these intermediate conditions and in exogenous Cushing's syndrome remains to be established, as well as its clinical relevance.
Here, we have analysed the global methylation state of blood DNA. Blood DNA is easy to obtain and DNA methylation marks are robust and convenient to investigate. In addition, DNA methylation is highly variable, enabling its use as a suitable biomarker.
Our analysis revealed a global DNA hypomethylation signature associated with endogenous Cushing's syndrome, demonstrating the direct impact of glucocorticoids on DNA methylation. DNA hypomethylation was already observed in an experimental model of mice treated with exogenous glucocorticoids (16). We could delineate a    (33,34). However, we could precisely infer white blood cell count from methylome profiles for each sample, and demonstrate that methylome prediction of glucocorticoid status remained significant after adjustment on white blood cell composition and therefore that methylome profiles variations do not only reflect blood composition variations. In addition, the methylome signature provided here in terms of differentially methylated CpG sites and regions is adjusted for white blood cell composition, thus focusing on differences not related to white blood cell composition. This suggests a global impact of glucocorticoids on methylation of some DNA regions, irrespective of their tissue of origin. After correction of glucocorticoid excess, longterm consequences have been reported (39,40,41,42). An easily measurable biomarker reflecting the dynamic of biological changes in time, such as blood DNA methylation, could help monitoring patients during follow-up. We observed that blood DNA hypomethylation progressively recovers in the years following remission. Similarly, a subtle DNA hypomethylation was observed by Glad et al. several years after Cushing's syndrome correction (43). The authors could correlate some blood methylation levels at specific genomic regions with longterm Cushing-associated neuropsychological sequels. Whether blood DNA methylation can properly help to monitor recovery after Cushing's syndrome correction remains to be specifically explored.
In case of eucortisolism and mild Cushing's syndrome, our methylation predictor showed intermediate classification, between overt Cushing's syndrome and adrenal insufficiency, sometimes discrepant from clinical appraisal. In such conditions, it was not possible to assess whether the methylation predictor was more accurate than clinical evaluation. In such cases, the most relevant judgmental criteria would be the correlation with longterm complications. Therefore, we explored the association between methylation and some well-appreciated glucocorticoid-related complications. Some selected CpG sites were identified as associated with hypertension and osteoporosis. New cohorts would be necessary to validate this association. In addition, whether these markers are specific to glucocorticoid excess remains to be established in prospective trials. Indeed, the inclusion of patients with no glucocorticoid excess would serve as negative controls, necessary for discarding markers of hypertension, diabetes, or osteoporosis not related to glucocorticoid excess.
Whole genome methylation profiling is not easily achievable in clinical routine, thus representing a limitation in using this new marker. A technology transfer to targeted methylation assays would be required, such as pyrosequencing, methylation specific-MLPA, or methylation-specific high-resolution melting analysis (44). An alternative would be the identification of surrogate DNA regions recapitulating this global information. One region could be represented by the FKBP5 gene locus, whose methylation and expression have been demonstrated to be modulated by glucocorticoids in different tissues (45,46,47,48,49,50), with hypomethylation of the promoter region associated with increased gene expression. In this study, we also identified the FKBP5 promoter region methylation as strongly associated with glucocorticoid excess and negatively correlated with FKBP5 gene expression. Particularly, one single CpG site from the FKBP5 gene locus could discriminate overt Cushing's syndrome samples.
Sample size is another limitation of this study, with 94 samples from 47 patients. However, samplings were performed at different times of the disease, corresponding to different glucocorticoid statuses. Additional samples from the ENSAT-HT cohort allowed to further validate the performance of our marker. Further extending the cohort would help to confirm this finding, helping to better characterize the association with specific complications of glucocorticoid excess.
In conclusion, glucocorticoids induce a dynamic whole blood DNA methylation signature. This signature could be used as a biomarker for assessing glucocorticoid action independently from hormone assays.

Supplementary materials
This is linked to the online version of the paper at https://doi.org/10.1530/ EJE-21-0907.

Declaration of interest
Guillaume Assié is on the editorial board of EJE. Guillaume Assié was not involved in the review or editorial process for this paper, on which he is listed as an author.

Funding
This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No633983, the Programme Hospitalier de Recherche Clinique 'CompliCushing' (PHRC AOM 12-002-0064), Agence Nationale pour la Recherche (ANR-18-CE14-0008-01), the Else Kröner-Fresenius Stiftung (2012_A103 and 2015_A228 to M R) and the Deutsche Forschungsgemeinschaft (DFG) within the CRC/Transregio 205/1 'The Adrenal: Central Relay in Health and Disease' (to M R, F B, A R).