Next Article in Journal
Genome-Wide Identification and Expression Analysis of the DMP and MTL Genes in Sweetpotato (Ipomoea batatas L.)
Previous Article in Journal
A New Cloud-Native Tool for Pharmacogenetic Analysis
Previous Article in Special Issue
Comprehensive Analysis of miRNA and mRNA Expression Profiles during Muscle Development of the Longissimus Dorsi Muscle in Gannan Yaks and Jeryaks
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Genome-Wide DNA Methylation Analysis and Functional Validation of Litter Size Traits in Jining Grey Goats

1
College of Animal Science, Xinjiang Agricultural University, Urumqi 830052, China
2
Institute of Animal Science and Veterinary Medicine, Shandong Academy of Agricultural Sciences, Jinan 250100, China
*
Authors to whom correspondence should be addressed.
Genes 2024, 15(3), 353; https://doi.org/10.3390/genes15030353
Submission received: 10 February 2024 / Revised: 7 March 2024 / Accepted: 7 March 2024 / Published: 12 March 2024
(This article belongs to the Special Issue Livestock: Genomics, Genetics and Breeding)

Abstract

:
DNA methylation (DNAm) is associated with the reproductive system. However, the genetic mechanism through which DNAm regulates gene expression and thus affects litter size in goats is unclear. Therefore, in the present work, genome-wide DNAm profiles of HP and LP Jining Grey goat ovary tissues were comprehensively analyzed via WGBS, and RNA-Seq data were combined to identify candidate genes associated with litter size traits in the Jining Grey goat. Finally, BSP and RT-qPCR were used to verify the sequencing results of the key genes. Notably, the DNMT genes were downregulated at the expression level in the HP group. Both groups exhibited comparable levels of methylation. A total of 976 differentially methylated regions (DMRs) (973 DMRs for CG and 3 DMRs for CHG) and 310 differentially methylated genes (DMGs) were identified in this study. Through integration of WGBS and RNA-Seq data, we identified 59 differentially methylated and differentially expressed genes (DEGs) and ultimately screened 8 key DMGs (9 DMRS) associated with litter size traits in Jining Grey goats (SERPINB2: chr24_62258801_62259000, NDRG4: chr18_27599201_27599400, CFAP43: chr26_27046601_27046800, LRP1B. chr2_79720201_79720400, EPHA6: chr1_40088601_40088800, TTC29: chr17_59385801_59386000, PDE11A: chr2_117418601_117418800 and PGF: chr10_ 16913801_16914000 and chr10_16916401_16916600). In summary, our research comprehensively analyzed the genome-wide DNAm profiles of HP and LP Jining Grey goat ovary tissues. The data findings suggest that DNAm in goat ovaries may play an important role in determining litter size.

1. Introduction

The Jining Grey goat, which originated in Shandong Province, China, is known for its high fecundity. Additionally, the Jining Grey goat represents the most famous local goat breed in the world and is an ideal animal model for studying the reproductive performance of goats [1]. Reproductive performance directly affects the economic benefits of goat farming, and the ovary, as an organ for oocyte production and reproductive hormone secretion, is directly related to litter size traits [2,3]. Existing studies have identified several genes associated with litter size traits in goats, such as DNAH1 [4], CFAP43 [5], CTSS [6], EPHA6 [7], and POU1F1 [8]. However, the key molecular regulatory mechanisms remain unclear.DNA methylation (DNAm) is a type of epigenetic modification that modulates gene expression. With increasing research on epigenetic modifications, it has been shown that DNAm is associated with the reproductive system. Genome-wide DNAm profiles used to study ovarian tissues from Jining Grey goats [9], Hu sheep [10], and pigs [11] with different litter sizes. In addition, Yao X et al. [2] performed DNAm and gene expression analyses via WGBS and RNA-seq techniques and reported that DNAm affects reproduction by influencing the expression of ITGB2 and LAPTM4B in Hu sheep’s ovaries. Stefano Frattini et al. [12] performed DNAm and gene expression analyses via methyl-CpG binding domain protein sequencing (MBD-seq) and RNA-seq combined with analysis and showed that specific differentially methylated genes in the ovary and hypothalamus of Saanen goats correlated with gene expression levels. However, studies on genome-wide DNAm and alterations between transcriptomes have not been conducted for HP or LP Jining Grey goats. Therefore, in the present study, WGBS and RNA-seq were combined to investigate the epigenetic regulatory mechanisms of DNAm in the ovaries of HPs and LP Jining Grey goats.
This research aimed to reveal the epigenetic regulation of DNAm at different levels of prolificacy in Jining Grey goats. We used WGBS to study the DNAm profiles of ovarian DNA from HPs and LP Jining Grey goats during the estrous period. Combining these findings with previous RNA-seq data, we explored the effects of DNAm on gene expression through joint analysis of DNAm and gene expression and screened key candidate genes that may be related to litter size traits in goats. In addition, our research provides a reference for exploring the potential regulatory mechanism of litter size traits in Jining Grey goats at the DNAm level.

2. Materials and Methods

2.1. Animals and Sample Collection

The experimental protocols and animal ethics were approved by the Ethics Committee of the Institute of Animal Science and Veterinary Medicine, Shandong Academy of Agricultural Sciences (IASVM-2022-003). The experimental animals were obtained from Jining Grey goats from the National Jining Grey Goat Breeding Farm (Jiaxiang, Jining, Shandong, China). Similarly, eight ewes (3–4 years of age) with the same feeding conditions who were healthy and disease-free were randomly selected and classified as exhibiting high prolificacy (HP, litter size ≥ 3 in three consecutive litters, n = 4) or low prolificacy (LP, litter size = 1 in three consecutive litters, n = 4). First, we used a vaginally embedded CIDR pessary to treat experimental animals for simultaneous estrus. The date of CIDR pessary implantation was recorded as Day 0 and 14 days after implantation in the vagina. The pessary was withdrawn, and 300 IU pregnant mare serum gonadotropin (PMSG) was injected intramuscularly into the neck to determine the estrus status of ewes using the method used to test the estrus of ewes and vaginal examination. The ewes were slaughtered within 12 h of estrus, and the samples of the ipsilateral ovary were collected immediately and kept at −80 °C for further experiments.

2.2. DNA and RNA Extractions

Total RNA was extracted from ovarian tissues using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), agarose gel electrophoresis was employed to determine RNA integrity, and the RNA’s purity was assessed using a NanoDropTM One spectrophotometer (Thermo Scientific, Waltham, MA, USA). The sample quality was satisfactory for the next step of the experiment. DNA samples were extracted from ovarian tissues using a TIANamp Genomic DNA Kit (Cat.#DP304; TianGen, Beijing, China). The purity of the RNA was determined using a NanoDropTM One spectrophotometer. Finally, DNA quality was determined by 2% gel electrophoresis, and the DNA was kept at −80 °C.

2.3. Library Preparation

Qualified genomic DNA was sheared into 200–300 bp fragments using an ultrasonic method. After DNA end repair, adenylation of the 3′ end and ligation of sequencing junctions were performed. The fragments were bisulfite-treated using an EZ ZYMO DNAm-Gold Kit, further desalted and recovered by gum cutting. Then, library fragment size selection was performed, followed by PCR amplification, to construct a sequencing library.

2.4. WGBS and Identification of DMRs

Qualified libraries were subjected to whole-genome bisulfite sequencing (WGBS) using the MGISEQ2000 platform (MGI Tech Co., Ltd., Shenzhen, China), and the raw reads were filtered through SOAPnuke V1.5.6. Reads containing splice contamination were first removed, after which reads with an N content greater than 1% or low-quality bases greater than 40% were filtered; the final filtered data were called clean reads, after which the Q30 and GC contents were computed.
Alignment of clean reads was conducted based on the goat reference genome (GCF_001704415. 1_ARS1) using Bismark v0.23.0 [13] software. Statistical metrics, including bisulfite conversion and comparison rates, were also computed for each sample. The deduplicate_bismark program was then used to eliminate duplication reads, enabling the bismark_methylation_extractor and coverage2cytosine programs to isolate the unit point methylation information. The parameters were as follows:
bismark_methylation_extractor --samtools_path --bedGraph --CX --gzip --parallel --buffer_size
coverage2cytosine --split_by_chromosome --genome_folder -CX --gzip
The DNAm level was calculated as follows [14,15]:
R m a v e r a g e = N m a l l N m a l l + N n m a l l 100 %
where N m a l l represents the number of methylated cytosine residues covering the site, and N n m a l l represents the number of nonmethylated reads.
This study used the R package DMRcaller [16] to calculate and analyze the DMRs. For different contexts, different minProportionDifference parameters were used. The default CG was 0.2, the CHG was 0.1, and the CHH was 0.05.

2.5. DMR-Associated Genes (DMGs) and Functional Enrichment Analysis

BEDTools [17] was employed to determine DMR-related genes or other genome elements according to DMR genome position. In this study, genes or gene elements that overlapped with DMRs by at least 1 bp were considered DMR-associated genes or DMR-related elements. Furthermore, the DMG was subjected to GO enrichment analysis and KEGG functional annotation using DAVID “https://david.ncifcrf.gov/ (accessed on 29 September 2023)” and was considered to be significantly enriched at p < 0.05.

2.6. RT-qPCR Analysis

cDNA synthesis was conducted with the PrimeScript RT Reagent Kit with gDNA Eraser (Perfect Real Time) (Takara, Dalian, China). The following cDNA was obtained by reverse transcription: 2 µL 5× gDNA Eraser Buffer, 1 µL gDNA Eraser, 1 µg total RNA and RNase free dH2O up to 10 µL. The mixture was subsequently placed in a PCR machine for the reaction, which was performed at 42 °C for 2 min. Afterward, the reaction mixture was removed, and 4 µL RT Primer Mix, 4 µL 5× PrimeScript Buffer 2, 1 µL PrimeScript RT Enzyme Mix 1, and 1 µL RNase free dH2O were added. The mixture was subsequently placed in the PCR instrument again for the reaction, after which the reaction protocol was 37 °C for 15 min, 85 °C for 5 s, and 4 °C, the cDNA was then diluted to 10 ng/µL and stored at −20 °C.
Primers pairs were designed with Primer Premier 5.0 to validate the expression levels of key genes by real-time quantitative PCR (RT-qPCR) using GAPDH as an internal reference gene, All primers were diluted to 10 µM. RT-qPCR was carried out on a Roche LightCycler®480 II with 10 µL TB Green Premix Ex Taq II (Tli RNaseH Plus) (2× Conc), 0.8 µL upstream and downstream primers each, 2 µL cDNA template and 6.4 µL RNase-free dH2O. The reaction procedure was as follows: 95 °C for 3 min, followed by 40 cycles of 95 °C for 5 s, 60 °C for 20 s, and extension at 72 °C for 5 s and 40 cycles. Three replicates were determined for each sample, and the change in target gene expression was determined via the 2−ΔΔCt method. The primer sequences are summarized in Table 1.

2.7. Bisulfite Sequencing PCR (BSP)

BSP was conducted to assess the DNAm level of the DMR region of key candidate genes. The genomic DNA was modified with sodium bisulfite using an EZ DNA Methylation-GoldTM Kit (Cat.#D5005; Zymo Research, Los Angeles, CA, USA) with a uniform dilution of 200 ng/µL and stored at −20 °C. We used MethPrimer 2.0 “http://www.urogene.org/methprimer/ (accessed on 29 October 2023)” online analysis software to design BSP primers for the DMR region, which were then delivered to Sangon Biotech (Shanghai, China) for synthesis. The DMR region was further amplified by PCR using bisulfite-converted DNA as a template in 12.5 µL Zymo TaqTM PreMix (Cat. #E2003; Zymo Research, Los Angeles, CA, USA), 1 µL upstream and downstream primers each, 2 µL gDNA, and 8.5 µL RNase-free dH2O. The reaction procedure was as follows: 95 °C for 10 min, followed by 40 cycles of 95 °C for 30 s, 52 °C for 30 s, 72 °C for 30 s, and 72 °C for 10 min. The PCR product was evaluated via 2% gel electrophoresis, and the qualified PCR products were placed in a PCR product purification kit (Cat.#BSC02M1A. Bioer Technology, Hangzhou, China) for purification, followed by ligation and cloning, and inserted into the pMD19-T vector (Cat. #6013; Takara, Osaka, Japan). Subsequently, the ligated product was transferred into DH5α receptor cells (Cat. G6016, AngYuBio, Shanghai, China), further coated on LA solid medium and cultured at 37 °C for 16 h. Afterwards, three white colonies were randomly picked from the medium, inoculated in LB liquid medium, and cultured for 3 h at 220 rpm/h. Afterward, PCR amplification was carried out in 10 µL 2× Taq PCR StarMix (Dye) (Cat# A012, GenStar, Beijing, China), 0.8 µL upstream and downstream primers each, 2 µL bacterial broth, and 6.4 µL RNase-free dH2O, with the conditions of the reaction described above. The PCR product was tested again via 2% gel electrophoresis, and the samples that passed quality control were then subsequently delivered to Sangon Biotech (Shanghai, China) for sequencing. The primer sequences are summarized in Table 2.

2.8. Statistical Analysis

The experimental data were preprocessed using Microsoft Excel, and independent sample t tests were conducted with SPSS 25.0. All values are shown as the mean ± SE, with a p value of 0.05 as the threshold.

3. Results

3.1. Genome-Wide Methylation Mapping

In this study, the ovaries of HP and LP Jining Grey goats were used as tissue test samples, and whole-genome DNAm profiles of the ovaries of HP and LP Jining Grey goats were constructed via WGBS sequencing and RNA-seq data to screen key candidate genes. Finally, the methylation and expression levels were verified via BSP and RT-qPCR (Figure 1A). After the sequencing data were filtered, an average of approximately 1.16 billion clean reads were produced for each sample. The conversion rate after bisulfite treatment was >99%, the Q30 of the clean reads ranged between 90.16% and 91.94%, and the Unique mapping ratio accounted for 74.17~79.17% of the total reads. The percentage of Total_mC to all C sites was 2.98~3.66%. The sequencing data’s quality is depicted in Table 3. Bismark statistics of methylation sites showed that the genome-wide mC levels of 98.539 ± 0.010% for CG, 0.362 ± 0.002% for CHG, and 1.099 ± 0.009% for CHH in HP. In LP, the levels were 98.511 ± 0.035% for CG, 0.368 ± 0.009% for CHG, and 1.121 ± 0.028% for CHH. There was no obvious difference between the groups (Figure 1B), and there were significantly more methylated cytosine sites in CG than in CHG and CHH (where H is T, C or A). The DNAm levels of diverse genomic elements in the HP and LP groups of Jining Gray goats showed that there were diverse DNAm levels in the three contexts. In the CG, the intergenic exhibited the highest DNAm level, and the exon had the lowest methylation level. However, there was almost no methylation in the CHH or CHG contexts (Figure 1C). The gene body region exhibited high methylation, followed by the UP3K and Down3K regions, and the lowest methylation level was found near the transcription start site (TSS) locus (Figure 1D).

3.2. Differential Methylation Region (DMR) Identification

In this study, a total of 973 CG DMRs (hyper: 460; hypo: 513) and 3 CHG DMRs (hypo: 3) were identified between the HP and LP groups of the Jining Grey goats, and no DMRs were identified in the CHH population (Figure 2A). The identified DMRs were compared to the goat reference genome (GCF_001704415.1_ARS1), and a total of 310 differentially methylated genes (DMGs) were identified, 309 of which were found in the CG (hyper:125; hypo:192), while only one differentially methylated gene was found in the CHG (Figure 2B). The majority of DMGs belonged to the CG type (>99%), and this study evaluated CG methylation for the next analysis. The distribution of DNAm levels in DMRs showed that the mean DNAm levels ranged from 0.51 to0.53, with no significant differences between groups (Figure 2C). The DMR gene region showed that DMRs were mostly distributed in the intergenic region, with only a small number of DMRs distributed in the DOWN3K, exon, and UP3K regions (Figure 2D).

3.3. DMG Functional Enrichment Analysis

To determine the regulatory mechanism of genes related to litter size in Jining Grey goats, we conducted GO and KEGG functional enrichment analyses of 309 CG-type DMGs detected in the ovarian tissues of HP and LP Jining Grey goats in the present study. GO analysis showed that the DMGs were enriched mainly in biological processes, including endochondral bone growth, cell adhesion, chondrocyte proliferation, and collagen fibrils. In addition, the GO term “cellular component” was enriched mainly in the integral component of the membrane, integral component of the plasma membrane, and collagen trimer, whereas the GO term “molecular function” was enriched mainly in the GO term “molecular function” and was predominantly enriched in ATP binding, metalloaminopeptidase activity, and polypeptide N-acetylgalactosaminyltransferase activity (Figure 3A). KEGG analysis demonstrated that the DMGs were predominantly enriched in the ABC transporter, cAMP signaling pathway, antifolate resistance, and bile secretion pathway (Figure 3B).

3.4. Screening for DMG Associated with Litter Size Traits

This study combined WGBS and RNA-Seq data to identify 59 differentially methylated and expressed genes (Figure 4A). Among these genes, 32 DMG methylation levels were positively correlated with expression levels and 29 DMGs were negatively correlated with expression levels (Figure 4B); the details are shown in Table 4. Finally, we identified eight key genes that might be related to litter size traits in Jining Grey goats through previous studies, and these have been previously reported in other studies. Eight key candidate genes were associated with SERPINB2, NDRG4, CFAP43, LRP1B, EPHA6, TTC29, PDE11A, and PGF. One DMR was identified in each of the SERPINB2 and NDRG4 genes located in the UP3K region, one DMR was identified in the HP group compared to the LP group, one DMR was identified in the CFAP43 region located in the intro, and one DMR was upregulated. Another DMR was identified in the LRP1B, EPHA6, TTC29, and PDE11A regions located in the intro, and the methylation level was downregulated. Two DMRs were identified in the PGF region located in the exon and intro; and the methylation level was downregulated (Table 5).

3.5. Key Candidate Gene Validation

The robustness of WGBS sequencing data was verified by BSP experiments, in which four DMR fragments were randomly selected from key candidate genes (NDRG4: chr18:27599410-27599622; SERPINB2: chr24:62258787-62258914; LRP1B: chr2:79444158-79444322; and EPHA6: chr1:39785850-39786019) and subjected to bisulfite sequencing PCR (BSP) (Figure 5A–D). The findings showed that the DNAm levels within the DMR in the HP subgroup were lower than those in the LP subgroup for NDRG4, SERPINB2, LRP1B, and EPHA6. In summary, the results of the BSP experiments were in good agreement with the WGBS sequencing results, and the sequencing results were robust. The mRNA levels of eight key candidate genes were validated through RT-qPCR, and those of DNMTs (DNMT1, DNMT3A, and DNMT3B) were detected. Notably, the expression levels of DNMTs (DNMT1, DNMT3A, and DNMT3B) in the HP group were lower compared to LP group (Figure 5E). The expression levels of the SERPINB2 and LRP1B genes were upregulated, while the expression levels of the PGF, EPHA6, NDRG4, CFAP43, TTC29 and PDE11A genes were downregulated (Figure 5F). In summary, DNMT expression was downregulated in the HP group, and the data from the BSP and RT-qPCR experiments were in good agreement with the sequencing results, meaning the findings are reliable.

4. Discussion

DNAm constitutes a prominent aspect of epigenetic regulation and plays a vital role in modulating gene expression [31]. Currently, genome-wide methylation profiles of pig [11,32], goat [9], and sheep [10] ovaries have been intensively studied. While genome-wide DNAm has been described in the ovaries of Jining Grey goats [9] and Saanen goats [12], there are still many unknowns about the epigenetic mechanisms of DNAm on litter size traits in goats, which have not been detected in studies on genome-wide DNAm and transcriptome-wide changes in litter size traits in goats. First, we assessed the genome-wide DNAm profile of the ovaries of Jining Grey goats using WGBS technology and found that methylated cytosine sites accounted for 2.98–3.66% of all cytosine sites and that the proportion of CG-type methylated sites in the genome was much greater than that of CHG and CHH, which is in agreement with the findings of other species [2,10,33,34,35,36]. It has been reported that the overall methylation level of the UP2K region in the ovary of Jining Grey goats is low and flat, the DNAm level increases after TSS, and the gene body region shows a high methylation level, which declines until after TTS [9]. This finding aligns closely with the results obtained in this study, in which the methylation level of cytosine sites near the TSS was significantly lower than that of the UP3K and DOWN3K regions in the ovaries of the Jining Grey goat, whereas the gene body region likewise showed a high degree of methylation, a result that is consistent with the DNAm patterns previously observed in the ovaries of Hu sheep [10] and pig [36] ovaries. A total of 976 DMRs and 310 DMGs were identified, and no obvious difference was observed in the DNAm level of DMRs between the two groups; most of the DMRs were detected within the intergenic region, and a small percentage of them were distributed in the DOWN3K, exon, and UP3K regions, a result that is similar to that of previous studies on the ovaries of Hu sheep [10] and pigs [36]. Pathway analysis revealed that DMGs were involved in cell division and proliferation through the p38MAPK and cAMP signaling pathways. No relationship was observed between tissue functions and specific pathways in the ovary, which may be attributed to limited gene annotation in the goat genome.
DNA methyltransferase (DNMT) genes are key genes involved in the construction of mammalian DNAm patterns [37]. Previous studies have shown that DNMT3a may have wide-ranging effects on reproductive functions [38], and we found that the mRNA expression of DNMTs was lower in the HP group compared to the LP group; therefore, we hypothesized that DNMTs could mediate the methylation levels of genes related to litter size traits in Jining Grey goats.
A total of eight key candidate genes that may be associated with litter size traits in Jining Grey goats were identified in this study (SERPINB2, NDRG4, CFAP43, LRP1B, EPHA6, TTC29, PDE11A, and PGF). Hypermethylation of promoters is usually associated with gene repression, and the DMRs of SERPINB2 and NDRG4 were found to be located within the promoter region in our research. SERPINB2 is a serpin family B member. Prior reports have shown that in porcine ovarian granulosa cells after lipopolysaccharide (LPS) stimulation, TNFα plays an essential role in cell proliferation by activating the Erk1/2 pathway, which mediates the upregulation of SERPINE1 and SERPINB2 expression [22]. In the present study, DMR hypomethylation in the promoter region of SERPINB2 promoted gene expression, which may have affected litter size traits in Jining Grey goats. The NDRG4 gene has been reported to affect the "adhesive switch" by producing epigenetic silencing through promoter hypermethylation, and this gene could be a potential mechanistic biomarker for breast cancer [27]. The NDRG4 gene, part of the N-myc downregulated gene family and belonging to the alpha/beta hydrolase superfamily, has garnered attention in recent studies for its potential role in estrogen-regulated embryo implantation in mice [26]. However, recent studies have shown that methylation of promoter regions can also activate gene expression [39]. Hypermethylation of the promoter region of the FoxA2 gene promotes its expression and thus regulates the development of endoderm in cells [40]. Herein, similar results were found for the DNAm level and expression level of DMR in NDRG4. Moreover, GO functional analysis indicated that NDRG4 was enriched in the vesicle docking term, and we further speculate that NDRG4 might affect its litter size traits through vesicle docking. The effect of DNAm in the gene body region on the regulation of gene expression is more complicated [12,41]. The CFAP43 gene, alternatively recognized as WDR96, exhibits widespread expression across gonadal tissues, and early studies have shown that mutations in the CFAP43 gene cause infertility in Trypanosoma and humans [19]. Meanwhile, indel mutations in CFAP43 in white cashmere goats from northern Shaanbei can affect their litter size [5] and body size traits [18]. A GWAS showed the LRP1B gene to be a key candidate gene for teat number in Luzhong meat sheep [20], following a study that showed a correlation between litter size and teat number [42]. The EPHA6 gene is also known as ephrin type-A receptor 6, and previous studies have shown that EPHA6 is a key candidate gene for prolificacy traits in Polish goats [7]. TTC29 is involved in cilium movement and cilium organization. Chunyu Liu et al. [43] showed that biallelic mutations in TTC29 can cause male infertility as well as asthenoteratospermia in humans and mice. The PDE11A gene encodes the PDE protein superfamily member and is associated with testicular germ cell tumors (TGCTs) [28,29]. Smoking can alter the overall methylation status of PDE11A in human spermatozoa and affect its transcription, which in turn affects male fertility [30]. Placental growth factor (PGF) enables growth factor activity and is closely related to fetal development. PGF has been identified as an important vascular growth factor in the placenta and contributes to placental and fetal development [24]. Apart from its roles in angiogenesis and modulating the growth and migration of various cell types, the presence of a CPG island hypermethylated on PGF exon 7 in the fetal growth restriction (FGR) placenta leads to the downregulation of PGF expression, which further regulates trophoblast cell proliferation and migration, affecting placental development and function [23]. Recent research has demonstrated that polymorphisms in the gene encoding placental growth factor (PGF) are strongly associated with stillbirths as well as calving ease in German dairy cows [25]. Herein, two DMRs were identified on the PGF, both of which showed hypomethylation and downregulated mRNA expression in the HP group, while GO analysis demonstrated that the PGF was enriched in the positive regulation of the cell division process. However, little has been reported on how DNAm of CFAP43, LRP1B, EPHA6, TTC29, PDE11A and PGF affects the litter size traits of female animals.
In this study, hypomethylation of the promoter region induced SERPINB2 expression, while repressing NDRG4 expression; hypermethylation within the gene body region repressed CFAP43 expression; and hypomethylation within the gene body region promoted LRP1B expression and repressed EPHA6, TTC29, PDE11A and PGF expression. Therefore, in this study, we identified eight genes, SERPINB2, NDRG4, CFAP43, LRP1B, EPHA6, TTC29, PDE11A, and PGF, as key candidate genes regulating litter size traits in Jining Grey goats, and the DMRs on these genes may directly or indirectly regulate phenotypic differences between HP and LP Jining Grey goats. However, the epigenetic mechanism through which these genes regulate litter size traits in Jining Grey goats still needs further study.

5. Conclusions

In this work, genome-wide DNAm profiles of ovarian tissues from HPs and LP Jining Grey goats were constructed, and it was concluded that DNAm modification is a potential factor affecting Jining Grey goats. Further studies showed that nine methylation modification sites (SERPINB2: chr24_62258801_62259000, NDRG4: chr18_27599201_27599400, CFAP43: chr26_27046601_27046800, LRP1B: chr2_79720201_79720400, EPHA6: chr1_40088601_40088800, TTC29: chr17_59385801_59386000, PDE11A: chr2_117418601_117418800 and PGF: chr10_ 16913801_16914000; chr10_16916401_16916600) may regulate lambing traits in Jining Grey goats by influencing gene expression to further affect ovarian development. A better understanding of the influence of epigenetics on litter size traits is likely to aid increased animal prolificacy.

Author Contributions

Conceptualization, C.Y., K.T. and X.H.; methodology, C.Y., J.H., J.M., G.L., Y.R., C.W. and G.Z.; software, C.Y. and J.H.; validation, C.Y.; formal analysis, C.Y. and K.T.; writing—original draft preparation, C.Y.; writing—review and editing, J.H., K.T. and X.H. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Agricultural Scientific and Technological Innovation Project of Shandong Academy of Agricultural Sciences (CXGC2023F10), the study of the Molecular mechanism of hair follicle development in high-quality fine-wool sheep revealed by single-cell transcriptome (32372851), and the Creation and Application of Healthy and Efficient Feed for Yak and Tibetan Sheep in Haibei Prefecture (YDZX2023012).

Institutional Review Board Statement

All animal experiments were approved by the Laboratory Animal Ethics Committee, Institute of Animal Science and Veterinary Medicine, Shandong Academy of Agricultural Sciences (Jinan, China, IASVM-2022-003).

Informed Consent Statement

Not applicable.

Data Availability Statement

All WGBS and RNA-seq data generated in this study were submitted to the NCBI SRA database under BioProject: NO.PRJNA1068426 (WGBS) and NO.PRJNA1068677 (RNA-seq).

Acknowledgments

Particular thanks are due to the veterinarians who helped collect the samples.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Wang, J.-J.; Zhang, T.; Chen, Q.-M.; Zhang, R.-Q.; Li, L.; Cheng, S.-F.; Shen, W.; Lei, C.-Z. Genomic Signatures of Selection Associated with Litter Size Trait in Jining Gray Goat. Front. Genet. 2020, 11, 286. [Google Scholar] [CrossRef]
  2. Yao, X.; Li, F.; Wei, Z.; Ei-Samahy, M.A.; Feng, X.; Yang, F.; Wang, F. Integrative Genome-Wide DNA Methylome and Transcriptome Analysis of Ovaries from Hu Sheep with High and Low Prolific. Front. Cell Dev. Biol. 2022, 10, 820558. [Google Scholar] [CrossRef]
  3. Liu, A.; Liu, M.; Li, Y.; Chen, X.; Zhang, L.; Tian, S. Differential expression and prediction of function of lncRNAs in the ovaries of low and high fecundity Hanper sheep. Reprod. Domest. Anim. 2021, 56, 604–620. [Google Scholar] [CrossRef]
  4. Wang, Z.; Wang, R.; Pan, C.; Chen, H.; Qu, L.; Wu, L.; Guo, Z.; Zhu, H.; Lan, X. Genetic Variations and mRNA Expression of Goat DNAH1 and Their Associations with Litter Size. Cells 2022, 11, 1371. [Google Scholar] [CrossRef] [PubMed]
  5. Wang, R.; Wang, Z.; Wang, X.; Li, Y.; Qu, L.; Lan, X. A novel 4-bp insertion within the goat CFAP43 gene and its association with litter size. Small Rumin. Res. 2021, 202, 106456. [Google Scholar] [CrossRef]
  6. Zhang, Y.; Chen, X.; Ruan, Y.; Guo, W.; Chen, J.; Tang, W.; Ji, Q.; Fu, K. Effect of CTSS non-synonymous mutations on litter size in Qianbei Ma goats. Front. Vet. Sci. 2023, 10, 1276673. [Google Scholar] [CrossRef] [PubMed]
  7. Smołucha, G.; Gurgul, A.; Jasielczuk, I.; Kawęcka, A.; Miksza-Cybulska, A. A genome-wide association study for prolificacy in three Polish sheep breeds. J. Appl. Genet. 2021, 62, 323–326. [Google Scholar] [CrossRef]
  8. Zhu, H.; Zhang, Y.; Bai, Y.; Yang, H.; Yan, H.; Liu, J.; Shi, L.; Song, X.; Li, L.; Dong, S.; et al. Relationship between SNPs of POU1F1 Gene and Litter Size and Growth Traits in Shaanbei White Cashmere Goats. Animals 2019, 9, 114. [Google Scholar] [CrossRef] [PubMed]
  9. Kang, B.; Wang, J.; Zhang, H.; Shen, W.; El-Mahdy Othman, O.; Zhao, Y.; Min, L. Genome-wide profile in DNA methylation in goat ovaries of two different litter size populations. J. Anim. Physiol. Anim. Nutr. 2021, 106, 239–249. [Google Scholar] [CrossRef]
  10. Zhang, Y.; Li, F.; Feng, X.; Yang, H.; Zhu, A.; Pang, J.; Han, L.; Zhang, T.; Yao, X.; Wang, F. Genome-wide analysis of DNA Methylation profiles on sheep ovaries associated with prolificacy using whole-genome Bisulfite sequencing. BMC Genom. 2017, 18, 759. [Google Scholar] [CrossRef]
  11. Zhou, X.; Yang, S.; Yan, F.; He, K.; Zhao, A. Genome-wide DNA methylation profiles of porcine ovaries in estrus and proestrus. Physiol. Genom. 2018, 50, 714–723. [Google Scholar] [CrossRef]
  12. Frattini, S.; Capra, E.; Lazzari, B.; McKay, S.D.; Coizet, B.; Talenti, A.; Groppetti, D.; Riccaboni, P.; Pecile, A.; Chessa, S.; et al. Genome-wide analysis of DNA methylation in hypothalamus and ovary of Capra hircus. BMC Genom. 2017, 18, 476. [Google Scholar] [CrossRef]
  13. Krueger, F.; Andrews, S.R. Bismark: A flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics 2011, 27, 1571–1572. [Google Scholar] [CrossRef]
  14. Lister, R.; Pelizzola, M.; Dowen, R.H.; Hawkins, R.D.; Hon, G.; Tonti-Filippini, J.; Nery, J.R.; Lee, L.; Ye, Z.; Ngo, Q.-M.; et al. Human DNA methylomes at base resolution show widespread epigenomic differences. Nature 2009, 462, 315–322. [Google Scholar] [CrossRef] [PubMed]
  15. Xiang, H.; Zhu, J.; Chen, Q.; Dai, F.; Li, X.; Li, M.; Zhang, H.; Zhang, G.; Li, D.; Dong, Y.; et al. Single base–resolution methylome of the silkworm reveals a sparse epigenomic map. Nat. Biotechnol. 2010, 28, 516–520. [Google Scholar] [CrossRef]
  16. Catoni, M.; Tsang, J.M.; Greco, A.P.; Zabet, N.R. DMRcaller: A versatile R/Bioconductor package for detection and visualization of differentially methylated regions in CpG and non-CpG contexts. Nucleic Acids Res. 2018, 46, e114. [Google Scholar] [CrossRef] [PubMed]
  17. Quinlan, A.R.; Hall, I.M. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics 2010, 26, 841–842. [Google Scholar] [CrossRef] [PubMed]
  18. Mi, F.; Wu, X.; Wang, Z.; Wang, R.; Lan, X. Relationships between the Mini-InDel Variants within the Goat CFAP43 Gene and Body Traits. Animals 2022, 12, 3447. [Google Scholar] [CrossRef] [PubMed]
  19. Coutton, C.; Vargas, A.S.; Amiri-Yekta, A.; Kherraf, Z.-E.; Ben Mustapha, S.F.; Le Tanno, P.; Wambergue-Legrand, C.; Karaouzène, T.; Martinez, G.; Crouzy, S.; et al. Mutations in CFAP43 and CFAP44 cause male infertility and flagellum defects in Trypanosoma and human. Nat. Commun. 2018, 9, 686. [Google Scholar] [CrossRef] [PubMed]
  20. Tao, L.; He, X.Y.; Wang, F.Y.; Pan, L.X.; Wang, X.Y.; Gan, S.Q.; Di, R.; Chu, M.X. Identification of genes associated with litter size combining genomic approaches in Luzhong mutton sheep. Anim. Genet. 2021, 52, 545–549. [Google Scholar] [CrossRef]
  21. Zhang, L.; Wang, F.; Gao, G.; Yan, X.; Liu, H.; Liu, Z.; Wang, Z.; He, L.; Lv, Q.; Wang, Z.; et al. Genome-Wide Association Study of Body Weight Traits in Inner Mongolia Cashmere Goats. Front. Vet. Sci. 2021, 8, 752746. [Google Scholar] [CrossRef]
  22. Qu, X.; Guo, S.; Yan, L.; Zhu, H.; Li, H.; Shi, Z. TNFα-Erk1/2 signaling pathway-regulated SerpinE1 and SerpinB2 are involved in lipopolysaccharide-induced porcine granulosa cell proliferation. Cell. Signal. 2020, 73, 109702. [Google Scholar] [CrossRef]
  23. Wu, W.B.; Xu, Y.Y.; Cheng, W.W.; Yuan, B.; Zhao, J.R.; Wang, Y.L.; Zhang, H.J. Decreased PGF May Contribute to Trophoblast Dysfunction in Fetal Growth Restriction. Reproduction 2017, 154, 319–329. [Google Scholar] [CrossRef]
  24. Vrachnis, N.; Kalampokas, E.; Sifakis, S.; Vitoratos, N.; Kalampokas, T.; Botsis, D.; Iliodromiti, Z. Placental growth factor (PlGF): A key to optimizing fetal growth. J. Matern.-Fetal Neonatal Med. 2013, 26, 995–1002. [Google Scholar] [CrossRef]
  25. Seidenspinner, T.; Tetens, J.; Habier, D.; Bennewitz, J.; Thaller, G. The placental growth factor (PGF)—A positional and functional candidate gene influencing calving ease and stillbirth in German dairy cattle. Anim. Genet. 2011, 42, 22–27. [Google Scholar] [CrossRef]
  26. Wang, H.; Yang, Q.; Gu, Y.; Zhang, X.; Wang, J.-M.; He, Y.-P.; Shi, Y.; Sun, Z.-G.; Shi, H.-J.; Wang, J. Uterine Expression of NDRG4 Is Induced by Estrogen and Up-Regulated during Embryo Implantation Process in Mice. PLoS ONE 2016, 11, e0155491. [Google Scholar] [CrossRef]
  27. Jandrey, E.H.F.; Moura, R.P.; Andrade, L.N.S.; Machado, C.L.; Campesato, L.F.; Leite, K.R.M.; Inoue, L.T.; Asprino, P.F.; da Silva, A.P.M.; de Barros, A.C.S.D.; et al. NDRG4 promoter hypermethylation is a mechanistic biomarker associated with metastatic progression in breast cancer patients. npj Breast Cancer 2019, 5, 11. [Google Scholar] [CrossRef] [PubMed]
  28. Faja, F.; Finocchi, F.; Carlini, T.; Rizzo, F.; Pallotti, F.; Spaziani, M.; Balercia, G.; Lenzi, A.; Paoli, D.; Lombardo, F. PDE11A gene polymorphism in testicular cancer: Sperm parameters and hormonal profile. J. Endocrinol. Investig. 2021, 44, 2273–2284. [Google Scholar] [CrossRef] [PubMed]
  29. Pathak, A.; Stewart, D.R.; Faucz, F.R.; Xekouki, P.; Bass, S.; Vogt, A.; Zhang, X.; Boland, J.; Yeager, M.; Loud, J.T.; et al. Rare inactivating PDE11A variants associated with testicular germ cell tumors. Endocr. Relat. Cancer 2015, 22, 909–917. [Google Scholar] [CrossRef] [PubMed]
  30. Laqqan, M.M.; Yassin, M.M. Influence of tobacco cigarette heavy smoking on DNA methylation patterns and transcription levels of MAPK8IP3, GAA, ANXA2, PRRC2A, and PDE11A genes in human spermatozoa. Middle East Fertil. Soc. J. 2021, 26, 41. [Google Scholar] [CrossRef]
  31. Zhang, W.; Zhang, S.; Xu, Y.; Ma, Y.; Zhang, D.; Li, X.; Zhao, S. The DNA Methylation Status of Wnt and Tgfβ Signals Is a Key Factor on Functional Regulation of Skeletal Muscle Satellite Cell Development. Front. Genet. 2019, 10, 220. [Google Scholar] [CrossRef]
  32. Yuan, X.-L.; Zhang, Z.; Li, B.; Gao, N.; Zhang, H.; Sangild, P.T.; Li, J.-Q. Genome-wide DNA methylation analysis of the porcine hypothalamus-pituitary-ovary axis. Sci. Rep. 2017, 7, 4277. [Google Scholar] [CrossRef] [PubMed]
  33. Fan, Y.; Liang, Y.; Deng, K.; Zhang, Z.; Zhang, G.; Zhang, Y.; Wang, F. Analysis of DNA methylation profiles during sheep skeletal muscle development using whole-genome bisulfite sequencing. BMC Genom. 2020, 21, 327. [Google Scholar] [CrossRef] [PubMed]
  34. Hao, Y.; Cui, Y.; Gu, X. Genome-wide DNA methylation profiles changes associated with constant heat stress in pigs as measured by bisulfite sequencing. Sci. Rep. 2016, 6, 27507. [Google Scholar] [CrossRef] [PubMed]
  35. Lister, R.; Pelizzola, M.; Kida, Y.S.; Hawkins, R.D.; Nery, J.R.; Hon, G.; Antosiewicz-Bourget, J.; O’Malley, R.; Castanon, R.; Klugman, S.; et al. Hotspots of aberrant epigenomic reprogramming in human induced pluripotent stem cells. Nature 2011, 471, 68–73. [Google Scholar] [CrossRef] [PubMed]
  36. Yuan, X.-L.; Gao, N.; Xing, Y.; Zhang, H.-B.; Zhang, A.-L.; Liu, J.; He, J.-L.; Xu, Y.; Lin, W.-M.; Chen, Z.-M.; et al. Profiling the genome-wide DNA methylation pattern of porcine ovaries using reduced representation bisulfite sequencing. Sci. Rep. 2016, 6, 22138. [Google Scholar] [CrossRef] [PubMed]
  37. Bourc’his, D.; Malik, H.S.; Molaro, A.; Satta, Y. Dynamic Evolution of De Novo DNA Methyltransferases in Rodent and Primate Genomes. Mol. Biol. Evol. 2020, 37, 1882–1892. [Google Scholar] [CrossRef]
  38. Stevenson, T.J.; Bowman, A.S.; Campbell, E.M.; Lorgen, M.; Coyle, C.S.; Lynch, E.W.J. Cyclical DNA Methyltransferase 3a Expression Is a Seasonal and Estrus Timer in Reproductive Tissues. Endocrinology 2016, 157, 2469–2478. [Google Scholar] [CrossRef]
  39. Harris, C.J.; Scheibe, M.; Wongpalee, S.P.; Liu, W.; Cornett, E.M.; Vaughan, R.M.; Li, X.; Chen, W.; Xue, Y.; Zhong, Z.; et al. A DNA methylation reader complex that enhances gene transcription. Science 2018, 362, 1182–1186. [Google Scholar] [CrossRef]
  40. Halpern, K.B.; Vana, T.; Walker, M.D. Paradoxical Role of DNA Methylation in Activation of FoxA2 Gene Expression during Endoderm Development. J. Biol. Chem. 2014, 289, 23882–23892. [Google Scholar] [CrossRef]
  41. Wan, J.; Oliver, V.F.; Wang, G.; Zhu, H.; Zack, D.J.; Merbs, S.L.; Qian, J. Characterization of tissue-specific differential DNA methylation suggests distinct modes of positive and negative gene expression regulation. BMC Genom. 2015, 16, 49. [Google Scholar] [CrossRef] [PubMed]
  42. Peng, W.F.; Xu, S.S.; Ren, X.; Lv, F.H.; Xie, X.L.; Zhao, Y.X.; Zhang, M.; Shen, Z.Q.; Ren, Y.L.; Gao, L.; et al. A genome-wide association study reveals candidate genes for the supernumerary nipple phenotype in sheep (Ovis aries). Anim. Genet. 2017, 48, 570–579. [Google Scholar] [CrossRef] [PubMed]
  43. Liu, C.; He, X.; Liu, W.; Yang, S.; Wang, L.; Li, W.; Wu, H.; Tang, S.; Ni, X.; Wang, J.; et al. Bi-allelic Mutations in TTC29 Cause Male Subfertility with Asthenoteratospermia in Humans and Mice. Am. J. Hum. Genet. 2019, 105, 1168–1181. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Genome-wide DNAm profiles of HP and LP Jining Grey goat ovaries. (A) Experimental design. (B) Mean proportion of diverse methylation contexts in the ovaries of Jining Grey goats. Blue, yellow and green represent methylated mCG, mCHH and mCHG, respectively. (C) Methylation levels in different genomic elements. Horizontal coordinates represent genomic elements, and vertical coordinates represent methylation levels; the value is the average methylation level within a gene element, with different colors representing different groups. (D) Distribution of DNAm level at 3K upstream/downstream. Horizontal coordinates show diverse regions, vertical coordinates represent methylation levels, with different colors representing different samples. Genebody: the active coding region of a gene from the beginning to the end of active coding (includes the exon and intron of the gene); UP3K: refers to 3000 bp before starting from the transcription start site (TSS); DOWN3k:refers to the transcription termination site (TTS) after 3000 BP.
Figure 1. Genome-wide DNAm profiles of HP and LP Jining Grey goat ovaries. (A) Experimental design. (B) Mean proportion of diverse methylation contexts in the ovaries of Jining Grey goats. Blue, yellow and green represent methylated mCG, mCHH and mCHG, respectively. (C) Methylation levels in different genomic elements. Horizontal coordinates represent genomic elements, and vertical coordinates represent methylation levels; the value is the average methylation level within a gene element, with different colors representing different groups. (D) Distribution of DNAm level at 3K upstream/downstream. Horizontal coordinates show diverse regions, vertical coordinates represent methylation levels, with different colors representing different samples. Genebody: the active coding region of a gene from the beginning to the end of active coding (includes the exon and intron of the gene); UP3K: refers to 3000 bp before starting from the transcription start site (TSS); DOWN3k:refers to the transcription termination site (TTS) after 3000 BP.
Genes 15 00353 g001
Figure 2. Identification of DMRs and DMGs. (A) Number of DMRs distributed in mCG, and mCHG, and number of hypo and hyper DMRs in CG. (B) Number of DMG distributions in mCG, and mCHG, and number of hypo and hyper DMG in CG. (C) Distribution of DMR methylation levels in CG. The horizontal axis indicates HP and LP groups, and the vertical coordinate indicates methylation level values. (D) DMR anchor region in CG. On the x-axis, each region’s type is indicated, while the number of hyper/hypo DMRs in each region is represented on the y-axis.
Figure 2. Identification of DMRs and DMGs. (A) Number of DMRs distributed in mCG, and mCHG, and number of hypo and hyper DMRs in CG. (B) Number of DMG distributions in mCG, and mCHG, and number of hypo and hyper DMG in CG. (C) Distribution of DMR methylation levels in CG. The horizontal axis indicates HP and LP groups, and the vertical coordinate indicates methylation level values. (D) DMR anchor region in CG. On the x-axis, each region’s type is indicated, while the number of hyper/hypo DMRs in each region is represented on the y-axis.
Genes 15 00353 g002
Figure 3. DMG functional enrichment analysis. (A) GO. horizontal coordinate represents −lg (p value), and vertical coordinate represents the corresponding GO terms. (B) KEGG: vertical and horizontal coordinate indicate enriched pathway and the corresponding pathway, respectively. The size of the dots indicates the amount of DMG contained in each pathway, and the color of the dots reflects the −log10 (p-value) of each pathway.
Figure 3. DMG functional enrichment analysis. (A) GO. horizontal coordinate represents −lg (p value), and vertical coordinate represents the corresponding GO terms. (B) KEGG: vertical and horizontal coordinate indicate enriched pathway and the corresponding pathway, respectively. The size of the dots indicates the amount of DMG contained in each pathway, and the color of the dots reflects the −log10 (p-value) of each pathway.
Genes 15 00353 g003
Figure 4. Combined DNAm and transcriptome analysis of HP and LP Jining Grey goat ovaries. (A) Venn diagram represents the intersection of DMG and DEG in WGBS and RNA-Seq; red represents DMG; and blue represents DEG. (B) DMG-DEG divided into groups according to hyper and hypo; vertical coordinate indicates the number of DMG-DEG in each group.
Figure 4. Combined DNAm and transcriptome analysis of HP and LP Jining Grey goat ovaries. (A) Venn diagram represents the intersection of DMG and DEG in WGBS and RNA-Seq; red represents DMG; and blue represents DEG. (B) DMG-DEG divided into groups according to hyper and hypo; vertical coordinate indicates the number of DMG-DEG in each group.
Genes 15 00353 g004
Figure 5. Key candidate gene validation. Methylation levels of (A) NDRG4, (B) SERPINB2, (C) EPHA6, and (D) LRP1B genes were validated by bisulfite sequencing PCR (BSP). Each circle represents a CpG dinucleotide, and black and white denote methylated and unmethylated sites, respectively. RT-qPCR was utilized to examine the mRNA expression of the (E) DNMTs in the ovaries and the (F) eight relative mRNA expressions of key candidate genes; mean ± SEM of 3 biological replicates of 2−ΔΔCt value.
Figure 5. Key candidate gene validation. Methylation levels of (A) NDRG4, (B) SERPINB2, (C) EPHA6, and (D) LRP1B genes were validated by bisulfite sequencing PCR (BSP). Each circle represents a CpG dinucleotide, and black and white denote methylated and unmethylated sites, respectively. RT-qPCR was utilized to examine the mRNA expression of the (E) DNMTs in the ovaries and the (F) eight relative mRNA expressions of key candidate genes; mean ± SEM of 3 biological replicates of 2−ΔΔCt value.
Genes 15 00353 g005
Table 1. Primer information for RT-qPCR.
Table 1. Primer information for RT-qPCR.
Primer NamesSequences (5′-3′)Product Length (bp)Gene Bank
ForwardReverse
GAPDHCGGCACAGTCAAGGCAGAGAACCACGTACTCAGCACCAGCATCAC115XM_005680968.3
DNMT1CCTGACTCCACCTACGAAGACCCTACTTGCTCCACCACGAACTG127XM_018051019.1
DNMT3ACCGCATTGTGTCTTGGTGGATGAGAACTTGCCGTCTCCGAACC86XM_018055548.1
DNMT3BTTGACTTGGTGATTGGTGGAAGCCGAGTGTAATTCAGCAGGTGGTAG124XM_018057711.1
PGFTGAGACTGTTCACTTGCTTCCTGGCTGCGGCTCCACACTTG131XM_005686088.3
EPHA6AGGATACAGTGGCTACAGTCAGAAGACGGCTGCGGTGGCTATC105XM_013972057.2
NFASCACGACAGCCTGGTGGACTATGGTCCGTCTCCTCCTTGTCCTTCTTG103XM_018060022.1
NDRG4CCAACCACCACGACCTTCCTGAGCCTTCGGTCCTTCAAGTATGC137XM_018062073.1
LRP1BGCACTGCGACTCTGATGATGACTGATCTGCCACTGGAACAACTGAACTG98XM_018062463.1
LOC102181552TGTTCCTGCTGCTTCCAGATGGGTACACCTCCACGTCATCCTCAC134XM_005697318.2
CFAP43GCAAGGAAGCAGGAGGAGAGGCCGCCAGTAAGTTGTCATAGTGTTC100XM_018041743.1
TTC29GCATTAGCAGTCCTCAACACTTACGCCGTGGCTTCATAGGCTCTCC85XM_005691202.2
PDE11ACTCTATGGAACCTCTGCCACCTTGATGTTGTGACCCTCGCTCTGAAG80XM_005675916.3
Table 2. Primers sequences of DNA-methylation-related genes.
Table 2. Primers sequences of DNA-methylation-related genes.
RegionSequences (5′-3′)Product Length (bp)Associate GeneGene Bank
chr18:27599410-27599622AAATATTTGTGTTTGGAGTTATTTT213NDRG4NC_030825.1
TACCCCAATCTTTATTTAAATTCCC
chr24:62258787_62258914TTTAGATTTAGTGATTAGATTAGTGGT128SERPINB2NC_030831.1
AACAATAAAACCTAACTCCATACC
chr2:79444158-79444322AAAATTTATGTAAATGTTAGAAGTTT165LRP1BNC_030809.1
AATAATTAAATAACATCACCAACTC
chr1:39785850-39786019TTATTTAGGTTAGGGGTTGAGGG170EPHA6NC_030808.1
ACACTAACAAAACCAAAAACAAAAC
Table 3. Whole-genome DNA bisulfite sequencing data.
Table 3. Whole-genome DNA bisulfite sequencing data.
GroupsSampleClean Reads 1Q30 1 (%)Total Mapping Reads 1 (%)Unique Mapping Ratio 1 (%)Bisulfite Conversion Rate 1 (%)Total_mC 1 (%)
HPHP-O1113947443691.3985.5479.1799.543.50
HP-O2114809557091.6283.577.5799.513.43
HP-O3112441907691.2184.5378.4099.513.55
HP-O4118541997491.3180.0574.1799.473.30
LPLP-O1123462894091.9483.8177.9099.533.24
LP-O2112093589890.9483.4877.4299.533.48
LP-O3134187365691.2184.1877.8199.512.98
LP-O4104224958490.1682.4976.8099.473.66
1 Clean base (Gb): total number of unambiguous base calls in the sequences; clean reads: total number of sequence reads after filtering for quality and contaminants; unique mapping ratio (%): proportion of clean reads that aligned uniquely to the reference genome out of the total clean sequences; bisulfite conversion rate (%): proportion of clean reads that underwent successful bisulfite conversion relative to the total amount of methylation observed in the clean sequences aligned to the reference genome; total mC (%): proportion of methylated cytosines within the clean sequences aligned to the reference genome relative to the total number of clean sequences aligned to the reference genome.
Table 4. Genes obtained from unite analysis of Whole-Genome Bisulfite Sequencing (WGBS) and RNA-Seq.
Table 4. Genes obtained from unite analysis of Whole-Genome Bisulfite Sequencing (WGBS) and RNA-Seq.
GroupGenesGene Counts
Hyper DMG-Hyper DEGNFASC,GLT1D1,DOCK10,JAKMIP3,GRM5,GRIK26
Hyper DMG-Hypo DEGCR2, CDH20, CFAP43 [5,18,19], DAGLA, GRM8, ZNF804B, RELN, TRHDE, LOC102172692, CLEC9A10
Hypo DMG-Hyper DEGISM1, MMP16, NCAM1, USH2A, GALNT9, VAT1L, UPK1A, LRP1B 1 [20,21], FMNL2, CTNND2, FSD2, MTCL1, SERPINB2 1 [22], ADAM12, DLG2, LOC102170378, GSTM3, GRIK2, TMEM176B19
Hypo DMG-Hypo DEGCADM2, EPHA6 1 [7], C1H3orf52, PGF 1 [23,24,25], ALDH1A2, LOC108637248, COL14A1, TTC29 1 [3], NDRG4 1 [26,27], HPN, ARHGAP27, PDE11A 1 [28,29,30], MFSD6, MYO7B, JAG2, CDH20, TRPM8, C3H1orf87, C3H1orf226, COL28A1, SUN3, CRACR2A, UNC5A, RFX2, FBXO16, LOC10863459426
1 Genes in bold are selected key candidate genes, and relevant references are included in [].
Table 5. DMR methylation levels of 8 critical genes of HP and LP Jining Grey goat.
Table 5. DMR methylation levels of 8 critical genes of HP and LP Jining Grey goat.
GeneChr 1StartEndElementMeth-DirectionMeth Diff 1Q-Value
SERPINB2246225880162259000UP3KHypo−0.3311.09 × 10−3
NDRG4182759920127599400UP3KHypo−0.3264.35 × 10−5
CFAP43262704660127046800introHyper0.4351.65 × 10−2
LRP1B27972020179720400introHypo−0.3122.81 × 10−2
EPHA614008860140088800introHypo−0.3091.19 × 10−4
TTC29175938580159386000introHypo−0.3281.36 × 10−5
PDE11A2117418601117418800introHypo−0.4033.33 × 10−5
PGF101691380116914000exonHypo−0.3098.42 × 10−4
1 Chr: chromosome, Meth diff: difference in DNAm levels between HP and LP.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Yang, C.; He, J.; Mao, J.; Ren, Y.; Liu, G.; Wei, C.; Zhang, G.; Tian, K.; Huang, X. Genome-Wide DNA Methylation Analysis and Functional Validation of Litter Size Traits in Jining Grey Goats. Genes 2024, 15, 353. https://doi.org/10.3390/genes15030353

AMA Style

Yang C, He J, Mao J, Ren Y, Liu G, Wei C, Zhang G, Tian K, Huang X. Genome-Wide DNA Methylation Analysis and Functional Validation of Litter Size Traits in Jining Grey Goats. Genes. 2024; 15(3):353. https://doi.org/10.3390/genes15030353

Chicago/Turabian Style

Yang, Cunming, Junmin He, Jingyi Mao, Yifan Ren, Guifen Liu, Chen Wei, Guoping Zhang, Kechuan Tian, and Xixia Huang. 2024. "Genome-Wide DNA Methylation Analysis and Functional Validation of Litter Size Traits in Jining Grey Goats" Genes 15, no. 3: 353. https://doi.org/10.3390/genes15030353

APA Style

Yang, C., He, J., Mao, J., Ren, Y., Liu, G., Wei, C., Zhang, G., Tian, K., & Huang, X. (2024). Genome-Wide DNA Methylation Analysis and Functional Validation of Litter Size Traits in Jining Grey Goats. Genes, 15(3), 353. https://doi.org/10.3390/genes15030353

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop