Genome-Wide DNA Methylation Analysis of Hypothalamus During the Onset of Puberty in Gilts

Although selection of the early age at puberty in gilts will make for a favorable effect on the reproductivity of sow, a large proportion of phenotypic variation in age at puberty of gilts cannot be explained by genetics. Previous studies have implicated hypothalamic DNA methylation in the onset of puberty in mammals. However, the underlying molecular mechanism regarding the regulation of the onset of puberty has remained largely unexplored in gilts. Herein, the genome-scale DNA methylation of hypothalamus was acquired, using the reduced representation bisulfite sequencing, to compare and describe the changes of DNA methylation across Pre-, In- and Post-pubertal gilts. In this study, the average methylation levels of CpGs and CpHs (where H = C, T, or A) in CpG islands- and gene-related regions were gradually decreased in hypothalamic methylomes during the pubertal transition. Comparisons of Pre- vs. In-, In- vs. Post-, and Pre- vs. Post-pubertal stage revealed that there were 85726, 92914, and 100421 differentially methylated CpGs and 5940, 14804, and 16893 differentially methylated CpHs (where H = C, T, or A) in the hypothalamic methylomes. The methylation changes of CpHs were more dynamic than that of CpGs, and methylation changes of CpGs and CpHs were likely to be, respectively, involved in the developmental processes of reproduction and the molecular processes of cellular communications in the hypothalamus. Moreover, methylation changes of CpHs were observed to overrepresent in the quantitative trait loci of age at puberty, and the biological function of these CpH methylation changes was enriched in the pancreas development in gilts. Furthermore, the mRNA levels of several differentially CpG or CpH methylated genes related to the transcription of RNA II polymerase, GnRH signaling pathway, Estrogen signaling pathway, PI3K-AKt signaling pathway, and Insulin signaling pathway, including MAX, MMP2, FGF11, IGF1R, FGF21, and GSK3B, were significantly changed across these pubertal stages in the hypothalamus. These results will help our understanding of how DNA methylation contributes to phenotypic variation of age at puberty.

Introduction

In female pigs, the onset of puberty is widely thought to be the processes of sexual maturations of gilts to be capable of reproduction. Many evidences suggest that gilts that have an earlier age at puberty are more likely to farrow multiple litters (Tummaruk et al., 2007), give birth to more piglets (Patterson et al., 2010) and have a longer reproductive longevity (Roongsitthichai et al., 2014). Nevertheless, with the intensive artificial selection for growth rate and lean meat, recent studies have suggested that approximately 30% of gilts are culled due to the delayed puberty (Stancic et al., 2009, 2011) with failure to display the first estrus by 240 days (Nonneman et al., 2014), which is consequent in a heavily economic burden on the modern commercial farms. However, the underlying molecular mechanism regarding the regulation of the onset of puberty has remained largely unexplored in pigs.

In gilts, age at puberty, a trait with moderate heritability (h² = 0.38 to 0.46) (Kuehn et al., 2009), is usually defined as the day of first standing estrus in modern commercial farms. Using 759 gilts genotyped by Illumina PorcineSNP 60 BeadChip, Nonneman et al. (2016) found two quantitative trait loci (QTLs) which explained 16.87% of the total genetic variance of age at puberty. Using the imputed whole genome sequences, Xin et al. (2018) identified 2339 SNPs which were associated with pubertal reproductive failure, explaining 16% of the phenotypic variances. These studies suggest that genetic factors only can explain and account for ∼16% of the phenotypic variances of age at puberty in pigs. Recently, changes of DNA methylation induced by 5-azacytidine have been demonstrated to disrupt the pulsatile release of gonadotropin-releasing hormone (GnRH) from the hypothalamus and delay the onset of puberty (Lomniczi et al., 2013, 2015), suggesting that DNA methylation plays a vitally regulatory role in the onset of puberty and may contribute to the phenotypic variation of age at puberty in pigs.

In female pigs, the onset of puberty is regulated and driven by GnRH released from the hypothalamus (Pelletier et al., 1981; Trout et al., 1984). We hypothesized that DNA methylation of hypothalamus might exhibit dramatic changes across the pubertal transition in gilts. Herein, the objectives of this study were to profile the genome-scale DNA methylation of hypothalamus for Pre-, In- and Post-pubertal gilts by using the reduced representation bisulfite sequencing (RRBS), and these methylation profiles were compared to describe the dynamics and changes across pubertal transition in both CpG and CpH (where H = C, T, or A) contexts. Then the association between methylation changes and QTLs of age at puberty in pigs were further explored, and the biological functions of these methylation changes were identified during the onset of puberty. These results are expected to provide useful insights into the epigenetic mechanism for the onset of puberty in gilts.

Materials and Methods

Ethics Statement

The pig cares and the experiments were conducted according to the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China) and were approved by the Animal care and Use Committer of South China Agricultural University, Guangzhou, China (approval number: SCAU#2013-10). The pigs were fed the same diet ad libitum and reared under the same conditions in the same environments.

Animals and Sample Preparation

The hypothalamus samples were collected from Landrace × Yorkshire crossbred gilts. The onset of pig puberty could be handily identified by the standing reflex with the back-pressure test and boar contact (Patterson et al., 2002). 25 Landrace × Yorkshire crossbred gilts aged at 160 days were selected and used in this study. Pubertal signs were checked and recorded twice daily at 09:00 and 15:30 by inspection of the vulva and assessment of the standing reflex for these 25 gilts. Three gilts aged 180 days without pubertal signs were selected as the Pre-pubertal gilts; three gilts at the day exhibiting the first estrus and standing reflex were selected as the In-pubertal gilts (about 205 days); and another three gilts in the dioestrus phase, 14 days after the day exhibiting the first estrus and standing reflex, were selected as the Post-pubertal gilts. Pigs were fed the same diet daily and reared in the same conditions and environments. The hypothalamus of each pubertal group was collected from and was frozen quickly in liquid nitrogen and then stored at −80°C for further using.

Constructions of RRBS Library

The library constructions and sequencing services were provided by RiboBio Co., Ltd., (Guangzhou, China) as previously described in our studies (Yuan et al., 2016, 2017a,b). The genomic DNA of these hypothalamus tissues was extracted using a DNeasy Blood and Tissue Kit (Qiagen, Beijing), and then, after checking on the quality of the extracted DNA, one library was built for each gilts based on previously published RRBS studies (Yuan et al., 2016, 2017a,b). The processes and procedures of RRBS libraries were briefed as follows. Firstly, the purified genomic DNA was digested overnight with MspI (New England Biolabs, United States). For the MspI digested segments, the sticky ends were filled with CG nucleotides and 3′ A overhangs were added. Secondly, methylated Illumina sequencing adapters with 3′ T overhangs were ligated to the digested segments, and the products obtained were purified. Then 110–220 bp fragments were selected (Yuan et al., 2017b) and converted by bisulfite using an EZ DNA Methylation Gold Kit (Zymo Research, United States). Lastly, libraries of 110–220 bp fragments were PCR amplified by 18 cycles, and each library was sequenced using one lane of an Illumina HiSeq 2500 and 100 bp paired-end reads.

Bioinformatic Analysis of RRBS Data

The first two nucleotides were trimmed from all the second read sequences to blunt-end the MspI site. All reads were trimmed using Trim Galore (v0.4.0) software (Babraham Bioinformatics¹) and a Phred quality score of 20 as the minimum. The adaptor pollution reads and multiple N reads (where N > 10% of one read) were removed to generate the clean reads. The quality control checks were performed by FastQC (v0.11.3) software (Babraham Bioinformatics). The clean reads were analyzed and performed by Bismark (v0.14.5) (Krueger and Andrews, 2011) using the default parameters. The clean reads were mapped to the pig reference genome (Sscrofa 11.1, downloaded from Ensembl²), by using “bismark” function, the PCR duplicates were removed by using “deduplicate_bismark” function, and then the DNA methylation calling was extracted by using “bismark_methylation_extraxtor” function.

After DNA methylation calling by Bismark for these nine RRBS data, 1320853 CpGs and 5884256 CpHs covered by at least five reads and coexisted across all tissues were remained for further analysis. The methylation level of the CpGs was calculated as the methylated reads divided the total covered reads. The average methylation level of each hypothalamic group was calculated by the average methylation level across the three replicates. For each specific region, the methylation level was measured as the average level of CpGs located in this region. To profile the DNA methylation patterns at the gene and CGI locations, the gene locations were divided into 20, 40, and 20 bins for 5 kb upstream region of the transcription start sites (TSSs), gene body and 5 kb downstream region of transcription end sites (TESs), respectively, and the CGI locations were divided into 20, 20, and 20 bins for 2 kb upstream region, CGIs and 2 kb downstream region, respectively. These analyses were performed by Perl and R scripts. This pipeline was carefully described by our previous studies (Yuan et al., 2016, 2017a).

Annotation of Detected CpGs and CpHs

The locations of porcine CGIs were downloaded from UCSC³. CGIs were denoted as the same to UCSC and described as regions >200 bp in length, with a C and G percentage >0.5, and a ratio of the observed CpG/expected CpG > 0.6. The expected CpG was calculated as the number of Cs multiplied by the number of Gs, divided by the length of the segment. The +/−2 kb regions outside of CGIs were defined as CGI shores, and the +/−2 kb regions outside of CGI shores were defined as CGI shelves.

The locations of genes were downloaded from Ensembl (see text footnote 2). Based on genic locations of Ensembl, the porcine genome was separated into five genic features, which were upstream, exonic, intronic, downstream and intergenic regions. The upstream region was 5 kb upstream region of the TSS. The exon was defined as the integration of 5′UTR, CDS, and 3′UTR arranging from the TSS to the TES. The intron was determined as the integration of introns arranged from the TSS to the TES. The downstream region was 5 kb downstream region of the TES. The intergenic region was denoted as the outside regions of upstream, exonic, intronic, and downstream regions.

Quantitative Real-Time PCR

The total RNA was extracted using TRIzol reagent (TaKaRa, Tokyo, Japan) and reverse-transcribed using a RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, United States) for mRNAs. The relative expression levels of the mRNAs were quantified using Maxima SYBR Green qRT-PCR Master Mix (2×) (Thermo Fisher Scientific) and THUNDERBIRD SYBR qPCR Mix (Toyobo) on a LightCycler Real-Time PCR system. The expression levels of GAPDH mRNAs were used as endogenous controls, and the fold changes were calculated using the 2^−ΔΔct method. The primer sequences are listed in Supplementary Table S1.

Statistics Analysis

The significant differences between two groups were tested using the Student’s test. The CGIs were identified as the differentially methylated CGIs (DMIs) if the average methylation levels of CGIs, which were covered with at least 20 CpGs or CpHs, were significantly different according to the P-value from the Student’s test corrected by the false discovery rate (P ≤ 0.05). The genes, including 5 kb upstream flanking region, gene body and 5 kb downstream flanking region, were identified as the differentially methylated genes (DMGs) if the average methylation levels of genes, which were covered with at least 20 CpGs or CpHs, were significantly different according to the P-value from the Student’s test corrected by the false discovery rate (P ≤ 0.05). The genes, including 5 kb upstream flanking region, gene body and 5 kb downstream flanking region, overlapped at least one DMIs were defined as the DMIs regarding genes.

The differentially methylated CpG (DMCs) and CpH sites (DMHs) were calculated by CGmapTools (Guo et al., 2018). The CpGs or CpHs, whose methylation levels changed more than 20%, were identified as DMCs or DMHs according to the P-value from a two-tail Fisher’s exact test corrected by the false discovery rate (P ≤ 0.05). For the comparison between two groups, CpGs and CpHs with higher methylation level in one group were considered as the hyper-methylated CpGs or CpHs in this group and were considered as the hypo-methylated CpGs or CpHs in another group, respectively. The genes, including 5 kb upstream flanking region, gene body and 5 kb downstream flanking region, contained at least one DMC were defined as the DMC regarding genes. The genes, including 5 kb upstream flanking region, gene body and 5 kb downstream flanking region, contained at least one DMH were defined as the DMH regarding genes. The enrichment for certain genomic regions was calculated by using a two-tail Fisher’s exact test. The Gene Ontology (GO) and Kyoto Encyclopedia of Gene and Genomes (KEGG) enrichment analysis were undertaken by R package “clusterProfiler” (Yu et al., 2012).

Results

Genome-Wide DNA Methylation of Hypothalamus During the Onset of Puberty

In the CpG context, the average methylation levels were 53.57% ± 0.07, 53.50% ± 0.41, and 53.01% ± 0.83% for Pre-, In-, and Post-puberty, respectively. The average methylation levels of CGI-related (Figure 1A) and gene-related regions (Figure 1B) were the highest in Pre- and the lowest in Post-pubertal methylome. Compared to other CGI- or gene-related regions (Figures 1A,B), the average methylation levels of CGIs and upstream regions displayed the lowest levels. In the CpH context, the average methylation levels were 0.75% ± 0.04, 0.73% ± 0.03, and 0.70% ± 0.05% for Pre-, In-, and Post-puberty, respectively. The CpH methylation levels of CGI-related (Figure 1C) and gene-related regions (Figure 1D) displayed the same pattern with that in CpG methylomes (Figures 1A,B). For Pre- vs. In-, In- vs. Post- and Pre- vs. Post-puberty, the Pearson’s correlation coefficients of methylation levels were all 0.99 (P < 2.22 × 10⁻¹⁶) (Figure 2A and Supplementary Table S2) for CpGs but were 0.82 (P < 2.22 × 10⁻¹⁶), 0.85 (P < 2.22 × 10⁻¹⁶), and 0.89 (P < 2.22 × 10⁻¹⁶) for CpHs (Figure 2B and Supplementary Table S2), respectively.

FIGURE 1

Figure 1. Genome-wide DNA methylation of hypothalamus tissues across pubertal transition in gilts. Average methylation levels of CpGs in CGI- (A) and gene-related regions (B). Average methylation levels of CpHs in CGI- (C) and gene-related regions (D).

FIGURE 2

Figure 2. Global DNA methylation of CpGs and CpH in hypothalamus tissues. The global CpG (A) and CpH (B) methylation levels in Pre- (track 1), In- (track 2), and Post-pubertal (track 3) methylomes of hypothalamus tissues, from outside to inside, were quantified per 1 Mb window. The densities of CGIs (track 4), genes (track 5) and CpGs (track 6 in A) or CpHs (track 6 in B) were also quantified per 1 Mb window. The labels outside of track 1 represent the chromosomes of the porcine genome.

DNA Methylation Patterns of CpGs and CpHs in CGI and Gene Locations

To explore the dynamics and similarities of DNA methylation at the locations of genes and CGIs during the onset of puberty, the methylation patterns of CpGs and CpHs were profiled for the hypothalamic methylomes of Pre-, In-, and Post-puberty (Figure 3), respectively. We found that the methylation levels of CpGs were likely to overlap together at the locations of genes (Figure 3A) and CGIs (Figure 3B). However, the methylation levels could be clearly distinguished from each other in the upstream or downstream regions around the CGIs (Figure 3B). Besides, the methylated patterns of CpHs (Figures 3C,D) exhibited the similarly methylated patterns with that of CpGs (Figures 3A,B).

FIGURE 3

Figure 3. Methylation patterns of CpGs and CpHs at CGI and gene locations. Methylation levels and patterns of CpGs at the locations of genes (A) and CGIs (B). Methylation levels and patterns of CpHs at the locations of genes (C) and CGIs (D). Up5k is denoted as the 5 kb upstream regions of transcription start sites. Down5k is denoted as the 5 kb downstream regions of transcription end sites. Up2k and Down2k are identified as the 2 kb upstream and downstream regions of the CGIs, respectively.

To further determine the DNA methylation dynamics at the locations of genes and CGIs, the DMIs and DMGs were calculated for Pre- vs. In-, In- vs. Post- and Pre- vs. Post-puberty, respectively. There were 221, 207, and 236 DMIs, and there were 290, 238, and 281 DMGs in the CpG context (Supplementary Figures S1A–C). There were 538, 568, and 571 DMIs, and there were 329, 300, and 321 DMGs in the CpH context for Pre vs. In-, In- vs. Post- and Pre- vs. Post-puberty (Supplementary Figures S1A–C), respectively. The GO and KEGG pathway analysis were performed on these DMGs and DMIs regarding genes. In the CpG context, the most significantly enriched GO and KEGG terms were Metanephric mesenchyme development, Kidney mesenchyme development, Signaling pathways regulating pluripotency of stem cells, mTOR signaling pathway, Wnt signaling pathway for DMIs regarding genes (Supplementary File S1); and were cardiovascular system development, and Signaling pathways regulating pluripotency of stem cells for DMGs (Supplementary File S1); In the CpH context, the most significantly enriched GO and KEGG terms in CpH were RNA polymerase II transcription factor binding, transcription factor binding, PI3K-Akt signaling pathway, MAPK signaling pathway and Oocyte meiosis for DMIs regarding genes (Supplementary File S2); and were phosphatidylinositol bisphosphate binding, phosphatidylinositol phosphate binding, and Lysine degradation for DMGs (Supplementary File S2).

Methylation Changes During the Onset of Puberty in Hypothalamic Genome

To investigate the methylation changes among hypothalamic methylomes across the Pre-, In-, and Post-puberty, the DMCs and DMHs were counted for Pre- vs. In-, In- vs. Post- and Pre- vs. Post-puberty. Respectively, 85726, 92914, and 100421 DMCs were identified for Pre- vs. In- (Figure 4A), In- vs. Post- (Figure 4B) and Pre- vs. Post-puberty (Figure 4C), representing 6.49, 7.03, and 7.60% of all detected CpGs (Table 1). These DMCs were observed to be hypo-methylated in Post-methylome but hyper-methylated in Pre- and In-methylomes (Figures 4A–C). Moreover, these DMCs were observed to overrepresent in CGI shores, CGI shelves, Introns, downstream and intergenic regions (relative enrichment = 1.11–1.49, P ≤ 9.29 × 10⁻³) (Table 1) but underrepresent in CGI, exons and upstream regions (relative enrichment = 0.46–0.68, P < 2.22 × 10⁻¹⁶) (Table 1).

FIGURE 4

Figure 4. Differentially methylated CpGs and CpHs in hypothalamic methylomes across pubertal transition. Distributions of differentially methylated CpGs between Pre- and In- (A), In- and Post- (B), and Pre- and Post-puberty (C). Distributions of differentially methylated CpHs between Pre- and In- (D), In- and Post- (E), and Pre- and Post-puberty (F).

TABLE 1

Table 1. Distribution of differentially methylated CpGs among these hypothalamic methylomes.

Then 15940, 14804, and 16893 DMHs were identified for Pre- vs. In- (Figure 4D), In- vs. Post- (Figure 4E), and Pre- vs. Post-puberty (Figure 4F), representing 0.27, 0.25, and 0.29% of all detected CpHs (Table 2), respectively. These DMHs were also observed to be hypo-methylated in Post-methylome but hyper-methylated in Pre- and In-methylomes (Figures 4D–F). Furthermore, we found that DMHs showed a significant underrepresentation in CGIs, upstream regions and exons (relative enrichment = 0.48–0.58, P < 2.22 × 10⁻¹⁶) (Table 2), but CGI shores, CGI shelves, introns and intergenic regions showed a significant overrepresentation of DMHs (relative enrichment = 1.07–1.41, P ≤ 6.53 × 10⁻³) (Table 2). Interestingly, DMHs of downstream regions enriched in a stage-specific manner among the comparisons of Pre- vs. In-, In- vs. Post-, and Pre- vs. Post-puberty (Table 2).

TABLE 2

Table 2. Distribution of differentially methylated CpHs among these hypothalamic methylomes.

Global Distributions of DMCs and DMHs in Hypothalamic Genomes

The distributions of DMCs and DMHs in the hypothalamic methylomes were profiled in Figure 5. We found that these DMCs and DMHs were likely to enrich at the ends of each chromosome (Figure 5). Whereas DMCs exhibited additional higher peaks in the bodies of chromosome 3, chromosome 5, and chromosome 6 (Figure 5A), and DMHs displayed additional higher peaks in the body of chromosome 6 (Figure 5B). Moreover, DMHs were likely to show a stage-specific distribution among the onset of puberty, e.g., the distributions of DMHs in the chromosome 3 (Figure 5B).

FIGURE 5

Figure 5. Global distribution of differentially methylated CpG and CpH sites in hypothalamic genomes during the pubertal transition. The global distributions of differentially methylated CpG (A) and CpH (B) sites in Pre- vs. In- (track 1), In- vs. Post- (track 2), and Pre- vs. Post-pubertal (track 3) methylomes of hypothalamus tissues, from outside to inside, were quantified per 1 Mb window. The densities of CGIs (track 4), genes (track 5), and CpGs (track 6 in A) or CpHs (track 6 in B) were also quantified per 1 Mb window. The labels outside of track 1 represent the chromosomes of the porcine genome.

The Pearson’s correlation coefficients among the distributions of DMCs were all 0.99 (P < 2.22 × 10⁻¹⁶) for Pre- vs. In-, In- vs. Post-, and Pre- vs. Post-puberty (Supplementary Table S3). The distributions of DMCs were observed to be highly positively correlated with the densities of CpGs (Pearson’s correlation coefficients = 0.86, P < 2.22 × 10⁻¹⁶) and CGIs (Pearson’s correlation coefficients = 0.94, P < 2.22 × 10⁻¹⁶) (Supplementary Table S3), and the distributions of DMCs were moderately linked with the densities of genes (Pearson’s correlation coefficients = 0.32, P < 2.22 × 10⁻¹⁶) (Supplementary Table S3). However, the Pearson’s correlation coefficients among the distributions of DMHs were 0.94 (P < 2.22 × 10⁻¹⁶), 0.93 (P < 2.22 × 10⁻¹⁶), and 0.95 (P < 2.22 × 10⁻¹⁶) for Pre- vs. In-, In- vs. Post-, and Pre- vs. Post-puberty (Supplementary Table S3), respectively. The distributions of DMHs were highly associated with the densities of CGIs (Pearson’s correlation coefficients = 0.85–0.86, P < 2.22 × 10⁻¹⁶) and CpHs (Pearson’s correlation coefficients = 0.54–0.57, P < 2.22 × 10⁻¹⁶) but were slightly correlated with the densities of genes in hypothalamic methylome (Pearson’s correlation coefficients = 0.16–0.18, P ≤ 3.94 × 10⁻¹⁰) (Supplementary Table S3).

Enrichments of DMCs and DMHs in the QTL of Age at Puberty in Pigs

To further explore the associations between DNA methylation changes and the QTLs of age at puberty in pigs, 237 QTLs of age at puberty were collected from PigQTLdb (Release 36). These QTLs covered 57332 CpGs and 251222 CpHs, respectively. The numbers of DMCs located at these QTLs were 3930, 4387, and 4691 for Pre- vs. In-, In- vs. Post-, and Pre- vs. Post-puberty, respectively (Table 1). Compared to the whole hypothalamic genome, these DMCs were not obviously enriched in QTLs (relative enrichment = 1.06–1.09, P ≤ 5.76 × 10⁻⁴) (Table 1). The number of DMHs located at QTLs were 796, 788, and 826 for Pre- vs. In-, In- vs. Post-, and Pre- vs. Post-puberty, respectively (Table 2). Moreover, these DMHs were observed to overrepresent in QTLs (relative enrichment = 1.15–1.25, P ≤ 1.83 × 10⁻⁴) (Table 2). To gain insight into the biological functions in which methylation changes of CpHs might be involved, we performed the biological processes of GO and KEGG pathway analysis on genes (Supplementary File S3) that associated with at least one DMH located at QTLs. We found no signaling pathway was enriched, and only one GO term, pancreas development, was enriched.

Biological Functions of DMCs and DMHs in Hypothalamic Genomes Across the Pubertal Transition

To further gain insight into the biological functions in which DMC or DMH genes might be involved, we performed the biological processes of GO and KEGG pathway analysis on genes who were associated with at least one DMC (Figure 6) and DMH (Figure 7). For DMCs regarding genes (Supplementary File S4), we found that the most significantly enriched GO terms were the development (head development, central nervous system development, brain development, gland development, etc.), developmental process involved in reproduction, regulation of transcription from RNA polymerase II promoter, animal organ morphogenesis, and so on (Figure 6A). Moreover, these genes were significantly enriched in Acute myeloid leukemia, Neurotrophin signaling pathway, Apoptosis, cAMP signaling pathway, Estrogen signaling pathway, and so on (Figure 6B).

FIGURE 6

Figure 6. Biological functions of DMCs in hypothalamic genomes across the pubertal transition. Significantly enriched GO terms (A) and KEGG pathways (B) of DMC genes.

FIGURE 7

Figure 7. Biological functions of DMHs in hypothalamic genomes across the pubertal transition. Significantly enriched GO terms (A) and KEGG pathways (B) of DMH genes.

For DMHs regarding genes (Supplementary File S4), the most significantly enriched terms were the ion transport (metal ion transport, regulation of metal ion transport, calcium ion transport, inorganic ion transmembrane transport, etc.), positive regulation of cell communication, positive regulation of signaling, import into cell, central nervous system neuron differentiation and so on (Figure 7A). Moreover, these genes were significantly enriched in PI3K-AKt signaling pathway, Proteoglycans in cancer, Focal adhesion, Axon guidance, Acute myeloid leukemia, Regulation of actin cytoskeleton, and so on (Figure 7B).

Furthermore, the mRNA expressions of several DMC or DMH regarding genes selected basing on literatures were detected across the pubertal transition in hypothalamus. We found that the mRNA levels of MAX (Figure 8A), MMP2 (Figure 8B), and FGF11 (Figure 8C), which were the DMC regarding genes and associated with the transcription of RNA II polymerase, GnRH signaling pathway and Estrogen signaling pathway, were significantly changed across Pre-, In-, and Post-puberty in hypothalamus. Also, the mRNA levels of IGF1R (Figure 8D), FGF21 (Figure 8E) and GSK3B (Figure 8F), which were DMH regarding genes and associated with PI3K-AKt signaling pathway and Insulin signaling pathway, were significantly changed across Pre-, In-, and Post-puberty in hypothalamus.

FIGURE 8

Figure 8. mRNA expression levels of several DMC or DMH regarding genes across the pubertal transition in hypothalamus. MAX (A), MMP2 (B), and FGF11 (C) were DMC regarding genes. IGF1R (D), FGF21 (E), and GSK3B (F) were DMH regarding genes. The significance of differences in means between two groups was analyzed by using Student’s t-test (two-tailed). Data were expressed as means ± standard error of biological replicates. ^∗∗ indicates P < 0.01.

Discussion

Although the early puberty is suggested to improve the reproductive performance and longevity (Patterson et al., 2010; Roongsitthichai et al., 2014), collecting age-at-puberty data for selection of the early puberty is tedious and labor intensive, which requires daily estrus detection in the presence of a boar. The appearance that approximately 30% of gilts are culled due to the delayed puberty (Stancic et al., 2009, 2011) markedly harms the financial stake of the modern commercial farms. Recent advance made in “omics” technologies are allowed for exploring the complicated mechanism underlying the onset of puberty and contributing to the selection and breeding practices of age at puberty in gilts.

In this study, we found that the methylation levels of CpGs were likely to overlap together at the genic regions (Figure 3A) and CGIs (Figure 3B). However, the methylation levels of CpHs were likely to be clearly distinguished from each other at gene bodies (Figure 3C) and at the locations of and CGIs (Figure 3D). Respectively, the Pearson’s correlation coefficients were all 0.99 (P < 2.22 × 10⁻¹⁶) (Figure 2A and Supplementary Table S2) for the methylation levels of CpGs but were 0.82–0.89 (P < 2.22 × 10⁻¹⁶) (Figure 2B and Supplementary Table S2) for the methylation levels of CpGs between Pre- vs. In-, In- vs. Post-, and Pre- vs. Post-puberty. Furthermore, during the pubertal transition, the numbers of DMIs and DMGs in the CpG context were more than that in the CpH context (Supplementary Figures S1A–C). These results indicated that the methylation levels of CpHs were more variable than that of CpGs among the hypothalamic methylomes of Pre-, In-, and Post-puberty.

The DNA methylation changes of CpGs were likely to play an important role in regulation of transcription from RNA polymerase II, Neurotrophin signaling pathway and Estrogen signaling pathway (Figure 6B), which were strongly associated with biological functions of hypothalamus for the onset of puberty in pigs (Trout et al., 1984; Ieda et al., 2014). Besides, compared to Pre- and Post-puberty, we found that MAX, which was MYC associated factor X and might stimulate the transcription of RNA II polymerase (Campbell and White, 2014), showed the lowest mRNA level in In-puberty (Figure 8A). MMP2, which was one of the key genes of the GnRH signaling pathway and Estrogen signaling pathway⁴, showed the highest mRNA level in In-puberty (Figure 8B). IGF1R, which was a key gene of Ovarian steroidogenesis pathway⁵ and played an important role in controlling the timing of puberty (Chang et al., 2018), showed the lowest mRNA level in In-puberty (Figure 8C). These observations suggested that the CpG methylation changes were likely to show a crucial regulatory role in controlling the transcription of genes related to the onset of puberty at hypothalamus level in gilts.

However, the DNA methylation changes of CpHs were mostly found to get involved in the ion metabolizing and the positive regulation of cell communication, import into cell, neuron differentiation (Figure 7A), which were likely to get involved in the maintenance of molecular function within cellular communications (Jang et al., 2017; Lee et al., 2017). Additionally, compared to Pre- and Post-puberty, FGF11, which was one of the key genes of PI3K-AKt signaling pathway, presented the lowest mRNA level in In-puberty (Figure 8D), but FGF21, which expressed in GnRH neurons and regulated GnRH secretion during the onset of puberty in mice (Misrahi, 2017), displayed the highest mRNA level in In-puberty (Figure 8E). GSK3B, which was one of genes of Insulin signaling pathway⁶ and associated with age of puberty in cattle (Fortes et al., 2013), displayed the lowest mRNA level in In-puberty (Figure 8F). These results indicated that methylation changes of CpGs were involved in the maintenance of molecular function within cellular communications during the initiation of puberty in gilts.

Recently, the genome-wide DNA methylation changes of hypothalamus have been investigated during the timing of puberty in rat (Luo et al., 2017) and goat (Yang et al., 2016), even the methylation variation between rat and goat (Yang et al., 2018). Results of rat (Luo et al., 2017) and goat (Yang et al., 2016) show that the global DNA methylation level is decreasing across Pre- to In-pubertal stages, which is in line with the present study. We found that the average CpG and CpH methylation levels of genomes, GCI- and gene-related regions in hypothalamic genomes were gradually decreased across Pre-, In-, and Post-puberty (Figure 1). It seems that the decreasing trend of global DNA methylation level from Pre- to In-pubertal stage happens whatever the mammalian species is, although the difference is not statistically significant. This decrease of the global DNA methylation, from Pre- to In-pubertal stage, may be related to the DNA methylation changes of genes which represses the initiating of puberty in Pre-puberty, and these genes may change from transcriptional repressions to activations. Moreover, many DMGs such as NLRC5, SMOC1, GRID1, ABAT, MAP3K4, and PTPRN reported in rat (Yang et al., 2018) and goat (Yang et al., 2018) across the pubertal transition were also identified in this study (Supplementary File S4). Furthermore, changes of the DNA methylation in the hypothalamus were all enriched in Neurotrophin signaling pathway and Estrogen signaling pathway for rat (Yang et al., 2018), goat (Yang et al., 2018)and pigs (Figure 6). These observations not only supported the accuracies of findings in the present study but also suggested the patterns of DNA methylation variations were like to be similar among goat, rat and pigs.

With the intensive artificial selection for growth rate and lean meat, recent studies have suggested that approximately 30% of gilts are culled due to the delayed puberty (Stancic et al., 2009, 2011) with failure to display the first estrus by 240 days (Nonneman et al., 2014). We found that, compared to DMCs (Table 1), DMHs were observed to overrepresent in QTLs of age at puberty in gilts (Table 1), and genes, which were associated with at least one DMH located at QTLs of age at puberty, were only significantly enriched in the GO term of pancreas development, which is highly associated with the onset of puberty in mammals (Bervoets et al., 2015). Previous studies have shown that the CpH methylation may play a repressive role in the gene’s transcription (Luo et al., 2016; Keown et al., 2017), and Insulin secreted from pancreases promotes the synthesis of glycogen, fat and protein (Bervoets et al., 2015; Singh et al., 2017), and the gilts with lean body composition were likely to show the delayed puberty (Gaughan et al., 1997). Moreover, compared to Pre- and Post-pubertal stages, IGF1R (Figure 8C) and GSK3B (Figure 8F), which are highly involved in synthesis and secretion of Insulin, show the lowest mRNA levels at In-pubertal stage. It is possible that the intensive artificial selection for growth rate and lean meat has damaged the development of pancreas and secretion of Insulin and then impacts the body composition and the timing of puberty in gilts. These observations recommend that DNA methylation may contribute to the phenotypic variation of age at puberty in pigs.

Conclusion

Collectively, during the pubertal transition in gilts, the methylation changes of CpHs were more dynamic than that of CpGs, and methylation changes of CpGs and CpHs were likely to be, respectively, involved in the developmental processes of reproduction and the molecular processes within cellular communications in the hypothalamus. Moreover, methylation changes of CpHs were observed to overrepresent in the QTLs of age at puberty, and the biological function of these CpG methylation changes were enriched in the pancreas development and might impact the body composition and the timing of puberty of gilts. These DNA methylation analyses will help our understanding of how DNA methylation contributes to phenotypic variation of age at puberty.

Data Availability

Publicly available datasets were analyzed in this study. This data can be found here: http://www.ebi.ac.uk/ena/data/view/PRJEB31293.

Author Contributions

XY, HZ, and JL conceived and designed the work. XY, XZ, YH, SY, and ZC acquired the biological samples and analyzed the data. XY, XZ, YK, and JL drafted the manuscript. NG, HZ, ZZ, and JL critically revised the manuscript. All authors reviewed and approved the final manuscript.

Funding

This work was supported by the earmarked fund for China Agriculture Research System (CARS-35), the Basic Work of Science and Technology Project (2014FY120800), Guangdong Sailing Program (2014YT02H042), and the Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.00228/full#supplementary-material

FIGURE S1 | Numbers of differentially methylated CpG islands and genes in hypothalamic methylomes across pubertal transition. The number of differentially methylated CpG islands (DMIs) and genes (DMGs) in Pre- vs. In-pubertal (A), In- vs. Post-pubertal (B), and Pre- vs. Post-pubertal (C) hypothalamus in the CpG or CpH context.

TABLE S1 | Primers used in the present study.

TABLE S2 | Correlation coefficients of CpG and CpH methylation with gene and CGI densities.

TABLE S3 | Correlation coefficients of DMCs and DMHs with gene and CGI densities.

FILE S1 | Significantly enriched GO terms and KEGG pathways of DMGs and DMIs regarding genes in the CpG context.

FILE S2 | Significantly enriched GO terms and KEGG pathways of DMGs and DMIs regarding genes in the CpH context.

FILE S3 | Regarding genes of DMCs and DMHs in QTLs.

FILE S4 | Regarding genes of DMCs and DMHs in hypothalamic genomes.

Footnotes

1. ^ http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/
2. ^ http://asia.ensembl.org/Sus_scrofa/Info/Index
3. ^ http://hgdownload.soe.ucsc.edu/goldenPath/susScr11/database/
4. ^ https://www.genome.jp/dbget-bin/www_bget?hsa:4313+hsa:4318
5. ^ https://www.genome.jp/kegg-bin/show_pathway?hsa04913+3480
6. ^ https://www.genome.jp/kegg-bin/show_pathway?ko04910+K03083

References

Bervoets, L., Mewis, A., and Massa, G. (2015). The shape of the plasma glucose curve during an oral glucose tolerance test as an indicator of beta cell function and insulin sensitivity in end-pubertal obese girls. Horm. Metab. Res. 47, 445–451. doi: 10.1055/s-0034-1395551

PubMed Abstract | CrossRef Full Text | Google Scholar

Campbell, K. J., and White, R. J. (2014). MYC regulation of cell growth through control of transcription by RNA polymerases I and III. Cold Spring Harb. Perspect. Med. 4:a018408. doi: 10.1101/cshperspect.a018408

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, H. P., Yang, S. F., Wang, S. L., and Su, P. H. (2018). Associations among IGF-1, IGF2, IGF-1R, IGF-2R, IGFBP-3, insulin genetic polymorphisms and central precocious puberty in girls. BMC Endocr. Disord. 18:66. doi: 10.1186/s12902-018-0271-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Fortes, M. R. S., Li, Y. T., Collis, E., Zhang, Y. D., and Hawken, R. J. (2013). The IGF1 pathway genes and their association with age of puberty in cattle. Anim. Genet. 44, 91–95. doi: 10.1111/j.1365-2052.2012.02367.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaughan, J. B., Cameron, R. D., Dryden, G. M., and Young, B. A. (1997). Effect of body composition at selection on reproductive development in large white gilts. J. Anim. Sci. 75, 1764–1772. doi: 10.2527/1997.7571764x

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, W., Zhu, P., Pellegrini, M., Zhang, M. Q., Wang, X., and Ni, Z. (2018). CGmapTools improves the precision of heterozygous SNV calls and supports allele-specific methylation detection and visualization in bisulfite-sequencing data. Bioinformatics 34, 381–387. doi: 10.1093/bioinformatics/btx595

PubMed Abstract | CrossRef Full Text | Google Scholar

Ieda, N., Uenoyama, Y., Tajima, Y., Nakata, T., Kano, M., Naniwa, Y., et al. (2014). KISS1 gene expression in the developing brain of female pigs in Pre- and peripubertal periods. J. Reprod. Dev. 60, 312–316. doi: 10.1262/jrd.2013-129

PubMed Abstract | CrossRef Full Text | Google Scholar

Jang, H. S., Shin, W. J., Lee, J. E., and Do, J. T. (2017). CpG and Non-CpG methylation in epigenetic gene regulation and brain function. Genes 8:E148. doi: 10.3390/genes8060148

PubMed Abstract | CrossRef Full Text | Google Scholar

Keown, C. L., Berletch, J. B., Castanon, R., Nery, J. R., Disteche, C. M., Ecker, J. R., et al. (2017). Allele-specific non-CG DNA methylation marks domains of active chromatin in female mouse brain. Proc. Natl. Acad. Sci. U.S.A. 114, E2882–E2890. doi: 10.1073/pnas.1611905114

PubMed Abstract | CrossRef Full Text | Google Scholar

Krueger, F., and Andrews, S. R. (2011). Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics 27, 1571–1572. doi: 10.1093/bioinformatics/btr167

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuehn, L. A., Nonneman, D. J., Klindt, J. M., and Wise, T. H. (2009). Genetic relationships of body composition, serum leptin, and age at puberty in gilts. J. Anim. Sci. 87, 477–483. doi: 10.2527/jas.2008-0936

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J. H., Park, S. J., and Nakai, K. (2017). Differential landscape of non-CpG methylation in embryonic stem cells and neurons caused by DNMT3s. Sci. Rep. 7:11295. doi: 10.1038/s41598-017-11800-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Lomniczi, A., Loche, A., Castellano, J. M., Ronnekleiv, O. K., Bosch, M., Kaidar, G., et al. (2013). Epigenetic control of female puberty. Nat. Neurosci. 16, 281–289. doi: 10.1038/nn.3319

PubMed Abstract | CrossRef Full Text | Google Scholar

Lomniczi, A., Wright, H., and Ojeda, S. R. (2015). Epigenetic regulation of female puberty. Front. Neuroendocrinol. 36:90–107. doi: 10.1016/j.yfrne.2014.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, C. Y., Lancaster, M. A., Castanon, R., Nery, J. R., Knoblich, J. A., and Ecker, J. R. (2016). Cerebral organoids recapitulate epigenomic signatures of the human fetal brain. Cell Rep. 17, 3369–3384. doi: 10.1016/j.celrep.2016.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Misrahi, M. (2017). beta-Klotho sustains postnatal GnRH biology and spins the thread of puberty. EMBO Mol. Med. 9, 1334–1337. doi: 10.15252/emmm.201708180

PubMed Abstract | CrossRef Full Text | Google Scholar

Nonneman, D., Lents, C., Rohrer, G., Rempel, L., and Vallet, J. (2014). Genome-wide association with delayed puberty in swine. Anim. Genet. 45, 130–132. doi: 10.1111/age.12087

PubMed Abstract | CrossRef Full Text | Google Scholar

Nonneman, D. J., Schneider, J. F., Lents, C. A., Wiedmann, R. T., Vallet, J. L., and Rohrer, G. A. (2016). Genome-wide association and identification of candidate genes for age at puberty in swine. BMC Genet. 17:50. doi: 10.1186/s12863-016-0352-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Patterson, J. L., Beltranena, E., and Foxcroft, G. R. (2010). 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. 88, 2500–2513. doi: 10.2527/jas.2008-1756

PubMed Abstract | CrossRef Full Text | Google Scholar

Patterson, J. L., Willis, H. J., Kirkwood, R. N., and Foxcroft, G. R. (2002). Impact of boar exposure on puberty attainment and breeding outcomes in gilts. Theriogenology 57, 2015–2025. doi: 10.1016/S0093-691X(02)00674-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Pelletier, J., Carrez-Camous, S., and Thiery, J. C. (1981). Basic neuroendocrine events before puberty in cattle, sheep and pigs. J. Reprod. Fertil. Suppl. 30, 91–102.

PubMed Abstract | Google Scholar

Roongsitthichai, A., Olanratmanee, E. O., and Tummaruk, P. (2014). Influence of growth rate and onset of boar contact on puberty attainment of replacement gilts raised in Thailand. Trop. Anim. Health Prod. 46, 1243–1248. doi: 10.1007/s11250-014-0634-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, R. G., Yoon, H. D., Poppitt, S. D., Plank, L. D., and Petrov, M. S. (2017). Ectopic fat accumulation in the pancreas and its biomarkers: a systematic review and meta-analysis. Diabetes Metab. Res. Rev. 33:e2918. doi: 10.1002/dmrr.2918

PubMed Abstract | CrossRef Full Text | Google Scholar

Stancic, I., Gagrcin, M., Jovanovic, S., and Stancic, B. (2009). Reproductive activity of gilts with prolonged preinsemination anoestrus after hormonal treatment. Acta Vet. Beograd 59, 525–533. doi: 10.2298/Avb0906525s

CrossRef Full Text | Google Scholar

Stancic, I., Stancic, B., Bozic, A., Anderson, R., Harvey, R., and Gvozdic, D. (2011). Ovarian activity and uterus organometry in delayed puberty gilts. Theriogenology 76, 1022–1026. doi: 10.1016/j.theriogenology.2011.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Trout, W. E., Diekman, M. A., Parfet, J. R., and Moss, G. E. (1984). Pituitary responsiveness to GnRH, hypothalamic content of GnRH and pituitary LH and FSH concentrations immediately preceding puberty in gilts. J. Anim. Sci. 58, 1423–1431. doi: 10.2527/jas1984.5861423x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tummaruk, P., Tantasuparuk, W., Techakumphu, A., and Kunavongkrit, A. (2007). Age, body weight and backfat thickness at first observed oestrus in crossbred Landrace x Yorkshire gilts, seasonal variations and their influence on subsequence reproductive performance. Anim. Reprod. Sci. 99, 167–181. doi: 10.1016/j.anireprosci.2006.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Xin, W. S., Zhang, F., Yan, G. R., Xu, W. W., Xiao, S. J., Zhang, Z. Y., et al. (2018). A whole genome sequence association study for puberty in a large Duroc x Erhualian F2 population. Anim. Genet. 49, 29–35. doi: 10.1111/age.12623

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, C., Ye, J., Li, X., Gao, X., Zhang, K., Luo, L., et al. (2016). DNA methylation patterns in the hypothalamus of female pubertal goats. PLoS One 11:e0165327. doi: 10.1371/journal.pone.0165327

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, C., Ye, J., Liu, Y., Ding, J., Liu, H., Gao, X., et al. (2018). Methylation pattern variation between goats and rats during the onset of puberty. Reprod. Domest. Anim. 53, 793–800. doi: 10.1111/rda.13172

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G. C., Wang, L. G., Han, Y. Y., and He, Q. Y. (2012). cluster Profiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, X. L., Gao, N., Xing, Y., Zhang, H. B., Zhang, A. L., Liu, J., et al. (2016). Profiling the genome-wide DNA methylation pattern of porcine ovaries using reduced representation bisulfite sequencing. Sci. Rep. 6:22138. doi: 10.1038/srep22138

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, X. L., Zhang, Z., Li, B., Gao, N., Zhang, H., Sangild, P. T., et al. (2017a). Genome-wide DNA methylation analysis of the porcine hypothalamus-pituitary-ovary axis. Sci. Rep. 7:4277. doi: 10.1038/s41598-017-04603-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, X. L., Zhang, Z., Pan, R. Y., Gao, N., Deng, X., Li, B., et al. (2017b). Performances of different fragment sizes for reduced representation bisulfite sequencing in pigs. Biol. Proced. Online 19:5. doi: 10.1186/s12575-017-0054-5

PubMed Abstract | CrossRef Full Text | Google Scholar