Abnormal DNA methylation within HPA-axis genes years after paediatric critical illness

Abstract

Background

Critically ill children suffer from impaired physical/neurocognitive development 2 years later. Glucocorticoid treatment alters DNA methylation within the hypothalamus–pituitary–adrenal (HPA) axis which may impair normal brain development, cognition and behaviour. We tested the hypothesis that paediatric-intensive-care-unit (PICU) patients, sex- and age-dependently, show long-term abnormal DNA methylation within the HPA-axis layers, possibly aggravated by glucocorticoid treatment in the PICU, which may contribute to the long-term developmental impairments.

Results

In a pre-planned secondary analysis of the multicentre PEPaNIC-RCT and its 2-year follow-up, we identified differentially methylated positions and differentially methylated regions within HPA-axis genes in buccal mucosa DNA from 818 former PICU patients 2 years after PICU admission (n = 608 no glucocorticoid treatment; n = 210 glucocorticoid treatment) versus 392 healthy children and assessed interaction with sex and age, role of glucocorticoid treatment in the PICU and associations with long-term developmental impairments. Adjusting for technical variation and baseline risk factors and correcting for multiple testing (false discovery rate < 0.05), former PICU patients showed abnormal DNA methylation of 26 CpG sites (within CRHR1, POMC, MC2R, NR3C1, FKBP5, HSD11B1, SRD5A1, AKR1D1, DUSP1, TSC22D3 and TNF) and three DNA regions (within AVP, TSC22D3 and TNF) that were mostly hypomethylated. These abnormalities were sex-independent and only partially age-dependent. Abnormal methylation of three CpG sites within FKBP5 and one CpG site within SRD5A1 and AKR1D1 was partly attributable to glucocorticoid treatment during PICU stay. Finally, abnormal methylation within FKBP5 and AKR1D1 was most robustly associated with long-term impaired development.

Conclusions

Two years after critical illness in children, abnormal methylation within HPA-axis genes was present, predominantly within FKBP5 and AKR1D1, partly attributable to glucocorticoid treatment in the PICU, and explaining part of the long-term developmental impairments. These data call for caution regarding liberal glucocorticoid use in the PICU.

Background

Critical illness in children, which requires treatment in a paediatric intensive care unit (PICU), represents a form of severe physical stress, hallmarked by a whole range of (neuro)endocrine abnormalities [1,2,3,4,5]. A typical response to the severe stress of critical illness is an acute activation of the hypothalamus–pituitary–adrenal (HPA) axis [1, 6, 7]. Indeed, critically ill patients show an acute rise in total and free cortisol, lasting shorter in children than in adults [1, 8,9,10,11]. Unlike long assumed, the rise in cortisol is not driven by a rise in ACTH, as ACTH levels have been shown to be normal or low during critical illness [1, 10, 11]. Instead, the elevated cortisol levels appeared explained by low levels of its binding proteins corticosteroid-binding globulin and albumin, combined with suppressed breakdown of cortisol [1, 11]. Suppressed cortisol breakdown has been revealed by reduced expression and/or activity of the cortisol metabolising enzymes 11β-hydroxysteroid dehydrogenase (11βHSD) 2 and the A-ring reductases, and by strongly reduced cortisol plasma clearance during tracer infusion and after administration of a hydrocortisone bolus [1, 11].

Follow-up studies of children who needed PICU admission for a wide variety of surgical or medical reasons have shown that these children reveal high risk for important impairments in physical and neurocognitive development and behavioural problems years after the acute illness, also in the absence of pre-existing conditions known to affect or possibly affect development [12,13,14,15,16,17]. Other severe early-life adverse events, such as abuse or neglect occurring during several sensitive neurodevelopmental windows, have been associated with similar neurocognitive impairments and behavioural problems and with risk of psychiatric and metabolic diseases later in life [18,19,20]. Some of these associations appear to be sex- and/or age-specific [20].

Research has suggested that stress affects development of brain regions crucial for many aspects of cognition and behaviour through increased glucocorticoid signalling that, when prolonged or excessive, can induce abnormal DNA methylation within different levels of the HPA-axis and glucocorticoid signalling [21,22,23,24]. More specifically, such DNA-methylation changes have been reported for genes encoding corticotropin-releasing hormone (CRH) and its receptor CRHR1, arginine vasopressin (AVP) and its receptor AVPR1b, pro-opiomelanocortin (POMC), corticotropin receptor (MC2R), glucocorticoid receptor GRα (NR3C1) and for the GRα-regulated chaperone protein FK506 binding protein 51 (FKBP5) [23, 25,26,27,28,29,30,31]. In addition, altered methylation of genes encoding the cortisol metabolising enzymes 11βHSD1 (encoded by HSD11B1), 11βHSD2 (HSD11B2) and the A-ring reductases (SRD5A1, SRD5A2, AKR1D1) as documented in brain regions or peripheral tissues was suggested to be involved [32,33,34,35,36,37,38]. Previous work further suggested that GR binding to glucocorticoid-responsive elements (GREs) can induce changes in methylation of GR-target genes, such as dual-specificity phosphatase-1 (DUSP1), annexin A1 (ANXA1), pro-protein convertase-1 (PCSK1), glucocorticoid-induced leucine zipper (GILZ, encoded by TSC22D3) and tumour necrosis factor (TNF)-α (TNF) [39, 40].

In a large genome-wide analysis, we identified de novo alterations in leukocyte DNA methylation in critically ill children that rapidly arose after PICU admission and remained present until PICU discharge [41, 42]. In another large genome-wide analysis on buccal mucosa DNA obtained from those patients 2 years later, we further documented abnormal methylation as compared with healthy children, in pathways known to be important for physical/neurocognitive development [43]. In view of the striking similarities in long-term impairments in physical and neurocognitive development and behavioural problems after critical illness and other severe early-life adverse events and the suggested involvement of abnormal DNA methylation within the levels of the HPA-axis and glucocorticoid signalling in the long-term developmental impairments after such other early-life adverse events, we aimed to look more in-depth into this axis in a pre-planned targeted hypothesis-driven sub-analysis. In the present study, we thus addressed the hypothesis that infants, children and adolescents who have been admitted to the PICU, sex- and age-dependently, show long-term abnormal buccal mucosa DNA methylation within the different levels of the HPA-axis and glucocorticoid signalling, possibly aggravated by treatment with glucocorticoids in the PICU, which may contribute to the long-term physical and neurocognitive impairments and behavioural problems [12, 13].

Methods

Study population and buccal mucosal swab sampling

This study is a pre-planned secondary analysis of the multicentre PEPaNIC-RCT (registered at ClinicalTrials.gov NCT01536275, 2012-2015) and its 2-year follow-up study (2014–2018) [12, 44]. The PEPaNIC-RCT included 1440 consecutive critically ill children aged 0–17 years admitted to the PICUs of Leuven (Belgium), Rotterdam (The Netherlands) or Edmonton (Canada), who had an expected PICU stay of at least 24 h, were at risk of malnutrition, and did not meet any of the exclusion criteria [44]. All patients, who survived and for whom written informed consent was obtained, were eligible for a follow-up study 2 years after PICU admission, to assess physical, neurocognitive and emotional/behavioural development [12]. Healthy children with comparable sex and age distribution as the former PICU patients were included as controls. These were either siblings and relatives of the patients or unrelated children from the same geographical area, attempting to adjust as much as possible for genetic and socio-economic/environmental background [12]. The number of healthy children to be recruited as control group was based on power to detect relevant differences in developmental outcomes [12].

At the 2-year follow-up of the PEPaNIC patients, upon which time also the healthy children were assessed, buccal mucosal swabs (Isohelix, Cell Projects, Harrietsham, Kent, England) were collected following a standardised collection procedure [43]. Swabs were stored in a DNA stabilising solution (DSK kit, Isohelix) at − 80 °C until further processing. All former patients and healthy children from whom a buccal mucosal swab was available were eligible for this DNA methylation study (Additional file 1: Fig. A1). These were patients and healthy children who had been recruited in Leuven and Rotterdam, as the Edmonton centre did not participate in swab collection.

Ethics

The institutional review boards at each participating site approved this study (Ethische Commissie Onderzoek UZ Leuven/KU Leuven: ML8052; Medische Ethische Toetsingscommissie Erasmus MC: NL49708.078). The study was performed in accordance with the 1964 Declaration of Helsinki and its amendments. Written informed consent was acquired from parents, legal guardians and/or the child if 18 years or older.

DNA extraction and DNA methylation data processing

As previously described, DNA was extracted from all available buccal mucosal swabs from patients and healthy children (DDK DNA isolation kit, Isohelix) [43]. DNA concentrations were quantified with the Qubit® 3.0 fluorometer (Thermo Fisher Scientific, Waltham, MA). Two-hundred ng DNA was subsequently subjected to bisulfite-conversion with use of the EZ-96 DNA-methylation-Direct® Kit (Zymo Research, Irvine, CA). Bisulfite-converted DNA was profiled using the Infinium® HumanMethylation EPIC BeadChip (Illumina Inc., San Diego, CA), which interrogates 865,859 CpG sites. Analyses were performed blinded for former PICU patient (subgroup) versus healthy child status.

Data were processed using R statistical software version 4.0.2 using the LICMEpigenetics package (version 0.1.0) and the ilm10b4.hg19 annotation file using human genome 19 (hg19) as a reference was used to annotate the CpG sites [43]. This package contains R functions to exclude low-quality samples (not showing the typical bi-peak curve of the methylation value distribution) and probes (detection p value greater than 0.01 in 50% of the samples, or probes spanning single nucleotide polymorphisms), normalise the methylation data, adjust for batch effect and find differentially methylated positions and regions, as described below [45, 46]. Non-biological or technical variation due to experimental conditions was corrected for by including the first 30 principal components (PCs) of the technical control probes located on the Infinium® HumanMethylation EPIC BeadChip, excluding the negative control probes, as covariates in all multivariable linear regression models, according to the method developed by Lehne et al. [47].

Selected genes

Figure 1 gives a schematic overview of the studied genes encoding key proteins within the different levels of the HPA-axis (CRH, CRHR1, AVP, AVPR1b, POMC, MC2R, NR3C1, FKBP5, HSD11B1, HSD11B2, SRD5A1, SRD5A2, and AKR1D1) and the selected GR-target genes (DUSP1, ANXA1, PCSK1, TSC22D3, and TNF). Methylation of these genes was studied in buccal mucosa as a surrogate for the tissues indicated in the figure, which are not accessible for obvious reasons. The CpG sites within these genes analysed in this study are listed in Additional file 1: Table A1.

Fig. 1

Schematic overview situating the studied genes within the different levels of the HPA-axis and glucocorticoid signalling. Please note that the DNA methylation status of the indicated genes was studied in buccal mucosa as a surrogate for the tissues in the figure. Created with BioRender.com. Abbreviations: GRE: glucocorticoid-responsive element

Statistical analyses

Patient demographics and medical characteristics are reported as median and interquartile range or as number and percentage. Group comparisons were performed with the Chi-square test for categorical variables and with the Wilcoxon signed rank test for continuous variables.

DNA methylation changes between former PICU patients and healthy children

Differences in DNA methylation between former PICU patients and healthy children were assessed, correcting for multiple testing by applying a false discovery rate (FDR) of less than or equal to 0.05 as determined with the Benjamini–Hochberg procedure [48]. Such differences, if any, are to be considered the sum of differences evoked by the paediatric critical illness and intensive medical care and those that may have been present in the former patients prior to PICU admission.

First, all differentially methylated positions (DMPs) in former PICU patients as compared with the healthy children were identified. For each CpG site, methylation status was compared between former patients and healthy children with use of a multivariable linear model, using the limma framework, adjusting for baseline risk factors [age, centre, race, sex, geographical origin, history of malignancy and predefined syndrome (rationale provided in Additional file 1: Method A1, definition of syndrome provided in Additional file1: Method A2)] and adjusting for technical variation (batch effect) [47, 49]. Next, all differentially methylated regions (DMRs), i.e. regions within the DNA where groups of neighbouring CpG sites are differentially methylated, in former PICU patients as compared with the healthy children were identified with the DMRcate package [50]. A stepwise explanation of the DMRcate method with an illustrative example is provided in Additional file 1: Method A3 [43].

Interaction with sex and age

We further assessed, for the above identified DMPs in former PICU patients as compared with healthy children, whether there was any interaction with sex and/or the age ‘at exposure’. These analyses allow to investigate whether any difference between former PICU patients and healthy children would depend on sex or ‘age at exposure’. For the former PICU patients, age ‘at exposure’ was the age at admission to the PICU. For the healthy children, the age ‘at exposure’ was imputed by subtracting 2 years from their age at 2-year follow-up assessment [12]. Assessment of interactions was done by adding an interaction term in the multivariable linear models described above.

Role of glucocorticoid treatment during PICU stay

We next investigated to what extent glucocorticoid treatment during PICU stay may have played a role in bringing about or aggravating any of the identified DMPs 2 years after PICU admission. To this end, we performed multivariable analyses with the limma framework, among former PICU patients, with use of the methylation status of each of the CpG sites identified above as the dependent variables, comparing glucocorticoid treatment versus no glucocorticoid treatment in the PICU. We adjusted for baseline risk factors and technical variation as described above, and for length of PICU stay, admission diagnosis, severity of illness [Paediatric Index of Mortality 3 (PIM3) score, Paediatric Logistic Organ Dysfunction (PeLOD) score], randomisation to one of two nutritional strategies and risk of malnutrition [STRONGkids score] (rationale provided in Additional file 1: Method A1).

Association with physical and neurocognitive/behavioural development

For each of the above identified DMPs present in former PICU patients as compared with healthy children, we next assessed with use of multivariable linear models, whether methylation status in patients was associated with the physical (weight, height and head circumference) and neurocognitive/behavioural outcomes (executive functioning, emotional/behavioural problems, intelligence, visual motor integration, alertness and memory) that were impaired as compared with healthy children at 2-year follow-up (a detailed description of these outcome measures is provided in Additional file 1: Method A4) [12], adjusting for all above listed baseline risk factors and technical variation, with the addition of ‘linguistic origin’, since this is a known confounder for the functional outcome assessment (rationale for each risk factor provided in Additional file 1: Method A1). To assess robustness of these linear models, we first performed a tenfold cross-validation and computed the p values of the 10 test folds using Fisher’s method [51]. This process was repeated in 100 iterations. The percentage of iterations with a significant p value across the tenfold of the cross-validation (α ≤ 0.05) was calculated and visualised in a heatmap. Additionally, with use of the mean coefficient from the multivariable models, we assessed whether the abnormal methylation that we observed in former PICU patients was associated with either better (B) or harmful (H) outcome.

Results

Sample collection and quality assessment

Buccal mucosal swabs were collected from 821 patients and from 392 healthy children (Additional file 1: Fig. A1). Three patients needed to be excluded for further analysis, as DNA yield was insufficient for 1 patient, and 2 samples from patients showed deviation from the typical bi-peak curve of the methylation value distribution. Of the 818 patients remaining for the analyses, 210 had received glucocorticoid treatment during PICU stay. Participants’ demographics and medical characteristics are shown in Table 1.

Table 1 Demographics and medical characteristics of participants included in the 2-year epigenetic follow-up analysis

DNA methylation differences between former PICU patients and healthy children

The adjusted comparison of the former PICU patients with the healthy children identified 26 DMPs within the HPA-axis genes. Their location within the genome, and the direction of change (hypo- or hypermethylated in former patients as compared with healthy children) are listed in Table 2. The absolute mean differences in DNA-methylation beta-values for the DMPs were 2.2% (SD 1.5%), ranging up to 5.5%, with mostly hypomethylation in patients (20/26 (76.9%), Table 2, Additional file 1: Table A1, Additional file 1: Fig. A2). The 26 DMPs in former PICU patients as compared with healthy children were located in 11 out of the 18 selected genes (61.1%), more specifically CRHR1 (1/41 CpG sites (2.4%)), POMC (1/27 CpG sites (3.7%)), MC2R (1/20 CpG sites (5.0%)), NR3C1 (3/89 CpG sites (3.4%)), FKBP5 (7/51 CpG sites (13.7%)), DUSP1 (1/33 CpG sites (3.0%)), TSC22D3 (2/50 CpG sites (4.0%)), TNF (3/27 CpG sites (11.1%)), TSC22D3 (1/50 CpG sites (2.0%)) and HSD11B1 (2/25 CpG sites (8.0%)), SRD5A1 (1/31 CpG sites (3.2%)), AKR1D1 (4/15 CpG sites, (26.7%)).

Table 2 Abnormal DNA methylation within the HPA-axis

Three DMRs were identified, located in the AVP gene (containing 4 CpG sites with a width of 665 bp), the TSC22D3 gene (containing 8 CpG sites with a width of 438 bp) and the TNF gene (containing 12 CpG sites with a width of 146 bp) (Additional file 1: Table A2, Additional file 1: Fig. A3). The DMRs in AVP and TNF were hypomethylated, whereas the DMR in TSC22D3 was hypermethylated in former PICU patients as compared with healthy children.

Interaction with sex and age

The identified differences in DNA methylation between former PICU patients and healthy children were not sex-dependent (Table 2, Additional file 1: Table A3). Age ‘at exposure’ affected the methylation status of only 4 of the differentially methylated CpG sites (Table 2, Fig. 2, Additional file 1: Table A4). These CpG sites were located in the 5’UTR region of FKBP5 (cg03591753 and cg03546163) and the promoter/5’UTR region (depending on transcription variant) of HSD11B1 (cg06571187 and cg14139038). The degree of hypomethylation of these CpG sites in former PICU patients as compared with healthy children increased with age.

Fig. 2

Interaction with age of abnormal DNA methylation in former PICU patients 2 years after critical illness vs healthy children. Methylation profiles are shown for each of the CpG sites that, adjusting for baseline risk factors and technical variation, showed a significant interaction between differential methylation in former PICU patients (n = 818) vs healthy children (n = 392) and ‘age at exposure’ (interaction p values obtained with multivariable linear regression analyses using the limma framework: cg03591753 p = 0.0060, cg03546163 p = 0.00063, cg06571187 p = 0.027, cg14139038 p = 0.0020). Univariate linear regression lines are drawn for methylation status of the CpG site in function of age at exposure, for former PICU patients (red) and healthy children (blue), where the shaded area represents the 95% confidence interval. Abbreviations: PICU: paediatric intensive care unit, UTR: untranslated region

Association between glucocorticoid treatment during PICU stay and DNA methylation differences

We observed that three of the identified DMPs within FKPB5 (cg23416081, cg15929276 and cg01839003), one DMP within SRD5A1 (cg22911074) and one DMP in AKR1D1 (cg02118020) appeared to be aggravated by glucocorticoid treatment given during PICU stay (Table 2, Fig. 3, Additional file 1: Table A5). The absolute differences in DNA methylation beta-values for these CpG sites, comparing former patients who had received glucocorticoid treatment with those who had not, were 3.1% (SD 2.3%), ranging up to 5.2%.

Fig. 3

Role of glucocorticoid treatment in the PICU as a contributor to abnormal DNA methylation within HPA-axis genes 2 years later. The boxplots show a univariate presentation of the methylation status (β value) of the CpG sites of which abnormal methylation was aggravated by glucocorticoid treatment, as identified in multivariable analyses adjusting for baseline risk factors and technical variation (p values obtained with multivariable linear regression analyses using the limma framework: cg23416081 p = 0.017, cg15929276 p = 0.0045, cg01839003 p = 0.036, cg22911074 p = 0.018, cg02118020 p = 0.017). Former PICU patients who had received glucocorticoid treatment in the PICU (n = 210) are depicted in orange and those who had not (n = 608) are depicted in green. The central lines of the boxplots depict the medians, the boxes the interquartile ranges, and the whiskers are drawn to the furthest point within 1.5 times the interquartile range from the box. Abbreviations: UTR: untranslated region

Association between altered DNA methylation and impaired physical and neurocognitive/behavioural development

The statistical association of DNA methylation differences between former PICU patients and healthy children with the impaired long-term physical growth and neurocognitive/behavioural functioning is visualised as a heatmap in Fig. 4. In the heatmap, we counted 325 significant associations of which 291/325 (89.5%) were found to be robust, being significant in 100% of the repetitions. Most associations were found for hypomethylation of CpG sites within the promoter/5’UTR of FKBP5 in former patients, which statistically associated with more impaired physical growth, intelligence, visual motor integration, alertness and memory, and to a lesser extent with executive functioning and behavioural problems. Associations with impairments in the same domains were also found for abnormal methylation within the gene body of AKR1D1. For the other genes, mostly more discrete associations of abnormal methylation with more developmental impairments were observed, except for abnormal methylation within CRHR1 and TSC22D3, which appeared to associate with less pronounced impairments in intelligence and/or VMI.

Fig. 4

Association of long-term abnormal methylation within HPA-axis genes in former PICU patients with impaired long-term developmental outcomes. This heatmap summarises the results of the multivariable linear regression analyses assessing associations between abnormal methylation within HPA-axis genes 2 years after critical illness and developmental outcomes of former PICU patients evaluated at that time (n = 713), adjusted for baseline risk factors and technical variation. Each box shows the percentage of the 100 iterations for which the p values across the 10 test folds (computed with Fisher’s method [51]) are significant (p < 0.05) for a given CpG site (X axis) and outcome (Y axis). Darker colour intensity indicates a larger percentage of iterations being significant. For each CpG site, the location within the respective gene is also indicated (due to multiple splice variants, multiple gene locations are possible). The labels within the boxes indicate whether the observed abnormal methylation correlates with a better score for the long-term outcome (benefit ‘B’), or with a worse score (harm ‘H’). Abbreviations: IQ: intelligence quotient, SD: standard deviation, UTR: untranslated region, VMI: visual motor integration

Discussion

Two years after admission to the PICU, children who had been critically ill revealed abnormal DNA methylation within several genes of the HPA-axis. This was the case for genes encoding the CRH receptor, the ACTH precursor POMC and the ACTH receptor, the glucocorticoid receptor GRα and its co-chaperone FKBP5, several cortisol metabolising enzymes (11ß-HSD1, 5α-reductase and 5β-reductase) and three GR-regulated proteins (DUSP1, GILZ and TNFα). The most extensive changes were observed for the FKBP5 and AKR1D1 genes which were also substantially aggravated by glucocorticoid treatment in the PICU. Most abnormal DNA methylation in former patients was hypomethylation, was sex- and age-independent, and found to be independently associated with impaired physical and neurocognitive development.

Several critical illness- and treatment-related factors could be potential drivers for altered DNA methylation. Examples are hyperglycaemia and insulin resistance [52,53,54,55,56], neuroendocrine abnormalities [1,2,3,4,5, 57,58,59] and exposure to endocrine-disrupting plasticisers leaching from indwelling medical devices [60,61,62]. Also, as critically ill patients are unable to eat, artificial feeding is often provided [63]. Feeding management of critically ill children has been shown to affect the DNA methylome up until PICU discharge [41, 42]. Since prolonged or excessive stress can induce abnormal DNA methylation within different levels of the HPA-axis and glucocorticoid signalling in the context of other adverse early-life exposures [21,22,23,24, 59], the development of such abnormalities after paediatric critical illness seemed plausible. Our present findings are indeed in line with those from studies on the impact of other forms of early-life stress which revealed mostly hypermethylation in NR3C1 and hypomethylation in FKBP5 within the gene body (intron 7) and the promotor [64, 65]. Via research performed in humans in combination with experiments in cells, it has previously been shown that stress reduces methylation in the FKBP5 gene, interestingly at a CpG site which we also found to be hypomethylated, whereby FKBP5 expression was found to be upregulated [66]. Of importance in the context of critical illness, a study in mice revealed that lipopolysaccharide injection increased FKBP5 expression in the hippocampus, driving increased neuroinflammation [67]. Whether the stress-induced rise in endogenous glucocorticoids mediated this effect is currently unclear. However, in patients with Cushing syndrome, it has been shown that long standing excessive endogenous hypercortisolism induces hypomethylation in FKBP5, explaining their long-lasting psychopathological sequelae [24]. In experimental models, also exogenous hypercortisolism has shown to reduce DNA methylation in FKBP5 in neuronal cells, and in vivo this was associated with increased FKBP5 expression across several brain regions [68]. These data provided a plausible explanation for our finding that hypomethylation in the FKBP5 gene was aggravated in patients who 2 years earlier had been treated with glucocorticoids in the PICU, taking into account potential confounders affecting the need of glucocorticoid treatment. It currently remains unclear what brings about the hypomethylation, though downregulated DNA methyltransferase-1 and cleavage of the DNA backbone directly by the GRα have been suggested to play a role [39, 69,70,71]. Some of the differentially methylated CpG sites that we found within NR3C1 and FKBP5 have also been described with other forms of stress. For example, cg15910486 within the promotor/5’UTR of NR3C1 was reported to be hypermethylated in children who experienced early-life adverse events; cg23416081 and cg15929276 located within the 5’UTR region of FKBP5 have previously been shown to alter cortisol reactivity and behaviour in children; cg20813374 located in the promotor/5’UTR of FKBP5 is a known stress-related epigenetic signature previously associated with myocardial infarction and inflammation; and cg03546163 in the 5’UTR of FKBP5 was shown to be differentially methylated in the context of Cushing’s syndrome [24, 65, 66, 72].

Also glucocorticoid metabolism can be altered by early-life stress [33]. We here found that genes encoding cortisol metabolising enzymes were differently methylated in former PICU patients, more specifically for 5α- and 5β-reductase and 11β-HSD1. These changes could again be partly explained by glucocorticoid treatment during PICU stay years before. The methylation differences within HSD11B1, SRD5A1 and AKR1D1 may result in altered systemic cortisol availability and hereby affect development as was suggested by our functional outcome analysis. The HSD11B1 hypomethylation we observed in former patients appeared to some extent protective against rather than contributing to some of the developmental impairments. However, a strong harmful association was found between most functional outcome measures and the abnormal DNA methylation in AKRD1, which encodes the 5β-reductase enzyme, an association that has not been reported before.

We also observed altered DNA methylation in former patients within the genes encoding CRH receptor, ACTH and its receptor MC2R, and three GR-regulated proteins, DUSP1, GILZ and TNFα, which has also been reported for other conditions of early-life stress [26, 28,29,30, 39, 40].

The observed differences in DNA methylation in former PICU patients as compared with healthy children were largely age- and sex-independent. However, we did find that the vulnerability to the illness-induced hypomethylation of 2 CpG sites within FKBP5 and 2 CpG sites within HSD11B1 was more pronounced with higher age at exposure. It has previously been shown that neonates and young infants are indeed relatively stress hyporesponsive, interpreted as a crucial protection against glucocorticoid-induced harmful effects on the developing brain [20, 73]. Higher vulnerability to illness-induced DNA-methylation changes in older than in younger children, in particular from age of adrenarche onwards and into puberty, has previously been shown by our group [74].

This study has some strengths and weaknesses to highlight. The large sample size and the multicentre, prospective study design with predefined long-term assessments of former PICU patients and healthy children were strengths. In addition, our methodology applying tenfold cross-validation over 100 iterations reduced the odds of findings by chance and reduced the impact of outliers. Our study also has some limitations. First, we have studied DNA methylation in buccal mucosa, whereas the HPA-axis links well-defined brain structures with the adrenal cortex via hormonal regulation. Studying the DNA-methylation markers of stress is obviously not possible in the physiological tissues and organs where it exerts its effects through the HPA-axis. It was thus for pragmatic reasons that we analysed buccal mucosa, as these epithelial cells are accessible to clinical research. These cells are also exposed to stress hormones, but we do not know whether the methylation changes observed in these cells reflect those supposed to occur in the HPA-axis itself. This difficulty has been encountered in other studies, where blood cells or buccal cells have been used as proxies to physiological tissues [23,24,25, 29, 31, 40,41,42,43, 66]. Second, due to lack of a sample before PICU admission, we cannot discriminate between pre-existing abnormal methylation and abnormal methylation induced by the critical illness. Third, although we preferentially recruited siblings and relatives of the patients to the control group, extended with unrelated children from the same geographical area and adjusted as much as possible for baseline risk factors, we cannot exclude residual confounding by genetic background and environment. Fourth, the treatment with glucocorticoids during PICU stay had not been randomised, arguing for caution when interpreting these results. Indeed, specific conditions triggering the need for glucocorticoid treatment theoretically may confound these results. However, with the extensive adjustment for risk factors which also included type, severity and duration of illness we aimed to reduce as much as possible the impact of this limitation. Nevertheless, we cannot exclude that there may be some residual unmeasured confounding in these and the other multivariable analyses. Finally, we were not able to assess what effect the abnormal DNA methylation might have on gene transcription, nor on cortisol or ACTH levels, as we did not have the samples for these analyses. However, differential methylation in the order of magnitude as observed in the present study was associated with differential gene expression in our earlier study on DNA methylation in muscle of adult critically ill patients and controls [75].

Conclusions

Two years after critical illness in children, buccal mucosa DNA revealed abnormal methylation of CpG sites within genes of the HPA-axis, most extensively within the FKBP5 and AKR1D1 genes, which occurred largely independent of sex and age. In addition, glucocorticoid treatment while in the PICU was found to be associated with an aggravation of the methylation changes in FKBP5 and AKR1D1 detected 2 years later. The observed abnormal DNA methylation within these genes in former PICU patients statistically explained part of the long-term physical and neurocognitive developmental impairments. These findings call for attention regarding safety of a liberal glucocorticoid use in the PICU in the absence of strong underlying evidence of benefit.