The chromatin mapping of H3K4me3 and H3K27me3 modication in ovaries during porcine puberty initiation

The biological structure and function of mammalian ovaries undergo important developmental changes during pre- puberty. Posttranslational modied histones dynamically regulate ovarian development, which also can respond to the onset of puberty with heritable changes in gene expression that are associated with nucleosome modications. This research used chromatin immunoprecipitation followed by sequencing (ChIP-seq) technology and bioinformatic analysis to conrm the histone-DNA interaction in pre-puberty ovarian follicles (PreOV) and in-puberty ovarian follicles (InOV). Then, employed qRT-PCR, and cell proliferation and apoptosis testing to further investigate the function of H3K4me3and H3K27me3 in gene transcription and cell growth. genome-wide modied prole of H3K4me3 and H3K27me3 showed that H3K27me3 modication was more active than H3K4me3 modication during the onset of puberty in porcine ovaries. H3K4me3 and H3K27me3 acted as activator and repressor to regulate gene transcription that are related to porcine puberty initiation process. Genes including FOXO3, BMPR1B, LEPTIN, RUNX1, involved in networks of female gamete generation pathway associating with ovarian follicle development, ovulation, ovulation cycle process and female gonad development, which showed the further direction to investigate the function of H3K4me3/H3K27me3-targets in puberty initiation.

Puberty is a critical physiological process in female pigs and involves the acquisition of reproduce capacity. The onset of puberty is marked by the rst estrus period and ovulation in response to boar stimulation [1]. In the pre-puberty, a period of follicle growth, atresia and degeneration, then the cycle repeats. After the onset of puberty, the functional development of the hypothalamic-pituitary-gonadal (HPG) axis improves, which ensure the higher gonadotropin-releasing hormone (GnRH) levels and prompt follicle maturation and ovulation [2]. The factors in uencing sow puberty include epigenetics, developmental programming, common diseases such as polycystic ovarian syndrome (PCOS), the cumulative effects of environmental factors, nutritional levels and boar stimulation [3]. However, internal and permanent changes in DNA dynamically regulate gene expression, and the programming, transcription and translation of the genetic machinery related to puberty remain incompletely elucidated. The onset of puberty is strongly determined by genetics [4], this study focused on the effect of epigenetics on sow puberty.
Epigenetics refers to a heritable phenomenon changes in gene expression without altering the DNA sequence; rather, the phenotype is changed without affecting genotype. There are many epigenetic phenomena, such as DNA methylation, histone posttranslational modi cation, genomic imprinting, maternal effects, gene silencing, and RNA editing [5][6][7][8]. Histone posttranslational modi cation including histone methylation [9], acetylation [10], ubiquitylation [11] and phosphorylation [12], is essential for gene transcription because it forms the chromatin on-off switch by affecting histone-histone, histone-DNA and histone-chaperone interactions [13]. Histone 3 lysine 4 tri-methylation (H3K4me3) and Histone 3 lysine 27 tri-methylation (H3K27me3) are typical chromatin marks of activating and repressing gene transcription respectively. In the promoter region, H3K4me3 loosens chromatin structure and is thus a positive regulator responsible for recruiting nucleosome remodeling enzymes and histone acetylases [14][15][16] to activate gene transcription, while H3K27me3 is the transcriptional inhibitory marks responsible for promoting a compact chromatin structure [17,18]. Epigenetic modi cations uctuate greatly from embryo implantation to early estrus in young females [19]. Mellon nds that a repressive histone modi cation H3K9me2 can elevate GnRH gene transcription [20]. H3K4me3 and H3K27me3 regulate development-related transcription factors (TFs) by forming bivalent domains (BDs), which can activate or inhibit the kisspeptin gene expression,and then affect the puberty initiation of female rats [21]. Generally, histone methylations are closely related to sexual maturity of rat, we speculate that H3K4me3 and H3K27me3 may be involved in the puberty initiation of sows.
In this current study, the main intent was to obtain genome-wide and time-speci c modi cation of H3K4me3 and H3K27me3 on the onset of puberty. Besides, the regulation of H3K4me3 and H3K27me3 in gene transcription was identify. Moreover, cell proliferation and apoptosis were analyzed after interfering with endogenous H3K4me3 or H3K27me3 in porcine granulosa cells. This research would expose the genome-wide H3K4me3-H3K27me3 bivalent modi cation during the porcine puberty.

Materials And Methods
Pre-puberty and in-puberty porcine ovaries The experimental groups conformed to a full sibling design. The six gilts came from 3 litters of full siblings and were separated into the pre-puberty ovarian (PreOV) group and the in-puberty ovarian (InOV) group. The three gilts in one group had one sibling gilt in the other group. Gilts in the PreOV group were slaughtered at day 180, before puberty initiation, whereas those in the InOV group were slaughtered on the day of puberty initiation. The follicles in one ovary were separated into three groups: big follicles (BF, > 5 mm), medium follicles (MF, 5 mm ≥ MF 3 mm) and small follicles (SF, ≥ 3 mm) [22]. The BF from the PreOV and InOV groups were xed in 1% paraformaldehyde. The follicles from the ovary were stored in liquid nitrogen at -80 ℃.

ChIP-Seq analysis of H3K4me3 and H3K27me3
The BF from each group were xed in 1% paraformaldehyde, resuspended in lysis buffer, and disrupted with a Branson 250 Soni er to obtain DNA in the size range of 200 to 1500 bp. Solubilized chromatin was diluted 10-fold in ChIP buffer; after removal of a control aliquot, the samples were incubated at 4 ℃ overnight with antibodies against H3K27me3 (#07-449, Millipore, German) and H3K4me3 (#ab8580, Abcam, Cambridge, UK). Immune complexes were precipitated with protein A-sepharose, washed sequentially with low-salt immune complex wash, LiCl immune complex wash, and TE, and then eluted in elution buffer. After cross-link reversal and proteinase K treatment, ChIP and control DNA samples were extracted with phenol-chloroform, precipitated in ethanol, treated with RNase and calf intestinal alkaline phosphatase, and puri ed with a MinElute Kit (Qiagen, Dusseldorf, German) [26].
The purity and concentration of DNA samples were determined with a Qubit® Fluorometer. DNA samples underwent end-repair, A tailing and adaptor ligation using a TruSeq Nano DNA Sample Prep Kit (#FC-121-4002, Illumina, San Diego, USA), following the manufacturer's instructions. Approximately 200-to 1500bp fragments were size selected using AMPure XP beads. The nal size of the library was con rmed by an Agilent 2100 Bioanalyzer. The samples were diluted to a nal concentration of 8 pM, and cluster generation was performed on the Illumina cBot using a HiSeq 3000/4000 PE Cluster Kit (#PE-410-1001, Illumina), following the manufacturer's instructions. Sequencing was performed on an Illumina HiSeq 4000 using a HiSeq 3000/4000 SBS Kit (300 cycles) (#FC-410-1003, Illumina), according to the manufacturer's introduction.

Data collection and analysis
Analysis of sequencing data. Sequence reads were generated by the Illumina HiSeq 4000 system, and image analysis and base calling were performed using off-line basecaller software (OLB V1.8.0). After passing the Illumina chastity quality lter, the clean reads were aligned to the pig genome (susScr3) using Peak annotation. The peaks in samples were annotated by the nearest gene (the nearest TSS to the centre of the peak region) using the newest UCSC RefSeq database (http://genome.ucsc.edu/). The online program BioVenn (http://www.biovenn.nl/index.php) was used to visualize the unique and common genes H3K4me3 or H3K27me3 modi cation in the PreOV and InOV groups, as well as to classify the modi cation as speci c to H3K4me3 or H3K27me3 or as bivalent modi cation of both histones.
Differentially modi ed regions. The signi cantly differentially enriched regions between the two groups (InOV vs. PreOV) were identi ed by diffReps (cut-off: FC = 1.5, p-value = 0.01). Then, the differentially modi ed genes were classi ed as upstream, promoter, exon, intron and intergenic enrichment of H3K4me3 or H3K27me3 as follows: 1) promoter region: 2,000 bp upstream and downstream of the TSS; 2) upstream region: > 2,000 bp but ≤ 20,000 bp upstream of the TSS; 3) intron region: de ned the same as UCSC RefSeq introns with a slightly different starting genomic coordinate of TSS + 2,000 bp; 4) exon region: de ned the same as UCSC RefSeq exons with a slightly different starting genomic coordinate of TSS + 2,000 bp; and 5) intergenic region: other genomic regions not classi ed as one of the above 4 regions. The enrichment centres of the differentially modi ed genes were de ned based on the classi cation of the peak. Various combinations of differential H3K4me3 or H3K27me3 modi cations were visualized via the online software Bioinformatics & Evolutionary Genomics (http://bioinformatics.psb.ugent.be/webtools/Venn/). The correlation of fold change and peak length for each differentially expressed gene was also included in the scope of investigation.
Pathway analysis. Based on the latest Kyoto Encyclopedia of Genes and Genomes database (KEGG) (https://www.kegg.jp), we performed a pathway analysis of the differentially expressed mRNAs. This analysis allows users to determine the biological pathway with signi cant enrichment of differentially expressed mRNAs. The P-value denotes the signi cance of the pathway; the lower the P-value is, the more signi cant the pathway (P-value ≤ 0.05). To better understand the biological function of genes with differentially modi ed promoters and of 88 genes with the other 4 classi cations of differential modi cation, GO and KEGG integrated analysis was conducted with ClueGO in Cytoscape (version 3.4.0) referring to the human GO and KEGG database [27,28]. Because the database of porcine genes is small, the identi ed genes were aligned with the human genome, and the human gene ID was used. Both GO terms and the KEGG pathways were ltered using a P-value < 0.05 as the selection criteria.
The biological function of H3K4me3 and H3K27me3 in porcine granulosa cells Granulosa cell proliferation assay. pGCs grown in 48-well plates were exposed to BCl-121, PBIT, GSK126, GSK-J4 and DMSO, and pGCs proliferation was analysed with EdU kit (#C10310-1, RiboBio, Guangzhou, China) according to the manufacturer's instructions and the literature [24]. In brief, 100 µL of 50 µM EdU solution was added to each well for 2 h. Then, the culture medium was discarded, 100 µL of cell xing solution (80% acetone) was added to each well for 30 min at room temperature, and the wells were washed twice with PBS. Next, 100 µL of penetrant (0.5% Triton X-100 (Sigma, St. Louis, MO, USA) in PBS) was added to permeabilize the cells, which were then washed once with PBS prior to the addition of 100 µL of 1 × Apollo Staining Solution for 30 min at room temperature in the dark. Subsequently, the staining solution was discarded, 100 µL of DAPI was added to each well for 30 min at room temperature in the dark, PBS was added, and pictures were taken under a microscope.
Granulosa cell apoptosis assay. pGCs grown in 6-well plates were exposed to BCl-121, PBIT, GSK126, GSK-J4 and DMSO. The cells were digested by trypsin, washed twice with ice-cold PBS and resuspended in 500 µL of binding buffer. Next, 1.25 µL of Annexin V-FITC (Sigma, St. Louis, MO, USA) was added, and the samples were incubated in the dark for 15 min at room temperature, followed by centrifugation at 1,000 × g for 5 min at room temperature to remove the supernatant. The cells were gently resuspended with 0.5 mL of pre-cooled 1 × solution, and 10 µL of propidium iodide (PI; 50 µg/mL) (MedChemExpress, NJ, USA) was added. Finally, the cells were analysed in a ow cytometer (Becton Dickinson Co., San Jose, CA, USA) using the FITC signal detector (FL1) and phycoerythrin emission signal detector (FL2). All experiments were performed at least three times. Cells in the lower right quadrant are annexin-positive/PInegative early apoptotic cells, and cells in the upper right quadrant are annexin-positive/PI-positive late apoptotic cells [24].
H3K4me3/H3K27me3-enriched gene analysis by qRT-PCR Different grade porcine follicles were used to measure gene expression. According to the size and characteristics, the follicles were classi ed as MF (3-5 mm), BF (> 5 mm) or Graa an follicles (GF, with an ovulation hole). qRT-PCR analysis was performed for the top three genes with differential histone modi cation and genes with the 4 different types of H3K4me3 and H3K27me3 modi cations, such as ENY2, SPATA24, NOBOX, MYC, LEPTIN, and BMPR1B. qRT-PCR primers are listed in Table S1. The relative gene expression was calculated with the 2 −ΔΔCt method.

H3K4me3 and H3K27me ChIP-Seq
The clean reads of H3K4me3 and H3K27me3 ChIP-Seq were aligned to the susScr3 reference genome to obtain the reads matched to the porcine genome. Peak calling (Table S2)  had H3K27me3 modi cations. Furthermore, there were 260 H3K4me3-modi ed genes unique to the PreOV group and 296 unique to the InOV group, and 3255 genes were modi ed in both phases (Fig. 3A). The number of phase-speci c H3K27me3-modi ed genes increased from 187 in the PreOV group to 1294 in the InOV group, and 1282 genes were common to these two groups (Fig. 3B).
The speci c and bivalent H3K4me3 and H3K27me3 modi cation analysis showed the H3K4me3-or H3K27me3-speci c modi cation of 2097 and 52 genes in the PreOV group (Fig. 4A) and of 1104 and 129 genes in the InOV group (Fig. 4B). There were 1418 and 2447 genes with bivalent modi cation in the PreOV and InOV groups, among which 172 and 1201 genes had speci c bivalent modi cations in PreOV and InOV, and 1246 genes had common bivalent modi cations in the two stages (Fig. 4C). The results presented almost double number of genes with bivalent modi cation of H3K4me3 and H3K27me3 in InOV vs. PreOV, was the results of recruiting H3K27me3 onto the genes with speci c H3K4me3 modi cation. In other word, some active genes in PreOV should be repressed in InOV to comply with relative biological processes in follicles during puberty initiation.

Differentially enriched genes in PreOV vs InOV
In total, H3K4me3 was signi cantly enriched at 946 genes and depleted at 916 genes upon comparison of InOV vs. PreOV. In addition, 1954 and 882 genes were signi cantly enriched (InOV vs. PreOV) or depleted (InOV vs. PreOV) in terms of H3K27me3 (Table S3). Differential peaks annotation showed that most of differentially peaks tittivated intergenic modi cation of H3K4me3 or H3K27me3 in PreOV and InOV, and the fewest peaks were promoter modi cation of H3K4me3 or H3K27me3 (Fig. 5).
The H3K4me3 and H3K27me3 modi cation in gene DNA sequences are usually in multi-sites (Fig. 6). The number of genes with unique differential H3K4me3 enrichment, H3K4me3 depletion, H3K27me3 enrichment and H3K27me3 depletion was 166, 131, 475 and 186, respectively. There were 88 genes even with four signi cantly differential modi ed sites. The number of genes with signi cantly differential H3K27me3 modi cation was much higher than that in other classi cations.
We found that more differentially enriched peaks with a higher enrichment fold change (FC) converged at 1-2 kb in all 4 cases of H3 posttranscriptional modi cation (Fig. 7A-D). The fold change generally ranged from 2-5 times. The differentially modi ed genes with a higher fold change caught our attention, such as ENY2, SPATA24, and MTRF in the H3K4me3 enriched group; LOC100520823, CPE and miR708 in the H3K4me3 depleted group; and COPS2, USP19 and MIR9811 in the H3K27me3 depleted group. To focus on 1-2 kb peaks, those longer than 2 kb were ignored, and more genes with shorter peaks and greater fold change were identi ed in the H3K4me3-enriched, H3K4me3-depleted and H3K27me3-depleted groups (Fig. 7E, F, H). Interestingly, the peaks in the H3K27me3-enriched group differed from those in the other three groups, with a more uniform distribution and a higher fold change focused for longer fragments, such as the top three genes ACTR2, CHBG and GPATCH2L (Fig. 7G).
qRT-PCR analysis of the top three differentially modi ed genes in mid-size follicles (3-5 mm, MF), bigsize follicles (5-6 mm, BF) and maturity follicles (with ovulation hole, GF) showed that the H3K4me3enriched genes ENY2 and SPATA24 were more highly expressed in development follicles (BF and GF) than in MF. CPE, which showed H3K4me3-deleted modi cation, had the highest expression in GF. The levels of COPS2, which showed increased H3K27me3 modi cation, were reduced in mature follicles. Genes depleted of H3K27me3, such as ACTR3, CHBG, and GPATCH2L showed a gradual increase in expression (Fig. 8A). The expression of genes related to female gamete generation, such as BMPR1B, FOXO3, and RUNX1, increased gradually, but LEPTIN and NOBOX expression decreased during the growth of follicles (Fig. 8B).

KEGG analyses of differentially enriched genes
The genes with a promoter signi cantly differentially enriched with H3K4me3, depleted of H3K4me3, enriched with H3K27me3 and depleted of H3K27me3 were involved in 3, 2, 12 and 0 biological pathways, respectively (Table 1). The citrate cycle, ether lipid metabolism, and ECM-receptor interaction were the three pathways involving the genes with a promoter differentially enriched in H3K4me3 in the InOV vs. PreOV comparison. The complement and coagulation cascades and Toll-like receptor signaling pathway were enriched in genes with a promoter differentially depleted of H3K4me3 in the InOV vs. PreOV comparison. Fatty acid elongation, microRNAs in cancer and prostate cancer were the top three most enriched pathways for genes with increased H3K27me3. Integrated analysis of the GO terms and KEGG pathways produced a simpler network of all genes with a promoter differentially modi ed by H3K4me3 or H3K27me3 (Fig. 9). Three prominent pathways were enriched: the BMP2-WNT4-FOXO1 pathway in human primary endometrial stromal cell differentiation, mitochondrial fatty acid beta-oxidation of unsaturated fatty acids and adenyl-nucleotide exchange factor activity. The BMP-WNT4-FOXO1 pathway in human primary endometrial stromal cell differentiation was related to the androgen receptor signaling pathway through DCN, FOXO1, CTNNB1, CREB1, and BAG1. Coincidentally, BAG1 linked adenyl-nucleotide exchange factor activity, ATPase regulator activity, and chaperone-mediated protein folding; this network was formed by BAG1, CRPEL2, FKBP3, and HSPH1. The last network of hydro-lyase activity, fatty acid beta-oxidation, fatty acid beta-oxidation and mitochondrial fatty acid beta-oxidation of unsaturated fatty acids was formed by HADHA, HADHB, UROS, DECR1, and ETFA.
The GO and KEGG con uence analysis of signi cantly differentially modi ed genes in terms of H3K4me3 and H3K27me3 modi cation showed that approximately 16% of the genes (14 of 88), including FOXO3, RUNX1, NOBOX, DMRT1, MYC, LEPTIN, RPL10, BDNF, TERT, FRZB, CD47, ACVR2, HSP90AB1 and BMPR1B, were involved in female gamete development associated with ovarian follicle development, ovulation, ovulation cycle process and female gonad development (Fig. 10). Other enriched pathways joined in female gamete generation pathway according to gene link, such as BMPR1B linking between positive regulation of bone and female gamete generation pathway, LEP linking between fatty acid transport pathway, and female gamete generation pathway, and RUNX1 linking between the hematopoietic stem cell differentiation pathway and female gamete generation pathway (Fig. 10).
The in uence of H3K4me3 or H3K27me3 on pGCs Compared to H3K4me3 antagonist NC, decreasing H3K4me3 signi cantly inhibited pGCs cell proliferation and induced cell apoptosis. Compared to H3K4me3 agonist NC, increasing H3K27me3 obviously promoted cell proliferation and anti-apoptosis (Fig. 11A, C). Compared to H3K27me3 antagonist NC/ agonist NC, both of decreasing or increasing H3K27me3 signi cantly restrained cell proliferation (Fig. 11B). however, there were the lowest cell proliferation in H3K27me3 agonist group. Compared to H3K27me3 antagonist NC, reducing H3K27me3 signi cantly repressed cell apoptosis, and accumulating H3K27me3 markedly induced cell apoptosis (Fig. 11D).

Discussion
H3K4me3 and H3K27me3 were trimethylated at sites of lysine residues 4 and 27 of histone 3, which performed activate and inactivate gene transcription to regulate organism development and growth [29,30]. It was con rmed that interfering endogenous H3K4me3 or H3K27me3 in uenced pGCs proliferation and apoptosis. Granulosa cells are critical components and functional units in ovarian follicles. The health of GCs determines follicle fate [31] and steroid biosynthesis and metabolism [32]. Synergism between granulosa cells and oocytes is the precursor to the initiation of puberty and maintains the normal oestrus cycle and reproduction.
The H3K4me3 and H3K27me3 modi ed pro les in the porcine follicular genome during the onset of puberty were determined in this study. The results of peak calling of H3K4me3 and H3K27me3 in PreOV and IOV indicated that H3K4me3 and H3K27me3 modi cations were dynamic during the transition from pre-puberty to puberty (Fig. 3). It seems like that the dynamical change of H3K27me3 modi cation was much drastic than H3K4me3 in PreOV and InOV, which comprised of the functional genes with different modi cation of H3K4me3 and H3K27me3 [15,18].
The analysis of speci c or bivalent modi cation in the PreOV and InOV groups revealed that H3K4me3 modi cation remained steady in the two phases but H3K4me3-speci c modi cation was reduced in puberty. Almost double number of genes with bivalent modi cation of H3K4me3 and H3K27me3 in InOV vs. PreOV, was the results of recruiting H3K27me3 onto the genes with speci c H3K4me3 modi cation. A total of 1246 of 2619 genes with bivalent modi cation were common to the PreOV and InOV groups (Fig. 4). In other word, some active genes in PreOV should be repressed in InOV under the regulation of H3K4me3 and H3K27me3 to comply with relative biological processes in follicles during puberty initiation.
Differential peaks annotation showed in Fig. 5 that most of differentially peaks showed intergenic modi cations of H3K4me3 or H3K27me3 in PreOV and InOV, and the fewest peaks were promoter modi cation of H3K4me3 or H3K27me3. The research indicated that the most HCNES, de ned as bivalent domains including a large H3K27me3-modi ed region encompassing a smaller H3K4me3modi ed region, encode key TFs that regulate embryonic development [29]. It pointed out that promoter modi cation of H3K4me3 and H3K27me3 related to DNA methylation [33] and chromatin compressing [34] to regulate gene transcription.
The different lengths of the H3K4me3-and H3K27me3-modi ed peaks were nteresting. Figure 7 presented the relationship between differentially enriched peak length and fold change (InOV vs. PreOV).
Higher fold changes were observed among H3K4me3-enriched peaks, H3K4me3-depleted peaks and H3K27me3-depleted peaks 1 kb in length, i.e., shorter peaks had a higher fold change in enrichment.
However, the fold change in H3K27me3-enriched peaks was dispersed over the 1-2 kb range, which indicates that most genes acquired H3K27me3 modi cation during the process of puberty. Zhang et al. noted that gene transcription and biological function were related to H3K4me3 peak width [35]. The results of qRT-PCR (Fig. 8) for the top three-fold change of H3K4me3-enriched, -deleted genes and H3K27me3-enriched, -deleted genes con rmed the function of H3K4me3-enriched modi cation was to upregulated gene transcription and H3K27me3-enriched modi cation was to repressed gene transcription [29].
According to the Venn of H3K4me3 and H3K27me3 differential modi cation (Fig. 6), it showed the regulation of histone posttranslational modi cation were in various sites. Not only different kinds of histone methylation [36], but also histone methylation and acetylation cooperated to regulated biology process [37]. The GO and KEGG conjoint analysis of the 88 genes with 4 kinds of differential modi cation of H3K4me3 and H3K27me3 participated in multifunction pathways relating to puberty. Approximately 16% of the genes (14 of 88), including FOXO3, RUNX1, NOBOX, DMRT1, MYC, LEPTIN, RPL10, BDNF, TERT, FRZB, CD47, ACVR2, HSP90AB1 and BMPR1B, were involved in female gamete development associated with ovarian follicle development, ovulation, ovulation cycle process and female gonad development (Fig. 10). FOXO3 is critical for oocyte reserve and female fertility [38,39]. RUXN1, a TF in periovulatory follicles, is induced by an LH surge to transactivate Ptgs2, which is essential for ovulation [40]. Oogenesis homeobox NOBOX encodes a TF that plays a role in folliculogenesis and the regulation of oocyte-speci c genes, among which the target Rsp2 is essential for follicular development [41].
DMRT1 is a testis-speci c TF associated with spermatogenesis and male infertility [34], and H3K4me3 and H3K27me3 modi cation might maximally inhibit the transcription of this gene in ovaries. The results clearly support future studies on all the above genes to investigate their function in the initiation of puberty in porcine ovaries.
The processes of female gamete generation and positive regulation of the bone mineralization pathway were correlation. The multifunctional nature of E2 leads to its involvement in different biological processes. The effect of E2 on the apoE receptor during osteoblast mineralization involves regulating the LDLR family, including LRP5, LRP6, LRP4, and ApoER2 [42]. In addition, BMPR1B is involved in fertility and follicular reserve and is strongly correlated with the FSH receptor in young women but not in old women. BMPR1B has a weak correlation with the LH receptor in young women but a strong correlation in old women [43]. ACVR2 knockdown in vivo inhibited its ability to maintain FSH homeostasis through the HPG axis in male mice [44].
Fatty acid biosynthesis, metabolism, and transport were enriched by 88 differentially expressed genes or differentially enriched promoters. We speculated that pathways related to fatty acid biosynthesis, metabolism and transport are essential for the initiation of puberty. Lipid droplets in steroidogenic tissues, such as ovaries enriched with PLIN1c, PLIN2, and PLIN3, have a cholesteryl ester (CE), which uses free cholesterol (FC) as the preferred cholesterol substrate for steroidogenesis, and HSL is the major neutral cholesterol esterase mediating the conversion of CE to FC [45]. The related genes include LEPTIN, FABP3, CROT, STARD4, PLIN2, FRZB, SULTIE1, KLF5, DECR1, ETFA, HADHA, and HADHB. LEPTIN is primarily expressed in fat and acts as a bridge between female gamete generation and the fatty acid transport pathway, which has been identi ed to be related to puberty [46]. STARD4 participates in both "Regulation of steroid biosynthesis process" and "Regulation of steroid metabolic process" in granulosa cells and is important for oogenesis, folliculogenesis, ovulation, and pregnancy [32]. SULT1E1 is an LHinducible gene [40] and is involves in the inactivation of oestrogen after exogenous GnRH stimulation [47]. Above all, these results broaden the role of fatty acid pathways in the onset of puberty, including the regulation of steroidogenesis, steroid metabolism, and oestrogen activity, and the associated gene expression was regulated at the chromatin level by histone posttranscriptional modi cation, i.e., H3K4me3 and H3K27me3.
Ovarian GCs are critical components and functional units in ovarian follicles. Synergism between granulosa cells and oocytes is the precursor to the initiation of puberty and maintains the normal oestrus cycle and reproduction. It was con rmed that interfering endogenous H3K4me3 or H3K27me3 in uenced pGCs proliferation and apoptosis (Fig. 11). The health of GCs determines follicle fate [31] and steroid biosynthesis and metabolism [32].

Conclusions
The genome-wide modi ed pro le of H3K4me3 and H3K27me3 showed that H3K27me3 modi cation was more active than H3K4me3 modi cation during the onset of puberty in porcine ovaries. H3K4me3 and H3K27me3 acted as activator and repressor to regulate gene transcription that are related to porcine puberty initiation process. Genes including FOXO3, BMPR1B, LEPTIN, RUNX1, involved in networks of female gamete generation pathway associating with ovarian follicle development, ovulation, ovulation cycle process and female gonad development, which showed the further direction to investigate the function of H3K4me3/H3K27me3-targets in puberty initiation.      qRT-PCR analysis of genes with signi cant differentially H3K4me3 or H3K27me3, or both modi cations.
Some of genes were the top three differentially modi ed genes in the H3K4me3-enriched, H3K4me3depleted, H3K27me3-enriched and H3K27me3-depleted groups, including ENY2, SPATA24, CPE, COPS2, ACTR3, CHBG, GPATCH2L. Some of genes were associated with the 4 classi cations of differential H3K4me3 and H3K27me3 modi cation in follicles of different grades. MF, BF and GF indicated 3-5 mm in diameter medium follicles, 5-8 mm big follicles and > 8 mm graa an follicles, respectively.

Figure 9
Network analysis of GO terms and KEGG pathways enriched in genes with differential H3K4me3 or H3K27me3 modi cation in promoter region. The three coloured nodes correspond to three pathways: BMP2-WNT4-FOXO1 pathway in human primary endometrial stromal cell differentiation, adenylnucleotide exchange factor activity and mitochondrial fatty acid beta-oxidation of unsaturated fatty acids. The 7 GO terms in purple were enriched in the genes in red.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.