Differential methylation pattern in pubertal girls associated with biochemical premature adrenarche

ABSTRACT

Biochemical premature adrenarche is defined by elevated serum DHEAS [≥40 μg/dL] before age 8 y in girls. This condition is receiving more attention due to its association with obesity, hyperinsulinemia, dyslipidemia, and polycystic ovary syndrome. Nevertheless, the link between early androgen excess and these risk factors remains unknown. Epigenetic modifications, and specifically DNA methylation, have been associated with the initiation and progression of numerous disorders, including obesity and insulin resistance. The aim of this study was to determine if prepubertal androgen exposure is associated with a different methylation profile in pubertal girls. Eighty-six healthy girls were studied. At age 7 y, anthropometric measurements were begun and DHEAS levels were determined. Girls were classified into Low DHEAS (LD) [<42 μg/dL] and High DHEAS (HD) [≥42 μg/dL] groups. At Tanner stages 2 and 4 a DNA methylation microarray was performed to identify differentially methylated CpG positions (DMPs) between HD and LD groups. We observed a differential methylation pattern between pubertal girls with and without biochemical PA. Moreover, a set of DNA methylation markers, selected by the LASSO method, successfully distinguished between HD and LD girls regardless of Tanner stage. Additionally, a subset of these markers were significantly associated with glucose-related measures such as insulin level, HOMA-IR, and glycaemia. This pilot study provides evidence consistent with the hypothesis that high DHEAS concentration, or its hormonally active metabolites, may induce a unique blood methylation signature in pubertal girls, and that this methylation pattern is associated with altered glucose metabolism.

KEYWORDS:

• Premature adrenarche
• DHEAS
• methylation
• glucose metabolism

Introduction

Adrenarche denotes the gradual growth and maturation of the innermost zone of the adrenal cortex (zona reticularis), leading to increased secretion of adrenal androgen precursors (AAP), including dehydroepiandrosterone (DHEA) and its sulphate ester (DHEAS) [1] [2]. Premature adrenarche (PA) can be defined biochemically by an increase of circulating AAP concentrations (above the age- and sex-specific reference range) [1,2]. Traditionally, prematurely elevated levels are indicated by serum DHEAS concentration ≥40 μg/dL before age 8 y in girls and 9 y in boys [2].

The incidence of PA is significantly higher in girls than in boys [3,4] and in children with higher BMI [1,2,4]. Premature adrenarche is receiving more attention due to its association with obesity and insulin resistance [1,2,5,6]; however, the mechanistic link between early androgen excess and these risk factors remains unknown.

Epigenetic modifications, such as DNA methylation, are important determinants of gene expression and have been associated with the pathogenesis of multiple diseases, including obesity and insulin resistance [7,8]. Moreover, epigenetic changes are themselves influenced by environmental and lifestyle factors, and do not remain stable throughout life [9]. Consequently, research on epigenetic modifications, their developmental changes, and their association with disease have become an important approach in understanding disease pathogenesis.

In this study, we investigated DNA methylation of peripheral blood cells during early and late puberty in girls who had exhibited prematurely elevated versus normal DHEAS concentrations prepubertally. Our observations identified methylation markers able to distinguish girls with and without a history of biochemically defined premature adrenarche. We also examined the association of these DNA methylation markers with glucose metabolic parameters.

Subjects, materials, and methods

Study population

Samples were originally collected to investigate the epigenetic age during pubertal development in girls. The study population for these girls has been described previously [10].

Our study population included a subset of 86 girls randomly selected from the longitudinal Growth and Obesity Cohort Study (GOCS) in Santiago, Chile, initiated in 2006, comprising children ages 2.6 to 4.0 y, recruited from public nursery schools of six counties in Santiago, who met the following inclusion criteria: 1) singletons born in 2002–2003 with birthweight between 2,500 and 4,500 g and 2) absence of physical (e.g., severe injury), medical (e.g., brain tumour), or endocrine diseases (e.g., diabetes, hypopituitarism, or thyroid disease) that could alter growth and/or onset of puberty.

All parents or guardians provided written informed consent, and the subjects gave their assent. The study protocol was approved by the Ethics Committee of the Institute of Nutrition and Food Technology, University of Chile (Approval number 028–2019).

Pubertal development

At age ~7 y, a single paediatric endocrinologist (VM) assessed breast development by inspection and palpation. The sexual maturity staging criteria and definitions were as described by Marshall and Tanner (Tanner staging) [11]. Thereafter, every 6 months, secondary sex characteristics were evaluated by a single female dietitian trained specifically for this purpose, with permanent supervision of a single paediatric endocrinologist (VM). Concordance for breast stages between dietitian and paediatric endocrinologist was 0.9 [12].

Anthropometric, hormonal, and biochemical measurements

Birth weight (kg) and length, as well as annual height and weight measurements prior to 2006, were obtained from medical charts. After study initiation, body fat percentage was estimated at each visit by Tanita-BC-418 MA bioelectrical impedance (Tanita-Corporation, Tokyo, Japan), according to the manufacturer’s guidelines and at measurement frequency of 50 kHz (accuracy ±0.1 kg) [13]. Weight and height were determined using standardized protocols (barefoot and light clothes). Sex- and age-adjusted height, weight, and BMI Z-scores at each visit were calculated based on World Health Organization (WHO) growth charts [14].

DHEAS concentration was determined in an early morning fasting blood sample at ~7 y of age and during pubertal progression at Tanner breast stages 2 (T2) and 4 (T4). Mothers were contacted the day before blood sampling to confirm the absence of fever (≥37.5°C) or symptoms of acute infection in the children. Samples were analysed at the Institute of Maternal and Child Research, University of Chile. Serum DHEAS at age ~7 y was determined as previously described by competitive specific binding RIA supplied by Diagnostic System Laboratories (Webster, TX, USA) [4]. Thereafter, DHEAS concentrations at T2 and T4 were analysed by liquid chromatography-mass spectrometry in a high-performance liquid chromatography (HPLC) Agilent system (Santa Clara, CA, USA) 1260 coupled to an AB Sciex 3200 Quantum ultratriple quadrupole mass spectrometer (Foster City, CA, USA).

Plasma levels of insulin (μU/mL) and glucose (mg/dL) were measured by electrochemiluminescence and colorimetric methods, respectively, at the Nutrition Laboratory of the Catholic University of Chile [15]. Insulin sensitivity was estimated from fasting insulin and glucose levels using the Homoeostasis model (HOMA-IR), which has been validated in non-diabetic children and adolescents (fasting insulin [μU/mL] × fasting glucose [mg/dL]/405).

DNA methylation measurements and preprocessing

Peripheral blood buffy coat DNA methylation (DNAm) at T2 and T4 was measured in bisulphite-treated DNA using the Illumina Infinium MethylationEPIC array at the USC Epigenome Center (Los Angeles, CA).

Array data analysis was performed using R version 4.0.2. The Minfi package v1.34.0 [16] was used for quality control checks at the sample level. Data were normalized using Noob for background correction and dye bias normalization from Minfi [16] and BMIQ (beta mixture quantile normalization) from WateRmelon package v1.32.0 [17]. We filtered out probes containing detection p-values >0.01 in one or more samples, on Y sex chromosome, cross-hybridizing probes, non-cg probes, and probes that lie near single nucleotide polymorphisms [18]. A total of 705,792 probes were available for further analysis.

A beta value was calculated for each locus as the ratio of methylated signal intensity to the sum of unmethylated and methylated signal intensities, with a range from 0 to 1 (proportion methylated). The beta values were log transformed into the corresponding M values, which were used for all analyses. Beta values only were reported in figures and tables for easy interpretation [19].

Estimated cell proportions (CD8+ T cells, CD4+ T cells, natural killer cells, B cells, monocytes, and granulocytes) were calculated using the Houseman approach [20] implemented in Minfi package. CpG sites were assigned to chromosomes, locations, and genes using the Illumina manifest for the EPIC array (hg19 genome reference).

Multidimensional scaling (MDS) was performed to examine sources of variation using the MDS plot function in the Minfi package. Since DNA methylation data were from two different plates, the M values were corrected for the effects of plate batch using removeBatchEffect function from Limma package [21].

Statistical analyses

A mixed model design was used to identify differentially methylated CpG positions (DMPs) between high and low DHEAS groups in paired samples from Tanner stages 2 and 4, adjusted for cell subtype composition (model 1). Age at DHEAS measurement was also included as a covariate in the model (model 2).

For both models, we applied a robust regression to down-weight the influence of outliers, in conjunction with an empirical Bayes method to moderate standard errors, both implemented in the R Limma package. To account for the Tanner 2 and Tanner 4 measurements from the same girl, the intra-sample correlation was first calculated by the duplicate Correlation function in the Limma package. We corrected for multiple testing using the False Discovery Rate (adjusted p-value <0.05) according to the Benjamini and Hochberg method [22].

To identify the strongest DNA methylation markers, a Least Absolute Shrinkage and Selection Operator (LASSO) regression was performed. LASSO regression, which involves a shrinkage process to identify variables strongly associated with variables corresponding to the target, adds a penalty equal to the absolute value of the magnitude of coefficients. Coefficients that become zero are eliminated from the model, which allows identification of CpGs that were significantly associated with antecedent premature adrenarche. The optimal parameter ‘lambda,’ obtained using 10-fold cross-validation, is a penalty parameter; cross validation indicates the variables and coefficients from the best model.

The CpG site selection was performed on the training set, considering the previously identified DMPs (FDR < 0.05) and anthropometric, hormonal, and biochemical variables, using the glmnet package v4.0–2. The training set of 64 paired samples (128 totals) was randomly selected, including 13 HD girls and 51 LD girls. Remaining paired samples (n = 22) were used as a test set. The performance of the 31-CpG methylation biomarkers was tested by sensitivity and specificity on the test set; the sensitivity was calculated as the number of HD girls that were correctly classified, divided by all HD girls, and the specificity was calculated as the number of LD correctly classified, divided by all LD individuals. The classification model, including 31 discriminating CpG sites, was assessed graphically using Heatmap based on the Manhattan metric. Finally, a classification prediction model was constructed using a logistic regression method.

The p-value distribution of all the CpG sites was compared to the theoretical distribution using a quantile-quantile plot and the genomic inflation factor lambda. The genomic inflation factor lambda was calculated by dividing the median of observed chi-square statistics by the median of theoretical chi-square statistics with 1 degree of freedom [23], which quantified the systematic bias of test statistic.

Correlations between differentially methylated CpGs and glucose metabolic variables at T1, T2, and T4 were tested using Spearman coefficients as implemented in R. Quantitative variable comparisons between the LD and HD groups were determined with the nonparametric Wilcoxon’s test.

Results

Cohort description

The DNA Methylation Microarray Data included 86 girls. To evaluate the relationship between serum DHEAS and DNA methylation changes through puberty, two blood samples were taken from each participant, one sample at T2 and the other at T4 (172 samples, Figure 1).

Figure 1. Flow chart of the study design. LD= Low DHEAS, HD= High DHEAS.

At age of recruitment, girls from the LD group were younger than those of the HD group (p-value = 0.043), and we therefore adjusted for age at recruitment in the Methylation Analysis. The groups did not differ in BMI SD or fat percentage at birth or at age of classification (Suplementary Table S1). The clinical data of the 86 girls at the pubertal stage assessments are shown in Table 1. Notably, the DHEAS difference at the Tanner stage 1 grouping into high and low DHEAS persisted throughout pubertal development (Tanner stages 2 and 4, Supplementary Figure S1).

Table 1. Clinical and anthropometric characteristics of 86 girls by DHEAS levels at Tanner 2 and 4.

	Tanner 2											Tanner 4
	Low DHEAS					High DHEAS						Low DHEAS					High DHEAS
	n	Median	IQR			n	Median	IQR			p-value	n	Median	IQR			n	Median	IQR			p-value
Age at DHEAS measurement, years	68	9.81	9.3	-	10.47	18	9.76	9.31	-	10.14	0.370	58	11.31	10.72	-	11.83	16	11.13	10.75	-	11.63	0.524
BMI-for-age Z-score	68	0.7	−0.38	-	1.62	18	0.7	0.19	-	1.54	0.730	58	0.71	−0.18	-	1.5	16	0.84	0.22	-	1.19	0.486
Waist circumference	68	65.45	59.05	-	71.65	18	63.62	60.12	-	68.01	0.903	58	69.03	63.4	-	73.99	16	67.6	64.84	-	75.16	0.670
Body fat percentage, %	68	25.55	22	-	29.85	18	25.35	22.98	-	27.27	0.844	58	25.4	22.15	-	27.73	16	25.6	22.62	-	27.82	0.679
DHEAS, μg/dL	66	58.45	41.67	-	71.83	18	90.45	87.33	-	116.85	<0.001	58	71.5	54.27	-	87.28	16	116.1	100.2	-	142.98	<0.001
Testosterone, ng/ml	49	0.06	0.05	-	0.09	14	0.07	0.06	-	0.09	0.467	58	0.17	0.09	-	0.23	15	0.15	0.11	-	0.25	0.806
Estradiol, pg/ml	68	14	9.75	-	20	18	14	11.25	-	18.5	1.000	52	28.5	15.75	-	41	12	36.5	12.25	-	50.75	0.816
Androstenedione, ng/ml	49	0.27	0.19	-	0.35	14	0.33	0.24	-	0.39	0.186	58	0.7	0.52	-	0.96	15	0.72	0.55	-	0.99	0.530
Fasting blood glucose, mg/dl	67	90	84	-	94	18	92	89	-	94	0.057	55	87	82.5	-	91.5	16	84.5	81	-	90	0.331
HOMA-IR	67	1.6	1.4	-	2.35	18	2	1.7	-	2.5	0.125	55	2	1.7	-	2.8	16	1.85	1.38	-	3.25	0.757
Insulin, mU/mL	67	7.3	6.42	-	10.4	18	8.75	7.9	-	10.12	0.197	55	9.7	7.85	-	13.1	16	8.9	7.05	-	15.65	0.978

Differences between DHEAS groups were analysed using Wilcoxon test. Values are presented as medians. IQR = interquartile ranges; LD = Low DHEAS; HD = High DHEAS.

Identifying differentially methylated CpGs at breast Tanner stage 2 and 4

To identify DNA methylation markers associated with DHEAS concentration that had a consistent profile across puberty, we performed a differential methylation analysis between LD and HD group paired samples. The study flowchart is presented in Figure 1. In model 1, we detected a total of 69 DMPs (λ inflation factor = 1.12) able to discriminate HD girls from LD girls in both Tanner stages (Figure 2a), of which 35 were up-methylated and 34 were down-methylated. Regarding the location of DMPs in relation to genes (Figure 2b), 33% were in gene body and 26% were in promoter regions (TSS200, TSS1500, and 5‘UTR). The percentage at promoter regions was higher for up-methylated sites (31%) than down-methylated sites (21%). DMP lists are available in Supplementary Table S2.

Figure 2. (a) Volcano plot showing differences in methylation between LD and HD groups (in 86 paired samples). Significant (FDR = 5%) DMPs are shown in red (down-methylated) and blue (up-methylated). The y axis shows –log10 P values and the x axis plots the mean difference in methylation M values. The horizontal line denotes significant CpGs at an FDR of 5%. (b) Location of DMPs relative to genomic location: 1st exon, 5’UTR (5’ untranslated region), gene body, TSS1500 (200–1500 bp upstream of the transcription start site), TSS200 (0–200 bp upstream of the transcription start site) and open sea (DNA methylation sites outside the CpG island regions). (c) Binomial Deviance of the LASSO model with different log, identification of the most appropriate log λ value (penalization coefficient) in the LASSO model, with 10-fold cross validation the optimal lambda was 0.00015. Numbers along the top x-axis represent the average number of predictors. (d) Heat map of unsupervised hierarchical clustering of the 31 sites of the 86 paired samples. Each column represents a girl DNA sample; each row represents a gCpg site; red and blue colors indicates DNA methylation levels as M value. The horizontal distances on each branch of the dendrogram represent the degree of similarity between DNA methylation profiles of the two DHEAS groups. logFC= Log fold change at base 2.

This list was further extended to 178 DMPs in model 2, when age at T1 was included as a covariate in the model (Supplementary Table S3), however the inflation factor was higher than in model 1 (λ inflation factor = 1.16) (Supplementary Fig S2A).

CpG candidate selection

The inflation factor lambda implied slight inflation of p-values. To improve the selection of CpGs and identify potential markers of high DHEAS, we performed a LASSO regression analysis. First, the paired samples were split randomly into training and testing groups at a 3:1 ratio. After splitting, we obtained a training group containing 64 samples and a testing group containing 22 samples. Then, we used the LASSO method in the trained set to further decrease the number of markers, considering the 69 DMPs resulting from model 1 (Figure 2c). CpGs with non-zero coefficients in the LASSO regression model were selected for further analysis. After LASSO, we obtained 31 methylation sites that could discriminate between LD and HD girls (Supplementary Table S4). Unsupervised hierarchical clustering of these 31 markers confirmed discrimination between LD and HD girls by the selected DMPs (Figure 2d).

To test these 31 sites as potential markers, we calculated the classification accuracy of the proposed model. As part of analysing the model performance, 75% of the data was taken for training, and 25% was taken for testing. Among the test group of 22 girls at Tanner 2 and 4 (44 samples), the classification sensitivity and specificity were 100% and 100% in HD girls, respectively. The confusion matrix generated out of the predicted variables is shown in Table 2. Next, we tested the performance of these 31 markers in the 86 girls (test and training groups, 172 samples), and the classification sensitivity and specificity reached 100% and 100%, respectively (Table 3); showing that these 31 DNA methylation markers could significantly distinguish T2 or T4 girls with a history of HD at age 7.

Table 2. Confusion table of the test group and performance of predictive CpG sets.

		Test group
		LD	HD
Group predicted	LD	34	0
	HD	0	10
Sensitivity		1
Specificity		1

Table 3. Confusion table of the 86 paired samples and performance of predictive CpG sets.

		Test group
		LD	HD
Group predicted	LD	136	0
	HD	0	36
Sensitivity		1
Specificity		1

We also performed this approach with 178 DMPs from model 2 and obtained 41 markers. However, using this subset of probes the samples did not correctly classify the respective groups by unsupervised hierarchical clustering as with model 1. Furthermore, the performance was inferior to model 1, yielding 80% sensitivity and 100% specificity in the test samples (Supplementary Figure S2). Therefore, model 2 results were not considered in further analyses.

Association between DMPs and glucose metabolic variables

Of the 31 DNA methylation markers, 16 were located within genes and 14 were close to genes (188 bp−145 Kbp). A literature search indicated that eight of these genes are associated with hyperandrogenism conditions or related to glucose metabolism (Table 4). The remaining DNA methylation markers, without mentioned literature link, were excluded from subsequent analysis.

Table 4. CpG sites related by literature review related to hormonal or metabolic parameters.

Gene	CpG site	chr	Position	Literature
DOCK1	cg19572051	chr10	128810484	Polycystic ovary syndrome [24]
SLC46A1	cg22436928	chr17	26733488	Cholesterol homeostasis [25]
NOM1	cg04969340	chr7	156740196	Obesity [26]
CACNA1H	cg09367545	chr16	1233640	Premature adrenarche [27–29]
CTGF	cg24684168	chr6	132271588	Polycystic ovary syndrome [30] Diabetes [22,31]
CPPED1	cg09756671	chr16	12968251	Glucose metabolism [32]
SORCS2	cg16283698	chr4	7610756	Diabetes [33]
PCDHB15	cg13337064	chr5	140625389	Diabetes biomarker [34]

Reference genome: hg19/GRCh37

Based on literature review indicating linkage with PA, we selected 8 of the 31 DNA methylation markers for exploratory analysis of whether methylation at these sites showed associations with glucose metabolic variables. We observed that methylation at four of the CpG sites at T2 and/or T4 in our cohort was significantly associated with glucose parameters (insulin, HOMA, or glycaemia, measured at age 7 y, T2, and/or T4) (Figure 3a, b, c, d). No association remained significant after Benjamini-Hochberg correction for multiple testing; all correlations for each time point are shown in Supplementary Table 5. To compare the methylation level at these four CpG sites between LD and HD groups, the nonparametric Wilcoxon test of methylation at these four CpG sites was conducted (Figure 4). Except for the cg16283698 site at Tanner 2, the other sites showed significant differences between the LD and HD groups at both T2 and T4.

Figure 3. A scatter plot of the Spearman correlation between methylation beta value at (a). cg13337064 (PCDHB15) at T2 and HOMA or insulin at T2. (b). cg24684168 (CTGF) at T2 and glycemia at T2. (c). cg16283698 (SORCS2) at T4 and HOMA at T1. d. cg09756671 (CPPED1) at T4 and glycemia at T4.

Figure 4. Boxplot comparison of methylation levels between DHEAS groups at Tanner stages 2 and 4. (a) cg13337064 (PCDHB15) (b) cg24684168 (CTGF) (c) cg16283698 (SORCS2) (d) cg09756671 (CPPED1). Differences were tested using Wilcoxon test.

CpG biomarkers and anthropometrics variables

To reduce the numbers of DNA methylation markers, we tested if anthropometric variables could contribute to the model. We repeated the LASSO analysis considering the 31 previously identified markers and birth weight, age, BMI-for-age, and body fat percentage at T1, T2, and T4. This reduced the number of CpG markers from 31 to 23 CpGs in the full model and the model including birth weight.

Discussion

This exploratory study showed a differential methylation pattern between pubertal girls with and without biochemical PA, and this effect was consistent across the early and late pubertal time-points examined. Indeed, a set of DNA methylation markers, selected by the LASSO method, successfully distinguished between HD and LD girls regardless of Tanner stage.

Altered DNA methylation patterns have previously been reported in animal models after exogenous DHEAS administration [35,36], and observations in humans also suggest that androgens may induce epigenetic modification [37–39]. Interestingly, DHEA treatment of addicts has been shown to reverse their altered DNA methylation pattern on peripheral white blood cells [40]. These observations suggest an impact of the androgen precursors DHEA and DHEAS, or of their active androgen and/or oestrogen metabolites, on DNA methylation. Further evidence suggests that androgens may induce epigenetic modifications by affecting the expression of DNA methyltransferase (DNMT) or by histone modifications [38,41]. Based on these earlier findings and our own observations, we hypothesize that biochemical PA may have induced the differential DNA methylation seen in the current study. However, our observational evidence of an association between biochemical PA and DNA methylation does not, by itself, allow assessment of causality.

We further explored potential associations between the differentially methylated CpG sites in our study and genes linked to disorders or markers of abnormal glucose metabolism. Literature review of genes associated with the 31 differentially methylated sites indicated eight such genes (Table 4). Among the exploratory association analyses, four CpGs (from the eight genes) showed a significant correlation between measures of glucose metabolism in our cohort and the level of DNA methylation.

Additionally, these four CpG sites showed significant differences between LD and HD groups at T2 and/or T4. Although there is overlap between the groups, the magnitude of the difference in methylation levels between groups can vary according to the tissue in which it was measured. This was observed by Crujeiras and co-workers, who found that the methylation map of circulating leukocytes mirrors the subcutaneous adipose DNA methylation profile between obese and non-obese individuals [42]. Therefore, the difference found in blood cells could be even higher in other important tissues.

The CpG site cg13337064, for which methylation status at Tanner 2 correlated with HOMA-IR and insulin level at Tanner 1 (Figure 3 A), is in exon 1 of gene PCDHB15 (Gene ID: 56121). It has recently been shown that PCDHB15 can be epigenetically silenced via promoter methylation [43].

The CpG site cg24684168, for which methylation status at Tanner 2 correlated inversely with glycaemia at Tanner 1 (Figure 3B), is located at exon 3 of gene CTGF (Gene ID: 100194970). Several studies have shown that DNA methylation regulates CTGF expression [44–46].

The CpG site cg16283698, for which methylation status at Tanner 4 correlated inversely with HOMA-IR at Tanner 1 (Figure 3 C), is located at intron 3 of gene SORCS2 (Gene ID: 57537). A specific methylation pattern of this gene was observed in children at risk for obesity associated with maternal BMI [33].

The CpG site cg09756671, for which methylation status at Tanner 4 correlated with glycaemia at Tanner 1 (Figure 3 D), is located upstream of gene CPPED1 (Gene ID: 553139). Downregulation of CPPED1 expression increased glucose uptake in adipose tissue [32].

Premature adrenarche in girls is associated with abnormalities of glucose metabolism, including hyperinsulinism in pre-puberty and puberty [5,6], increased prevalence of hyperglycaemia during puberty [47], and higher HOMA-IR and prevalence of insulin resistance in young adulthood [48], but the causality for these associations is controversial. Increased BMI has been hypothesized to cause both insulin resistance and premature adrenarche [6,48]. Interestingly, when we include girls’ anthropometric variables with the DNA methylation markers, the LASSO method selected birth weight as the only informative anthropometric variable. Studies in two large UK cohorts [49,50] also showed that low birth weight was associated with higher adrenal androgen levels at 8 (347 girls) and 9 (101 girls) y old, respectively [51]. Moreover, those children with rapid postnatal catch-up weight gain exhibited higher adrenal androgen levels. To these observations, we add the association with differential methylation of CpG sites in girls with biochemical PA. Because DNA methylation patterns have highly variable spatio-temporal dynamics, longitudinal evaluation at multiple time points would be required to determine the temporal relationship between the observed epigenetic pattern and elevated DHEAS levels. Such studies may indicate whether the observed epigenetic patterns are a likely cause or a consequence of elevated DHEAS levels and the associated abnormalities of glucose metabolism.

Limitations of this study include its relatively small sample size, exploratory nature, and limited sampling points, and the need for validation in an external, independent, and larger cohort. However, the study benefits from a prospective and unique paediatric cohort in which data on a large array of variables have been acquired.

In summary, data from this pilot study are consistent with the hypothesis that DHEAS concentration may induce a unique blood methylation signature in girls with premature adrenarche – and that this signature remains consistent from Tanner stage 2 to 4. Furthermore, the observed associations with measures of glucose metabolism suggest that DNA methylation markers may be a causal link between premature adrenarche and its reported association with insulin resistance. Further studies will be necessary to confirm these hypotheses and their potential pathophysiologic significance.

Acknowledgments

We thank all the children participating in the GOCS cohort. The genotyping of GOCS was performed at the Spanish National Cancer Research Centre, in the Human Genotyping lab, a member of CeGen, PRB3, and is supported by grant PT17/0019, of the PE I+D+i 2013-2016, funded by ISCIII and ERDF.

Disclosure statement

The authors declare that there are no actual or perceived conflicts of interest that could prejudice the impartiality of the research reported or any financial or other potential conflict of interest.

Data availability statement

The data that support the findings of this study are available from the corresponding author (VM) upon reasonable request.https://repositorio.uchile.cl/handle/2250/183290

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15592294.2023.2200366.

Additional information

Funding

This work was supported by the National Research and Development Agency through the Regular Fondecyt Project N° 1140447 and 1190346 to Dr. Verónica Mericq; the Regular Fondecyt Project N° 1120326 to Dr. Camila Corvalán; and the Regular Fondecyt Project N° 1150416 to José L. Santos and the R01 grant (R01CA158313) to Karin B. Michels.