Whole transcriptome analysis in oviduct provides insight into microRNAs and ceRNA regulative networks that targeted reproduction of goat (Capra hircus)

Abstract

Background

Reproduction traits are crucial for livestock breeding and represent key economic indicators in the domestic goat (Capra hircus) industry. The oviduct, a critical organ in female mammals, plays a pivotal role in several reproductive processes; however, its molecular mechanisms remain largely unknown. Non-coding RNA and mRNAs are essential regulatory elements in reproductive processes; yet their specific roles and regulatory networks in goat oviducts remain unclear.

Results

In this study, we conducted small RNA sequencing of the oviduct in high- and low-fecundity goats during the follicular (FH and FL groups, n = 10) and luteal (LH and LL groups, n = 10) phase, profiling 20 tissue samples. Combinatorial whole-transcriptome expression profiles were applied to the same samples to uncover the competitive endogenous RNA (ceRNA) regulation network associated with goat fecundity. RT-qPCR was employed to validate the miRNA profiling results, and ceRNA regulatory networks were analyzed through luciferase assay. Gene set enrichment analysis (GSEA) confirmed that the cytokine-cytokine receptor interaction and TGF-β signaling pathway, both related to embryonic development, were enriched in DEM target genes. Additionally, miR-328-3p, a core miRNA, targets SMAD3 and BOP1, which are involved in the negative regulation of cell growth and embryonic development. TOB1 and TOB2, targeted by miR-204-3p, regulate cell proliferation via the protein kinase C-activating G-protein coupled receptor signaling pathway. Analyses of ceRNA regulatory networks revealed that LNC_005981 − miR-328-3p − SMAD3 and circ_0021923 − miR-204-3p − DOT1L may affect goats’ reproduction, findings that were validated using luciferase assay.

Conclusion

Analysis of whole-transcriptome profiling of goat oviducts identified several key miRNAs and ceRNAs that may regulate oocyte maturation, embryo development, and the interactions between the oviduct and gametes/early embryos, providing insights into the molecular mechanisms of reproductive regulatory networks.

Introduction

Kidding number is a key economic trait in goats that directly influences productivity and is regulated by genetic and non-genetic factors, such as nutrition and environment [1, 2]. Therefore, investigating the genetic determinants of goat prolificacy remains one of the most effective approaches to improving productivity. The oviduct, a vital reproductive organ in mammals, links the ovary and the uterus, facilitating key processes such as gamete maturation, sperm storage, fertilization, early embryo development, and embryo transport [3]. Recent studies have highlighted the critical roles of non-coding RNAs (ncRNAs) in germ cell biogenesis and regulating reproductive organ activity [4], indicating their potential involvement in oviductal function and goat fecundity.

Advances in high-throughput sequencing technology have gradually unraveled the mysteries of genome transcription [5, 6], showing that protein-coding genes constitute less than 2% of the entire genome, while non-coding RNAs (ncRNAs) comprise over 98% [7, 8]. ncRNAs are primarily classified into three categories based on their size and structure: microRNA (miRNA), long non-coding RNA (lncRNA), and circular RNA (circRNA) [9]. By binding to the 3’-untranslated region (3’-UTR) of target messenger RNA (mRNA), miRNAs, a small class of endogenous ncRNA approximately 21 nucleotides in length, regulate gene expression at the post-transcriptional level, leading mRNAs degradation or silencing [10]. Accumulating evidence indicates that miRNAs regulate the growth of tissues and cells, including breast tissue, skeletal muscle, and follicles, by inhibiting the translation of their target mRNAs [11]. Furthermore, miRNAs are integral to various reproduction processes, including sex differentiation [12], oogenesis [13], fertilization [14], zygotic genome activation (ZGA) [15], early embryonic development, implantation [16,17,18], and pregnancy [19]. For instance, miR-338-3p can affect early follicular development and normal ovary functions by interfering with the proliferation and estradiol production of GCs [20]. Similarly, miR-92b-3p is a novel regulator of primordial follicle assembly by negatively regulating TSC1 in mTOR/Rps6 signaling [21]. Another study revealed that inhibiting miRNA-155 improves the preimplantation developmental competence of porcine embryos by upregulating ZEB2 and downregulating ATF4 [22].

Competitive endogenous RNA (ceRNA) represents a paradigm wherein coding and non-coding RNAs (circRNA and lncRNA) regulate gene expression by competing for shared miRNA targets [23]. In comparison to the microRNA regulatory network, the ceRNA network represents a novel and more intricate model of gene expression regulation, involving a greater variety of RNA molecules [11]. Circ0001470, functioning as a ceRNA, has been shown to facilitate embryonic development by sponging miR-140-3p, regulating downstream PTGFR expression, activating the MAPK reproductive pathway, and sustaining pregnancy maintenance [24]. Another study revealed that miR-34a and miR-34c target circ-8073 and CEP55, modulating the RAS/RAF/MEK/ERK and PI3K/AKT/mTOR pathways to trigger apoptotic in caprine endometrial epithelial cell (CEEC) [25]. There have also been studies on lncRNAs acting as ceRNAs. For example, lncRNA366.2 enhances endometrial epithelial cell (EEC) proliferation, migration, and growth factor expression by sponging miR-1576, upregulating WNT6 expression, and activating the Wnt/β-catenin pathway, thereby identifying a candidate gene related to fertility [26]. Hu et al. discovered that lncRNA TCONS_00814106 is upregulated in the ovarian tissue of high-fecundity sows and functions as a ceRNA, regulating TGFBR1 by sponging miR-1343. This promotes the proliferation of porcine granulosa cells (GC) and inhibits apoptosis [27]. As mentioned, above, these findings underscore the critical role of the ceRNA network in reproductive regulation. However, the specific ceRNA regulatory network comprising lncRNAs, circRNAs, miRNAs, and mRNAs that influences goat fertility in the oviduct remains unclear.

In this study, we aimed to investigate the roles of miRNAs and ceRNA networks in prolificacy by collecting oviduct samples from high- and low-fecundity Yunshang black goats. First, we examined the regulatory functions and molecular mechanisms of miRNAs, which were validated using RT-qPCR. miRNA − mRNA interaction pairs were identified by integrating whole transcriptome sequencing data from the same samples. Subsequently, GSEA was conducted to identify key DEMs. Furthermore, we confirmed lncRNA − DEM and circRNA − DEM interaction pairs and constructed the ceRNA networks (lncRNA − DEM − mRNA and circRNA − DEM − mRNA) in the oviduct at different phases to elucidate the complex genetic architecture underlying prolificacy. Finally, the LNC_005981 − miR-328-3p − SMAD3 and circ_0021923 − miR-204-3p − DOT1L networks were experimentally validated using luciferase assay. These findings offer novel insights into fecundity mechanisms and contribute to the identification of embryonic development-related candidate biomarkers.

Materials and methods

Ethical statement

All animals were treated under protocol No. IAS2019-63 approved by the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences, Beijing, China. Ethical oversight was given by the Animal Ethics Committee of the IAS-CAAS.

Experimental design and workflow

The experimental design and workflow are illustrated in Fig. 1. Including i: Oviduct samples were collected from Yunshang Black goats with high- and low-fecundity during the follicular and luteal phases, with 5 replicates for each group. ii: miRNA sequencing was performed using the Illumina NovaSeq platform, followed by differential miRNA and target gene analyses. The miRNA sequencing results were validated using RT-qPCR. iii: Based on our previously published data, ceRNA interaction networks for lncRNAs and circRNAs were constructed and subjected to biological analysis; Finally, the ceRNA-target interaction relationships were validated using luciferase reporter assays.

Fig. 1

Experimental design and workflow for the whole transcriptomics-based comparison of high- and low-fecundity oviduct from Yunshang Black goats during the follicular and luteal phases

Animals and sample collection

Twenty Yunshang black goats with no significant differences in age (3 years old), weight (52.22 ± 0.43 kg), height (63.2 ± 0.58 cm), or genetic background were selected. Based on kidding number records (p < 0.05), the goats were divided into a high fecundity group (H, n = 10, kidding number = 3.3 ± 0.11) and a low fecundity group (L, n = 10, kidding number = 1.75 ± 0.08). The goats were raised under standard conditions with free access to food and water and exposed to natural temperature and lighting. Before the experiment, all goats were implanted with a controlled internal drug release device (CIDR, Inter Ag Co., Ltd. New Zealand; progesterone 300 mg) to synchronize estrus. The CIDR was removed after 16 days, marking 0 h as the starting time point. Subsequently, 10 goats (five high-fecundity and five low-fecundity) were humanely euthanized at 48 h (follicular phase; FH and FL) via intravenous injection of sodium pentobarbital at a dose of 100 mg/kg body weight (It was selected due to its rapid onset of action and minimal side effects), and another 10 goats (five high-fecundity and five low-fecundity) were euthanized at 168 h (luteal phase; LH and LL). Oviduct tissues were immediately frozen in liquid nitrogen for short-term storage, transported to the laboratory, and stored at -80 °C for long-term preservation until RNA extraction. All procedures were conducted with proper consideration for animal welfare.

MicroRNA library construction and sequencing

Total RNA for miRNA sequencing (miRNA-seq) was extracted from 20 oviduct tissue samples of Yunshang black goats using TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA). The RNA concentration and purity (OD 260/280) were measured using the Nanodrop 2000 instrument, and RNA integrity was assessed by agarose gel electrophoresis. RNA integrity was further evaluated using the Agilent Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, USA) with an RNA Nano 6000 Assay Kit. Library construction and sequencing followed our previously published methods [28]. In brief, 20 cDNA libraries were generated from oviduct tissues collected during the follicular (FH and FL) and luteal (LH and LL) phases. The Small RNA (sRNA) Sample Pre Kit was employed for miRNA library preparation from qualified RNA samples. Total RNA served as the input material, leveraging the unique structures of sRNA 3’ and 5’ ends. The sRNA ends were ligated directly and reverse-transcribed into cDNA. Following PCR amplification, DNA fragments were separated by PAGE gel electrophoresis and cDNA libraries were retrieved via gel excision. The quality and quantity of the cDNA libraries were assessed using Qubit 2.0 and Agilent Bioanalyzer 2100, resulting in the final miRNA library. SE50 (single-end, 50 bp) sequencing was conducted on the Illumina Novaseq platform for miRNAs.

Quality control, mapping, and assembly

Raw data in fastq format were initially processed using in-house scripts. Q20, Q30, and GC contents metrics of the clean data were calculated. Low-quality reads, reads contaminated with 5’ connector, reads lacking 3’ adaptor or insert fragment, reads containing continuous A/T/C/G sequences, and reads with abnormal lengths were removed. Subsequent analyses utilized clean, high-quality reads with lengths ranging from 18 to 35nt in length. sRNA tags were aligned to Bowtie (v.2.2.3) reference sequences to identify known miRNAs. Potential miRNAs and their secondary structures were identified using modified miRdeep (v.2.0.0.5) and sRNA-tools-cli, with miRBase (v.22.0) [29] as the reference. Small RNA tags were compared against the RepeatMasker, Rfam (14.2), and species-specific databases to exclude tags derived from protein-coding genes, repetitive sequences, ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), small nucleolar RNA (snoRNA). Novel miRNAs were predicted using miREvo (v.2.0.0.5) and miRdeep2 (v2.0.0.5) through analysis of secondary structures, Dicer cleavage sites, and minimum free energy of unannotated sRNA tags. The expression levels of known and novel miRNAs were quantified and normalized to TPM (Transcripts Per Kilobase of exon model per Million mapped reads) using the formula: TPM = (read count × 1,000,000) / libsize, where libsize represents the total miRNA read count in the sample.

Noncoding RNA and mRNA expression profiles

Details of RNA extraction, RNA-seq library preparation, sequencing, and bioinformatics analyses of mRNAs, lncRNAs, and circRNAs from the same samples were provided in our previous studies [30,31,32]. In this study, a p-value threshold of < 0.05 was applied to identify differentially expressed mRNAs (DEGs), lncRNAs (DELs), and circRNAs (DECs).

Differential expression MiRNAs analysis and target genes prediction

Raw sequencing reads were first aligned to the reference genome using Bowtie. Gene-level raw counts were then quantified using HTSeq. These raw counts were used as input for the differential expression analysis conducted with DESeq2 (v3.18.1) [33], which is designed to process count-based data for modeling using a negative binomial distribution. In this study, DESeq2 was employed to identify differentially expressed miRNAs (DEMs) in the FL vs. FH and LL vs. LH comparisons using raw read count data as input. The analysis provided fold change (FC) values and adjusted p-values (padj/q-values) for assessing differential expression. For visualization and downstream analysis, normalized expression values TPM (transcripts per kilobase per million reads) were calculated separately. A p-value < 0.05 and|log₂FC|>1 was used as default thresholds for defining significant DEMs. MiRNAs exert their biological functions primarily by suppressing or silencing target gene expression through post-transcriptional regulation. Potential miRNA target genes were predicted using MiRanda (v3.3a) [34] and qTar (https://github.com/zhuqianhua/ qTar.git) to assess their biological functions. To improve the accuracy of target gene prediction, putative DEMs target genes identified by miRanda were cross-referenced with DEGs from the same samples in our previous research [30, 31], and miRNA-mRNA correlations were analyzed.

Gene set enrichment analysis

GSEA is a computational method for identifying whether a gene set exhibit coordinated up- or down-regulation under specific conditions. The overlap between predicted DEM target genes and transcript genes was analyzed using GSEA and implemented via GSEA software (version 4.2.1). Gene sets analyzed in this study were derived from the Hallmark (H), KEGG (C2.CP. KEGG), and Gene Ontology (C5.GO. BP) collections within the Molecular Signature Database (MSigDB, v2022.1. Mm) database [35]. GSEA calculated the enrichment score (ES), normalized enrichment score (NES), and nominal p-value (Nom p-value) for each gene set. Gene set with nominal (Nom) p-value < 0.05,|NES| > 1, and false discovery rate (FDR) < 0.25 was considered statistically significant.

Prediction of competing endogenous RNA (ceRNA) network

The ceRNA regulatory network was constructed by integrating all previously identified competing triplets, and the key ceRNA network was visualized using a Sankey diagram. Additionally, the lncRNA-circRNA-miRNA-mRNA ceRNA network was constructed and visualized using Cytoscape software. To identify reproduction-related miRNAs and explore the molecular mechanisms regulating goat fecundity in the oviduct, ceRNA networks were constructed based on miRNA functions. The ceRNA network was constructed according to ceRNA theory. Pearson correlation analysis was applied to calculate miRNA-mRNA interactions using miRNA and transcriptome expression profiles, focusing on DEM-DEG pairs with negatively correlated expression. Similarly, MiRanda and qTar (https://github.com/zhuqianhua/qTar.git) were used to predict miRNA target circRNAs and lncRNAs. Pearson correlation analysis was also used to identify negatively correlated miRNA-lncRNA and miRNA-circRNA pairs. Negatively co-expressed lncRNA–miRNA, circRNA–miRNA, and miRNA–mRNA pairs were identified based on a Pearson correlation coefficient (cor) < − 0.4 and p-value of < 0.05. The ceRNA regulatory network was constructed by integrating all previously identified competing triplets, and the key ceRNA network was visualized using a Sankey diagram. Additionally, the lncRNA-circRNA-miRNA-mRNA ceRNA network was constructed and visualized using Cytoscape software.

Bioinformatics analysis of ceRNA network

The target genes of DEMs in the ceRNA network were analyzed using Gene Ontology Biological Processes (GO BP) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis through the Database for Annotation, Visualization and Integrated Discover (DAVID, version 2021, https://david-d.ncifcrf.gov/) [36]. Each list of DEMs target DEGs was submitted separately to DAVID for analysis, with Capra hircus used as the background. GO terms and KEGG pathways with a p-value < 0.05 were considered significantly enriched.

Data validation by RT-qPCR

We have adhered to the Minimum Information for Publication of Quantitative Real-Time PCR Experiments (MIQE) guidelines for real-time quantitative polymerase chain reaction (RT-qPCR) of miRNA. Briefly, DEMs were selected to validate the accuracy of miRNA sequencing data. Primers of miRNA were designed using Primer Premier 6 software. The primer information is provided in Table S1. All primers were synthesized by Sangon Biotech (Shanghai, China). For reverse transcription, total RNA was isolated from the samples using TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA). The RNA concentration and purity (OD 260/280) were measured using the Nanodrop 2000 instrument, and RNA integrity was assessed by agarose gel electrophoresis. The reverse transcription reaction was performed using the MiRcute Plus miRNAs First-Strand cDNA Kit (TIANGEN, Beijing, China). The reverse transcription reaction was carried out as follows: incubation at 42 °C for 60 min, heat inactivation at 95 °C for 5 min. Immediately, miRNA quantification was performed using the miRcute Plus miRNA qPCR Kit (TIANGEN, Beijing, China). The RT-qPCR reaction was performed on a RocheLight Cycler^®480 II system (Roche Applied Science, Mannheim, Germany). All RT-qPCR reactions were carried out in triplicate to ensure the reliability and reproducibility of the results. The RT-qPCR thermal cycling conditions was conducted in the following procedure: initial denaturation at 95 °C for 15 min, followed by 40 cycles of denaturation at 94 °C for 20 s, then annealing at 60 °C for 34 s. For relative quantification of miRNA expression, the 2^−ΔΔCt method was used, with U6 as the reference gene. All the RT-qPCR results were presented as the mean ± SEM. The relative expression level of each miRNA was calculated using the following formula:

$$\:\text{R}\text{e}\text{l}\text{a}\text{t}\text{i}\text{v}\text{e}\:\text{E}\text{x}\text{p}\text{r}\text{e}\text{s}\text{s}\text{i}\text{o}\text{n}\hspace{0.17em}=\hspace{0.17em}2^{-\left(\Delta\Delta\text{C}\text{t}\right)}$$

Note

ΔCt is the difference between the ΔCt value of the target miRNA and the reference gene (U6); ΔΔCt is the difference between the ΔCt of the experimental sample and the ΔCt of the control group.

Vector construction and luciferase reporter assay

Forty-eight hours post-transfection, HEK293T cells were harvested, and luciferase activity was measured using the Dual-Luciferase Reporter Assay System (TECAN, Switzerland) following the manufacturer’s instructions. Firefly luciferase activity was normalized to Renilla luciferase activity as an internal control. The wild-type (WT) and mutant (MUT) sequences of LNC_005981 (SMAD3-3’UTR) and circ_0021923 (DOT1L-3’UTR), containing potential chi-miR-328-3p and chi-miR-204-3p MRE motifs, were synthesized and cloned into the SacI/XhoI sites of the pmirGLO luciferase reporter vector (GenePharma, Shanghai, China). MiR-328-3p and miR-204-3p mimics, as well as the mimics negative control (NC), were synthesized by GenePharma (Shanghai, China). To overexpress LNC_005981 and circ_0021923, their full-length sequences were amplified, BspeI/EcoRI digestion sites were introduced at both ends, and the fragments were ligated into the pEGFP-C1 overexpression vector (Sangon Biotech, Shanghai, China), resulting in pEGFP-LNC_005981 and pEGFP-circ_0021923. All vectors were verified by Sanger sequencing, and luciferase reporter and overexpression vectors, along with mimics and NC, were transfected into cells according to the manufacturer’s protocol. HEK293T cells were seeded in 24-well plates and cultured in DMEM supplemented with 10% fetal bovine serum and 1% penicillin − streptomycin. Forty-eight hours post-transfection, HEK293T cells were harvested, and luciferase activity was measured using the Dual-Luciferase Reporter Assay System (TECAN, Switzerland) following the manufacturer’s instructions. Firefly luciferase activity was normalized to Renilla luciferase activity as an internal control.

Statistics

To assess the statistical significance of differences between groups, Student’s t test was conducted at the kidding number significance level and the result of luciferase activities using SPSS 25.0 (SPSS, USA). The relative miRNA expression levels were determined using RT‒qPCR data and the 2^−ΔΔCt method. The data are presented as the mean ± SEM and were derived from a minimum of three separate experiments. GraphPad Prism version 9.3.10 (GraphPad Software, USA) was used for graphical representation of the data. Differences between groups were denoted as *p < 0.05, **p < 0.01, and ***p < 0.001.

Results

MicroRNA profiling analysis

To identify and characterize miRNAs associated with prolificacy traits during the follicular and luteal phases in oviducts, we constructed 20 sRNA libraries from prolific and non-prolific Yunshang black goats, comprising 20 individuals across the two phases. Using the Illumina NovaSeq platform, expression profiles of sRNAs were analyzed in oviduct samples. High-throughput sequencing yielded 263 million and 249 million raw reads from oviduct libraries during the follicular and luteal phases, respectively (Supplementary Data 1). After removing low-quality reads, 252 million (95.81% of the total) and 240 million (96.08% of the total) clean reads were retained. On average, 97.44% and 97.56% of these clean reads mapped to the reference genome (Supplementary Data 1). The Q30 quality score for the clean reads ranged from 96.81 to 97.64% in the follicular phase and 96.46 to 97.84% in the luteal phase, while the GC content exceeded 49.09% and 48.69%, respectively (Supplementary Data 1). Length distributions of sRNA (18-35nt) were consistent across the four groups, with most sRNA molecules measuring 21–22 nt, consistent with typical miRNA length characteristics (Fig. 2A). After filtering repeated reads, unique reads in the follicular and luteal phases exceeded 223,460 and 204,779, respectively (Supplementary Data 1). The clean and unique reads were categorized into known miRNA, rRNA, tRNA, snRNA, snoRNA, repeat, novel miRNA, exon, intron, and other categories. RNA class compositions for each group are presented in Fig. 2B (Supplementary Data 1). In this study, the total rRNA proportions in the FH, FL, LH, and LL groups were 3.02,% 2.73%, 2.88%, and 2.88%, respectively, confirming the high quality of the samples (Supplementary Data 1). In the four groups, 71.76%, 70.88%, 76.75%, and 79.43% of the clean reads and 1.28%, 1.46%, 1.54%, and 1.87% of the unique reads were identified as known miRNAs (Fig. 2B). Although known miRNAs accounted for a small portion of all identified ncRNAs, novel miRNAs comprised 0.33%, 0.38%, 0.30%, and 0.41% of the total clean reads and 0.40%, 0.46%, 0.47%, and 0.62% of the unique reads across the FH, FL, LH, and LL groups, respectively (Fig. 2B).

The dataset was compared to miRbase (v22.0) to identify known miRNAs in goat oviducts, detecting 425 mature miRNAs (Supplementary Data 1). Additionally, novel miRNAs were predicted using the miREvo and miRDeep2 algorithms, which map precursor sequences to the genome. A total of 691 potential novel miRNAs were identified (Supplementary Data 1). Predictions were further refined by analyzing secondary structures, predicting Dicer cleavage sites, and calculating binding energy.

Differential expression of oviduct MiRNAs in follicular and luteal phases

As depicted in Fig. 2C, differentially expressed miRNAs were identified using DESeq2 based on raw counts, and results were reported as log₂ fold changes and adjusted p-values. In addition, the miRNA expression levels in each group were quantified and normalized to TPM, showing comparable distributions across the four groups (Supplementary Data 1). TPM values were calculated for visualization and downstream analyses but not used for differential expression analysis. A total of 25 and 66 DEMs were identified in the FL vs. FH (9 upregulated and 16 downregulated) and LL vs. LH (41 upregulated and 25 downregulated) comparisons, respectively. as illustrated in the pie chart (Fig. 2D, Supplementary Data 2), Among these, 4 DEMs (miR-216b, miR-204-5p, miR-204-3p, and novel_92) were co-expressed in both comparisons, while 21 DEMs and 62 DEMs were uniquely expressed in the follicular phase and luteal phase, respectively (Fig. 2E). Cluster analysis of the four groups revealed distinct expression patterns for the DEMs (Fig. 2F, Supplementary Data 2).

Fig. 2

Characterization of microRNA (miRNA) profiling. (A) Length distribution of clean reads from identified miRNA fragments. (B) Categories of identified noncoding RNAs (ncRNAs) via sequencing in FH, FL, LH, and LL groups. (C) The distribution of miRNA expression levels in each group. The expression of miRNAs was normalized by TPM. (D) Differentially expressed miRNAs (DEMs) in FL vs. FH and LL vs. LH comparisons. (E) The Venn diagram shows the intersection of DEMs between two comparisons. (F) Heatmap with expression patterns of DEMs in FH, FL, LH, and LL groups. FL: Follicular phase with low fecundity; FH: Follicular phase with high fecundity; LL: Luteal phase with low fecundity goats; LH: Luteal phase with high fecundity goats

GSEA analysis of the DEMs-targeted genes

Using MiRanda, the prediction of miRNA target genes identified 37,457 transcripts from 18,730 target genes across 1116 miRNAs (Supplementary Data 3). To investigate the critical miRNAs associated with prolificacy and oviduct function, we focused on DEMs identified in the FL vs. FH and LL vs. LH comparisons. By overlapping the DEMs’ target genes with the transcriptome data from our previous study [30, 31] via Venn diagrams, we identified 4,132 and 4,209 target genes for 25 and 66 DEMs, respectively (Fig. 3A and B, Supplementary Data 3). GSEA was performed on these target genes. The enrichment plot of significant KEGG pathways are presented in Fig. 3C and D, while GOBP-based enrichment is shown in Fig. 3E and F (Supplementary Data 3). Several genes sets associated with reproduction and oviduct function were significantly enriched (p-value < 0.05), in the follicular phase, including cytokine-cytokine receptor interaction, cell cycle, hematopoietic cell lineage, tight junction, cell adhesion, cell motility, carbohydrate derivative metabolic process, motile cilium assembly, extracellular transport, and insulin response (Fig. 3C and E). In the luteal phase, pathways such as ABC transporters, valine leucine and isoleucine degradation, cilium organization, carbohydrate derivative metabolic process, TGF-β signaling pathway, cell adhesion, cell activation, regulation of growth, positive regulation of the developmental process, positive regulation of MAPK cascade, and cell motility were enriched (Fig. 3D and F). Notably, the TGF-β signaling pathway, associated with embryo development, showed significant enrichment (NES = − 1.67, Nom p-value = 0.037; Supplementary Data 3).

Analysis of the targeting relationships of miRNA-mRNA co-expression network

In the FL vs. FH comparison, 1,996 DEM-DEG pairs displayed negatively correlated expression (cor < − 0.4, p-value < 0.05; Supplementary Data 4). These pairs included 1,773 DEGs targeted by 15 DEMs, comprising 6 known miRNAs and 9 novel miRNAs. To explore the interactions between DEMs and their target DEGs, a mini-DEM-DEG regulatory network was constructed using Cytoscape software (Fig. 4A). Within this network, miR-204-3p emerged as a core miRNA targeting FGFR2, AKT2, and EPHA2, etc., which are strongly linked to the Ras signaling pathway. In the LL vs. LH comparison, 4325 DEM-DEG pairs with negatively correlated expression (cor < − 0.4, p-value < 0.05) were identified, encompassing 35 DEMs and 2,902 DEGs (Supplementary Data 4). A mini-DEM-DEG interactive network (Fig. 4B) highlighted the putative target genes of miR-328-3p, including SMAD3, and BOP1, etc., which are involved in the negative regulation of cell growth and embryonic development (Supplementary Data 4). Additionally, TOB1, TOB2, FES, LMO7, DGKD, DGKE, DGKH, and CISH were identified as targets of miR-204-3p, playing roles in the negative regulation of cell proliferation, cell adhesion, and the protein kinase C-activating G-protein coupled receptor signaling pathway.

Fig. 3

The gene set enrichment analysis (GSEA) of DEMs-targeted genes. Overlapped genes in FL vs. FH (A) and LL vs. LH (B) comparisons between transcripts and miRNA-targeted genes. GSEA was performed according to the expression of DEMs-targeted genes based on the KEGG in the FL vs. FH comparisons (C) and LL vs. LH comparisons (D). GSEA was performed according to the expression of DEMs-targeted genes based on the Gene Ontology biological process (GOBP) in FL vs. FH (E) and LL vs. LH comparisons (F). NES: Normalized enrichment score

Fig. 4

Overview of miRNA and target genes networks. (A) The DEM-DEG interaction network in the FL vs. FH comparison. (B) The DEM-DEG interaction network in the LL vs. LH comparison. Triangles and circles denote miRNAs and target genes, respectively

Validation of candidate MiRNAs by RT‒qPCR

To confirm the accuracy of the sequencing data, four DEMs were randomly selected from the oviduct tissues during the follicular and luteal phases for expression validation through RT-qPCR analysis. The results demonstrated that the expression levels of miR-204-3p, miR-204-5p, and miR-3955-5p were downregulated (p < 0.05), whereas the expression levels of novel_996 showed an increasing trend (p < 0.05) in high fecundity goats during the follicular phase (Fig. 5A). Similarly, during the luteal phase, miR-10b-5p and miR-328-3p were upregulated (p < 0.05), while miR-195-5p, miR-33a-3p, and miR-34a exhibited decreasing trends (p < 0.05) in high fecundity goats (Fig. 5B). These RT-qPCR results were consistent with the RNA-seq findings, confirming the reliability and accuracy of the miRNA expression level observed in the sequencing data.

Fig. 5

Validation of differentially expressed miRNAs (DEMs) by RT-qPCR in the FL vs. FH comparison (A) and LL vs. LH comparison (B). *P < 0.05, **P < 0.01

LncRNA–ceRNA network analysis and functional annotation

Based on this, we hypothesized that lncRNAs might contribute to the regulation of prolificacy via the ceRNA mechanism. To investigate this, we conducted a comprehensive prediction of lncRNA-miRNA interactions using whole-transcriptome data, identifying a total of 295 and 434 DEL-DEM pairs during the follicular and luteal phases, respectively, based on Pearson correlation coefficients (< − 0.4) and p-value < 0.05 (Supplementary Data 5). Subsequently, we constructed lncRNA − miRNA − mRNA co-regulated networks by combining the DEL − DEM pairs and DEM-DEG pairs under the same criteria (Supplementary Data 5).

During the follicular phase, the ceRNA network comprised 347 known DEGs, 6 DEMs (miR-105a, miR-204-3p, miR-204-5p, miR-3955-5p, miR-767, novel_294), and 136 DELs. To explore the potential biological functions modulated by these target DEGs, GOBP and KEGG pathway enrichment analyses were performed (Supplementary Data 5). Key target genes, such as ATP2B1, ATP2B2, FGFR2, MAPK7, TAB 1, CDH5, ACVRL1, and APC, were enriched in pathways related to cGMP-PKG, MAPK, TGF-β, positive regulation of BMP signaling pathways, and negative regulation of Wnt signaling pathways, all of which are involved in reproductive processes (Fig. 6A). Given the complexity of the ceRNA network, a mini-ceRNA network highlighting critical interactions was constructed. This mini-network included 129 DELs, 5 DEMs, and 48 target genes (Fig. 7A). Similarly, during the luteal phase, the lncRNA-miRNA-mRNA network consisted of 1529 DEGs, 22 DEMs, and 309 DELs. Enrichment analyses identified significant pathways associated with reproductive processes, such as growth hormone synthesis, secretion, and action, ECM-receptor interaction, oxytocin signaling pathway, MAPK signaling pathway, regulation of cell adhesion, TGF-β receptor signaling pathway, and embryonic development (Supplementary Data 5). Notable DEGs involved in embryonic development (PCGF2, PRKCSH, SMAD3, MBD3, and ACVRL1) were also linked to the TGF-β receptor signaling pathway and the negative regulation of cell growth (ACVRL1, SMAD3, CDH5, and ITGB8) (Fig. 6B). To further highlight critical interactions, a mini-ceRNA network was constructed for the luteal phase, comprising 332 lncRNA–miRNA–mRNA interactions, which included 37 DELs, 14 DEMs, and 46 target DEGs (Fig. 7B).

Fig. 6

The ceRNA interaction network analysis. Sankey plot showcasing target gene to the lncRNA–miRNA pairs were involved in each of the enriched pathways or biological processes in the follicular (A) and luteal (B) phases

Fig. 7

The lncRNA–miRNA–mRNA interaction network in the follicular (A) and luteal (B) phases. Hexagon, V, and round rectangle denote circRNA, miRNAs, and target genes, respectively

CircRNA–ceRNA network construction and functional annotation

In this study, we further screened DECs and DEGs regulated by the same DEMs. Ultimately, we identified 19,177 circRNA-miRNA-mRNA interactions in oviduct samples, involving 59 DECs, 50 DEMs, and 7,990 DEMs (Supplementary Data 6). Based on Pearson correlation coefficients (< − 0.4) and p-value < 0.05, 3 DEC-DEM pairs were identified in both follicular and luteal phases (Supplementary Data 6).

During the follicular phase, we identified 1,107 circRNA − miRNA − mRNA interactions, comprising 3 DECs (chi_circ_0021923, chi_circ_0021540, and chi_circ_0020309), 3 DEMs (miR-204-3p, novel_294, and novel_565), and 1036 DEGs (Supplementary Data 6). GOBP and KEGG pathway enrichment analyses revealed several pathways and processes associated with reproduction, including MAPK, TGF-β receptor, negative regulation of Wnt signaling pathways, negative regulation of cell proliferation, and biological processes termed negative regulation of cilium assembly and cell proliferation (Fig. 8A). We constructed a ceRNA regulatory network that included 3 DECs (chi_circ_0021923, chi_circ_0020309, and chi_circ_0021540), 3 DEMs (miR-204-3p, novel-565, and novel-294), and 61 DEGs (including FGFR2, MAPK7, MAPK8IP3, and MAP4K2), which are implicated in pathways critical for reproductive function and cellular regulation (Fig. 8C).

During the luteal phase, we identified 1,253 circRNA–miRNA–mRNA interactions, involving 3 DECs (chi_circ_0014438, chi_circ_0021923, and chi_circ_0005898), 3 DEMs (miR-10b-5p, miR-204-3p, and novel_88), and 1,107 DEGs (Supplementary Data 6). Enrichment analyses highlighted 73 target genes significantly associated with reproductive and embryonic development processes. These genes were enriched in the progesterone receptor signaling pathway, TGF-β receptor signaling, and pathways related to embryonic development, suggesting their potential roles in regulating fecundity (Fig. 8B, Supplementary Data 6). Key genes identified in this phase include TRIM28, ZMIZ2, LMO7, GATA6, PATZ1, ESRRA, SETDB1, XRCC3, TCF3, KDM5C, ZC3H4, DENND1A, CEP164, ADAR, MECP2, DNMT3A, STK40, DOT1L, FOXJ2, and PER1. These genes were regulated by 3 DEC–DEM pairs (chi_circ_0005898, chi_circ_0014438, and chi_circ_0021923, paired with miR-204-3p, miR-10b-5p, and novel_88) and were components of a reconstructed ceRNA regulatory network (Fig. 8D).

Fig. 8

The ceRNA interaction network analysis. Sankey plot showcasing target gene to the circRNA–miRNA pairs were involved in each of the enriched pathways or biological processes in the follicular (A) and luteal (B) phases. The circRNA–miRNA–mRNA interaction network in the follicular (C) and luteal (D) phases; Circles, V, and round rectangle denote circRNA, miRNAs, and target genes, respectively

LNC_005981 and circ_0021923 serves as sponges for miR-328-3p and miR-204-3p, respectively, attenuating their Inhibition of SMAD3 and DOT1L

To further investigate the candidate ceRNA subnetwork contributing to the regulation of goat fecundity traits, specific DECs, and DELs were selected to construct a ceRNA regulatory network hypothesized to influence reproduction. Ultimately, a lncRNA-circRNA–miRNA-mRNA regulative network model was established, comprising 7 miRNAs, 13 lncRNAs, 2 circRNAs, and 14 mRNAs, with a total of 66 interaction relationships (Fig. 9A).

To elucidate the ceRNA mechanism, we selected LNC_005981 − miR-328-3p − SMAD3 and circ_0021923 − miR-204-3p − DOT1L interactions for luciferase assay, as these pathways might impact mammalian embryo development. Using miRanda V3.3a and RNAhybrid, we confirmed the binding relationships between: LNC_005981 and miR-328-3p, miR-328-3p and SMAD3, circ_0021923 and miR-204-3p, and miR-204-3p and DOT1L. To validate these interactions, wild-type (WT) and mutant (MUT) pmirGLO vectors for LNC_005981 and SMAD3 were constructed (Fig. 9B) and co-transfected with miR-328-3p mimics or negative control (NC) into HEK293T cells. The assay revealed that miR-328-3p significantly reduced the luciferase activity of pmirGLO-LNC_005981-WT and pmirGLO-SMAD3-WT, whereas no changes were observed in the luciferase activity of pmirGLO-LNC_005981-MUT and pmirGLO-SMAD3-MUT (Fig. 9D). These results demonstrated that LNC_005981 functions as a sponge for miR-328-3p, and SMAD3 is a direct target of miR-328-3p. To further confirm this ceRNA mechanism, we co-transfected pEGFP-LNC_005981 with miR-328-3p mimics and pmirGLO-SMAD3-WT, observing that LNC_005981 restored luciferase activity in a dose-dependent manner (Fig. 9E). Similarly, WT and MUT pmirGLO vectors for circ_0021923 and DOT1L were also constituted, respectively (Fig. 9C), and co-transfected with miR-204-3p mimics/NC into HEK293T cells. Luciferase activity was significantly reduced in constructs but remained unchanged in their mutant counterparts (pmirGLO-circ_0021923-MUT and pmirGLO-DOT1L-MUT) (Fig. 9F). This demonstrated that circ_0021923 acts as a sponge for miR-204-3p, and DOT1L is a direct target of miR-204-3p. Moreover, overexpression of circ_0021923 restored luciferase activity in a dose-dependent manner (Fig. 8G). Taken together, these results collectively suggest that LNC_005981 and circ_0021923 act as a sponge for miR-328-3p and miR-204-3p, respectively, thereby attenuating the inhibitory effects of these miRNAs on their target genes, SMAD3 and DOT1L. This regulatory mechanism highlights their potential roles in controlling reproductive processes and mammalian embryonic development.

Fig. 9

Verification of the competing endogenous RNA (ceRNA) network. (A) The Sankey diagram shows the ceRNA regulatory network. (B) and (C) Luciferase reporter vectors. (D) The miR-328-3p mimics/NC was co-transfected with pmirGLO-LNC_005981-WT/MUT and pmirGLO-SMAD3-WT/MUT, and luciferase activities were measured after transfection, respectively. (E) The pmirGLO-SMAD3-WT was co-transfected with the miR-328-3p mimic and 1× or 2× pEGFP-LNC_005981 into HEK293T cells, and luciferase activities were measured after transfection, respectively. (F) The miR-204-3p mimics/NC was co-transfected with pmirGLO-circ_0021923-WT/MUT and pmirGLO-DOT1L-WT/MUT, and luciferase activities were measured after transfection, respectively. (G) The pmirGLO-DOT1L-WT was co-transfected with the miR-204-3p mimic and 1× or 2× pEGFP-circ_0021923 into HEK293T cells, and luciferase activities were measured after transfection, respectively. Values are means ± SEM for three individuals. *p < 0.05, **p < 0.01, and **p < 0.001

Discussion

The mammalian oviductal presents a complex internal environment that plays a critical role in gamete maturation, sperm storage, fertilization creation, early embryo development, and embryo transport [3], all of which are key determinants of pregnancy success. Therefore, gaining a comprehensive understanding of oviductal molecular dynamics and regulatory networks is essential for improving the fertility of female goats. MiRNAs, a small class of non-coding RNA in the eukaryotic transcriptome, serve as crucial regulators of various biological processes. A single miRNA that can regulate multiple mRNAs post-transcriptionally, while individual the mRNA can also be targeted and regulated by multiple miRNAs [37]. Reza et al. systematically elucidated the roles of miRNAs in mammalian reproduction [38]. Our findings suggest that the expression patterns of oviductal miRNAs are correlated with prolificacy traits. However, a detailed investigation of oviductal miRNAs to understand their influence on reproductive competence remains limited in goats. In this study, we characterized the miRNA expression profiles at different estrous stages in goat oviducts using miRNA sequencing and bioinformatics analysis. Moreover, we established the whole transcriptome of the oviducts during the follicular and luteal phases in Yunshang black goats, building on our previous studies [30,31,32]. Additionally, we constructed ceRNA regulatory networks, systematically elucidating the molecular mechanisms underlying the influence of the oviduct on goat fecundity. These findings enhance our understanding of the molecular pathways and biological processes that enable the oviduct to sustain high kidding rates.

Several researchers have increasingly focused on the role of miRNAs in regulating gene expression during fertilization and early embryo development [38, 39]. To identify key miRNAs associated with female reproduction in goats, intact oviducts were collected from high- and low-fecundity goats during the follicular and luteal phases. These samples were homogenized for RNA extraction, followed by small RNA (sRNAs) sequencing using the Illumina Novaseq 6000 platform. The sequencing results showed that sRNA lengrhs in both libraries were primarily distributed between 20 and 24 nucleotides. The dominant sRNAs size in the oviduct was 22 nucleotides, which is consistent with the classical length of mammalian miRNAs, and like observations in sheep [40], bovines [41], and pigs [42]. During the follicular phase, miR-124a and miR-376a were significantly differentially expressed in the oviduct, indicating their potential roles in reproductive regulation. Maternal endometrial miRNAs act as transcriptomic modifiers of the peri-implantation endometrium and embryo development, with overexpression of miR-124-3p negatively affecting the migration and proliferation of endometrial cells, thereby impairing endometrial receptivity formation and negatively impacting embryo implantation [16]. Marson et al. [43] noted that miR-124a expression promotes the differentiation of mouse embryonic stem cells (ESCs). Furthermore, the expression of miR-376a is downregulated in mouse blastocysts compared to the morula stage, suggesting its involvement in trophectoderm development [44]. Additionally, miR-376a modulates proliferating cell nuclear antigen (PCNA), thereby reducing oocyte apoptosis, increases the number of primordial follicles, and regulating follicle assembly in mice [45]. MiR-122 was significantly upregulated in the oviducts of high-fecundity goats during the luteal phase. MiR-122 has been reported to play a role in female reproduction by inhibiting testosterone biosynthesis [46]. Reza et al. [38] demonstrated that miR-18a-5p and miR-582-5p are involved in the proliferation and growth of granulosa cells (GCs) and thecal cells (TCs). MiR-21-3p, mediated via extracellular vesicles, has been shown to significantly influence the proliferation and migration of vascular smooth muscle cells [47]. Furthermore, miR-21-3p regulates the apoptosis of ovarian granulosa cells in polycystic ovarian syndrome-related aberrant follicular formation and restores their normal proliferation [48]. These findings indicate that miR-18a-5p, miR-582-5p, miR-34a, and miR-21-3p may play roles in goat oocyte maturation, embryo development, and crosstalk between the oviductal environment and gametes/early embryos, particularly in the oviduct. Although these miRNAs were expressed at low levels in the oviducts of high-fecundity goats, their low abundance may promote processes such as enhancing embryo survival or improving oviduct conditions by reducing the inhibitory effects on target genes, which warrants further investigation.

To better understand the regulation of prolificacy by miRNAs and their target genes, we analyzed the interactions between miRNAs and mRNAs. GSEA of the target genes revealed significant enrichment in pathways such as cytokine-cytokine receptor interaction, cell cycle, tight junction, cell adhesion, cell motility, motile cilium assembly, and insulin response during the follicular phase. Numerous target genes involved in these pathways have been associated with reproduction, including FGFR2 and AKT2, which are regulated by miR-204-3p. Paracrine interactions through FGFR1 and FGFR2 receptors regulate the development of preimplantation mouse chimeric embryos, where downregulation of FGFR1 and FGFR2 expression in 8-cell embryos disrupts intercellular interactions between the components, leading to an inverse proportion of primitive endoderm and epiblast in the inner cell mass (ICM), resulting in abnormal embryo development [49]. Moreover, the maintenance of germ cells in mouse ovaries depends on somatic FGFR2 [50]. Knockdown of FGFR2 in the trophoblast ectoderm of sheep embryos has been shown to retard embryonic development via the MAPK pathway [51]. During mouse embryo development, AKT2, a member of the AKT family, was found to specifically regulate blastomere proliferation, which is indispensable for early blastomere proliferation and embryonic development [52]. These findings suggest that miR-204-3p regulates the processes by targeting FGFR2 and AKT2. Interestingly, during the luteal phase, miR-204-3p also targets genes involved in pathways such as valine leucine, and isoleucine degradation, cilium organization, TGF-β signaling pathway, cell adhesion, cell activation, positive regulation of developmental process, and the MAPK cascade. These genes include TOB1, TOB2, FES, LMO7, DGKD, DGKE, DGKH, and CISH. In mammals, the TOB family comprises antiproliferative proteins, including TOB1 and TOB2, characterized by a conserved C-terminal region. Overexpression of TOT proteins has been reported to inhibit G1-to S-phase cell-cycle progression and reduce cell proliferation in various cell types. TOB1 and TOB2 are involved in early mouse embryonic development and mouse ESCs [53]. Therefore, we hypothesize that miR-204-3p may be involved in embryonic development; however, further studies are needed to elucidate its exact mechanism of action.

The ceRNA theory has been applied to provide a deeper understanding of diverse genetic systems [54]. In this study, we constructed circRNA-miRNA-mRNA and lncRNA-miRNA-mRNA regulated networks for goat oviduct tissues. To explore the potential functions of the ceRNA networks, we conducted functional enrichment analysis on all targeted genes involved in these networks. Notably, these target genes are extensively involved in signaling pathways associated with reproduction, including the MAPK signaling pathway, regulation of cell adhesion, TGF-β receptor signaling pathway, negative regulation of Wnt signaling pathways, and embryonic development. These findings indicate that multiple signaling pathways from a complex regulatory network governing reproduction. From the constructed ceRNA networks, we identified LNC_005981-miR-328-3p-SMAD3 and circ_0021923-miR-204-3p-DOT1L as potential regulators of goat fertility. Subsequently, dual-luciferase report gene assays confirmed the targeted regulatory relationship among these genes. The transforming growth factor-β (TGF-β)-SMAD signaling pathway is critical for regulating multiple aspects of female reproduction. SMAD3 is essential for normal embryonic developmental processes, and both SMAD2 and SMAD3 are required for early embryonic development and for mediating the stimulatory effects of follicular inhibitors on bovine embryos at the 8- to 16-cells stage. of the bovine embryo and that the rate of blastocysts is attenuated when Inhibition of SMAD2/3 signaling has been shown to reduce the blastocyst formation rate [55]. Additionally, DOT1L is predicted to be regulated by miR-204-3p, which has been associated with in high-glucose-induced proapoptotic and dysfunction via the downregulation of Bdkrb2 [56]. Liao et al. [57] demonstrated that DOT1L is critical for mouse embryo development. Homozygous mutations in DOT1L result in various developmental abnormalities and embryonic lethality before 11.5 days of gestation. Maternal DOT1L deposition in the uterus and embryos at the preimplantation stage, as well as H3K79 methylation, is indispensable for mouse development. Furthermore, a developmental study on porcine somatic cell nuclear transfer (SCNT) embryos found that DOT1L inhibitor significantly increased the blastocyst rate and markedly reduced H3K79me2 levels during SCNT 1-cell embryos development, suggesting that DOT1L-mediated H3K79me2 act as a barrier to early development in porcine SCNT embryos [58]. This mechanism may also play a key role in regulating goat embryo development. In this study, LNC_005981 and circ_0021923 were validated as regulators of miR-328-3p and miR-204-3p, respectively. However, the molecular functions of LNC_005981 and circ_0021923 remain unclear. We hypothesize that they may be involved in the regualting goat reproduction. These mechanisms of action need to be further investigated.

Conclusion

In summary, this study explores the regulatory relationship between transcriptome alterations and the prolificacy trait in goat oviduct tissues. By analyzing miRNA expression profiles during the follicular and luteal phases of Yunshang black goats, we identified key miRNAs (e.g., miR-124a, miR-376a, miR-122) and candidate genes (e.g., FGFR2, AKT2, TOB1, TOB2) critical for embryonic development, particularly regulated by miR-204-3p. Additionally, LNC_005981 − miR-328-3p − SMAD3 and circ_0021923 − miR-204-3p − DOT1L were validated, which may play important roles in the reproduction of goats. Collectively, the whole-transcriptome profiling provides several key candidate miRNAs and ceRNAs (DE lncRNAs/DE circRNAs-DE miRNAs-DE genes) that may regulate oocyte maturation, embryo development, and the interactions between the oviduct and gametes/early embryos, providing insights into the molecular mechanisms of reproductive regulatory networks from the perspective of coding and non-coding RNAs.