Skip to main content

Identification of differential genomic DNA Methylation in the hypothalamus of pubertal rat using reduced representation Bisulfite sequencing

Abstract

Background

There are many variables affecting the onset of puberty in animals, including genetic, nutritional, and environmental factors. Recent studies suggest that epigenetic regulation, especially DNA methylation, plays a majorrole in the regulation of puberty. However, there have been no reports on DNA methylation of the pubertal genome.

Methods

We investigated DNA methylation in the female rat hypothalamus at prepuberty and puberty using reduced representation bisulfite sequencing technology. The identified genes and signaling pathways exhibiting changes to DNA methylation in pubertal rats were determined by Gene Ontogeny and Kyoto Encyclopedia of Genes and Genomes analysis.

Results

The distribution of the three types of methylated C bases in promoter and CpG island (CGI) regions in the hypothalamus was as follows: 87.79% CG, 3.05% CHG, 9.16% CHH for promoters, and 88.35% CG, 3.21% CHG, 88.35% CHH for CGI in prepubertal rats; and 90.78% CG, 2.13% CHG, 7.09% CHH for promoters, and 88.59% CG, 88.59% CHG, 8.35% CHH for CGI in pubertal animals. CG showed the highest percentage of methylation, and was the highest methylation state in CGI. Compared to prepubertal hyoyhalamus samples, we identified ten genes with altered methylation in promoter regions in the pubertal hypothalamus samples, and 43 genes with altered methylation in the CGI. Changes in DNA methylation were found in gonadotropin-releasing hormone signaling pathways, and the oocyte meiosis pathway.

Conclusion

Our results demonstrate changes in DNA methylation occur in female rats from prepuberty to puberty suggestng DNA methylation may play a crucial role in the regulation of puberty onset. This study provides essential information for future studies on the role of epigenetics in the regulation of puberty.

Background

Puberty is defined as the time of transition from childhood to adulthood [1]. Physiologically, puberty is the point at which an animal is first capable of reproducing sexually; in the female, this is the first occurrence of menstruation and in the male the testes begin to undergo spermatogenesis. In humans, early and rapidly progressive puberty is harmful [2,3,4,5]. In contrast, domesticated animals are bred for early puberty [6, 7]. In animal production, breeding of early puberty livestock can save feeding costs and improve the utilization of animal, while it can also shorten the generation interval of fine animals and accelerate the genetic breeding process [8, 9].

Pubertal timing is influenced by complex interactions involving genetic [10], nutritional [11], environmental [12], and (in humans) socioeconomic factors [13]. The traditional view is that hypothalamic gonadotropin-releasing hormone (GnRH) neurons are sensitive to negative feedback by the gonadal hormones before puberty. Therefore, due to inhibition by gonadal hormones, a low amplitude and frequency, pulsed secretion of GnRH is observed from birth to puberty onset. The onset of puberty begins with the release of this inhibitory mechanism in the hypothalamus, and re-activation of a high amplitude and frequency GnRH secretion pulse. In turn, GnRH acts on GnRH receptors in anterior pituitary follicle stimulating hormone (FSH) and luteinizing hormone (LH) secretory cells, thereby promoting secretion of LH and FSH. In turn, LH and FSH act on the gonads, stimulating the secretion of testosterone in males or estradiol secretion in females, thereby promoting the rapid development of reproductive organs and sexual characteristics.

A specific switch may exist to control hypothalamic GnRH secretion, and induce puberty onset; however, its precise mechanism remains to be determined. Evidence suggests the initiation of puberty requires the coordinated activity of genes organized into functional networks [14,15,16]. A large number of genes involved in neuroendocine regulation control of the onset of puberty [17, 18]. Elks et al. reported that more than 30 gene variants were associated with age at menarche in humans [19], supporting the claim to that no one gene controls the onset of puberty [20]. At present, important genes in the hypothalamus related to the onset of puberty include the Kiss-1/GPR54 system [21, 22], NPY [23, 24], Leptin [25, 26], LIN28B/let-7 [27] and NKB [28,29,30]. Polygenic inheritance is the basis of the onset of puberty [20, 31, 32], and is closely related to hypothalamic correlation factors, transcription and changes in the expression levels of receptor genes. However, how these puberty-related genes are activated, and how the dynamic regulation of gene expression is coordinated during the onset of puberty, remains unclear.

Epigenetic regulation of gene expression is now recognized as an important player in the complexity of DNA and histone modification patterns, leading to huge differences in cell phenotype and function [32]. Epigenetic regulation is sensitive to stimuli, and cells are able to adapt to their environment through these modifications. However, the modification of these complex patterns has not been fully explained. The neuroendocrine system is known to have a high degree of sensitivity to environmental changes, so in recent years, epigenetic control of neuroendocrine function has become a focus of much research. It has been confirmed that DNA methylation plays an important role in the puberty onset in rats [32,33,34]. In particular, hypothalamic expression of two key genes (Eed and Cbx7) decreased and methylation of their promoters increased before puberty; inhibiting DNA methylation blocked both events and resulted in pubertal failure [32].

Despite these data, little information is available on DNA methylation profiling of the rat pubertal genome. In this research, we performed reduced representation bisulfite sequencing (RRBS), a cost-efficient method for genome-scale DNA methylation analysis [35]. The aim of the current study was to investigate the DNA methylation profile of the hypothalamic genome in the pubertal rat at a single base resolution.

Methods

Animals

For the first 19 days after pairing, litters were assessed daily, and the day of birth was considered postnatal day 1. The animals were weaned on day 21. The rats were kept under standard conditions (12:12 h light-dark cycle with lights on between 06:00 and 18:00 h; temperature, 22 ± 1 °C; rat chow and water provided ad libitum; cage size: L × W × H,46.5 × 34.5 × 20 cm). Prepubertal samples were collected at postnatal day 25 and pubertal samples were collected on the day of vaginal opening (postnatal day 36.2 ± 0.58) at the stage of estrus.

Sample preparation and genomic DNA isolation

Female rats were sacrificed using cervical dislocation. The hypothalamus was dissected using the tuber cinereum as the ventral landmark for cuts to remove the frontal lobe and lateral and posterior portions of the brain. The cortex was peeled away from the remaining ventral brain piece containing the hypothalamus [36, 37]. The position of the hypothalamus was indicated by HE staining (Fig. 1a,b). Samples were frozen in liquid nitrogen immediately and stored at −80 °C until use.

Fig. 1
figure 1

The location of the hypothalamus in the brain.Representative horizontal (a) and vertical (b) images of histological section of a rat brain. The black circle indicates the position of the hypothalamus

The hypothalami from five rats each in the prepubertal of pubertal groups were mixed together and genomic DNA was extracted. Genomic DNA was isolated using AxyPrep™ Multisource Genomic DNA Miniprep Kit (Corning, AP-MN-MS-GDNA-50G), according to the manufacturer’s protocol. The quality of the genomic DNA was checked by the NANODROP 2000 UV DNA Analyzer of which the minimum detection sample volume is 0.5 μl and the wavelength range is 190-840 nm. (Thermo scientific).

Library construction and sequencing

The library construction and sequencing was performed at BGI-Shenzhen (Shenzhen, China). Briefly, at least 3 μg of total genomic DNA from each sample was digested with MspI (NEB, R0106T). The digestion was followed by DNA end-repair, the addition of a single A nucleotide, and multiplexed adapter ligation. The DNA samples were purified by agarose gel electrophoresis to isolate 40–220 bp long fragments. We then used the ZYMO EZ DNA Methylation-Gold™ kit (ZYMO RESEARCH, D5006) to get bisulfite conversion, and amplified the products by PCR. The quality of library was analyzed by the Agilent 2100 Bioanalyzer (Agilent Technologies), of which the DNA size analysis range is 25–12,000 bp and quantified by the ABI StepOnePlus Real time-PCR System (Thermo Fisher). The resulting fragments were sequenced on an Illumina Hiseq2000 analyzer (Illumina). To determine the location of each clean read before doing methylation analysis, we used Bsmap 2.74 to compare with the reference genome (R. norvegicus, assembly Rno-6.0, http://www.ncbi.nlm.nih.gov/assembly/GCA_0000018 95.4).

RRBS data analysis

RRBS utilizes the restriction enzyme Mspl to digest the entire genome, and allows for enrichment of sequences with relatively high CpG content (i.e., promoter and CpG islands [CGI] regions). These fragments are sequenced after bisulfite treatment (to convert unmethylated cytosine), and analysis of the methylation status of each base site is performed [38, 39]. Although only restriction fragments are sequenced, they mainly cover CpG-rich regions, making it possible to judge the state of methylation of the whole genome from RRBS results.

We analyzed the numbers of promoters and CGI (CG, CHG, CHH), cytosine site coverage analysis, the cumulative distribution of cytosine effective sequencing depth and compared the theoretical coverage and the actual depth of coverage on CG sites. Methylation analysis of promoter and CGI regions was done through the distribution ratio of CG, CHG and CHH in methylated C bases, the average methylation level of C bases and the distribution of different types of methylation levels. Regional methylation profiling was calculated by the association between methylation level and CpG density in a specific area, and the CpG density was defined for each CpG site within a window of 200 bp. Analysis of differentially methylated regions (DMRs) using sliding window method was performed to calculate the DMRs between prepubertal and pubertal samples.

Functional annotation of DMRs

DMRs were used for functional annotation, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. GO analysis is the international standard for gene function classification, which fully describes the properties of genes and gene products in organisms [40]. In brief, all the identified genes are mapped to each term in GO database, gene numbers for each term are calculated, and then a hyper geometric test is used to find significantly enriched GO terms in DMRs compared to the reference gene background and obtained a p-value. The Bonferroni Correction for the p-value was used to obtain a corrected p-value. Differential expression is determined and the corrected p-value (<0.05) is the threshold for significant differences in gene expression within each GO term, and allows for identification of differentially expressed genes in pathways of major biological functions. The gene sequences were then compared to the NCBI non-redundant protein database by using BLAST, and the results were annotated to GO using Blast2 GO.

KEGG is the major public databases analysis of pathway [41] that can be used to better understand the biological function of genes. Significant enrichment of pathways is determined by hypergeometric test identifying the significantly altered pathways in differentially expressed genes after comparing to the entire genome. The Q Value is set at <0.05 as the threshold for significance, and gene sequences were annotated to KEGG databases by BLAST. Rich Factor analysis was then used to determine the ratio of the number of the genes located in pathways of differentially expressed genes and in all annotated genes; the greater the Rich Factor, the greater the degree of enrichment.

Isolation of total RNA and quantitative real-time PCR

Total RNA was isolated using a Transzol Up Plus RNA Kit (Transgen Biotech, Beijing, China) following the manufacturer’s instructions. RNA degradation and contamination was monitored on 1% agaroes gels. RNA purity was detected using the NanoDrop 2000c spectrophotometer (Thermo Scientific, IL, USA). Reverse transcription was carried by HiScript® QRT SuperMix for qPCR(+gDNA wiper)(R123–1, Vazyme, Nanjing, China). Real-time PCR was performed using a StepOnePlus Real-time PCR Instrument (Applied Biosystems) in a 15 reaction system as described in the instructions of AceQTM qPCR SYBR® Green Master Mix (Q121–02, Vazyme, Nanjing, China). Reaction were incubated in a 96-well optical plate (Applied Biosystems) at 95 °C for 10 min, followed by 40 cycles of 95 °C for 30s, 60°Cfor 30s, and 72 °C for 30s. Each sample was run in triplicate. Primer sequences were synthesized by Sangon Biotech (Shanghai, China) on the basis of the gene sequences obtained from the NCBI database. Primer sequences are lised in Table 1. Gene expression was quantified relative to β-actin expression using the 2-ΔΔCT method.

Table 1 Primers used in real-time PCR

Results

Coverage analysis of promoter and CGI regions

After library construction, alignments were compared to the reference genome sequences. We had 248,727,316 and 244,131,868 clean reads from the prepubertal and pubertal samples, respectively (Additional file 1: Table S1). The enzyme rate which is the the unique reads containing enzyme loci percentage of total unique reads was 77.19% and 75.64%, and the bisulfite conversion rate was 99.32% and 99.39%, respectively. Our RRBS covered most promoters, although not all of the genome; we list the information covered, including the total number covered in the genome, the theoretical value in enzyme fragments and the actual value covered by at least 5 CGs in Additional file 1: Table S2. The total number of promoters covered was 18,845, the theoretical value (target region) was 174, and the actual coverage value was 66 in prepubertal and 64 in pubertal samples (Additional file 1: Table S2).

For CGI, the theoretical value indicates location of cytosine in the theoretical enzyme cutting region, and the experimental value is the actual number of cytosines covered by our sequencing reads. We investigated the numbers of each cytosine type identified (CG, CHG and CHH; H represents non-G base, hereinafter inclusive). In the promoter regions, the actual coverage rate of C sites (CG, CHG, CHH), was 35.69% CG, 18.61% CHG and 13.56% CHH in prepubertal samples, but was 33.66% CG, 16.96% CHG and 12.14% CHH in pubertal samples (Table 2). In CGI in prepubertal samples, the actual distribution of C sites (CG, CHG, CHH) was 52.17% CG, 53.76% CHG and 50.17% of CHH; in pubertal samples, the distribution was 51.36% CG, 52.92% CHG and 49.32% CHH (Table 2). The number of C sites in the genome, the number of theoretical fragments of C sites, and the number of the actual coverage sites are shown in Table 2.

Table 2 Cytosine distribution

Methylation analysis of promoter and CGI regions

We calculated the ratio of the three types of methylated cytosine (mCG, mCHG, mCHH) in promoters and CGI. We found that compared with the prepubertal samples, in promoter regions from pubertal samples, the proportion of mCG increased 2.99%, while mCHG and mCHH decreased 0.92% and 2.07% (Fig. 2a and c). The same result was seen in CGI, with the proportion of mCG increased 0.24%, and mCHG and mCHH decreased 0.16% and 0.19% (Fig. 2b and d).

Fig. 2
figure 2

Different methylation at C bases in the hypothalamus of prepubertal and pubertal rats. Methylation at CG, CHG, and CHH in promoter regions in (a) prepubertal samples and (b) pubertal samples, as well as in CGI in (c) prepubertal samples and (d) pubertal samples. mCG, mCHG, and mCHH represent the methylated forms of CG, CHG and CHH

We then analyzed the average methylation level of all C bases. CG had the highest average level of methylation in both prepubertal and pubertal samples, whether in the promoter or CGI. Investigating further, we found that mCHH had the highest average level of methylation in promoters in both groups. However, in CGI, mCHH had the highest average level of methylation in the prepubertal samples, while mCG as the highest in the pubertal samples group (Table 3). Finally, we analyzed the percentage of methylation of methylated cytosines in each sequence context. In promoter regions, mCG was distributed evenly in pubertal samples. The distribution of mCHH was almost zero in low methylation state (<20% of the distribution), but was more than 40% of mCHH in high methylation states (>90%). mCHG were mainly distributed at approximately 60% and approximately 100% high methylation levels in prepubertal samples, and in pubertal samples were mainly distributed in 10%, 50% and >90% of the three methylation status (Fig. 3a and c). The distribution of methylation level in CGI was different from the distribution in promoters: mCG, mCHG, mCHH methylation levels were highly distributed (>90%), and especially for mCHH, nearly 45% were hypermethylated (Fig. 3b and d).

Table 3 Average methylation levels of C and mC
Fig. 3
figure 3

Distribution of mC bases in the hypothalamus of prepubertal and pubertal rats. Distribution of methylation at mCG, mCHG, and mCHH in promoter regions in (a) prepubertal samples and (b) pubertal samples, as well as in CGI in (c) prepubertal samples and (d) pubertal samples. The x-axis indicates methylation level, and the y-axis shows the percentage of all mC at a certain level of methylation. mCG, mCHG, and mCHH are represented by the red, blue, and green lines, respectively

Regional methylation profiling

We used heat mapping to represent methylation distribution, CpG density distribution and the relationship between methylation and density. In prepubertal samples, there were 108,083 restriction fragments, 10 fragments of promoters, and 37,003 fragment of CGI were identified in this analysis (Fig. 4a). In all enzyme fragments, both in prepubertal and pubertal samples, methylation levels were dichotomized: high methylation levels with low CpG density and low methylation levels with high CpG density. CpG density was generally low in promoters and high in CGIs; both regions showed low levels of methylation (Fig. 4a and b). Hypermethylated and hypomethylated C in enzyme fragments indicated a higher proportion than others while the majority proportion C was in low methylation levels in promoter and CGI. The heat map of CHG and CHH showed similar results.

Fig. 4
figure 4

Heat maps showing methylation distribution and CpG density distribution. Heat maps of (a) Prepubertal and (b) pubertal samples. n refers to the number of analyzed CpGs (per-strand depth ≥ 4) within that feature. CpG density (x-axis) is defined as the number of CpG dinucleotides within the 200 bp window. Methylation level (y-axis) is defined as average methylation level of cytosines in CpGs. The thin black lines within each heat map denote the median methylation level of CpGs at the given local density. The red gradient indicates abundance of CpGs that fall into bins of given methylation levels and CpG densities. The blue bar charts above each heat map show the distribution of CpG densities, projected onto the x-axis of the heat maps. The green bar charts to the right of the heat maps show the distribution of methylation levels, projected onto the y-axis of the heat maps

Identification of differentially methylated regions and bioinformatics analysis

We identified DMRs between prepubertal and pubertal samples, as shown in Table 4. The 23 chromosomes investigated showed a total of 152 DMRs. The identity of the differentially methylated genes are shown in Table 5; if the log2 ratio value is positive, that indicates upregulation, while a negative log2 ratio value indicates downregulation. In promoter regions, there were eight differentially methylated genes that were upregulated, and two were downregulated. In CGI regions, 42 differentially methylated genes were upregulated, and 10 were downregulated.

Table 4 Methylation regional analysis
Table 5 Gene with Differentiation Methylated Regions across puberty

Statistical significance of GO and KEGG enrichment analysis to identify differentially methylated genes was evaluated by the hypergeometric tes. GO functional analysis indicated that differentially methylated genes in promoter regions were enriched in 106 major functional groups. In promoter regions, six differentially methylated genes were found involved in the 14 cellular component groups (Additional file 1: Table S3), six differentially methylated gene were involved in 13 molecular function groups (Additional file 1: Table S4), and seven differentially methylated genes were involved in 79 biological process groups (Additional file 1: Table S5).

The numbers of differentially methylated genes were identified in the functional groups of single-organism process, biological regulation, cellular process in biological process and binding functional groups in molecular function (Fig. 5a). In CGI, differentially methylated genes were enriched in 348 functional groups; there were 31 differentially methylated genes involved in the 47 cellular component groups (Additional file 1: Table S6), 28 differentially methylated genes involved in 63 molecular function groups (Additional file 1: Table S7), and 28 differentially methylated genes involved in 238 biological processes (Additional file 1: Table S8). Cellular process, single-organism process, cell, cell part and binding functional groups occupied the dominant positions (Fig. 5b). This evidence suggests that differentially methylated genes are involved in a variety of biological processes.

Fig. 5
figure 5

The Gene Ontology functional analysis of prepubertal and pubertal rat hypothalamus. Gene Ontology (GO) analysis of differentially methylated genes between the hypothalamic genome of prepubertal and pubertal rats in (a) the promoter and (b) the CGI regions. The y-axes on the left represents the percentage of genes and the one on the right is the number of genes involved in different functional groups. The x-axis is the name of the functional groups within the three GO terms, cellular component, molecular function and biological process

Using KEGG Pathway analysis, we found that differentially methylated genes were involved in 25 (Additional file 1: Table S9) and 57(Additional file 1: Table S10) pathways, respectively in promoter and CGI. The pathways with the largest number of differentially methylated genes in the promoter region were the melanogenesis pathways, and the pathway with the largest Rich Factor value was the vasopressin-regulated water reabsorption pathway (Fig. 5a). Rich factor relates to the ratio between the number of genes enriched in pathway and annotated in DMR. Rich factor and enrichment are positively correlated. The GnRH signaling pathway also had a high degree of enrichment (Fig. 6a). The pathway with the largest number of differentially methylated genes in CGI was the vasopressin-regulated water reabsorption pathway, and the highest Rich Factor value pathway was the glycosaminoglycan biosynthesis-keratin sulfate (Fig. 6b). This analysis identified the GnRH signaling pathway and the oocyte meiosis pathway (Fig 5b). Thus, differentially methylated genes were shown to be extensively involved in many aspects of cell metabolism. Interestingly, the expression of Gnas and Itpr2 whcih were enriched in the GnRH signaling pathways and the oocyte meiosis pathways were detected and showed no variation from prepuberty to puberty (Fig. 7).

Fig. 6
figure 6

Top 20 enriched pathways in the hypothalamic genome of prepubertal and pubertal rats. KEGG pathway analysis showing the top 20 enriched pathways in the (a) promoter and (b) CGI regions of differentially methylated genes from the prepubertal and pubertal rats. Rich factor relates to the ratio between the number of genes enriched in the pathway and annotated in differentially methylated regions (DMRs). Rich factor and enrichment are positively correlated. The Q value is the P value that was corrected after multiple hypothesis testing. Q value and enrichment are in a negatively correlated

Fig. 7
figure 7

Expression of the Itpr2 and Gnas in prepubertal and pubertal samples by quantitative PCR

Discussion

In this study, we found methylation predominated at CG sites, with differential methylation levels across the CGI versus promoter sites. We also found an increase in the proportion of methylated Cs in CG sites (over CHG and CHH sites) in promoter regions in pubertal hypothalamus samples, and a reduction in the methylation levels at mCHG sites in promoter regions. Genes that were differentially regulated during puberty were involved in the vasopressin-regulated water reabsorption and GnRH signaling pathways. The results presented here provide the initial targets for future studies of epigenetic regulation in the onset of puberty, and enhance our understanding of the epigenetic mechanisms involved.

Our data showed that the highest proportion of methylation occurs on CG, which conforms to the known characteristics of mammal methylation. There are differences in DNA methylation between different species, and in different tissues from the same species. For example, transposon elements and promoter methylation level is higher than in the gene and the coding sequence [39, 42]. In plants, Ziller et al. confirmed that DNA methylation states differ in different plant organs and different developmental stages of the same organ [43]. Different types of C site methylation (mCG, mCHG, mCHH) occur in the genome in different proportions; the ratio of each type of mC site in all the mC sites reflects the characteristics of methylation in a particular species and period to a certain extent.

In this study, we examined the proportion of the three types of C bases (mCG, mCHG, mCHH) in the promoter and CGI regions. Comparing pubertal to prepubertal samples, the ratio of CG increased 2.99%, while CHG and CHH decreased 0.29% and 2.07% respectively, in promoter regions. The ratio of CG increase 0.24% and the ratio of CHG and CHH decrease 0.09% and 0.16% respectively in CGI regions. Recent Recent studies have shown non-CpG methylation (mCHG, mCHH) plays a role in transcriptional repression [44, 45]. In this study, the ratio of CG increased, while CHG and CHH decreased in both the promoter regions and the CGI regions. Those changes may be related to the initiation of transcription of some genes during puberty onset, although whether these changes can be used to monitor puberty onset requires further research.

In addition, this study determined the levels of each type of mC base. We found that in the promoter regions, the methylation level was distributed more evenly; however, high methylation levels (> 90%) were more apparent in the CGI regions (about 25%). In a 2008 study, Cokus et al. found that methylation levels polarize; if an mC is not in a low methylation state, it is in a high methylation state [46]. We consider such a difference is mainly due to different species; our study used rats, and different species have different characteristics of DNA methylation.

The effects of DNA methylation on the onset of puberty have been reported in previous studies. In particular, the hypothalamic expression of two key Polycomb group (PcG) genes, Eed and Cbx7, decreased, and methylation of their promoters increased before puberty [32]. Before the onset of puberty, Eed represses of Kiss1 expression (which encodes the protein kisspeptin) by binding to the Kiss1 promoter. Epigenetic silencing of Eed leads to its dissociation from the Kiss1 promoter, resulting in an increase in Kiss1 mRNA expression at the puberty onset. Indeed, the typicap increase in pulsatile GnRH release that occurs at puberty onset was disrupted when Eed remained bound to the Kiss1 promoter, perhaps due to inhibition of DNA methylation [32]. Unfortunately, we did not find methylation changes in these two genes in our study. However, undetected promoter methylation changes may have occurred due to the absence of restriction enzyme sites used in the RRBS technology.

Pervious studies have also shown an abundance of H3K4me3 and H3K9ac14 at the Kiss1 promoter in prepuberty is associated with the loss of PcG inhibition [32, 47]. In particular, the TrxG (Trithorax group) complex, which antagonizes PcG proteins, catalyzes H3K4me3 and facilitates H3 acetylation. Some evidence indicates the TrxG complex is involved in the control of puberty, for example, CHD7 mutations result in hypothalamic hypogonadism in humans [48]. In particular, the TrxG complex may activate puberty-associated genes by post-translationally modifying histones during the onset of puberty [47]. However, we found no link between the TrxG complex and DNA methylation in our study. The involvement of post-translational modification of histones in the control of puberty onset is a future direction of our research.

We applied GO and KEGG pathway analysis of the DMRs in prepuberty and puberty. We found the vasopressin-regulated water reabsorption pathway was the most enriched in these pathways. Vasopressin is released due to high salt could,and could activate vasopressin receptors expressed on neurons that release kisspeptin in the anteroventral periventricular nucleus. This could stimulated GnRH and LH secretion, which affect puberty onset [49,50,51,52,53]. Regardless of whether the methylated regions were in the promoter or CGI, they are related to the GnRH signaling pathway (Gnas and Itpr2), and the key role of GnRH in puberty is well established [54]. In CGI, the methylated regions were linked to the oocyte meiosis signaling pathway (Itpr2).Although the methylation of Gnas and Itpr2 were different, the expression level showed no significance from prepuberty to puberty, indicating the correclation between methylation levels yaries, which is consistent with previous studies [55,56,57]. Indeed, the role of Gnas and Itpr2 in the onset of puberty requires further study.

Except for Gnas and Itpr2, we found that more than one differentially methylated gene was involved in a wide range of signaling pathways, including cellular processes, biological processes, biological regulation, cells and cell components; these are all closely related to growth and development. These results also indicate that the onset of puberty involves multiple signaling pathways simultaneously, and these signaling pathways could have crosstalk interactions. This further confirms the complexity and multifaceted nature of the regulatory mechanisms of puberty. Exactly which gene silencing or activation pathway plays a regulatory role in puberty onset, and their regulatory mechanisms, requires further screening and verification.

The puberty onset of the male is different from the female, and only female rats were examined in the manuscript. The differences between males and females will be solved in a future study. In addition, there was DNA methylation change in prepubertal and pubertal female samples. However, our data is not sufficient to prove that DNA methylation regulation onset of puberty, because so many factors can influence the epigenetic change, such as estrogen [58]. We will be exploring what specific factors cause this methylation changes and whether this changes can be used as a marker to indicate the puberty onset in the future research.

Conclusions

Our results demonstrate changes in DNA methylation occur from prepubertal to puberty in rats, suggesting that DNA methylation may play a crucial role in the regulation of puberty onset. The study provides essential information for future studies of the role of epigenetics in puberty.

Abbreviations

CGI:

Cytosine-phosphate-guanine island

CpG:

Cytosine-phosphate-guanine

DMRs:

Differentially methylated regions

DNA:

Deoxynucleic acid

GnRH:

Gonadotropin-releasing hormone

GO:

Gene ontology

KEGG:

Kyoto encyclopedia of genes and genomes

Kiss-1/GPR54:

Gene encoding kisspeptin/G-protein receptor

NKB:

Neurokinin B

NPY:

Neuropeptide Y

RRBS:

Reduced Representation Bisulfite Sequencing

References

  1. Banerjee I, Clayton P. The genetic basis for the timing of human puberty. J Neuroendocrinol. 2007;19:831–8.

    Article  CAS  PubMed  Google Scholar 

  2. Kaltiala-Heino R, Rimpela M, Rissanen A, Rantanen P. Early puberty and early sexual activity are associated with bulimic-type eating pathology in middle adolescence. J Adolesc Health. 2001;28:346–52.

    Article  CAS  PubMed  Google Scholar 

  3. Kaltiala-Heino R, Kosunen E, Rimpela M. Pubertal timing, sexual behaviour and self-reported depression in middle adolescence. J Adolesc. 2003;26:531–45.

    Article  PubMed  Google Scholar 

  4. Golub MS, Collman GW, Foster PM, Kimmel CA, Rajpert-De Meyts E, Reiter EO, Sharpe RM, Skakkebaek NE, Toppari J. Public health implications of altered puberty timing. Pediatrics. 2008;121(Suppl 3):S218–30.

    Article  PubMed  Google Scholar 

  5. Elks CE, Ong KK, Scott RA, van der Schouw YT, Brand JS, Wark PA, Amiano P, Balkau B, Barricarte A, Boeing H, et al. Age at Menarche and Type 2 Diabetes Risk The EPIC-InterAct study. Diabetes Care. 2013;36:3526–34.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Patterson JL, Beltranena E, Foxcroft GR. The effect of gilt age at first estrus and breeding on third estrus on sow body weight changes and long-term reproductive performance. J Anim Sci. 2010;88:2500–13.

    Article  CAS  PubMed  Google Scholar 

  7. Rosales Nieto CA, Ferguson MB, Macleay CA, Briegel JR, Martin GB, Thompson AN. Selection for superior growth advances the onset of puberty and increases reproductive performance in ewe lambs. Animal. 2013;7:990–7.

    Article  CAS  PubMed  Google Scholar 

  8. Saito H, Sasaki Y, Koketsu Y. Associations between Age of Gilts at First Mating and Lifetime Performance or Culling Risk in Commercial Herds. J Vet Med Sci. 2011;73:555–9.

    Article  PubMed  Google Scholar 

  9. Nieto CAR, Thompson AN, Macleay CA, Briegel JR, Hedger MP, Ferguson MB, Martin GB. Relationships among body composition, circulating concentrations of leptin and follistatin, and the onset of puberty and fertility in young female sheep. Anim Reprod Sci. 2014;151:148–56.

    Article  Google Scholar 

  10. Macedo DB, Brito VN, Latronico AC. New Causes of Central Precocious Puberty: The Role of Genetic Factors. Neuroendocrinology. 2014;100:1–8.

    Article  CAS  Google Scholar 

  11. Soliman A, De Sanctis V, Elalaily R. Nutrition and pubertal development. Indian J Endocrinol Metab. 2014;18:S39–47.

    PubMed  PubMed Central  Google Scholar 

  12. Soliman A, De Sanctis V, Elalaily R, Bedair S. Advances in pubertal growth and factors influencing it: Can we increase pubertal growth? Indian J Endocrinol Metab. 2014;18:S53–62.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Lomniczi A, Wright H, Ojeda SR. Epigenetic regulation of female puberty. Front Neuroendocrinol. 2015;36:90–107.

    Article  CAS  PubMed  Google Scholar 

  14. Krewson TD, Supelak PJ, Hill AE, Singer JB, Lander ES, Nadeau JH, Palmert MR. Chromosomes 6 and 13 harbor genes that regulate pubertal timing in mouse chromosome substitution strains. Endocrinology. 2004;145:4447–51.

    Article  CAS  PubMed  Google Scholar 

  15. Ojeda SR, Lomniczi A, Mastronardi C, Heger S, Roth C, Parent AS, Matagne V, Mungenast AE. Minireview: the neuroendocrine regulation of puberty: is the time ripe for a systems biology approach? Endocrinology. 2006;147:1166–74.

    Article  CAS  PubMed  Google Scholar 

  16. Ojeda SR, Lomniczi A, Loche A, Matagne V, Kaidar G, Sandau US, Dissen GA. The transcriptional control of female puberty. Brain Res. 2010;1364:164–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Teles MG, Bianco SD, Brito VN, Trarbach EB, Kuohung W, Xu S, Seminara SB, Mendonca BB, Kaiser UB, Latronico AC. A GPR54-activating mutation in a patient with central precocious puberty. N Engl J Med. 2008;358:709–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Freitas KC, Ryan G, Brito VN, Tao YX, Costa EM, Mendonca BB, Segaloff D, Latronico AC. Molecular analysis of the neuropeptide Y1 receptor gene in human idiopathic gonadotropin-dependent precocious puberty and isolated hypogonadotropic hypogonadism. Fertil Steril. 2007;87:627–34.

    Article  CAS  PubMed  Google Scholar 

  19. Elks CE, Ong KK. Whole genome associated studies for age at menarche. Brief Funct Genomics. 2011;10:91–7.

    Article  CAS  PubMed  Google Scholar 

  20. Ojeda SR, Dubay C, Lomniczi A, Kaidar G, Matagne V, Sandau US, Dissen GA. Gene networks and the neuroendocrine regulation of puberty. Mol Cell Endocrinol. 2010;324:3–11.

    Article  CAS  PubMed  Google Scholar 

  21. d'Anglemont de Tassigny X, Colledge WH. The role of kisspeptin signaling in reproduction. Physiology (Bethesda). 2010;25:207–17.

    Article  Google Scholar 

  22. Navarro VM. New insights into the control of pulsatile GnRH release: the role of Kiss1/neurokinin B neurons. Front Endocrinol. 2012;3:48.

    Article  Google Scholar 

  23. Tillet Y, Picard S, Bruneau G, Ciofi P, Wankowska M, Wojcik-Gladysz A, Polkowska J. Hypothalamic arcuate neuropeptide Y-neurons decrease periventricular somatostatin-neuronal activity before puberty in the female lamb: morphological arguments. J Chem Neuroanat. 2010;40:265–71.

    Article  CAS  PubMed  Google Scholar 

  24. Kamegai J, Minami S, Sugihara H, Hasegawa O, Higuchi H, Wakabayashi I. Growth hormone receptor gene is expressed in neuropeptide Y neurons in hypothalamic arcuate nucleus of rats. Endocrinology. 1996;137:2109–12.

    Article  CAS  PubMed  Google Scholar 

  25. Atala A. Re: leptin signaling in GABA neurons, but not glutamate neurons, is required for reproductive function. J Urol. 2014;191:1928.

    Article  CAS  PubMed  Google Scholar 

  26. Hausman GJ, Barb CR, Lents CA. Leptin and reproductive function. Biochimie. 2012;94:2075–81.

    Article  CAS  PubMed  Google Scholar 

  27. Park SW, Lee ST, Sohn YB, Cho SY, Kim SH, Kim SJ, Kim CH, Ko AR, Paik KH, Kim JW, Jin DK. LIN28B polymorphisms are associated with central precocious puberty and early puberty in girls. Korean J Pediatr. 2012;55:388–92.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Topaloglu AK, Reimann F, Guclu M, Yalin AS, Kotan LD, Porter KM, Serin A, Mungan NO, Cook JR, Ozbek MN, et al. TAC3 and TACR3 mutations in familial hypogonadotropic hypogonadism reveal a key role for Neurokinin B in the central control of reproduction. Nat Genet. 2009;41:354–8.

    Article  CAS  PubMed  Google Scholar 

  29. Ramaswamy S, Seminara SB, Ali B, Ciofi P, Amin NA, Plant TM. Neurokinin B stimulates GnRH release in the male monkey (Macaca mulatta) and is colocalized with kisspeptin in the arcuate nucleus. Endocrinology. 2010;151:4494–503.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Tusset C, Noel SD, Trarbach EB, Silveira LF, Jorge AA, Brito VN, Cukier P, Seminara SB, Mendonca BB, Kaiser UB, Latronico AC. Mutational analysis of TAC3 and TACR3 genes in patients with idiopathic central pubertal disorders. Arq Bras Endocrinol Metabol. 2012;56:646–52.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Ojeda SR, Lomniczi A. Puberty in 2013: Unravelling the mystery of puberty. Nat Rev Endocrinol. 2014;10:67–9.

    Article  PubMed  Google Scholar 

  32. Lomniczi A, Loche A, Castellano JM, Ronnekleiv OK, Bosch M, Kaidar G, Knoll JG, Wright H, Pfeifer GP, Ojeda SR. Epigenetic control of female puberty. Nat Neurosci. 2013;16:281–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Wyatt AK, Zavodna M, Viljoen JL, Stanton JA, Gemmell NJ, Jasoni CL. Changes in methylation patterns of kiss1 and kiss1r gene promoters across puberty. Genet Epigenet. 2013;5:51–62.

    CAS  PubMed  PubMed Central  Google Scholar 

  34. Semaan SJ, Dhamija S, Kim J, Ku EC, Kauffman AS. Assessment of Epigenetic Contributions to Sexually-Dimorphic Kiss1 Expression in the Anteroventral Periventricular Nucleus of Mice. Endocrinology. 2012;153:1875–86.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Gu H, Bock C, Mikkelsen TS, Jager N, Smith ZD, Tomazou E, Gnirke A, Lander ES, Meissner A. Genome-scale DNA methylation mapping of clinical samples at single-nucleotide resolution. Nat Methods. 2010;7:133–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Clegg DJ, Brown LM, Zigman JM, Kemp CJ, Strader AD, Benoit SC, Woods SC, Mangiaracina M, Geary N. Estradiol-dependent decrease in the orexigenic potency of ghrelin in female rats. Diabetes. 2007;56:1051–8.

    Article  CAS  PubMed  Google Scholar 

  37. Gao H, Sisley S, Yallampalli C. Blunted hypothalamic ghrelin signaling reduces diet intake in rats fed a low-protein diet in late pregnancy. Physiol Rep. 2015;3

  38. Gu H, Smith ZD, Bock C, Boyle P, Gnirke A, Meissner A. Preparation of reduced representation bisulfite sequencing libraries for genome-scale DNA methylation profiling. Nat Protoc. 2011;6:468–81.

    Article  CAS  PubMed  Google Scholar 

  39. Zhang X, Yazaki J, Sundaresan A, Cokus S, Chan SW, Chen H, Henderson IR, Shinn P, Pellegrini M, Jacobsen SE, Ecker JR. Genome-wide high-resolution mapping and functional analysis of DNA methylation in arabidopsis. Cell. 2006;126:1189–201.

    Article  CAS  PubMed  Google Scholar 

  40. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25:25–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, Katayama T, Kawashima S, Okuda S, Tokimatsu T, Yamanishi Y. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36:D480–4.

    Article  CAS  PubMed  Google Scholar 

  42. Zilberman D, Gehring M, Tran RK, Ballinger T, Henikoff S. Genome-wide analysis of Arabidopsis thaliana DNA methylation uncovers an interdependence between methylation and transcription. Nat Genet. 2007;39:61–9.

    Article  CAS  PubMed  Google Scholar 

  43. Ziller MJ, Gu H, Muller F, Donaghey J, Tsai LT, Kohlbacher O, De Jager PL, Rosen ED, Bennett DA, Bernstein BE, et al. Charting a dynamic DNA methylation landscape of the human genome. Nature. 2013;500:477–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Dyachenko OV, Shevchuk T, Kretzner L, Buryanov YI, Smith SS. Human non-CG methylation Are human stem cells plant-like? Epigenetics. 2010;5:569–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Yan J, Zierath JR, Barres R. Evidence for non-CpG methylation in mammals. Exp Cell Res. 2011;317:2555–61.

    Article  CAS  PubMed  Google Scholar 

  46. Cokus SJ, Feng S, Zhang X, Chen Z, Merriman B, Haudenschild CD, Pradhan S, Nelson SF, Pellegrini M, Jacobsen SE. Shotgun bisulphite sequencing of the Arabidopsis genome reveals DNA methylation patterning. Nature. 2008;452:215–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Lomniczi A, Ojeda SR. The Emerging Role of Epigenetics in the Regulation of Female Puberty. Endocr Dev. 2016;29:1–16.

    PubMed  Google Scholar 

  48. Kim HG, Kurth I, Lan F, Meliciani I, Wenzel W, Eom SH, Kang GB, Rosenberger G, Tekin M, Ozata M, et al. Mutations in CHD7, Encoding a Chromatin-Remodeling Protein, Cause Idiopathic Hypogonadotropic Hypogonadism and Kallmann Syndrome. Am J Hum Genet. 2008;83:511–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Watson RE Jr, Langub MC Jr, Engle MG, Maley BE. Estrogen-receptive neurons in the anteroventral periventricular nucleus are synaptic targets of the suprachiasmatic nucleus and peri-suprachiasmatic region. Brain Res. 1995;689:254–64.

    Article  CAS  PubMed  Google Scholar 

  50. Williams WP 3rd, Jarjisian SG, Mikkelsen JD, Kriegsfeld LJ. Circadian control of kisspeptin and a gated GnRH response mediate the preovulatory luteinizing hormone surge. Endocrinology. 2011;152:595–606.

    Article  CAS  PubMed  Google Scholar 

  51. Vida B, Deli L, Hrabovszky E, Kalamatianos T, Caraty A, Coen CW, Liposits Z, Kallo I. Evidence for suprachiasmatic vasopressin neurones innervating kisspeptin neurones in the rostral periventricular area of the mouse brain: regulation by oestrogen. J Neuroendocrinol. 2010;22:1032–9.

    Article  CAS  PubMed  Google Scholar 

  52. Pitynski D, Flynn FW, DC S: Does salt have a permissive role in the induction of puberty? Med Hypotheses 2015, 85:463-467.

  53. Kriegsfeld LJ, Leak RK, Yackulic CB, LeSauter J, Silver R. Organization of suprachiasmatic nucleus projections in Syrian hamsters (Mesocricetus auratus): An anterograde and retrograde analysis. J Comp Neurol. 2004;468:361–79.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Plant TM. Neuroendocrine control of the onset of puberty. Front Neuroendocrinol. 2015;38:73–88.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Rakyan VK, Down TA, Thorne NP, Flicek P, Kulesha E, Graf S, Tomazou EM, Backdahl L, Johnson N, Herberth M, et al. An integrated resource for genome-wide identification and analysis of human tissue-specific differentially methylated regions (tDMRs). Genome Res. 2008;18:1518–29.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Fan S, Zhang X. CpG island methylation pattern in different human tissues and its correlation with gene expression. Biochem Biophys Res Commun. 2009;383:421–5.

    Article  CAS  PubMed  Google Scholar 

  57. Yang C, Ye J, Li XM, Gao XX, Zhang KF, Luo L, Ding JP, Zhang YH, Li YS, Cao HG, et al. DNA Methylation Patterns in the Hypothalamus of Female Pubertal Goats. PLoS One. 2016;11

  58. Dumasia K, Kumar A, Deshpande S, Balasinor NH. Estrogen signaling, through estrogen receptor beta, regulates DNA methylation and its machinery in male germ line in adult rats. Epigenetics. 2017;12:476–83.

    Article  PubMed  Google Scholar 

Download references

Acknowledgements

We are grateful to BGI Genomics Co.,Ltd. for sequencing works. We would like to thank the members of Anhui Provincial Laboratory of Animal Genetic Resources Protection and Breeding, Anhui Provincial Laboratory for Local Livestock and Poultry Genetic Resource Conservation and Bio-Breeding for their fruitful discussions and valuable sggestions.

Funding

This work was supported by the National Natural Science Foundation of China [grant number31472096] and the Anhui Provincial Natural Science Foundation [grant number 1408085MKL40].

Availability of data and materials

All data generated or analyzed during this study are included in this published article (and its Additional file 1).

Author’s contributions

LL and ZQY conceived of the study, participated in its design and coordination and drafted the manuscript; YJ, TY and YC conducted qRT-PCR validation and statistical analysis; GXX and SM performed the statistical analysis; LY, LYH, ZXR and ZYH carried out the analysis of data, animal experiments, and surgical processes; FFG and ZYH participated in the design and coordination and helped with manuscript writing. All authors read and approved the final manuscript.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Fugui Fang.

Ethics declarations

Ethics approval and consent to participate

All experimental procedures were approved by the ethics committee of Anhui Agricultural University, Anhui, China.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional file

Additional file 1:

Table S1. Comparison of reads within reference sequence. Table S2. The number of promoters and CGI covered by RRBS. Table S3. The cellulara component term with differentially methylated genes in promoter. Table S4. The molecular function term with differentially methylated genes in promoter. Table S5. The biological process term with differentially methylated genes in promoter. Table S6. The cellulara component term with differentially methylated genes in CGI. Table S7. The molecular function term with differentially methylated genes inCGI. Table S8. The biological process term with differentially methylated genes in CGI. Table S9. The pathways with differentially methylated genes in promoter. Table S10. The pathways with differentially methylated genes in CGI. (DOC 830 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Luo, L., Yao, Z., Ye, J. et al. Identification of differential genomic DNA Methylation in the hypothalamus of pubertal rat using reduced representation Bisulfite sequencing. Reprod Biol Endocrinol 15, 81 (2017). https://doi.org/10.1186/s12958-017-0301-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12958-017-0301-2

Keywords