Introduction

Approximately 10% of infants are born prematurely (<37 weeks’ gestation) in the United States (U.S.).1 Globally, prematurity is the leading cause of death in children under 5 years of age.2 Bronchopulmonary dysplasia (BPD), the most common neonatal morbidity among very preterm infants, is a respiratory disorder in which infant lung development is disrupted by mechanisms that include oxidative stress and inflammation, necessitating treatment with supplemental oxygen, assisted ventilation, and/or pharmacological interventions, such as anti-inflammatory corticosteroids, to improve gas exchange. However, these interventions can result in oxidative stress and perturb the hypothalamic-pituitary-adrenal (HPA) axis, an innate human stress response system, in patients with BPD increasing the physical and physiological stress. Infants born very preterm (<30 weeks’ gestational age) are at a particularly high risk for developing poor neuromotor, cognitive, and behavioral outcomes that can persist through adulthood.3,4,5,6,7,8,9,10,11,12,13,14 Additionally, very preterm infants with BPD are at increased risk for respiratory conditions and neurodevelopmental impairments.15,16 Very preterm infants also have longer stays in the neonatal intensive care unit (NICU) and more hospital encounters following NICU discharge compared with infants without BPD.17 Despite the public health burden of preterm birth and BPD, there is a dearth of tools that can be easily used at NICU discharge to understand the ramifications of prematurity, and existing clinical indicators are limited in their ability to predict longer-term developmental impairments. Thus, there is a need to better understand the links between acute neonatal morbidities, such as BPD, and long-term negative health outcomes.

DNA methylation (DNAm) is a biological process in which methyl groups are added to DNA at cytosine (C) bases that are followed by guanine (G) bases, referred to as CpG sites. This process influences the expression potential of genes but does not change the sequence of DNA. Preterm infants have different methylation patterns than those born at term, and these differences can be detected in multiple tissues, with special emphasis on minimally invasive sampling from tissues such as the placenta, cord blood, and saliva.18,19,20 Additionally, children’s DNAm profiles capture information about their early life environmental exposures and health conditions.21 Our group has shown that neonatal neurobehavioral responses22 and an increasing number of neonatal morbidities, including BPD,23,24 are associated with differences in DNAm throughout the genome, specifically, at NR3C1 which encodes a primary glucocorticoid (GC) receptor.25 Another study observed a relationship between BPD and an epigenetic marker of biological aging.26 Additionally, several epigenome-wide studies of peripheral tissues collected at birth, such as cord blood and placenta, have identified differential DNAm associated with BPD.27,28,29 Thus, early life exposures and morbidities, including BPD, appear to leave a fingerprint on the early life epigenome, even in peripheral tissues. These epigenetic associations may provide information on important developmental systems whose epigenetic patterns contribute to BPD, and/or have been affected by these early exposures and experiences.

Environmental, biological, and psychosocial stressors during critical windows of development have been associated with changes in DNAm and linked to poorer health and disease.30,31,32 Epigenetic alterations to the stress response systems may play a critical role in explaining some of these long-term effects. For example, a cross-tissue polyepigenetic GC score was developed by Provençal et al. with the goal of creating a tool to study how early life stress may impact psychiatric outcomes later in life, and in subsequent analyses the GC score was associated with exogenous GC exposure, maternal stress, and a higher risk of children’s negative mental and behavioral outcomes.32 In humans, exposure to dexamethasone and hydrocortisone, GC agonists, alter DNAm in peripheral tissues, such as blood33,34 and buccal cells.35 Similarly, in hippocampal progenitor cells, dexamethasone induces persistent changes in DNAm and gene expression when the exposure occurs during cellular proliferation and differentiation,32 which happens between gestational weeks 10 and 25 in humans. Additionally, DNAm in human buccal tissue was found to be associated with high-dose GC exposure, specifically in the FKBP5 and NR3C1 genes, which are vital for HPA axis function.35

Among neonates born very preterm, infants with BPD are among those with the greatest exposure to exogenous GC,36 suggesting that BPD may be associated with altered DNAm at stress response genes. We also hypothesized that epigenetic differences in the stress responses system would be detectable in peripheral tissues like buccal cells, since stress response is a systemic process and can impact the activity of cells throughout the body.37 However, the epigenetic effects of BPD and steroid treatment on HPA axis genes in infants born preterm has not been well studied, although evidence suggests that the early environment can permanently influence the genome through epigenetic mechanisms and modify endocrine function, metabolism, and behavior of offspring.38,39

This exploratory and hypothesis-generating study aimed to demonstrate whether an innovative epigenetic biomarker for early life endogenous and exogenous GC exposure or CpG-specific methylation of HPA axis genes were associated with BPD severity or with antenatal steroid treatment in very preterm infants. Exogenous steroid treatment is an exposure that could directly interact with the HPA axis system, while BPD itself exerts physiological stress changes that potentially triggers the HPA system and impacts the methylome of infants born very preterm. We hypothesized that the polyepigenetic GC score and methylation levels of 8 HPA axis genes would be altered in association with antenatal steroid exposure and increasing BPD severity.

Methods

Characterization of the study population

The goal of Environmental influences on Child Health Outcomes (ECHO) program is to investigate the effects of environmental exposures on child health and development.40,41,42 The two ECHO cohorts involved in this study are prospective studies of early life influences on child health outcomes among infants born very preterm for which epigenetics data were available.

Demographics of the Neonatal Neurobehavior Outcomes in Very Preterm Iinfants (ECHO NOVI) cohort

A total of 538 infants in the ECHO NOVI study had DNAm, BPD, antenatal steroid, and post-menstrual age (PMA) data available. In this cohort, gestational age (GA) at birth ranged from 22 to 29.86 weeks (mean = 26.99) and PMA ranged from 32.14 to 51.43 weeks (mean = 39.15 weeks) (Table 1). ECHO NOVI is a diverse cohort based on maternal self-reported race and ethnicity, where 52% of infants were White (n = 279); 23% were Black (n = 121); 9.5% were Other Race (n = 51); 7.6% were Asian (n = 41); 7.1% were Hawaiian/Pacific Islander (n = 38); >1% were Native American (n < 5); and <1% did not report race (n < 5). Overall, 22% of infants were of Hispanic ethnicity. Approximately 21% of infants were outborn (N = 112), with 74 of these infants having a BPD diagnosis (66%). In this cohort, 7% of infants were classified as having fetal growth restriction, based on birthweight-for-gestational-age percentile. Over 51% of infants had BPD, with 130 (24.2%), 122 (22.7%), and 25 (4.2%) classified as mild, moderate, and severe BPD, respectively. Almost 89% of the infants’ mothers were treated with antenatal steroids during pregnancy (N = 477); 15% of infants were treated with steroids for BPD (N = 81); and 12% of infants were exposed to both antenatal steroids and steroids for BPD (N = 67).

Table 1 Distribution of demographic characteristics, neonatal morbidities, and maternal/fetal characteristics of the ECHO NOVI and ELGAN ECHO cohort populations.
Full size table

Demographics of the Extremely Low Gestational Age Newborn (ELGAN ECHO) cohort

DNAm, GA, and antenatal steroid data were available for 365 infants in the ELGAN ECHO cohort (Table 1). The average GA of this cohort was 26.07 weeks (range = 23.00–27.86 weeks). Approximately 63% of infants were White (n = 230); 27% were Black (n = 97); 4.3% were Other Race (n = 16); >2% were Asian (n < 5); >1% were Native American (n < 5); 3% were Multiracial (n = 11); and >1% did not report race. In this cohort, 10% of infants were of Hispanic ethnicity (n = 35). All infants were born in the hospital where NICU care was provided (no outborn). Overall, 5% of infants were classified with fetal growth restriction. Approximately 54% of infants had BPD, with 142 (38.9%), 31 (8.5%), and 23 (6.3%) classified as mild, moderate, and severe BPD, respectively. Over 90% of the infants’ mothers were treated with antenatal steroids during pregnancy (N = 331).

Description of study samples

This study was implemented using data from two cohorts of very preterm infants who were participating in the ECHO Program: the ECHO NOVI study and the ELGAN ECHO study. The timeline of sample collection for each cohort is shown in Fig. 1. The ECHO NOVI cohort includes 704 infants who were born between 2014–2016 at <30 weeks’ gestation and were cared for in 9 NICUs in the U.S. Prenatal demographic, socioeconomic, physical and mental health, and substance use data were collected during hospital stays and parent interviews. We obtained inform consent from the mothers of 651 study participants for buccal swab collection prior to discharge from the NICU, and sufficient high-quality DNA was extracted and profiled for DNAm from 542 of these buccal samples. In the ECHO NOVI cohort, BPD was defined as requiring supplemental oxygen at or before 36 weeks PMA and was classified into three levels of severity: mild (oxygen alone or high-flow nasal cannula without supplemental oxygen), moderate (high-flow nasal cannula with supplemental oxygen or continuous positive airway pressure or nasal intermittent mandatory ventilation), and severe (mechanical ventilation or high frequency ventilation through an endotracheal tube).43

Fig. 1: Timeline of ECHO cohort sample collection.
figure 1

ELGAN blood spot samples were collected at birth (gestational age) before BPD diagnosis and treatment. NOVI buccal swabs were collected at NICU discharge (post-menstrual age) after BPD diagnosis and treatment.

Full size image

The ELGAN ECHO study44 includes 1506 neonates born between 2002–2004 at <28 weeks’ gestation in 14 U.S. hospitals. A total of 1222 infants survived to discharge from neonatal intensive care, and 1198 survived to 2 years. Approximately 1190 are presumed alive and currently range in age from 18 to 20 years. Prenatal demographic, socioeconomic, physical health, pregnancy complication, and neonatal medical complication data were collected from review of medical records during NICU stays and maternal interviews around the time of birth of the ELGAN study participant. Neonatal dried blood spots were collected during the first several days of life.

Inclusion/exclusion criteria for this study

Among the participants in ECHO NOVI and ELGAN ECHO, those with missing DNAm, BPD, and antenatal steroid data, and age at sample collection (GA and/or PMA) were excluded from this study. In the ECHO NOVI cohort, 4 of the 542 participants were excluded due to missing antenatal steroid data; therefore, 538 infants were included in this study. In the ELGAN ECHO cohort, 857 of the 1222 infants were excluded, primarily due to missing DNAm, leaving 365 infants included in the study.

DNA methylation (DNAm) measurement, quality control, and preprocessing

Both cohorts used their pre-developed preprocessing and quality control pipelines that aligned with their previously published work. To assess the methylation profiles of the ECHO NOVI study cohort, we analyzed the DNAm data from buccal tissue at NICU discharge as measured using the Illumina MethylationEPIC microarray. These data have undergone established quality control and processing pipelines. We excluded samples with >5% of probes yielding detection p values > 1.0E−5, with mismatch between reported and predicted sex, or incomplete phenotype data.23 Functional normalization and beta-mixture quantile normalization were performed,45 and we excluded probes located on the X and Y chromosomes, those with single-nucleotide polymorphisms (SNPs) within the binding region, and those able to cross-hybridize to other regions of the genome.46,47 After these exclusions, a total of 706,278 probes were available from 538 ECHO NOVI samples.

The ELGAN ECHO cohort used similar quality control procedures, such as the exclusion of blood spot samples with a detection p value > 0.01 for greater than 10% of samples, the removal of sex mismatches, and functional normalization using minfi48 and ShinyMethyl49 R packages. Final beta values were generated following batch effect correction using PCA checks to identify batches and ComBat50 to remove these effects. After exclusions, a total of 786,267 probes were available from 365 ELGAN samples.

Data harmonization

The definition of BPD severity (4-level factor) was harmonized between the ECHO NOVI and ELGAN ECHO cohorts, with a 4-level factor for BPD severity: no BPD (no oxygen requirement), mild BPD (supplemental oxygen or high-flow nasal cannula), moderate BPD (supplemental oxygen via high-flow nasal cannula, nasal continuous positive airway pressure or nasal positive airway pressure ventilation), and severe BPD (mechanical ventilation via endotracheal tube) at ≥36 weeks PMA.51,52 GA was calculated in weeks, and PMA at NICU discharge was defined as the number of days from the last day of the mother’s menstrual cycle to the day of discharge from the NICU, also calculated in weeks. Outborn status,53 maternal smoking during pregnancy, and steroid variables were binary “yes” versus “no.” Fetal growth restriction was also binary with “yes” defined as a birthweight for GA more than 2 standard deviations below the population-specific growth curve. Antenatal steroids were given to mothers who were at risk of preterm delivery to accelerate fetal organ maturation and the variable was defined as a binary “yes” or “no.”

GC score calculation

The neonatal GC score was calculated by multiplying the DNAm beta values by the 24 CpG weights from the cross-tissue polyepigenetic GC score established by Provençal et al.32 and summing those weighted methylation levels for each infant in both cohorts. The weights represent the coefficients from an elastic net regression model with lower GC score indicating greater GC exposure.32 Additionally, the GC score was approximately normally distributed, but with a minimal right skew. The estimated polyepigenetic GC scores from the neonatal DNAm data in the ECHO NOVI cohort yielded scores that ranged from −1.78 to −1.10 with a median score of −1.61. While the GC scores from the neonatal data in the ELGAN ECHO cohort ranged from −0.375 to −2.09 with a median score of −1.19.

HPA axis genes

Genes previously associated with the HPA axis were identified via literature review.54 The HPA system involves complex interactions across many genes. We focused on a subset of genes that have been previously studied in relation to prospective outcomes, and we identified 8 genes for this analysis: FKBP5, NR3C1, NR3C2, CRH, CRHR1, CRHR2, POMC, and HSP90AA1. To capture CpGs annotated to these genes, the UCSC Genome Browser (hg19)55 was used to find regions ±2500 bp beyond the start and stop of each gene. The regions for each gene are as follows: FKBP5 (chr6: 35538862–35698897), NR3C1 (chr5: 142654996–142817577), NR3C2 (chr4: 148076262–148447198), CRH (chr8: 67086112–67093346), CRHR1 (chr17: 45781820–45838328), CRHR2 (chr7: 30649442–30684956), POMC (chr2: 25381222–25394059), and HSP90AA1 (chr14: 102544575–102608586). A total of 344 CpGs were annotated to these genes. Due to differences in quality control and processing steps between the two cohorts, the ECHO NOVI cohort captured 317 of these CpGs while the ELGAN ECHO cohort captured 294 CpGs.

Data analysis

For both cohorts, DNAm was profiled using the Illumina MethylationEPIC BeadArray, which measures over 850,000 CpG loci at a single-nucleotide resolution. After data cleaning and harmonization, our analysis included 538 infants from the ECHO NOVI cohort and 365 infants from the ELGAN ECHO cohort with complete data. Since DNA was profiled using two different tissues (neonatal blood spot samples and buccal cells) and two different sample collection time points (near birth or at NICU discharge), we performed and interpreted all analyses separately within the ECHO NOVI and ELGAN ECHO cohorts while adjusting for the same covariates consistent with the sample collection timeline and relevant exposures (Fig. 1). We used complete cases for children with DNAm and phenotype data (BPD and steroid treatment). We used directed acyclic graphs to identify potential confounders (Supplementary Fig. 1). First, we tested whether the rank-normalized neonatal GC scores were associated with BPD and increasing BPD severity while adjusting for potential confounders, including weeks of PMA at NICU discharge, GA at birth, sex, study site, outborn status, fetal growth restriction, maternal age, maternal smoking during pregnancy, cell type proportions (estimated from reference methylomes),56,57,58,59,60 and maternal antenatal steroid exposure. Buccal tissue from the ECHO NOVI cohort consisted primarily of epithelial cells (>90%), and the blood spot samples from the ELGAN ECHO cohort consisted of CD8 T cells, CD4 T cells, natural killer cells, B cells, monocytes, granulocytes, and nucleated RBCs. Next, we explored whether BPD severity and steroid exposure were associated with differential methylation of the 8 HPA genes of interest while adjusting for the same confounders. Lastly, we stratified the data for each cohort by sex to investigate whether there were sex-specific associations between HPA DNAm, BPD severity, and steroid treatment. All statistical analyses were performed in the R statistical environment, and all regression models were fit with generalized estimating equations. We generated coefficients, confidence intervals, and p values. We used a significance threshold of α = 0.05 to assess statistical significance for the GC score models, whereas we used a 5% FDR threshold to account for multiple testing in the CpG-specific analyses of HPA genes; all tests were two-sided.

Results

GC score and BPD

After regressing the GC score on BPD and BPD severity while adjusting for antenatal steroids and steroids for BPD, we observed that the GC score decreased with increasing BPD severity in the ECHO NOVI cohort. After adjusting for additional confounders, the inverse relationship between BPD severity and GC score was not statistically significant in this cohort (Fig. 2a). The proportion of epithelial cells was strongly associated with GC score and likely confounded our initial findings in the ECHO NOVI cohort. In the ELGAN ECHO cohort, BPD severity was not associated with the GC score in any of the models, although GC score estimates decreased with increasing BPD severity (Fig. 2b), similar to what we observed in the ECHO NOVI cohort. Cell proportion for both monocytes (p = 4.97 × 10−13) and granulocytes (p = 0.0001) and GA (p = 0.02) were associated with the GC score in the ELGAN ECHO cohort. Overall, while we found weak evidence that GC score decreases with increasing BPD severity, the strong correlations between cell composition and GC score appeared to explain much of this relationship.

Fig. 2: Forest plots of the associations between BPD severity and the polyepigenetic glucocorticoid (GC) score in ECHO NOVI and ELGAN ECHO.
figure 2

The established polyepigenetic GC score was not a significant measure of (a) risk in the ECHO NOVI cohort or (b) a predictor of BPD severity in the ELGAN ECHO cohort. No significant association was observed between GC score and BPD severity when adjusting for cell proportion, fetal growth restriction, maternal age, antenatal steroids, and/or steroids for BPD. However, GC score decreased as BPD severity increased, which is indicative of exposure to increased glucocorticoids (endogenous and exogeneous).

Full size image

HPA axis genes and BPD severity

Of the 706,323 CpGs captured in the ECHO NOVI cohort, 317 CpGs were annotated to the 8 HPA genes (Supplementary Table 1). BPD severity was not significantly associated with differential methylation of the HPA genes after accounting for multiple testing. However, we did observe nominally significant associations (raw p value < 0.05): mild BPD (11 CpGs), moderate BPD (27 CpGs), and severe BPD (2 CpGs) (Supplementary Table 2). Interestingly, we observed significant associations between antenatal steroid exposure and the methylation of 10 CpGs within 5 HPA genes (FKBP5, NR3C1, NR3C2, CRHR1, POMC, and HSP90AA1) (false discovery rate [FDR] <0.05) (Table 2). The most significant differentially methylated CpGs were cg06798479 (HSP90AA1; coef = 0.003 p value = 2.16 × 10−5), cg05857597 (NR3C1; coef = −0.003; p value = 7.8 × 10−5), and cg04281268 (HSP90AA1; coef = 0.001; p value = 8.59 × 10−5). One of the top-most significant CpGs associated with antenatal steroid exposure, cg20598211 annotated to NR3C1 (coef = −0.044; p value = 4.85 × 10−4), was also found to be nominally significant with mild BPD (coef = −0.025; p value = 0.019) and moderate BPD (coef = −0.038; p value = 0.005).

Table 2 CpGs significantly associated with antenatal steroid exposure in the ECHO NOVI cohort after adjustment for potential confoundersa.
Full size table

In the ELGAN cohort, 294 CpGs were annotated to HPA genes (Supplementary Table 3). As was found in the ECHO NOVI cohort, BPD severity was not significantly associated with the methylation of the 8 HPA genes at the FDR threshold, although several CpGs did yield raw p values < 0.05 (Supplementary Table 4). Similarly, antenatal steroid exposure was not significantly associated with HPA gene methylation, although 10 CpGs produced nominally significant raw p values < 0.05 (Supplementary Table 5).

Sex-specific HPA DNAm

Studies have shown that antenatal steroids affect males and females differently;61 therefore, we conducted a secondary analysis that was stratified by newborn sex. After sex-stratification, we observed sex-specific associations between antenatal steroid treatment and altered DNAm of HPA genes in the ECHO NOVI cohort. Males were observed to have differential methylation at 10 CpGs annotated to 4 HPA genes: CRHR1, HSP90AA1, NR3C1, and NR3C2 (Table 3). Additionally, 6 CpGs were uniquely altered in males, with cg09775582 (CRHR1; coef = 0.005; p value = 7.55 × 10−5) being the most significantly associated with antenatal steroid treatment. The other 4 CpGs were also observed in the unstratified model, with cg03396557 (CRHR1; coef = 0.006; p value = 6.28 × 10−5) and cg04281268 (HSP90AA1; coef = 0.002; p value = 1.31 × 10−4) being the most strongly associated with antenatal steroid treatment in males. We conducted a Spearman correlation analysis of the sex-specific regression coefficients and observed that males and females had a positive correlation for the differential methylation that was associated with antenatal steroid treatment (rho = 0.32; p value = 4.12 × 10−9). Similarly, after comparing males and females in the ECHO NOVI cohort, regression coefficients were correlated for moderate BPD (rho = 0.21; p value = 2 × 10−4) and severe BPD (rho = 0.125; p value = 0.025). These findings suggest that the methylomes of both males and females tend to be impacted similarly by antenatal steroids and moderate to severe BPD, but males may experience greater differential methylation at specific CpGs.

Table 3 Sex-specific associations between HPA axis CpGs and antenatal steroid exposure in ECHO NOVI males (N = 298) after adjustment for potential confoundersa.
Full size table

In the ELGAN ECHO cohort, we observed significant sex-specific associations with severe BPD (FDR < 5%) at two CpGs: cg26981809 (FKBP5, coef = −0.0200; p value = 2.49 × 10−6) and cg02757179 (POMC; coef = 0.0791; p value = 3.43 × 10−5) (Table 4). However, the Spearman correlation analysis did not identify significant correlations between the sex-specific regression coefficients of BPD severity or antenatal steroids between males and females in this cohort.

Table 4 Sex-specific associations between HPA-axis CpGs and antenatal steroid exposure in ELGAN ECHO males (N = 194) after adjustment for potential confoundersa.
Full size table

Discussion

Our study focused on BPD- and steroid-associated perturbations to the stress response system using an established polyepigenetic GC score and CpG-specific DNA methylation of 8 HPA genes. We examined whether the GC score was associated with BPD severity (in the ELGAN ECHO cohort), whether this score was responsive to BPD severity (in the ECHO NOVI cohort), and whether BPD and/or antenatal steroids were associated with the differential methylation of known HPA genes.

The polyepigenetic GC score was not associated with BPD severity, neither in the ELGAN ECHO cohort or the ECHO NOVI cohort. Tissue specificity and time of collection are likely crucial when utilizing the polyepigenetic GC score,32 as the differences between our samples and those used to train the GC score likely impacted its utility in the current study. The original GC score was calculated based on overlapping epigenetic responses to GC in human hippocampal progenitor cell lines and adult blood samples then it was tested in newborn cord blood, whereas our DNA samples were collected from buccal tissue and neonatal blood spot from infants born very and extremely preterm. Our minimally adjusted models appeared to reveal associations between increasing BPD severity and lower GC scores, which would have aligned with our expectations as higher GC exposure was associated with a lower GC score in the original investigations of the score. However, the adjustment for cellular proportions and age metrics substantially attenuated these initial findings. Thus, age and cell type may be strong confounders of these relationships. These two factors should be considered when utilizing this polyepigenetic GC score. Buccal tissue primarily consists of epithelial cells with very small proportions of immune cells and fibroblasts, whereas blood consists of several cell types, including nucleated red blood cells (nRBCs) and white blood cells (WBCs). Thus, our samples in the ECHO NOVI cohort (mostly epithelial) and in the ELGAN ECHO cohort (nRBCs and WBCs) potentially provide different methylome profiles from those described by Suarez et al.62 and Provençal et al.32 Lastly, the demographics of the participants in the Provençal et al. study were less racially diverse than our study, and thus the polyepigenetic GC score may not be as generalizable to our population.32

Several studies have found that adversity during pregnancy is associated with altered HPA function in offspring,63,64 and a combination of epidemiologic and experimental studies have recently shown that synthetic corticosteroids induce changes in DNAm33 and has lasting effects.32 Altered programming of the HPA axis has been implicated as a potential mechanism linking early life stress exposure to neurodevelopmental impairments. BPD is a substantial physiologic stressor in infants born very preterm, and its treatments could constitute additional stress. Contrary to our hypothesis, BPD severity did not show a strong impact on the methylome of 8 HPA genes. However, we did observe differential DNAm at two CpGs among males that developed severe BPD in the ELGAN cohort. As hypothesized, antenatal steroid exposure was associated with differential methylation at some CpG sites, but only in buccal tissue, suggesting that this association was tissue- or age-specific. These findings suggest that CpG methylation in HPA genes in some tissues are responsive to prenatal steroid exposure while other tissues are not. FKBP5, CRHR1, HSP90AA1, NR3C1, NR3C2, and POMC were the most significantly associated with antenatal steroid exposure in buccal tissue (FDR < 0.05). FKBP5 is a co-chaperone in the GC receptor complex and is often used as a proxy for HPA axis function.65 CRHR1 is the receptor of corticotropin-releasing hormone (CRH), which stimulates the secretion of adrenocorticotropin hormone (ACTH) when bound.66 HSP90AA1 is a molecular chaperone protein that regulates hormone signaling in response to stress.67 NR3C1 is a glucocorticoid receptor that when bound by cortisol negatively regulates CRH and ACTH production.68 NR3C2 is a mineralocorticoid receptor that also influences cortisol secretion.69 POMC precedes ACTH and endorphin production in the pituitary and is responsible for cortisol secretion from the adrenal cortex.70 Alterations to these genes have been associated with panic disorder development, depression, and anxiety.54

We also examined whether there were sex-specific associations between antenatal steroid exposure and HPA gene methylation, since prior work, including some from our group, has shown that males and females are impacted differently by early life exposures, including those introduced by the external environment, from maternal behavior, and from an adverse intrauterine environment.71,72,73,74,75,76,77 CpGs within CRHR1, HSP90AA1, NR3C1, and NR3C2 were significantly associated with antenatal steroid exposure in buccal tissue from the ECHO NOVI cohort; with 6 CpGs uniquely observed in the sex-stratified analysis for males only and 4 CpGs appearing in the full analysis. Thus, it is possible that some individual CpGs exhibit sex-specific responses to antenatal steroids, but overall, we observed relatively similar findings between males and females when we compared the regression coefficients from both models.

Although we did not observe significant association between BPD severity and the GC score or a strong significant association with DNAm of HPA genes, this does not rule out the possibility that BPD influences the methylome of other genes and biological systems. In previous work, our group observed that out of 4 neonatal morbidities (BPD, brain injury, retinopathy, and serious infection), BPD had the strongest impact on the neonatal methylome of preterm infants.23,24 Additionally, Martin et al. found that BPD was associated with more adverse neurodevelopmental and behavioral outcomes in NOVI infants at 2 years adjusted age.51 While this project specifically focused on potential impacts within the HPA system, it is certainly possible that other genes or biological pathways are differentially methylated among infants diagnosed with BPD. Additionally, observational studies78 and clinical trials,79 suggest that adrenal insufficiency, which might decrease cumulative secretion of adrenal glucocorticoids, is associated with a higher risk of BPD. It is plausible that the relationship between adrenal insufficiency and BPD may have confounded our analyses, if this condition limited endogenous GC secretion among our neonates with BPD. Such confounding would have biased our results toward the null. These limitations highlight the complexity of studying this medical condition as well as antenatal treatments, and all of its associated factors. Despite these limitations, our findings provide valuable insights into some epigenetic features that appear to not be responsive to BPD and its treatments, while some individual genes within the HPA axis warrant further research. Future studies should use untargeted epigenome-wide approaches to explore which genes and systems are most differentially methylated after BPD diagnosis and treatment, and to identify the impact of BPD severity on methylation.

While our findings are interesting and align with our hypothesis that steroid treatment and BPD may impact the methylation of some HPA axis genes, these findings should be considered along with the limitations of this study.80 First, a limitation of this study is that the steroid variables used in this analysis were limited to a binary “yes” or “no,” which does not tell us about timing, dose, or steroid type. Therefore, our findings are hypothesis-generating to further explore the impact of steroids on HPA methylation in cohorts with richer data on prenatal steroid exposure. Another limitation of this study is that there were differences in the timing and tissue type for DNA methylation samples. Thus, we could not perform a pooled analysis that would have improved the statistical power to detect small effects. Additionally, because we studied epigenetic associations in buccal cells, our study is unable to explore the DNAm landscape of the lung tissue that is directly impacted by BPD, and where the epigenetic responses and indicators may be the strongest. Overall, there was limited evidence of association between BPD severity and HPA axis methylation in infants born very preterm, and there is some evidence of HPA methylation differences with antenatal steroid exposure, but these findings require follow-up studies for confirmation and further exploration.

Conclusions

In this study, the GC score was not associated with BPD, and BPD severity did not strongly influence the methylome at CpGs within the stress response system. BPD and its treatments do not appear to be associated with the methylation of HPA genes in buccal tissue at NICU discharge; however, blood spot samples collected shortly after birth from males with severe BPD were found to have differential methylation within two HPA genes: FKBP5 and POMC. Antenatal steroids may have a lasting impact on the methylome of neonatal buccal cells at CpGs within FKBP5, CRHR1, HSP90AA1, NR3C1, NR3C2, and POMC. These findings were tissue-specific, as no differential methylation was observed in neonatal blood spots. Overall, the epigenome of the HPA system appears to be modestly impacted by antenatal steroid exposure and not BPD severity. Despite not observing a strong association between BPD severity and the methylation of HPA genes, methylation patterns of other biological pathways may still be impacted by BPD and its severity. Additionally, future studies should explore the effects of antenatal steroid exposure on the methylome in similar cohorts with richer prenatal and postnatal steroid exposure data and examine whether any identified epigenetic responses are risk factors for developmental outcomes.

Disclaimer

The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.