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 fixed in 1% paraformaldehyde. The follicles from the ovary were stored in liquid nitrogen at -80 ℃.
Porcine granulosa cell culture
The ovaries used for pGCs culture were collected from a slaughterhouse in Guangzhou. Ovaries that lacked a corpus luteum but had clustered follicles full of follicular fluid were chosen. The ovaries were stored in phosphate-buffered saline (PBS) (Hyclone, Logan, UT, USA) containing 1% penicillin-streptomycin (Thermo Fisher, Waltham, CF, USA) on ice and were used for pGCs culture under sterile conditions. The pGCs culture method was modified from published procedures [23]. The pGCs were separated from 3- to 5-mm antral follicles with a 1 mL disposable sterile syringe (with a 4.5-gauge needle), washed twice with Dulbecco’s Modified Essential Medium (Hyclone, Logan, UT, USA) and cultured in medium supplemented with 10% foetal bovine serum (Gibco, USA) and 1% penicillin-streptomycin (Thermo Fisher, Waltham, CF, USA) in an incubator with 5% CO2 at 37 °C [24].
H3K4me3 and H3K27me3 inhibitors and agonists
H3K4me3 levels were decreased by 150 nM antagonist BCl-121 (#SML1817, Sigma, St. Louis, MO, USA) and increased by2 nM agonist PBIT (#SML1058, Sigma, St. Louis, MO, USA). The antagonist NC and agonist NC for H3K4me3 were 150 nM and 2 nM dimethyl sulfoxide (DMSO) (Sigma, St. Louis, MO, USA). H3K27me3 levels were decreased by 6 nM antagonist GSK126 (#1346574-57-9, MedChemExpress, NJ, USA) and increased by 2 nM agonist GSK-J4 (#1373423-53-0, MedChemExpress, NJ, USA) [25]. The antagonist NC and agonist NC for H3K27me3 were 6 nM and 2 nM DMSO (Sigma, St. Louis, MO, USA).
ChIP-Seq analysis of H3K4me3 and H3K27me3
The BF from each group were fixed in 1% paraformaldehyde, resuspended in lysis buffer, and disrupted with a Branson 250 Sonifier 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 purified 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 1500-bp fragments were size selected using AMPure XP beads. The final size of the library was confirmed by an Agilent 2100 Bioanalyzer. The samples were diluted to a final 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 filter, the clean reads were aligned to the pig genome (susScr3) using BOWTIE software (v2.1.0).
Peak calling. MACS v1.4 (Model-based Analysis of ChIP-seq) software was used to detect the peaks from ChIP-sEq. Statistically significant ChIP-enriched regions (peaks) detected by MACS were identified by comparison of IP to Input or to a Poisson background model using a p-value threshold fold of 10− 4. The proportions of different peak lengths were calculated and mapped.
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 modification in the PreOV and InOV groups, as well as to classify the modification as specific to H3K4me3 or H3K27me3 or as bivalent modification of both histones.
Differentially modified regions. The significantly differentially enriched regions between the two groups (InOV vs. PreOV) were identified by diffReps (cut-off: FC = 1.5, p-value = 0.01). Then, the differentially modified genes were classified 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: defined the same as UCSC RefSeq introns with a slightly different starting genomic coordinate of TSS + 2,000 bp; 4) exon region: defined 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 classified as one of the above 4 regions. The enrichment centres of the differentially modified genes were defined based on the classification of the peak. Various combinations of differential H3K4me3 or H3K27me3 modifications 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 significant enrichment of differentially expressed mRNAs. The P-value denotes the significance of the pathway; the lower the P-value is, the more significant the pathway (P-value ≤ 0.05). To better understand the biological function of genes with differentially modified promoters and of 88 genes with the other 4 classifications of differential modification, 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 identified genes were aligned with the human genome, and the human gene ID was used. Both GO terms and the KEGG pathways were filtered 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 fixing 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 flow 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/PI-negative 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 classified as MF (3–5 mm), BF (> 5 mm) or Graafian follicles (GF, with an ovulation hole). qRT-PCR analysis was performed for the top three genes with differential histone modification and genes with the 4 different types of H3K4me3 and H3K27me3 modifications, 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.