Differentially methylated CpGs in response to growth hormone administration in children with idiopathic short stature

Background Recombinant human growth hormone (rhGH) has shown a great growth-promoting potential in children with idiopathic short stature (ISS). However, the response to rhGH differs across individuals, largely due to genetic and epigenetic heterogeneity. Since epigenetic marks on the methylome can be dynamically influenced by GH, we performed a comprehensive pharmacoepigenomics analysis of DNA methylation changes associated with long-term rhGH administration in children with ISS. Results We measured DNA methylation profiles before and after GH treatment (with a duration of ~ 18 months in average) on 47 healthy children using customized methylC-seq capture sequencing. Their changes were compared and associated with changes in plasma IGF1 by adjusting sex, age, treatment duration and estimated blood proportions. We observed a considerable inter-individual heterogeneity of DNA methylation changes responding to GH treatment. We identified 267 response-associated differentially methylated cytosines (DMCs) that were enriched in promoter regions, CpG islands and blood cell-type-specific regulatory elements. Furthermore, the genes associated with these DMCs were enriched in the biology process of “cell development,” “neuron differentiation” and “developmental growth,” and in the TGF-beta signaling pathway, PPAR Alpha pathway, endoderm differentiation pathway, adipocytokine signaling pathway as well as PI3K-Akt signaling pathway, and cAMP signaling pathway. Conclusion Our study provides a first insight in DNA methylation changes associated with rhGH administration, which may help understand mechanisms of epigenetic regulation on GH-responsive genes. Supplementary Information The online version contains supplementary material available at 10.1186/s13148-022-01281-z.


Introduction
Growth hormone (GH) acts directly on intracellular pathways downstream of the GH receptor or via the stimulation and action of IGF1 [1]. Physiological effects of GH are mainly on growth, body composition and metabolism. In the epiphyseal growth plate, GH effects are largely mediated by promoting IGF1 production to stimulate skeletal growth. In addition, GH is known to increase muscle mass. GH also contributes to the acute metabolic response to stressful situations, such as fasting [2,3], exercise [4], injuries, critical Open Access *Correspondence: Xiaojian.shao@nrc-cnrc.gc.ca; tpastinen@cmh.edu; pierre. bougneres@inserm.fr illnesses [5] and infectious diseases [6] by acting on the liver, pancreatic beta cells, adipose tissue and muscle. Overall, most of these GH effects are mediated by activating the transcription of numerous genes.
Non-physiological exposure to GH occurs in the human clinic when children with short stature (SS) receive supra-physiological doses of recombinant human GH (rhGH) to stimulate their skeletal growth. These cases differ widely in the level and duration of exposure to rhGH which may leave lasting marks in different tissues that could conceivably be based on epigenetic mechanisms.
Changes in CpG methylation of the genome are one of the most studied epigenetic mechanisms. They can facilitate or reduce the transcription of certain genes when they occur in the regulatory zones of these genes. Once installed, some of them cannot be erased. However, they can show a certain plasticity and vary with age or environmental factors [7][8][9][10][11]. Most of them are variable depending on the tissue [12][13][14].
Modifications of CpG methylation have been associated with cancer, metabolic and cardiovascular diseases [15][16][17]. Following binding to its receptor, GH activates GHR-associated Src family kinases, acting via other intracellular pathways [18], such as ERK and Jun, which are known to affect CpG methylation [19].
It was therefore interesting to investigate whether GH was capable of producing its own epigenetic marks on the methylome. The only opportunity to detect them in the human clinic is to study the methylome before and after exposure to the hormone. This is made possible in one particular circumstance, the application of rhGH treatment to normal small children. Given the considerable individual variability of methylation marks, it would have been difficult to compare rhGH treated with untreated children. Therefore, we searched for acquired CpG modifications by studying the same children before and during rhGH exposure. As for most DNA methylation studies in humans, it was only possible to study blood cells in these children, with the expectation that changes in blood cells may reflect those occurring in the physiological target tissues of GH action. It is probably important to clarify that our study does not address at all the epigenetic variations that could contribute to the mechanisms of short stature, which have been the subject of a number of previous studies. Such methylation marks are not expected to vary upon GH administration and are therefore outside the scope of our approach devoted to the epigenetic consequences of exposure to GH action.

Cohort-sample collection
DNA samples were obtained from 47 children with nonpathological idiopathic short stature (ISS) recruited among 317 children who served as controls in a study of type 1 diabetes (T1D) epigenetics [20]. The 47 studied children with short stature all had ISS, as defined [21]. Briefly, this diagnosis requires that child's height is < − 2SD (standard deviation) from the population average, with a birth length > − 2SD for gestational age (none of them had intrauterine growth retardation) and appropriate for parents' height, then linear and normal growth rate, and no detectable etiology for the short stature, such as chromosomal, endocrine, or skeletal diseases. Classically, testing of such short children includes clinical examination, bone X-rays, IGF1measurement and/or GH stimulation tests (to exclude GH deficiency). The 47 studied children all went through these screening analyses, and none of them showed any abnormality in the cited parameters. This study was supported by the Programme Hospitalier de Recherche Clinique of the French Ministry of Health according to the French bioethics law with the objective of studying gene-environment factors in young patients with T1D and age-matched controls. Families were carefully informed about the investigational nature of the study and signed an informed consent agreed by CPP (number DC-2008-693; NI 2620, Comité de Protection des Personnes). The studied children had received rhGH treatment for a duration of 6-38.4 months.
None had deficiencies of GH or other hormones, chromosomal disorders or syndromes, skeletal dysplasia or metabolic disease. All children received rhGH as daily sc injections 6 days/wk, starting with a 40 µg/kg day rhGH dose and following a target-to-treat rhGH dosing protocol based on the growth response to treatment, so that individual average dose ranged 40-113 µg/kg day across the studied children. Children gave two blood samples, one before onset of rhGH treatment and the other after 6-38.4 months of rhGH treatment. All were seen as outpatients every 6 months for clinical examination and measurements of serum IGF-I. Children were healthy at time of study, with no sign of viral or other intercurrent infection. Blood samples were collected in the course of routine medical evaluation of patients. DNA methylation was measured in peripheral blood mononuclear cells (PBMC) samples from these patients before and during treatment.

Isolation of genomic DNA
Peripheral blood mononuclear cells (PBMC) were isolated from fresh blood using a density gradient. Five milliliters of fresh blood were mixed with 5 ml of NaCl 154 mM, and 5 ml of Lymphoprep solution (Eurobio, Paris, France) was added to the diluted blood and centrifuged for 20 min at room temperature at 800 g. After centrifugation, the interphase containing PBMC was carefully aspirated and the cells were mixed with NaCl. The cell suspension was centrifuged at 300 g and the cell pellet washed with PBS. PBMC were frozen at − 80 °C. Nucleic acids were extracted from PBMC using Gentra Puregene blood kit (Qiagen, Hilden, Germany).

Measurement of IGF1 concentrations
IGF1 concentrations were measured in serum samples around 07.00 am to 08.00 pm before breakfast using a chemiluminescent immunometric assay after pre-treatment with acid using Immulite ® 2000 (Siemens Healthcare Diagnostics Products Llanberis, UK). IGF1 values under treatment were averaged for analysis.

DNA methylation capture sequencing
Methylation capture sequencing (MCC-Seq) was performed as previously described [22][23][24]. In brief, the MCC-Seq protocol was carried out using the SeqCap Epi Enrichment System protocol (Roche NimbleGen) [22][23][24]. Specifically, a whole-genome sequencing library was prepared and bisulfite converted, amplified and then a capture enriching for targeted bisulfite-converted DNA fragments was carried out. Equal amounts of multiplexed libraries (12 samples per capture) were combined and were further amplified. Lastly, the MCC-Seq libraries were sequenced on the Illumina HiSeq4000 or NovaSeq 6000 system using 100 bp paired-end sequencing. More specifically, whole-genome sequencing libraries were generated from 700 to 1000 ng of genomic DNA spiked with 0.1% (w/w) unmethylated λ DNA (Promega) previously fragmented to 300-400 bp peak sizes using the Covaris focused-ultrasonicator E210. Fragment size was controlled on a Bioanalyzer DNA 1000 Chip (Agilent), and the KAPA High Throughput Library Preparation Kit (KAPA Biosystems) was applied. End repair of the generated dsDNA with 3′-or 5′-overhangs, adenylation of 3′-ends, adaptor ligation and clean-up steps were carried out as per KAPA Biosystems' recommendations. The cleaned-up ligation product was then analyzed on a Bioanalyzer High Sensitivity DNA Chip (Agilent) and quantified by PicoGreen (Life Technologies). Samples were then bisulfite converted using the Epitect Fast DNA Bisulfite Kit (Qiagen), according to the manufacturer's protocol. Bisulfite-converted DNA was quantified using OliGreen (Life Technologies) and, based on quantity, amplified by 9-12 cycles of PCR using the Kapa Hifi Uracil + DNA polymerase (KAPA Biosystems), according to the manufacturer's protocol. The amplified libraries were purified using Ampure Beads and validated on Bioanalyzer High Sensitivity DNA Chips, and quantified by PicoGreen.
The hybridization procedure of the amplified bisulfiteconverted library was performed as described by the manufacturer, using 1 μg of total input of library, which was evenly divided by the libraries to be multiplexed, and incubated at 47 °C for 72 h. Washing and recovering of the captured library, as well as PCR amplification and final purification, were carried out as recommended by the manufacturer. Quality, concentration and size distribution of the captured library were determined by Bioanalyzer High Sensitivity DNA Chips.
This blood MCC-Seq panel covers (1) the majority of human gene promoters, blood-cell-lineage-specific enhancer regions and methylation footprint regions [25] observed in blood, (2) CpGs from Illumina Human Methylation 450 Bead Chips and (3) published autoimmune-related SNPs as well as SNPs in their LD regions with r2 > 0.8. Overall, it covers 4,861,805 CpGs which have been applied to multiple blood-based epigenomewide association studies [24,26].

MCC-Seq data process
Targeted MCC-Seq HiSeq and NovaSeq reads were aligned using the Epigenome Pipeline available from the DRAGEN Bio-IT platform (Edico Genomics/Illumina). Specifically, the MCC-Seq paired-end raw reads were first demultiplexed into FASTQ files using Illumina's bcl2Fastq2-2.19.1 software. Reads were then trimmed for quality (phred33 ≥ 20) and Illumina adapters using trimgalore v.0.4.2 (https:// www. bioin forma tics. babra ham. ac. uk/ proje cts/ trim_ galore/), a wrapper tool around Cutadapt [27] and FastQC (https:// www. bioin forma tics. babra ham. ac. uk/ proje cts/ fastqc/). Then, the trimmed reads were aligned, per sequencing lane, to the bisulfiteconverted GRCh37 reference genome using DRAGEN EP v2.6.3 or later in paired-end mode using the directional/Lister methylation protocol presets; alignments were calculated for both strands, and the unique alignment with highest quality was retained. Lane bam files were merged and then de-duplicated using Picard (version 2.9). A genome-wide cytosine methylation report was generated by DRAGEN to record counts of methylated and unmethylated cytosines at each cytosine position in the genome. Methylation counts are provided for the CpG, CHG and CHH cytosine contexts. DNA methylation level of each CpG was calculated by the number of methylated reads over the total number of sequenced reads. CpGs that were located within sex chromosomes, overlapping with SNPs (dbSNP 137), the DAC Blacklisted Regions or Duke Excluded Regions (generated by the ENCODE project) were removed. CpG sites with less than 20X read coverage were also discarded.
Variants/SNPs (including homozygous alternate and heterozygous genotypes) were inferred using Bis-SNP (version 0.82.2) [28] on the de-duplicated bam files. The homozygous reference genotypes of individuals on these SNPs were extracted from the aligned bam files by requiring > = 10X read coverage aligned to the reference allele. Hierarchical clustering was performed based on the genotype profiles of SNPs on chromosome 1 where genotypes were inferred from all samples.

Statistical analysis
A linear regression model (LM) was built to investigate the association between the changes of DNA methylation level before and after treatment (delta beta) and the changes of IGF1 level (delta IGF1) by correcting age onset, the treatment duration and the estimated blood cell proportions. We used the R function lm() to fit the model and calculated p values for variables of interest. Due to limited sample size, the nominal p value < 1e−4 was used as a threshold of statistical significance to determine differentially methylated CpGs (DMCs). To reduce the impact of zero-inflated methylation differences on inferring the regression, we limited our analysis on CpGs that have > = 30% of pairs with nonzero methylation differences. Blood deconvolution was done using constrained linear projection [29] via the projectMix function of the RefFreeEWAS package, using a custom panel of 30,455 cell-type-specific hypo-methylated and hyper-methylated CpGs. The blood reference epigenome profiles include neutrophil, monocyte, B cell and T cell.

Genome features and function enrichment analysis
We downloaded the genome feature annotation tables, including transcription start sites (TSSs), 3'UTRs, 5'UTRs, first exons, exons, introns and transcription end sites (TESs), from the UCSC genome browser with the hg19 build version (https:// genome. ucsc. edu/). We considered both TSS200 (200 bp from TSSs) and TSS1500 (1500 bp from TSSs) for the promoter regions. We also downloaded the CpG islands (CGI) annotation table from the UCSC genome browser. Furthermore, CGI north and south shores were defined as the 2-kb flanking sequences on upstream and downstream of CGIs, respectively, and north and south shelves were defined as the 2-kb flanking sequences beyond the shores. Genome feature enrichment analyses of GH response DMCs were performed using Fisher's exact test for significance where the background set was the all testable CpGs. The gene ontology enrichment analyses were performed using homer [30] (version 4.11) with gene sets detected from the immune panel as the background set.

Descriptions of the GH cohort
Main characteristics of the 47 subjects (pairs) included in this study are shown in Table 1. All had ISS as defined by [31]. Briefly, ISS is a condition in which the height of the individual is more than 2 standard deviations (SD) below the corresponding mean height for a given age, sex and population, in whom no identifiable disorder is present.
The sex-matched participants (24 females vs. 23 males) were aged 4.9 to 13.1 years at time of starting rhGH treatment. They received a mean supra-physiological dose of rhGH of 70.8 μg/kg per day for a duration range from 6 to 38.4 months. During rhGH therapy, the mean increase in plasma IGF1 across children ranged from 4 to 658 ng/ml.

Differentially methylated CpGs in response to GH treatment
The average sequence genome coverage in targeted regions for these 94 samples was 25-fold (Additional file 1: Table S1). Close to 5.3 million CpGs (including those out of the targeted panels) were captured at autosomes with more than 20-fold read coverage and in at least one sample. When restricting attention to CpGs with good sample coverage in at least 10 before and after treatment pairs, 3,342,494 CpGs at autosomes remained for downstream analysis. In addition, the paired samples before and after GH treatment for 47 children were clustered well according to their genotype profiles which were inferred from the methylation sequencing data (Additional file 2: Figure S1).
We first compared the DNA methylation difference of each CpG per individual pair. For these ~ 3.3 M CpGs across 47 individual pairs, 68.1% of all comparisons were measured, see the sample (pair) coverage distribution of CpGs in Additional file 3: Figure S2A. A majority of the CpGs showed low mean absolute methylation difference (e.g., > 95% of these CpGs are showing mean methylation differences ≤ 5%) and low standard deviation (Fig. 1A, B). Specifically, it showed majority of them have DNA methylation difference of zero as indicated in the middle peak of Fig. 1A. Meanwhile, it also showed two lower sub-peaks which indicated the potential epigenetic responses to the GH treatment and represent the average methylation changes (i.e., ~ ± 5%). When considering each of the measured CpG pairs across 47 children (n = 107,019,585), the mean absolute value of the methylation changes between the baseline and rhGH stimulated samples is 4.8 ± 6.6% (SD). There is no difference between male and female groups (Additional file 3: Figure S2B).  Figure S2C demonstrates the methylation changes pattern for the top 5% most variable CpGs whose methylation profiles were measured for all individuals.
We observed that roughly 14.64%, 3.76% and 1.07% of the CpGs showed methylation difference > 10%, 20% and 30% comparing the baseline and after treatment samples, respectively (Fig. 1C). We then correlated the ratio of CpGs showing methylation difference > 10% of individuals with other phenotypes such as onsite age, difference of blood proportions and changes in plasma IGF1 concentration under treatment. We observed that the proportion changes of T cell and neutrophil are significantly corrected with the ratio of CpGs showing > = 10% differences (r = 0.39 and − 0.32, with p = 0.006 and 0.03, respectively) (Fig. 1D, E), but not for others parameters. This correlation trends are further increased when checking the ratio of CpGs showing > = 20% or 30% methylation level changes (Additional file 1: Table S2). Furthermore, we also observed that the proportion changes of B cell are showing significant correlation with the ratio CpGs showing 30% methylation level changes. Additional file 1: Table S2 presents the detailed correlations between CpG ratios at different levels of methylation changes and various phenotypic features. This implied that the methylation level changes might be confounded by the blood proportions changes. To identify the differentially methylated CpGs upon long-term exposure to rhGH, we fit a linear regression model on the difference between the DNA methylation values of all baseline and post-treatment PBMC samples and the difference in IGF1 values (Delta-IGF1) by adjusting the age at onset of treatment, treatment duration and differences in blood cell proportions before and after treatment. Figure 2A, B illustrates the QQ-plot and the Manhattan plot for this analysis. We did not observe any significant associations between Delta-IGF1 and the methylation level changes of CpGs (response DMCs) between baseline and post-treatment samples after multiple test corrections at a FDR of 0.05. However, we identified 2599 CpGs showing response DMCs (p < 1e−3) among baseline and post-treatment samples where 1317 response DMCs are negatively correlated and 1282 DMCs are positively correlated. With a more stringent p value threshold, we identified 267 DMCs at p value < 1e−4 with 123 negatively correlated and 144 positively correlated response DMCs. The DNA methylation change's pattern of the DMCs at p < 1e−3, which were also measured by all the individuals, is illustrated in Fig. 2C. Table 2 Fig. 2D-G, and Additional file 1: Table S3 lists the full list of these 267 response DMCs.

Genome feature and function enrichment analysis of response DMCs
We first performed a genomic feature enrichment analysis of response DMCs with the p value < 1e−4 and p value < 1e−3. We observed that the response DMCs with p value < 1e−3 were slightly enriched in the first exon, TES200 and CGIs (Fig. 3A) but no significant regions for the response DMCs with p value < 1e−4 (Fig. 3A). Furthermore, we also observed slightly enrichments on the blood-specific regulatory elements-DNase I hypersensitive sites (DHSs) regions of CD19 and CD20 for the response DMCs with p value < 1e−4 but not for the response DMCs with p value < 1e−3 (Fig. 3A).
We compared our identified DMCs with previously reported GH-related DMCs (n = 239, [32]). We did not find that any reported significant DMCs were replicated at our study, but when expanded to 1 kb distance-based neighboring CpGs, we observed a slight enrichment (fold-change of 1.73) of these reported DMCs at our response DMCs (p value < 0.001) compared with our non-significant CpGs (p value > 0.001) although this Table 2 The top 20 response DMCs list with p value < 1e−4. The response DMCs were sorted by the p value. CpG chromosome and position, regression p value, beta value (coefficient) and the annotated closest gene information (including genomic Annotation, Distance to TSS, Gene Name, Gene Type, and Gene Description of the closest gene) were provided enrichment was not significant (p value = 0.15) due to the small number of overlapped CpGs being counted. Interestingly, one response DMC (with p value = 1e−5) at chr11:68608981, which is located within the CTP1A gene, was within 200 bp to the reported GH DMC at chr11:68609166 (p value = 1e−5, [32]). It was reported that the effect of GH treatment on growth is potentially influenced by individual's SNPs [33]. We first linked our response DMCs to 37 potential SNPs which were reported to be GH-response-associated SNPs. We observed that only 14 SNPs have their surrounding CpGs in our tested CpG sets within 5 kb distance. Interestingly, three of these 14 potential SNPs were showing significant response DMCs (p value < 0.01) within a 5 kb distance to these SNPs. Particularly, one SNP (rs6600230, chr16:738477) is overlapping with the gene WDR24 where multiple CpGs were showing p value < 0.05 with the significant one (chr16:739598) located at the first exon region (p value = 0.01).

Discussion
The emerging field of pharmacoepigenomics will provide promising insights into the role drugs play in modulating the host epigenome and in addressing inter-individual variability in drug response and adverse effects. Although there is growing evidence that pharmacoepigenetics has the potential to become an important element of personalized medicine, we know of no study that has evaluated the changes in the individual methylome of the same group of patients undergoing a treatment, as performed in the current study. An additional advantage of our pharmacoepigenetic study is that clinical (growth rate, height) and biological outcomes (IGF1) can be quantified in response to precise rhGH dosing carefully injected by parents. More specifically, our data provide the first comprehensive pharmacoepigenomics analysis on rhGH treatment in children with ISS by comparing DNA methylation marks before and after several months of rhGH treatment. We identified 267 response DMCs which are associated with 265 genes and these genes were enriched in the biological process of cell differentiation, system development and different growth-related pathways such as endoderm differentiation, adipocytokine signaling, PPAR alpha and TGF-beta signaling pathways. This pilot study thus supports the existence of dynamic epigenetic changes in response to rhGH treatment. Again, it should be recalled that these are methylation changes induced by prolonged GH administration, and not epigenetic marks associated with short stature, an example of which can be found in one of our previous studies [34].
Our customized methylation sequencing panel captured more than 5 million CpGs, which is much larger than the previously used 450 K array data (i.e., > 10 folder larger), representing an unprecedented level of resolution. After quality control, more than 3.3 million CpGs remained for the response association analysis, providing the potential to discover novel signals. Indeed, of identified 267 response DMCs (with p value < 1e−4), only 114 (43%) of them are located within 100 bp distance to the known 450 K loci and only 25 (9.4%) of them are exactly located at the 450 K loci. Among the top 20 response DMCs listed in Table 2, half of them does not have neighboring CpGs (with 100 bp distance) in 450 K loci.
Previous studies had shown that the TGF-beta signaling pathway plays an important role in regulating osteoblast differentiation and could inhabit IGF-1/Akt signaling pathway [35]. The adipocytokine pathway and cAMP-signaling pathway are downstream signaling pathways upon activation of IGF1 receptor and contribute to the signal transduction of insulin-like growth factors on growth [36]. Another study reported that the GH modulates EGFR expression and signaling and further activates PI3K-Akt signaling, which was enriched in our response DMCs (p value < 1e−3) [37]. Moreover, our response DMC-associated genes were enriched in DNA-binding transcription factors as well as proteins with the common domain families of "WH/WH-like DNA" and "TIMP. " Particularly, TIMP3 is known to modulate GHR abundance and GH sensitivity [38], and NFKB1 is a known gene associated with short stature [39] and the growthpromoting effects of the transcription factor family of NFKB seems to be facilitated by GH and IGF-1 [40], while FOXA1, FOXN1 are regulators for GH activation [41]. Here, our identified genes with rhGH-associated methylation changes were enriched in these pathways, supporting the biological relevance of our findings. The genes involved in these pathways include CDKN2B, LEFTY2, PPP2R1B, CPT1A, RXRA, NFKB1, KCNMA1, BORCS8-MEF2B, MRVI1, PPIF and GATA4 (See the full response DMC list at p value < 1e−3, Additional file 1: Table S4).
The current study identified marked intra-individual responses of DNA-methylation to long-term rhGH treatment. A study by Kolarova et al. investigated 24 patients at baseline and after only 4 days of rhGH administration [32]. The studied patients had various forms of GH deficiency (N = 13) or other pathological conditions, which could influence the epigenetic responses to rhGH and complicate the interpretation. In comparison, only healthy children were selected for the current study and were either prepubertal or with minimal manifestations of puberty in order to avoid epigenetic changes that are associated in blood cells with advancing puberty [42].
Array-based DNA-methylation profiling of paired peripheral blood mononuclear cell samples in the Kolarova et al. 's study revealed clustering according to individuals rather than treatment [32]. Supervised analysis identified 239 CpGs as significantly differentially methylated between baseline and acutely GHstimulated samples, which nevertheless did not retain significance after adjustment for multivariate analysis. In a companion study, Kolarova et al. investigated the long-term effects of prolonged rhGH treatment on the DNA-methylome and analyzed peripheral blood cells from an independent cohort of 36 rhGH-treated children born small for gestational age (SGA) compared to 18 untreated controls. These were not paired samples which had to face major unwanted inter-individual variance of children methylome. No differentially methylated targets reached the level of significance in this long-term rhGH-treated cohort [32]. Our study did not replicate any of these 239 DMCs but observed a slight enrichment if considering significant response DMCs in 1 kb distance to them. The lack of high replicates might due to different etiologies of short stature (intrauterine growth retardation may influence epigenetic marks in the Kolarova et al. 's study) design, different treatment durations and dosing, different ages and more importantly the considerable inter-individual heterogeneity, while our study investigated paired intraindividual changes in methylation.
Of interest, MEX3C, the top gene in our response association analysis, was reported to be a translational regulator of IGF expression in mice [43]. IGF1 protein expression in bone cells was decreased upon MEX3C deficiency in Mex3c homozygous mutant mice. Given that MEX3C is highly conserved among mammalian species, the observation in mice might be relevant to the human IGF1 regulation and warrants further investigation. Among the top 20 signals, the response DMC at CPT1A (chr11:68608981) was close to the reported locus (cg20228509, chr11:68609166) within 1 kb in the Kolarova et al. 's short-term rhGH treatment study [32]. CPT1A was observed to be a genetic regulation of fatty acid metabolism, and missense mutation reduces height [44]. Although this evidence was not revealed in epigenetic studies, the potential pathway CPT1A involved (such as the adipocytokine signaling pathway, an important pathway related to IGF signaling) might indicate its indirect association with GH. In the same study, more than 3 CpGs in SLC15A4 were identified as differentially methylated loci and two of them were further validated with bisulfite pyrosequencing [32]. Our data showed the top signal was at the downstream of SLC15A4, and a few CpGs with nominal significance at p value < 0.01 were located in the intron of SLC15A4 gene. In addition, we also observed a couple of non-coding RNAs at our top signals list. Their functions related to GH are currently unknown and need further exploration.