Differential age-related transcriptomic analysis of ovarian granulosa cells in Kazakh horses

Introduction: The Kazakh horse, renowned for its excellence as a breed, exhibits distinctive reproductive traits characterized by early maturity and seasonal estrus. While normal reproductive function is crucial for ensuring the breeding and expansion of the Kazakh horse population, a noteworthy decline in reproductive capabilities is observed after reaching 14 years of age.

Methods: In this study, ovarian granulosa cells (GCs) were meticulously collected from Kazakh horses aged 1, 2, 7, and above 15 years old (excluding 15 years old) for whole transcriptome sequencing.

Results: The analysis identified and selected differentially expressed mRNAs, lncRNAs, miRNAs, and circRNAs for each age group, followed by a thorough examination through GO enrichment analysis. The study uncovered significant variations in the expression profiles of mRNAs, lncRNAs, miRNAs, and circRNAs within GCs at different stages of maturity. Notably, eca-miR-486-3p and miR-486-y exhibited the highest degree of connectivity. Subsequent GO, KEGG, PPI, and ceRNA network analyses elucidated that the differentially expressed target genes actively participate in signaling pathways associated with cell proliferation, apoptosis, and hormonal regulation. These pathways include but are not limited to the MAPK signaling pathway, Hippo signaling pathway, Wnt signaling pathway, Calcium signaling pathway, Aldosterone synthesis and secretion, Cellular senescence, and NF-kappa B signaling pathway—essentially encompassing signal transduction pathways crucial to reproductive processes.

Discussion: This research significantly contributes to unraveling the molecular mechanisms governing follicular development in Kazakh horses. It establishes and preliminarily validates a differential regulatory network involving lncRNA-miRNA-mRNA, intricately associated with processes such as cell proliferation, differentiation, and apoptosis and integral to the developmental intricacies of stromal follicles. The findings of this study provide a solid theoretical foundation for delving deeper into the realm of reproductive aging in Kazakh mares, presenting itself as a pivotal regulatory pathway in the context of horse ovarian development.

1 Introduction

The Kazakh horse, a vital constituent of the equine family, plays an important role in human activities spanning production, daily life, competitive sports, and recreational pursuits. The distinctive seasonal patterns of estrus in Kazakh horses also impact their reproductive capabilities. Throughout the female reproductive cycle, follicles undergo continuous development, initially growing at a similar rate during each estrous period. Follicles, fundamental units of mammalian ovaries, require normal development to produce oocytes capable of fertilization (1). When a follicle, typically the largest, attains a critical size (approximately 35 millimeters in horses), a phenomenon known as diameter deviation occurs. This involves the continued growth of the dominant follicle (occasionally two or three follicles) while all other follicles (subordinate follicles) cease growth and enter a quiescent state (2). GCs, constituting a primary component of the follicle wall (3, 4), play a pivotal role in regulating follicular development. They not only provide nourishment to oocytes but also exert significant control over follicle growth and development through the secretion of hormones, cytokines, and protein regulation. Furthermore, they regulate the meiotic division process of oocytes and orchestrate follicular atresia through apoptosis. Post-maturation, follicular atresia is primarily induced by GCs (5). To discern the potential connection between granulosa cell proliferation and age (6), a growing body of research suggests that follicular atresia is mainly induced by apoptosis in GCs (7, 8). In horses, GCs stand out as the most crucial cell type within follicles, and their quiescence results in a diminishing reproductive function. Studies affirm that follicular atresia predominantly constitutes a hormone-regulated process involving GC apoptosis (9).

In females, oocytes remain in a state of meiotic arrest for several decades, with cycle lengths being similar (22 and 28 days, each featuring a 14-day luteal phase). Only one follicle ovulates per cycle, and the decline in oocyte quality has been identified as a major age-related factor in infertility (10, 11). Aging initiates a cascade of functional declines and diseases, leading to the cessation of reproductive function in both humans and animals. Given the substantial impact of aging on livestock production, the identification of molecular targets associated with aging holds crucial significance. Mares, due to their ovarian dynamics and reproductive capacity resembling humans, serve as crucial animal models for exploring ovarian aging (12). Research on reproductive lifespan remains an underexplored area. Follicles, as functional units of the reproductive lifespan, comprise oocytes (female gametes) and a specialized group of supporting cells. Unraveling the biology of follicles requires more dedicated efforts. In animal models, biological interventions aimed at manipulating pathways related to aging to maintain the quality and quantity of follicles may extend female reproductive and overall health span. Studies indicate intervention strategies targeting the molecular mechanisms driving ovarian aging and menopause (13). Therefore, investigating the mechanisms of follicular development plays a pivotal regulatory role in the reproductive and developmental processes of horses.

In recent years, with the advancements in biotechnology and functional genomics, whole transcriptomics has emerged as a crucial tool for delving into the characteristics and functions of biological cells. By scrutinizing the complete transcriptome of cells, distinctive gene expression patterns specific to each cell type can be discerned, thereby revealing information intricately linked to cell functions. Consequently, the exploration of the whole transcriptome characteristics and functions of ovarian GCs in Kazakh horses holds paramount significance for gaining a profound understanding of the cellular regulatory mechanisms governing the reproductive and developmental processes in this breed. Presently, there is a notable dearth of research in this specific domain. Thus, the primary objective of this study is to investigate the whole transcriptome of ovarian GCs in Kazakh horses, aiming to unravel their transcriptomic characteristics and the functions associated with them. This research endeavors not only to provide crucial support and guidance for the comprehensive study of reproduction and development in equine species but also holds the promise of offering novel insights and methodologies for genetic enhancement and the advancement of reproductive technologies in equine animals.

2 Materials and methods

2.1 Experimental animals and locations

Experimental Animals: From March 21, 2022, to April 5, 2022, a total of 37 mares were selected from a group of 200 horses. Twenty Kazakh mares were chosen for the study and divided into four groups as follows: Group A - 1 year old (14), five individuals; Group B - 2 years old (15), five individuals; Group C - seven years old (16), five individuals; Group D - above 15 years old (16) (excluding 15 years old), five individuals. The Kazakh horses were sourced from the Tacheng region in Xinjiang, China, and served as the subjects for the research, with experimental samples collected. Throughout the experiment, the mares were housed in individual enclosures in Qiaxia Town, Tacheng City, Xinjiang, China. They were provided with a diet consisting of high-quality dry alfalfa and corn kernels, along with free access to water.

2.2 Collection and preservation of ovarian GCs

From May 8, 2022, to August 8, 2022, ovarian follicles larger than 35mm (Grade III) were identified based on the ovarian follicular development of mares using a B-mode ultrasonic diagnostic instrument. After a 12-hour fasting period, surgical procedures were carried out the next day. The surgery was performed by a licensed veterinarian, assisted by another practicing veterinarian and four assistants. Following the incision of the abdominal wall, follicular fluid was aspirated using a 50 mL syringe and promptly transported to the laboratory for centrifugation. Subsequently, the follicular fluid was discarded, and 2 mL of PBS was added for rinsing before undergoing another round of centrifugation. This process was repeated, and after the final centrifugation, TRIzol was added for preservation. The samples were then stored in liquid nitrogen for future use (See Figure 1).

Figure 1

Figure 1 Schematic diagram of sample collection (All mares subjected to surgery for follicular fluid collection received a 14-day intravenous treatment. Daily disinfection and dressing changes were performed at the surgical site. Please refer to the attached document for information on the drugs used.).

2.3 RNA Extraction, library construction, and sequencing

Total RNA was extracted using a Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s protocol. The quality of the RNA was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and verified through RNase-free agarose gel electrophoresis. Upon RNA extraction, eukaryotic mRNA was enriched using Oligo(dT) beads [For prokaryotes: After total RNA was extracted, prokaryotic mRNA was enriched by removing rRNA by Ribo-ZeroTM Magnetic Kit (Epicentre, Madison, WI, USA)]. Subsequently, the enriched mRNA was fragmented into short fragments using fragmentation buffer and was reverse transcribed into cDNA by using NEBNext Ultra RNA Library Prep Kit for Illumina (NEB #7530, New England Biolabs, Ipswich, MA, USA). The resulting double-stranded cDNA fragments underwent end repair, had an A base added, and were ligated to Illumina sequencing adapters. The ligation reaction was purified with the AMPure XP Beads (1.0X) and subjected to polymerase chain reaction (PCR) for amplification. The resulting cDNA library was sequenced using Illumina Novaseq6000 by Gene Denovo Biotechnology Co. (Guangzhou, China).

2.4 Bioinformatics analysis

2.4.1 Filtering of clean reads

Reads obtained from the sequencing machines include raw reads containing adapters or low-quality bases, which can adversely impact subsequent assembly and analysis. To obtain high-quality clean reads, we further filtered the reads using fastp (17) (version 0. 18.0).

2.4.2 Alignment with ribosome RNA

The short-read alignment tool Bowtie2 (18) (version 2.2.8) was employed to map reads to the ribosome RNA (rRNA) database. The reads mapped to rRNA were subsequently removed. The remaining clean reads were further used for assembly and gene abundance calculation.

2.4.3 Alignment with reference genome

An index of the reference genome was constructed, and paired-end clean reads were mapped to the reference genome using HISAT2. 2.4 (19), with other parameters set to default.

2.4.4 Quantification of gene abundance

The mapped reads of each sample were assembled by using StringTie v1.3. 1 (20, 21) in a reference-based approach. For each transcription region, an FPKM (fragment per kilobase of transcript per million mapped reads) value was calculated to quantify its expression abundance and variations using RSEM (22) software.

The FPKM formula is shown as follows:

$$\text{FPKM}=\frac{{10}^6\text{C}}{\text{NL}/{10}^3}$$

Where FPKM(A) represents the expression of gene A, C represents the number of fragments mapped to gene A, N represents the total number of fragments mapped to reference genes, and L represents the number of bases on gene A. The FPKM method effectively eliminates the influence of different gene lengths and sequencing data amounts on the calculation of gene expression. Consequently, the calculated gene expression can be directly used to compare differences in gene expression among samples.

2.5 Relationship analysis of samples

2.5.1 Correlation analysis of replicas

Correlation analysis was performed using R to assess the reliability and operational stability of two parallel experiments. The correlation coefficient between replica samples was calculated to evaluate repeatability. A correlation coefficient approaching 1 indicates higher repeatability between the two parallel experiments.

2.5.2 Principal component analysis

Principal component analysis (PCA) was performed in this study using the R package gmodels (http://www.rproject.org/). PCA is a statistical procedure that converts hundreds of thousands of correlated variables (gene expression) into a set of values of linearly uncorrelated variables known as principal components. PCA is widely employed to elucidate the structure/relationship among the samples/data.

2.6 Differentially expressed genes

Differential expression analysis of RNAs was performed using DESeq2 (23) software between different groups (and by edgeR (24) between two samples). Genes/transcripts with a false discover rate FDR parameter below 0.05 and an absolute fold change≥2 were considered differentially expressed genes/transcripts.

2.6.1 GO enrichment analysis

Firstly, all DEGs were mapped to GO terms in the Gene Ontology database (http://www.geneontology.org/), and gene numbers were calculated for every term. Significantly enriched GO terms in DEGs, compared to the genome background, were determined using the hypergeometric test (25). The calculating formula for the P-value is:

$$P=1-{\displaystyle \sum _{i=0}^{m-1}\frac{\left(\begin{array}{c}M\\ i\end{array}\right)\left(\begin{array}{c}N-M\\ n-i\end{array}\right)}{\left(\begin{array}{c}N\\ n\end{array}\right)}}$$

Where N represents the total number of genes with GO annotation; n represents the number of DEGs in N; M represents the total number of genes annotated to certain GO terms; and m represents the number of DEGs in M. The calculated p-value underwent False Discovery Rate (FDR) correction, with an FDR ≤ 0.05 threshold. GO terms meeting this condition were defined as significantly enriched GO terms in DEGs. This analysis helped identify the primary biological functions exercised by DEGs.

2.6.2 Pathway enrichment analysis

KEGG, a prominent public pathway-related database, was employed for pathway enrichment analysis, identifying significantly enriched metabolic or signal transduction pathways in DEGs compared with the whole genome background (26). The calculating formula is the same as that in GO analysis.

$$P=1-{\displaystyle \sum _{i=0}^{m-1}\frac{\left(\begin{array}{c}M\\ i\end{array}\right)\left(\begin{array}{c}N-M\\ n-i\end{array}\right)}{\left(\begin{array}{c}N\\ n\end{array}\right)}}$$

Where N represents the total number of genes with KEGG annotation, n represents the number of DEGs in N, M represents the number of all genes annotated to specific pathways, and m represents the number of DEGs in M. The calculated p-value underwent FDR correction, with an FDR ≤ 0.05 threshold. Pathways meeting this condition were defined as significantly enriched pathways in DEGs.

2.7 Protein-protein interaction

The Protein-Protein Interaction network was constructed using String v10 (27), designating genes as nodes and interaction as lines. The resulting network file was visualized using Cytoscape (v3.7.1) (28) software to depict core and hub gene biological interactions.

2.8 Single-nucleotide polymorphism analysis

The GATK (29) (version 3.4-46) was used for variant calling of transcripts, and ANNOVAR was used for SNP/InDel annotation. The function, genome site, and type of variation of SNPs were also analyzed.

2.9 Validation through RT-qPCR

To confirm the accuracy of the sequencing data, we randomly chose a set of six differentially expressed mRNAs for fluorescent quantitative testing using RT-qPCR. The reference genes employed for this validation were GADPH and β-Actin. For detailed primer sequences, please refer to Supplementary Materials (see Supplementary Table S1 for information on differentially expressed mRNA primers, materials and methods).

3 Results

3.1 Comparison of clean reads to the reference genome

This study reveals that, through transcriptome sequencing of GCs from 20 Kazakh mares, 50.73% to 92.82% of clean reads were successfully aligned to the reference genome (Ensembl_release100). Based on the alignment of transcriptome data with the reference genome, a total of 291 new genes were annotated, and the SNP count ranged from 1,053,699 to 2,281,283 (please refer to Supplementary Table S2).

3.2 Identification and analysis of differentially expressed mRNAs

In the Group A vs Group B comparison, three mRNAs were identified (one upregulated and two downregulated): SNRPN, RPS17, and LRP2. For the Group A vs Group C comparison, twelve mRNAs were identified (eight upregulated and four downregulated), including NID1, SNRPN, and COL1A1, among others. In the Group A vs Group D comparison, a total of 186 mRNAs were identified (84 upregulated and 102 downregulated), such as LARGE2, CDH24, TENM4, and others. In the Group B vs Group C comparison, one mRNA was identified (no upregulation and one downregulation), specifically MMP12. For the Group B vs Group D comparison, a total of 80 mRNAs were identified (18 upregulated and 62 downregulated), including IHH, TENM4, and CDH24, among others. Lastly, in the Group C vs Group D comparison, 992 mRNAs were identified (567 upregulated and 425 downregulated), such as TENM4, SLC39A10, CTNNAL1, and others. Group C vs Group D exhibited the highest number of differentially expressed mRNAs, as depicted in Figure 2.

Figure 2

Figure 2 Statistics of differentially expressed mRNAs in each comparison group. (A) mRNAs Venn diagram of differences between groups; (B) mRNAs statistical histogram of differences among groups.

3.3 Identification and classification of lncRNA

The new LncRNAs have been categorized into five major classes based on their relative positions in the genome compared to protein-coding genes: Intergenic Long Non-Coding RNAs (Intergenic LncRNAs): 12,553; Bidirectional Long Non-Coding RNAs (Bidirectional LncRNAs): 465; Intronic Long Non-Coding RNAs (Intronic LncRNAs): 54; Antisense Long Non-Coding RNAs (Antisense LncRNAs): 935; Sense Overlapping Long Non-Coding RNAs (Sense Overlapping LncRNAs): 595. Transcript reconstruction using stringtie and prediction of transcript coding potential through CPC2, CNCI, and Feelnc software resulted in reliable predictions, as depicted in Figure 3.

Figure 3

Figure 3 LncRNA identification and analysis. (A) Venn diagram illustrating lncRNA coding potential; (B) Bar Chart displaying lncRNA types.

In Group A vs Group B, one lncRNA (MSTRG.14266.1) was identified, showing no upregulation and one downregulation. In Group A vs Group C, two lncRNAs (ENSECAT00000071442 and MSTRG.4012.8) were identified, both exhibiting upregulation and no downregulation. In Group A vs Group D, three lncRNAs (MSTRG.40650.1, ENSECAT00000049786, MSTRG.49929.4) were identified, with two upregulated and one downregulated. In Group B vs Group C, one lncRNA (ENSECAT00000071442) was identified, showing upregulation and no downregulation. Regarding Group B vs Group D, a total of 17 lncRNAs (e.g., MSTRG.24142.1, MSTRG.44109.1, MSTRG.14314.1) were identified, with 14 upregulated and three downregulated. In Group C vs Group D, a total of 81 lncRNAs (e.g., MSTRG.3596.12, MSTRG.39723.1, MSTRG.10102.1) were identified, with 74 upregulated and seven downregulated. The comparison of Group C vs Group D exhibited the highest number of differentially expressed lncRNAs, as depicted in Figure 4.

Figure 4

Figure 4 Statistics of differentially expressed LncRNAs in each comparison group. (A) LncRNAs Venn diagram of differences between groups; (B) LncRNAs statistical histogram of differences among groups.

3.4 Identification and analysis of differentially expressed miRNAs

In the Group A vs B group comparison, a total of 168 miRNAs were identified (92 upregulated and 76 downregulated), including MSTRG.10102.1, eca-miR-146a, eca-miR-143, and others. For Group A vs Group C, 115 miRNAs were identified (57 upregulated and 58 downregulated), such as eca-miR-146a, novel-m0007-3p, eca-miR-450b-5p, and others. In Group A vs Group D, 152 miRNAs were identified (101 upregulated and 51 downregulated), including eca-miR-199a-5p, eca-miR-143, eca-miR-199a-3p, and others. In Group B vs Group C, 76 miRNAs were identified (24 upregulated and 52 downregulated), including eca-miR-451, eca-miR-199a-3p, eca-miR-199b-3p, and others. In Group B vs Group D, a total of 143 miRNAs were identified (82 upregulated and 61 downregulated), such as miR-99-z, eca-miR-142-5p, miR-486-x, and others. Lastly, in Group C vs Group D, 134 miRNAs were identified (84 upregulated and 50 downregulated), including eca-miR-199a-5p, eca-miR-20a, and miR-2478-z, as depicted in Figure 5.

Figure 5

Figure 5 Statistics of differentially expressed miRNAs in each comparison group. (A) miRNAs Venn diagram of differences between groups; (B) miRNAs statistical histogram of differences among groups.

3.5 Identification and analysis of differentially expressed circRNAs

In the Group A vs Group B comparison, a total of 104 circRNAs were identified (45 upregulated and 59 downregulated), including TEX14, PCMT1, SLAMF6, and others. For Group A vs Group C, 134 circRNAs were identified (68 upregulated and 66 downregulated), such as MERTK, EIF4B, PSD3, and others. In Group A vs Group D, 141 circRNAs were identified (87 upregulated and 54 downregulated), including SEC24B, ACTB, THBS1, and others. Regarding the Group B vs Group C comparison, 85 circRNAs were identified (40 upregulated and 45 downregulated), including PPM1D, MYO6, USP7, and others. In Group B vs Group D, a total of 155 circRNAs were identified (95 upregulated and 60 downregulated), including CPSF6, MDN1, GNS, and others. Lastly, in Group C vs Group D, 176 circRNAs were identified (96 upregulated and 80 downregulated), featuring VCAN, MYO1D, and STAG1, as depicted in Figure 6.

Figure 6

Figure 6 Statistics of differentially expressed circRNAs in each comparison group. (A) circRNAs Venn diagram of differences between groups; (B) circRNAs statistical histogram of differences among groups.

It is noteworthy that there are 23 differential genes for mRNA in Group A vs Group D, Group B vs Group D, and Group C vs Group D, respectively. For lncRNA, there is one differential gene in each of the Group A vs Group D, Group B vs Group D, and Group C vs Group D comparisons. Additionally, miRNA exhibits 23 differential genes in Group A vs Group D, Group B vs Group D, and Group C vs Group D, respectively. As for circRNA, there are 19 differential genes in Group A vs Group D, Group B vs Group D, and Group C vs Group D, respectively. These shared differentially expressed mRNA, lncRNA, miRNA, and circRNA can be further investigated as key regulators in the aging process of mare reproduction, influencing age-related changes in ovarian function and follicular development, as illustrated in Figure 7.

Figure 7

Figure 7 Venn diagrams of: (A) the shared differentially expressed mRNAs in Group A vs Group D, Group B vs Group D, and Group C vs Group D; (B) the shared differentially expressed lncRNAs in Group A vs Group D, Group B vs Group D, and Group C vs Group D; (C) the shared differentially expressed miRNAs in Group A vs Group D, Group B vs Group D, and Group C vs Group D; (D) the shared differentially expressed circRNAs in Group A vs Group D, Group B vs Group D, and Group C vs Group D.

3.6 Enrichment analysis of differentially expressed mRNAs

In the Group A vs Group B, Group A vs Group C, Group A vs Group D, Group B vs Group C, Group B vs Group D, and Group C vs Group D comparisons, 2, 9, 146, 1, 55, and 808 differentially expressed mRNAs were respectively annotated in the GO database. In Group A vs Group B, significant enrichments were observed in U2-type prespliceosome (Cellular Component - CC), U5 snRNP (CC), chaperone binding (Molecular Function - MF), positive regulation of oligodendrocyte progenitor proliferation (Biological Process - BP), and regulation of oligodendrocyte progenitor proliferation (BP). In Group A vs Group C, major enrichments were identified in proteinaceous extracellular matrix (CC), Wnt-protein binding (MF), acetyltransferase activator activity (MF), laminin-1 binding (MF), and Wnt signaling pathway (BP). In Group A vs Group D, enrichments were prominent in cell junction (CC), germinal vesicle (CC), protein homodimerization activity (MF), protein binding (MF), cellular component organization (BP), and cellular component organization or biogenesis (BP). In Group B vs Group C, notable enrichments were found in extracellular matrix (CC), metalloendopeptidase activity (MF), negative regulation of hydrogen peroxide metabolic process (BP), and regulation of hydrogen peroxide metabolic process (BP). In Group B vs Group D, major enrichments included cell surface (CC), cell junction (CC), calcium ion binding (MF), vascular endothelial growth factor receptor 2 binding (MF), and cell adhesion molecule binding (MF). In Group C vs Group D, enrichments were predominantly observed in cell surface (CC), membrane (CC), endomembrane system (CC), protein binding (MF), binding (MF), cytoskeletal protein binding (MF), positive regulation of response to stimulus (BP), and regulation of response to stimulus (BP) (see Figure 8).

Figure 8

Figure 8 GO enrichment analysis of differentially expressed mRNAs. (A) Group A vs Group B; (B) Group A vs Group C; (C) Group A vs Group D; (D) Group B vs Group C; (E) Group B vs Group D; (F) Group C vs Group D.

Bar chart of GO enrichment classification: (The horizontal axis represents secondary GO terms, and the vertical axis represents the number of differentially expressed genes in that term, with red indicating upregulation and green indicating downregulation; the same applies below).

In the Group A vs Group B, Group A vs Group C, Group A vs Group D, Group B vs Group C, Group B vs Group D, and Group C vs Group D comparisons, differentially expressed mRNAs were significantly enriched in 4, 0, 1, 5, 0, and 15 signaling pathways, respectively (q-value < 0.05). In Group A vs Group B, the differentially enriched pathways include Hedgehog signaling pathway, Thyroid hormone synthesis, Cholesterol metabolism, and Ribosome. In Group A vs Group D, the differentially enriched pathway is Pathways in cancer. In Group B vs Group C, differentially enriched pathways include the IL-17 signaling pathway, Prostate cancer, TNF signaling pathway, Rheumatoid arthritis, and Transcriptional misregulation in cancer. In Group C vs Group D, differentially enriched pathways include Hematopoietic cell lineage, Phagosome, Leishmaniasis, Hippo signaling pathway - fly, Tight junction, Leukocyte transendothelial migration, Hippo signaling pathway, Hippo signaling pathway - multiple species, Toll-like receptor signaling pathway, Hedgehog signaling pathway, Platelet activation, Longevity regulating pathway, Pathways in cancer, Regulation of lipolysis in adipocytes, and MAPK signaling pathway (refer to Figure 9).

Figure 9

Figure 9 KEGG enrichment analysis of differentially expressed mRNAs. (A) Group A vs Group B; (B) Group A vs Group C; (C) Group A vs Group D; (D) Group B vs Group C; (E) Group B vs Group D; (F) Group C vs Group D.

KEGG enrichment bubble chart: (Constructed using the 20 pathways with the lowest Q-values for visualization. The vertical axis represents pathways, the horizontal axis represents the enrichment factor (the ratio of differentially expressed genes in that pathway to the total number in that pathway), the bubble size indicates quantity, and the color becomes deeper red as the Q-value decreases; the same principles apply below).

3.7 Enrichment analysis of differentially expressed lncRNA target genes

Predictions for target genes of differentially expressed lncRNAs between groups were made. Subsequently, GO functional and KEGG pathway enrichment analyses were performed on these target genes to delve into the roles of the identified differentially expressed lncRNAs during the ovarian follicular development of Kazakh horses at various ages. Cis-target gene predictions were based on the positional relationship between lncRNAs and mRNAs, considering mRNAs detected within a specified range upstream and downstream of lncRNAs as their target genes. Results showed that within the 10 kb upstream and downstream range of lncRNAs, a total of 3067 lncRNAs corresponded to 3006 protein-coding genes, forming 3972 lncRNA-mRNA pairs. Additionally, trans-expression target gene prediction, based on the correlation between lncRNAs and mRNAs (with a screening criterion of an absolute value of Pearson correlation coefficient greater than 0.95), identified 1235 potential regulatory relationships between 50 lncRNAs and 252 mRNA targets.

To gain a more comprehensive understanding of the functional roles of these differentially expressed lncRNAs, GO functional enrichment and KEGG pathway analyses were separately conducted for cis and trans target genes. The results revealed that under the cis-regulatory mode, the target genes of differentially expressed lncRNAs were significantly enriched in 138 GO terms (P<0.05, Figures 2–9). These terms included 123 BP, 9 MF, and 6 CC. Notably, some of the top-ranking GO terms included “positive regulation of interferon-alpha biosynthetic process” (BP), “negative regulation of B cell receptor signaling pathway” (BP), “cellular response to external stimulus” (BP), “regulation of interferon-alpha biosynthetic process” (BP), “double-stranded RNA binding” (MF), “eukaryotic translation initiation factor 2alpha kinase activity” (MF), “siRNA binding” (MF), “pattern recognition receptor activity” (MF), “kinase activity” (MF), “early phagosome” (CC), “endosome” (CC), “phagocytic vesicle” (CC), “recycling endosome” (CC). In addition, a total of KEGG pathways showed significant enrichment (P<0.05, Figure 10). Among the top-ranking KEGG pathways were “Influenza A,” “Measles,” and “Glycerolipid metabolism” signaling pathways. Notably, enrichment was also observed in some common signaling pathways associated with testicular tissue development or spermatogenesis, such as “Cell adhesion molecules” and “Calcium signaling pathway” (see Figure 10).

Figure 10

Figure 10 (A) Antisense GO enrichment analysis of differentially expressed lncRNAs target gene. (B) Antisense KEGG enrichment analysis of differentially expressed lncRNAs target gene. (C) CIS GO enrichment analysis of differentially expressed lncRNAs target gene. (D) CIS KEGG enrichment analysis of differentially expressed lncRNAs target gene. (E) Trans GO enrichment analysis of differentially expressed lncRNAs target gene (F) Trans KEGG enrichment analysis of differentially expressed lncRNAs target gene.

Under the trans-regulatory mode, the target genes of differentially expressed lncRNAs were notably enriched in 138 GO terms (P<0.05, Figure 10). These included 123 BP, 9 MF, and 6 CC. Important GO terms consisted of “positive regulation of interferon-alpha biosynthetic process” (BP), “negative regulation of B cell receptor signaling pathway” (BP), “cellular response to external stimulus” (BP), “interferon-alpha biosynthetic process” (BP), “regulation of interferon-alpha biosynthetic process” (BP), “double-stranded RNA binding” (MF), “eukaryotic translation initiation factor 2alpha kinase activity” (MF), “siRNA binding” (MF), “pattern recognition receptor activity” (MF), “transferase activity, transferring phosphorus-containing groups” (MF), “early phagosome” (CC), “endosome” (CC), “phagocytic vesicle” (CC), “recycling endosome” (CC), “endocytic vesicle” (CC), and others. Moreover, a total of 23 KEGG signaling pathways were identified, including “Influenza A,” “Measles,” and “Glycerolipid metabolism,” among others (refer to Figure 10).

3.8 Enrichment analysis of differentially expressed miRNA target genes

For Group A vs Group B, Group A vs Group C, Group A vs Group D, Group B vs Group C, Group B vs Group D, and Group C vs Group D, there were 13,363, 12,551, 13,056, 11,287, 12,854, and 12,730 annotated differentially expressed miRNA target genes in the GO database, respectively. In Group A vs Group B, differentially expressed miRNA target genes were enriched in “cell cortex” (CC), “cytoplasm” (CC), “endomembrane system” (CC), “binding” (MF), “protein binding” (MF), “carbohydrate derivative binding” (MF), “single-multicellular organism process” (BP), “regulation of biological quality” (BP), “enzyme-linked receptor protein signaling pathway” (BP). In Group A vs Group C, differentially expressed miRNA target genes were enriched in “binding” (MF), “protein binding” (MF), “carbohydrate derivative binding” (MF), “single-multicellular organism process” (BP), “regulation of biological quality” (BP), “enzyme-linked receptor protein signaling pathway” (BP). In Group A vs Group D, differentially expressed miRNA target genes were enriched in “cytoplasm” (CC), “endomembrane system” (CC), “cell cortex” (CC), “binding” (MF), “protein binding” (MF), “Ras GTPase binding” (MF), “single-multicellular organism process” (BP), “regulation of biological quality” (BP), “cellular response to organic substance” (BP). In Group B vs Group C, differentially expressed miRNA target genes were enriched in “intracellular” (CC), “intracellular part” (CC), “cytoplasm” (CC), “binding” (MF), “protein binding” (MF), “growth factor binding” (MF), “single-multicellular organism process” (BP), “enzyme-linked receptor protein signaling pathway” (BP), “regulation of biological quality” (BP). In Group B vs Group D, differentially expressed miRNA target genes were enriched in “cell cortex” (CC), “focal adhesion” (CC), “cytoplasm” (CC), “protein binding” (MF), “binding” (MF), “Ras GTPase binding” (MF), “single-multicellular organism process” (BP), “system development” (BP), “multicellular organism development” (BP). In Group C vs Group D, differentially expressed miRNA target genes were enriched in “cell cortex” (CC), “focal adhesion” (CC), “cytoplasm” (CC), “protein binding” (MF), “binding” (MF), “Ras GTPase binding” (MF), “single-multicellular organism process” (BP), “anatomical structure development” (BP), “system development” (BP) (refer to Figure 11).

Figure 11

Figure 11 GO enrichment analysis of differentially expressed miRNAs target gene. (A) Group A vs Group B; (B) Group A vs Group C; (C) Group A vs Group D; (D) Group B vs Group C; (E) Group B vs Group D; (F) Group C vs Group D.

In comparisons between Group A vs Group B, Group A vs Group C, Group A vs Group D, Group B vs Group C, Group B vs Group D, and Group C vs Group D, the target genes of differentially expressed miRNAs are significantly enriched in 137, 135, 140, 118, 113, and 113 pathways, respectively (q-value < 0.05). Notably, the enriched pathways vary among the comparisons: In Group A vs Group B, enriched pathways include the Ras signaling pathway, Regulation of actin cytoskeleton, and Graft-versus-host disease. In Group A vs Group C, enriched pathways involve the Ras signaling pathway, Regulation of actin cytoskeleton, and Graft-versus-host disease. In Group A vs Group D, enriched pathways encompass the Ras signaling pathway, Regulation of actin cytoskeleton, and Graft-versus-host disease. In Group B vs Group C, enriched pathways consist of the Ras signaling pathway, Regulation of actin cytoskeleton, and Proteoglycans in cancer. In Group B vs Group D, enriched pathways include Regulation of actin cytoskeleton, Graft-versus-host disease, Viral myocarditis, and ECM-receptor interaction. In Group C vs Group D, enriched pathways involve the Ras signaling pathway, Regulation of actin cytoskeleton, Cell adhesion molecules, PI3K-Akt signaling pathway, Calcium signaling pathway, AMPK signaling pathway, MAPK signaling pathway, and GnRH signaling pathway. (refer to Figure 12).

Figure 12

Figure 12 KEGG enrichment analysis of differentially expressed miRNAs target gene. (A) Group A vs Group B; (B) Group A vs Group C; (C) Group A vs Group D; (D) Group B vs Group C; (E) Group B vs Group D; (F) Group C vs Group D.

3.9 Enrichment analysis of differentially expressed circRNA target genes

The circRNAs are generated through the reverse splicing of exons from protein-coding genes. Annotating the source genes can offer insights into the functions of circRNAs to some extent. In the comparisons between Group A vs Group B, Group A vs Group C, Group A vs Group D, Group B vs Group C, Group B vs Group D, and Group C vs Group D, 83, 105, 103, 65, 109, and 133 differentially expressed circRNA source genes, respectively, were annotated in the GO database. In Group A vs Group B, differentially expressed circRNA source genes are primarily annotated for intracellular (CC), intracellular part (CC), cytoplasmic ribonucleoprotein granule (CC), enzyme binding (MF), high-density lipoprotein particle binding (MF), protein binding (MF), and endothelial cell migration (BP). In Group A vs Group C, differentially expressed circRNA source genes are mainly annotated for macromolecular complex (CC), myosin complex (CC), protein complex (CC), motor activity (MF), transforming growth factor beta-activated receptor activity (MF), protein binding (MF), antigen processing and presentation via MHC class Ib (BP), and antigen processing and presentation of peptide antigen via MHC class Ib (BP). In Group A vs Group D, differentially expressed circRNA source genes are primarily annotated for intracellular (CC), intracellular part (CC), protein complex (CC), binding (MF), lipoprotein particle binding (MF), protein-lipid complex binding (MF), and regulation of catabolic process (BP). In Group B vs Group C, differentially expressed circRNA source genes are mainly annotated for nucleoplasm part (CC), ion binding (MF), binding (MF), kinase binding (MF), antigen processing, and presentation of endogenous peptide antigen via MHC class Ib (BP), and antigen processing and presentation of peptide antigen via MHC class Ib (BP). In Group B vs Group D, differentially expressed circRNA source genes are primarily annotated for intracellular (CC), intracellular part (CC), nuclear inclusion body (CC), clathrin heavy chain binding (MF), ion binding (MF), macromolecular complex binding (MF), and membrane organization (BP). In Group C vs Group D, differentially expressed circRNA source genes are mainly annotated for intracellular (CC), intracellular part (CC), protein complex (CC), ATP binding (MF), adenyl ribonucleotide binding (MF), antigen processing and presentation of endogenous peptide antigen via MHC class Ib (BP), and regulation of metabolic process (BP) (refer to Figure 13).

Figure 13

Figure 13 GO enrichment analysis of differentially expressed circRNAs target gene. (A) Group A vs Group B; (B) Group A vs Group C; (C) Group A vs Group D; (D) Group B vs Group C; (E) Group B vs Group D; (F) Group C vs Group D.

The KEGG analysis results for differentially expressed circRNA host genes indicate enrichment in 1, 14, and 2 pathways for Group A vs Group B, Group B vs Group C, and Group C vs Group D, respectively (q-value<0.05). No enrichment was observed in the other groups. In Group A vs Group B, the differentially expressed host genes were enriched in the Viral carcinogenesis pathway. In Group B vs Group C, the enriched pathways included Viral carcinogenesis, Phagosome, Human T-cell leukemia virus type 1 infection, Epstein-Barr virus infection, Graft-versus-host disease, and Cellular senescence, among others. In Group C vs Group D, the enriched pathways comprised Phagosome and Endocytosis (refer to Figure 14).

Figure 14

Figure 14 KEGG enrichment analysis of differentially expressed circRNAs target gene. (A) Group A vs Group B; (B) Group A vs Group C; (C) Group A vs Group D; (D) Group B vs Group C; (E) Group B vs Group D; (F) Group C vs Group D.

3.10 Protein-protein interaction

1, 1, 105, 0, 35, and 765 pairs of protein interactions were identified in Group A vs Group B, Group A vs Group C, Group A vs Group D, Group B vs Group D, and Group C vs Group D, respectively. In the protein interaction network of Group A vs Group D, MCMDC2, KIT, ABL1, ACTA2, RHOJ, RASL11A, PNPLA6, AURKA, MMP2, CALM2 were pinpointed as core nodes. For Group B vs Group D, IHH, KALRN, NCALD, VWCE, ALDH3A2, CORIN, IRS2, TNF, GLIS2, and STK32A emerged as core nodes in the protein interaction network. In Group C vs Group D, MCMDC2, RPS6KA3, MAPK4, PTPRC, BTK, FLT3, RRAD, TNF, KALRN, RAP1B, and others were recognized as core nodes in the protein interaction network (please refer to Supplementary Figures S1–S3).

3.11 ceRNA network

Finally, a predicted ceRNA network comprised 58 mRNA pairs, 12 lncRNA pairs, and 8 circRNA pairs. For the top 10 connected mRNA/lncRNA/circRNA, a Sankey diagram illustrating their regulatory relationships with miRNAs was constructed. From the diagram, it is apparent that eca-miR-486-3p and miR-486-y exhibit the highest connectivity (please refer to Supplementary Table S3, Figures S4, S5).

3.12 RT-qPCR

The consistent trends in the fluorescence quantification results of the six mRNAs with the RNA-Seq data, as depicted in Figure 15, affirm the high reliability of the sequencing results.

Figure 15

Figure 15 RT-qPCR verification of differently expressed mRNAs.

4 Discussion

Ovaries and uterus, as essential reproductive organs, undergo the development and ovulation of dominant follicles on the ovaries, fertilization in the oviduct, and subsequent development in the uterus, leading to parturition under the influence of reproductive hormones. Transcriptomic studies related to ovaries and follicles have been applied to various species, such as humans (30), cows (17, 31, 32), pigs (33, 34), sheep (35), chickens (36), and geese (37). However, research on the ovaries of equine animals, especially horses, is relatively scarce. Horses, as typical seasonal breeders, experience estrus and ovulation influenced by various factors, including age, environment, and nutritional status. Reproductive aging in mares starts around 15 years old, and exceptional mares are essential for breeding high-quality racehorses, exemplified by horses like Dahlia, Allez France, and Oh So Sharp. Therefore, investigating the reproductive mechanism of mares using molecular biology techniques becomes crucial. The Kazakh horse exhibits evident seasonal estrus, occurring 2-3 times per year, with each estrus period lasting for 21 days. According to previous research by our team, mares aged 7-14 years in moderate-to-good body condition and exposed to a temperature range of 18-24°C show optimal estrus conditions. After 15 years of age, reproductive functions begin to noticeably decline, with a distinct pattern of follicular development and higher fertility compared to other age groups. This characteristic makes them an ideal animal model for studying mare follicular development. In-depth research has been conducted globally on the endocrine regulation of ovarian follicular development in equine animals, including mares. However, the specific molecules and cellular processes involved in the regulation of follicular growth and development at different ages in mares remain to be further explored.

In the comparison between groups A, B, C, and D in this study, 22 differentially expressed mRNAs exhibited significant downregulation, while FCER2 (ENSECAG00000009646) demonstrated a notable upregulation. Although existing research on FCER2 predominantly focuses on cancer, it is hypothesized that FCER2 plays a role in the regulation of ovarian aging. The varied expression states of these differentially expressed genes may regulate the proliferation and apoptosis of GCs, leading to follicular atresia or ovarian dysfunction. Among the 23 shared differentially expressed genes, AUTS2, ZFHX2, MAST4, SLC24A3, B3GNT6, and PPP1R13L showed a downregulation, while in Group A vs Group B, Group A vs Group C, and Group B vs Group C, they were upregulated. This suggests that AUTS2 is involved in regulating RNA transcription, splicing, localization, and stability. Follicular development and ovulation are highly complex biological processes, and lncRNA has emerged as a major translational regulatory factor in reproductive processes for humans and animals, such as gonadal development (38) and sex hormone secretion (39). GCs play a highly complex role in oocyte function, involving the regulation of multiple factors. This study suggests that lncRNA plays an important regulatory role in the follicular development and hormonal regulation of mares at different stages before sexual maturity, physical maturity, youth, and old age. By analyzing the lncRNA expression profiles, we identified differentially expressed genes in GCs at different ages in Kazakh horses. It is noteworthy that MSTRG.40650.1 showed significant upregulation, and its expression changed with age, suggesting that MSTRG.40650.1 may be involved in the regulation of aging. In Group C vs Group D, there were the most differentially expressed genes, with 74 upregulated and seven downregulated. MiRNA is a class of endogenous RNAs, approximately 18–24 nucleotides in length, which act by specifically binding to the 3’-UTR region of mammalian mRNAs to inhibit or reduce their translation, thereby post-transcriptionally regulating gene expression (40). Previous studies have shown that miRNAs play crucial roles in cell proliferation, differentiation, growth, migration, and apoptosis (41–44). MiRNAs may inhibit the expression of proteins involved in the proliferation and development of reproductive cells at the post-transcriptional level (45). Eca-miR-34b-3p, miR-12245-y, and miR-21-z show distinct expression patterns when compared under the same conditions, providing substantial evidence for their crucial regulatory roles in the sexual gland development of Kazakh horses. The CircRNA results indicate significant age-related features in the expression profile of GCs in mare ovarian follicle. Notably, novel_circ_001327 showed a downregulation in Group B vs Group C while displaying an upregulation in the other groups. Conversely, novel_circ_009685 exhibited upregulation in the B group compared to the C group but demonstrated downregulation in the other groups. novel_circ_008283 displayed a consistently decreasing trend across all stages. novel_circ_018383 and novel_circ_020541 exhibited an initial increase followed by a rapid decline, reaching their peak expression in the 7-year-old group (Group C) and then rapidly decreasing in mares aged 15 years and above (Group D). This suggests that novel_circ_001327, novel_circ_009685, novel_circ_018383, and novel_circ_020541 may be involved in signaling pathways related to proliferation, differentiation, and reproductive aging. Currently, there is limited research about the involvement of circRNAs, including novel_circ_001327, in ovarian and follicular development, and further studies are needed for validation.

Organisms achieve biological processes through the coordinated action of different genes and their products. The results of the GO enrichment analysis in this study show that differentially expressed genes in ovarian GCs of Kazakh mares at different ages are significantly enriched in biological processes and cellular components. The Mitogen-Activated Protein Kinase (MAPK) cascade is a highly conserved module involved in various cellular functions, including cell proliferation, differentiation, and migration (46). Mammals express at least four distinctly regulated MAPKs: Extracellular Signal-Regulated Kinase (ERK)-1/2, c-Jun N-terminal Kinase (JNK1/2/3), p38 proteins (p38 α/β/γ/δ), and ERK5. These MAPKs are activated by specific MAPKKs: MEK1/2 for ERK1/2, MKK3/6 for p38, MKK4/7 (JNKK1/2) for JNK, and MEK5 for ERK5 (47). However, each MAPKK can be activated by multiple MAPKKKs, increasing the complexity and diversity of MAPK signaling (48). It is speculated that each MAPKKK is responsive to different stimuli. For example, the activation of ERK1/2 by growth factors depends on the MAPKKK c-Raf, while other MAPKKKs may activate ERK1/2 under pro-inflammatory stimuli. The MAPK signaling pathway is one of the important pathways in eukaryotic signal transduction. MAPKs are activated by a series of extracellular stimuli, mediating signals from the cell membrane to the nucleus and regulating various physiological processes such as cell proliferation, differentiation, apoptosis, and death (49, 50). The Hippo signaling pathway is an evolutionarily conserved pathway that controls organ size from flies to humans. In humans and mice, this pathway consists of MST1 and MST2 kinases, their cofactors Salvador, and LATS1 and LATS2 (51, 52). When activated at high cell density, phosphorylated LATS1/2 activates the transcriptional co-activators YAP and TAZ, promoting their cytoplasmic localization, leading to apoptosis and restricting excessive organ growth (53). When the Hippo pathway is inactive at low cell density, YAP/TAZ translocates into the nucleus, binds to the transcription enhancer factor (TEAD/TEF) family of transcription factors, and promotes cell growth and proliferation (54). YAP/TAZ also interact with other transcription factors or signaling molecules, linking processes mediated by the Hippo pathway with other key signaling cascades, such as those mediated by TGF-β and Wnt growth factors (54). The Wnt signaling pathway is essential for the fundamental development of many different species and organs, regulating cell fate determination, progenitor cell proliferation, and controlling asymmetric cell division (55, 56). There are at least three different Wnt pathways: the canonical pathway, the Planar Cell Polarity (PCP) pathway, and the Wnt/Ca²⁺ pathway. In the canonical Wnt pathway, the primary role of Wnt ligands binding to their receptors is to stabilize cytoplasmic β-catenin by inhibiting the degradation complex, thereby stabilizing cytoplasmic β-catenin. Subsequently, β-catenin enters the nucleus freely, interacts with the T-cell factor (TCF) family of transcription factors and associated coactivators, and activates Wnt-regulated genes (57). This indicates that these differentially expressed mRNAs may participate in mare follicular development and the proliferation and apoptosis of GCs through the pathways mentioned above.

Mammalian TLRs are expressed on innate immune cells, such as macrophages and dendritic cells, and respond to membrane components of both Gram-positive and Gram-negative bacteria (58). TLRs recognize pathogens, inducing the production of pro-inflammatory cytokines and upregulating co-stimulatory molecules, thereby rapidly activating the innate immune response (59). The TLR signaling pathway is divided into two groups: the MyD88-dependent pathway, leading to the production of pro-inflammatory cytokines and the rapid activation of NF-Kappa B and MAPK, and the MyD88-independent pathway, inducing INF-β and INF-induced genes and slowly activating NF-Kappa B and MAPK. The MAPK signaling pathway is inhibited by prolactin and is crucial for normal follicular development (59, 60). Cell adhesion molecules are proteins on the cell surface that regulate physiological processes such as cell morphology, migration, proliferation, and differentiation (61). These findings suggest that lncRNAs play a crucial role in regulating protein-coding genes through both cis and trans actions, influencing follicular development in mares during estrus and follicular development. This, in turn, impacts the follicular development of mares at different ages.

Notably, differentially expressed miRNA target genes also exhibit enrichment in the Longevity regulating pathway - multiple species. GCs regulate the development and blockade of follicles through various factors, including gonadotropin receptors, steroid hormones, as well as various growth factors, and cytokines (62–64). Steroid hormones, including estrogen, progestogen, and androgen, directly regulate the maturation of follicles (65). Ras protein, as a GTPase and a molecular switch in signaling pathways, regulates cell proliferation, survival, growth, migration, differentiation, or cytoskeletal dynamics (66). It includes 627 genes, and the Ras signaling pathway is closely linked to 168 pathways, including the MAPK signaling pathway, Calcium signaling pathway, PI3K-Akt signaling pathway, and Regulation of actin cytoskeleton. The FOXO transcription factor family regulates gene expression in cellular physiological events, including apoptosis, cell cycle control, glucose metabolism, oxidative stress resistance, and lifespan (67). Besides PKB, JNK, and AMPK, FOXOs undergo various post-translational modifications, including phosphorylation, acetylation, methylation, and ubiquitination, involving 191,623 genes and 102 pathways, including the MAPK signaling pathway, Cell cycle, and PI3K-Akt signaling pathway. The PI3K-Akt signaling pathway, activated by various stimuli or toxic damage, regulates fundamental cellular functions such as transcription, translation, proliferation, growth, and survival (68). Once activated, Akt controls crucial cellular processes through the phosphorylation of substrates involved in apoptosis, protein synthesis, metabolism, and the cell cycle (69). Gonadotropin-releasing hormone (GnRH) is secreted by the hypothalamus, acts on receptors in the anterior pituitary, and regulates the production and release of gonadotropins, LH, and FSH (70). Downstream signaling of protein kinase C (PKC) leads to the deactivation of epidermal growth factor (EGF) receptors and the activation of MAPKs, including ERK, JNK, and p38 MAPK. Active MAPKs move into the cell nucleus, leading to the activation of transcription factors and the rapid induction of early genes (71). ECM-receptor interaction: The extracellular matrix (ECM) is a complex mixture of structural and functional macromolecules that play a crucial role in tissue and organ morphogenesis, as well as in maintaining the structure and function of cells and tissues (72). Specific interactions between cells and ECM are mediated by transmembrane molecules, primarily integrins, as well as possibly proteoglycans, CD36, or other cell surface-related components (73). These interactions directly or indirectly control cellular activities such as adhesion, migration, differentiation, proliferation, and apoptosis (74). Cellular senescence is an irreversible state of cell growth arrest triggered by various factors such as telomere shortening, oncogene activation, radiation, DNA damage, and oxidative stress (75). Its characteristics include enlarged and flattened morphology, increased activity of senescence-associated β-galactosidase (SA-β-gal), and the secretion of inflammatory cytokines, growth factors, and matrix metalloproteinases, forming part of the senescence-associated secretory phenotype (SASP) (76). Cellular senescence is functionally linked to various biological processes, including aging, tumor suppression, placental biology, and embryonic development. Ovarian Steroidogenesis is critical for normal uterine function, the establishment and maintenance of pregnancy, as well as mammary gland development. The ovarian steroids, 17-β estradiol (E2), and progesterone (P4) play pivotal roles in these processes (77).

This study established networks of protein interaction for differentially expressed genes in GCs of Kazakh horses across various ages, identifying core node genes within each group. Notably, TNF (Tumor Necrosis Factor) emerged as a core node in both Group B vs Group D and Group C vs Group D, exhibiting a significant upregulation in expression across both comparisons. TNF participates in 134 signaling pathways, including the MAPK signaling pathway, NF-kappa B signaling pathway, C-type lectin receptor signaling pathway, and TNF signaling pathway. Previous research indicates that TNF plays a dual role in follicular development across species like pigs and mice, promoting both cell survival and apoptosis. TNF stands out as a crucial cytokine capable of triggering diverse intracellular signaling pathways, encompassing apoptosis, cell survival, inflammation, and immunity. Activated TNF forms homotrimers and binds to its receptors (TNFR1, TNFR2), leading to the trimerization of TNFR1 or TNFR2. TNFR1, expressed in nearly all cells, serves as the primary receptor for TNF (also known as TNF-α). Conversely, TNFR2 is expressed in specific cell types, including CD4 and CD8 T lymphocytes, endothelial cells, microglia, oligodendrocytes, neuronal subtypes, cardiomyocytes, thymocytes, and human mesenchymal stem cells, acting as a receptor for both TNF and LTA (also known as TNF-β). Upon ligand binding, TNFR mediates the recruitment of adaptor proteins (such as TRADD or TRAF2), initiating signal transduction. TNFR1 signaling triggers the activation of numerous genes, predominantly regulated by two distinct pathways: the NF-kB pathway and the MAPK cascade, or cell apoptosis and necroptosis. TNFR2 signaling activates the NF-kappa B pathway, including the PI3K-dependent NF-kappa B pathway and the JNK pathway, promoting cell survival. TNF participates in GC apoptosis through the MAPK signaling pathway, and the apoptosis of GCs is directly linked to follicle development and the decline in ovarian function. The findings of this study suggest that TNF is similarly involved in regulating aging in large mammals beyond mice and humans (78), underscoring its high conservation. Additionally, TNF plays a critical role in multiple signaling pathways, including the mTOR signaling pathway, TGF-beta signaling pathway, and Cytokine-cytokine receptor interaction. Thus, the core nodes regulated by the protein interaction network present potential regulatory factors influencing the follicular development of Kazakh horses at different ages. However, a more in-depth investigation is required to elucidate the specific regulatory mechanisms. Another core node, KALRN, was identified in both Group B vs Group D and Group C vs Group D, exhibiting a significant downregulation in expression across both comparisons. This suggests that KALRN may regulate the proliferation and apoptosis of follicular GCs, thereby influencing follicular development through downregulation.

The ceRNA theory serves as a supplement to the traditional miRNA→RNA paradigm, shedding light on a reverse regulatory mechanism in RNA→miRNA interactions. In this context, ceRNAs competitively engage with the same miRNA recognition elements (MREs), thereby influencing the silencing of target genes triggered by miRNA (79). This theory unveils a novel mechanism of RNA interaction, garnering widespread attention in recent years across various domains, including tumor research, pathological studies, and investigations into growth and development. Currently, research on ceRNA mechanisms is predominantly focused on animals such as cows, pigs, sheep, chickens, and geese, with limited studies on horses. Particularly scarce are investigations into the regulatory network of ceRNA involved in regulating ovarian and follicular development. Through the analysis of mRNA differences in GCs of Kazakh horses at various ages, certain critical candidate genes related to follicular development have become the focal points of our study. Cell growth and differentiation, embryonic development, and the onset and progression of tumors are vital aspects of cellular physiology. The findings of this study reveal that within the ceRNA network, eca-miR-486-3p and miR-486-y have the highest number of nodes. Their target genes participate in multiple signaling pathways associated with reproduction, such as cGMP-PKG signaling, Basal transcription factors, Hippo signaling pathway, Longevity regulating pathway, Ovarian steroidogenesis, Endocrine and other factor-regulated calcium reabsorption, Cortisol synthesis and secretion, and more. Notably, ovarian steroids, 17-β estradiol (E2), and progesterone (P4) play a crucial role in normal uterine function, the establishment and maintenance of pregnancy, and mammary gland development (80). Furthermore, the locally essential actions for normal ovarian physiology rely on the endocrine, paracrine, and autocrine effects of E2, P4, and androgens (81). In most mammals, including humans and mice, ovarian steroidogenesis occurs based on the two-cell/two-gonadotropin theory. This theory elucidates how GCs and follicular cells collaborate to produce ovarian steroids (82, 83). Follicular membrane cells respond to LH signals by increasing the expression of enzymes needed to convert cholesterol into androgens (such as androstenedione and testosterone). GCs respond to FSH signals by increasing the expression of enzymes needed to convert androgens derived from the follicular membrane into estrogens (E2 and estradiol) (84). Cortisol is the primary endogenous glucocorticoid that influences various physiological functions, including lipid and glucose metabolism, metabolic homeostasis, and stress adaptation. Cortisol production is primarily regulated by adrenocorticotropic hormone (ACTH) from the pituitary gland (85). The stimulatory effect of ACTH on cortisol synthesis depends on cAMP-dependent signal transduction but also involves membrane depolarization and cytoplasmic Ca²⁺ increase. Both cAMP and Ca²⁺ induce the expression of StAR, stimulating cholesterol transfer within mitochondria the steroidogenic enzymes (such as CHE, CYP17A1, and CYP11B1) in the cholesterol-to-cortisol pathway (86, 87). While we have constructed a regulatory network related to the mechanisms of action in GCs at different ages and validated the accuracy of the data through RT-qPCR, the intricate mechanisms within the regulatory network still require further in-depth research and discussion.

5 Conclusion

This study employed whole transcriptome sequencing to analyze ovarian GCs from a total of 20 Kazakh horses across four different age groups. The results revealed that differentially expressed mRNAs and miRNA target genes were significantly enriched in pathways such as MAPK signaling, Hippo signaling, Wnt signaling, Calcium signaling, Aldosterone synthesis and secretion, Cellular senescence, and NF-kappa B signaling. By constructing protein interaction networks for mRNAs and miRNA target genes, potential regulatory genes associated with reproductive aging were identified. Noteworthy genes for further investigation in the context of horse reproductive aging and regulation include MCMDC2, TNF, and the core regulatory factors eca-miR-486-3p and miR-486-y. Functional validation at the cellular level will be conducted in the future to elucidate the specific mechanisms underlying differential gene expression in our experimental results.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

The animal studies were approved by Laboratory Animal Care and Use Ethics Committee at Xinjiang Agricultural University. The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent was obtained from the owners for the participation of their animals in this study.

Author contributions

WR: Conceptualization, Methodology, Writing – review & editing, Data curation, Investigation, Validation, Writing – original draft. JW: Investigation, Validation, Writing – original draft, Software. YZ: Visualization, Writing – review & editing. TW: Writing – review & editing, Formal analysis, Validation. JM: Project administration, Resources, Supervision, Writing – review & editing. XY: Writing – review & editing, Conceptualization, Funding acquisition, Methodology, Project administration, Supervision.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This research was supported by the Xinjiang Uygur Autonomous Region Major Science and Technology Project (2022A02007-1), Autonomous Region Innovation Environment (Talent, Base) Construction project (PT2311) and the Graduate Student Research Innovation Project of Xinjiang Uygur Autonomous Region (XJ2022G127).

Acknowledgments

We extend our gratitude to Guangzhou Genedenovo Biotechnology Co., Ltd. for their support in sequencing (http://www.genedenovo.com/) and to KetengEdit (http://www.ketengedit.com) for their assistance with language editing for the manuscript.

Conflict of interest

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

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

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

References

1. Fan H, He J, Bai Y, He Q, Zhang T, Zhang J, et al. Baicalin improves the functions of granulosa cells and the ovary in aged mice through the mTOR signaling pathway. J Ovarian Res (2022) 15((1):34. doi: 10.1186/s13048-022-00965-7

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Ginther OJ, Beg MA, Donadeu FX, Bergfelt DR. Mechanism of follicle deviation in monovular farm species. Anim Reprod Sci (2003) 78(3-4):239–57. doi: 10.1016/S0378-4320(03)00093-9

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Matsuda F, Inoue N, Manabe N, Ohkura S. Follicular growth and atresia in mammalian ovaries: regulation by survival and death of granulosa cells. J Reprod Dev (2012) 58(1):44–50. doi: 10.1262/jrd.2011-012

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Li Q, Du X, Wang L, Shi K, Li Q. TGF-β1 controls porcine granulosa cell states: A miRNA-mRNA network view. THERIOGENOLOGY (2021) 160:50–60. doi: 10.1016/j.theriogenology.2020.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Zhang X, Yu T, Guo X, Zhang R, Lin P. Ufmylation regulates granulosa cell apoptosis via ER stress but not oxidative stress during goat follicular atresia. THERIOGENOLOGY (2021) 169(1):47–55. doi: 10.1016/j.theriogenology.2021.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Terenina E, Bonnet A, Monniaux D, Tosser-Klopp G. Differentially expressed genes and gene networks involved in pig ovarian follicular atresia. Physiol Genomics (2016) 49(2):67–80. doi: 10.1152/physiolgenomics.00069.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Manabe N, Goto Y, Matsuda-Minehata F, Inoue N, Maeda A, Sakamaki K, et al. Regulation mechanism of selective atresia in porcine follicles: regulation of granulosa cell apoptosis during atresia. J Reprod Dev (2004) 50(5):493–514. doi: 10.1262/jrd.50.493

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Hughes FM, Gorospe WC. Biochemical identification of apoptosis (programmed cell death) in granulosa cells: evidence for a potential mechanism underlying follicular atresia. ENDOCRINOLOGY (1991) 129(5):2415–22. doi: 10.1210/endo-129-5-2415

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Ma M, Wang H, Wang M, Yao W, Zhang J, Pan Z. Identification of circINHBB during the process of porcine follicular atresia and its impact on granulosa cell apoptosis. Chinese Agricultural Science (2021) 54(18):10. doi: 10.3864/j.issn.0578-1752.2021.18.017

CrossRef Full Text | Google Scholar

10. Carnevale EM, Ginther OJ. Defective oocytes as a cause of subfertility in old mares1. Biol Reprod (1995) (monograph_series1):209–14. doi: 10.1093/biolreprod/52.monograph_series1.209

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Fitzgerald C, Zimon AE, Jones EE. Aging and reproductive potential in women. Yale J Biol Med (1998) 71(5):367–81.

PubMed Abstract | Google Scholar

12. Grady ST, Watts AE, Thompson JA, Penedo MCT, Konganti K, Hinrichs K. Effect of intra-ovarian injection of mesenchymal stem cells in aged mares. J Assist Reprod Gen (2019) 36(3):543–56. doi: 10.1007/s10815-018-1371-6

CrossRef Full Text | Google Scholar

13. Dong L, Teh DBL, Kennedy BK, Huang Z. Unraveling female reproductive senescence to enhance healthy longevity. Cell Res (2023) 33(1):19. doi: 10.1038/s41422-022-00718-7

CrossRef Full Text | Google Scholar

14. Aurich C. Reproductive cycles of horses. Anim Reprod Sci (2011) 124(3-4):220–8. doi: 10.1016/j.anireprosci.2011.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Marinone AI, Losinno L, Fumuso E, Rodríguez EM, Redolatti C, Cantatore S, et al. The effect of mare's age on multiple ovulation rate, embryo recovery, post-transfer pregnancy rate, and interovulatory interval in a commercial embryo transfer program in Argentina. Anim Reprod Sci (2015) 158:53–9. doi: 10.1016/j.anireprosci.2015.04.007

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Rizzo M, Ducheyne KD. Advanced mare age impairs the ability of invitro-matured oocytes to correctly align chromosomes on the metaphase plate. EQUINE Vet J (2019) 51(2):252–7. doi: 10.1111/evj.12995

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Hatzirodos N, Irving-Rodgers HF, Hummitzsch K, Harland ML, Rodgers RJ. Transcriptome profiling of granulosa cells of bovine ovarian follicles during growth from small to large antral sizes. BMC Genomics (2014) 15(1):24. doi: 10.1186/1471-2164-15-24

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Langmead B, Salzberg SL, Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods 9: 357-359. Nat Methods (2012) 9(4):357–9. doi: 10.1038/nmeth.1923

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods (2015) 12((4):357. doi: 10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol (2015) 33(3):290–5. doi: 10.1038/nbt.3122

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc (2016) 11(9):1650–67. doi: 10.1038/nprot.2016.095

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Dewey CN, Bo L. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinf (2011) 12(1):323. doi: 10.1186/1471-2105-12-323

CrossRef Full Text | Google Scholar

23. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol (2014) 15(12):550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. BIOINFORMATICS (2010) 26(1):139. doi: 10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Ashburner M, Ball CA, Blake JA, Botstein D. Gene Ontology: tool for the unification of biology. Nat Genet (2000) 25((1):25–9. doi: 10.1038/75556

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Kanehisa G. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res (2000) 27((1):29–34. doi: 10.1093/nar/28.1.27

CrossRef Full Text | Google Scholar

27. Damian S, Andrea F, Stefan W, Kristoffer F, Davide H, Jaime HC, et al. STRING v10: protein–protein interaction networks, integrated over the tree of life. Nucleic Acids Res (2015) 43(D1):D447–52. doi: 10.1093/nar/gku1003

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Shannon P. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res (2003) 13(11):2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Der-Auwera GAV, Carneiro MO, Hartl C, Poplin R, Thibault J. From FastQ data to high confidence variant calls: the Genome Analysis Toolkit best practices pipeline. Curr Protoc Bioinf (2013) 43(1110):10–1. doi: 10.1002/0471250953.bi1110s43

CrossRef Full Text | Google Scholar

30. Zhang Y, Yan Z, Qin Q, Nisenblat V, Yan L. Transcriptome landscape of human oocytes and granulosa cells throughout folliculogenesis. bioRxiv (2018). doi: 10.1101/285445

CrossRef Full Text | Google Scholar

31. Lan DL, Xiong XR, Wei YL, Xu T, Zhong JC, Zhi XD, et al. RNA-Seq analysis of yak ovary: improving yak gene structure information and mining reproduction-related genes. Sci China Life Sci (2014) 57((9):925–35. doi: 10.1007/s11427-014-4678-2

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Bahrami A, Miraie-Ashtiani SR, Sadeghi M, Najafi A. reproduction mirna-mrna network involved in folliculogenesis interactome: systems biology approach. REPRODUCTION (2017) 154((1):51–65. doi: 10.1530/REP-17-0049

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Knapczyk-Stwora K, Nynca A, Ciereszko RE, Paukszto L, Jastrzebski JP, Czaja E, et al. Flutamide-induced alterations in transcriptional profiling of neonatal porcine ovaries. J Anim Sci BIOTECHNO (2019) 3:10–35. doi: 10.1186/s40104-019-0340-y

CrossRef Full Text | Google Scholar

34. Liang-Ting T, Xue-Qin R, Ning M, Fu-Ping Z, Xi N, Yi-Qi R, et al. Analysis of alternative splicing events by RNA sequencing in the ovaries of Xiang pig at estrous and diestrous. THERIOGENOLOGY (2018) 119:60. doi: 10.1016/j.theriogenology.2018.06.022

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Bao Y, Yao X, Li X, Ei-Samahy M, Yang H, Liang Y, et al. INHBA transfection regulates proliferation, apoptosis and hormone synthesis in sheep granulosa cells. THERIOGENOLOGY (2021) 175:111–22. doi: 10.1016/j.theriogenology.2021.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Liu X, Lin X, Zhang S, Guo C, Li J. Lycopene ameliorates oxidative stress in the aging chicken ovary via activation of Nrf2/HO-1 pathway. Aging (2018) 10((8):2016–36. doi: 10.18632/aging.101526

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Zhao X, Wu Y, Li H, Li J, Yao Y, Cao Y, et al. Comprehensive analysis of differentially expressed profiles of mRNA, lncRNA, and miRNA of Yili geese ovary at different egg-laying stages. BMC Genomics (2022) 23(1):1–20. doi: 10.1186/s12864-022-08774-4

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Huang T, Gu W, Liu E, Shi X, Wang B, Wu W, et al. Comprehensive analysis of miRNA-mRNA/lncRNA during gonadal development of triploid female rainbow trout (Oncorhynchus mykiss). GENOMICS (2021) 113(6):3533–43. doi: 10.1016/j.ygeno.2021.08.018

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Xing B, Bai X, Guo H, Chen J, Hua L, Zhang J, et al. Long non-coding RNA analysis of muscular responses to testosterone deficiency in Huainan male pigs. Anim Sci J (2017) 88. doi: 10.1111/asj.12777

CrossRef Full Text | Google Scholar

40. Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. CELL (2004) 116. doi: 10.1016/S0092-8674(04)00045-5

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Ambros V. The function of animal MicroRNAs. NATURE (2004) 431(7006):350–5. doi: 10.1038/nature02871

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Zamore PD, Du T. Beginning to understand microRNA function. Cell Research (2007) 17(8):661–3. doi: 10.1038/cr.2007.67

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Tomé M, López-Romero P, Albo C, Sepúlveda JC, Fernández-Gutiérrez B, Dopazo A, et al. miR-335 orchestrates cell proliferation, migration and differentiation in human mesenchymal stem cells. Cell Death Differ (2011) 18:985–95. doi: 10.1038/cdd.2010.167

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Najafi A, Tavallaei M, Hosseini SM. A systems biology approach for miRNA-mRNA expression patterns analysis in non-small cell lung cancer. Cancer biomark (2016) 16((1):31. doi: 10.3233/CBM-150538

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Vasileva A, Tiedau D, Firooznia A, Muller-Reichert T, Jessberger R. Tdrd6 is required for spermiogenesis, chromatoid body architecture, and regulation of miRNA expression. Curr Biol (2009) 19((8):630–9. doi: 10.1016/j.cub.2009.02.047

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Tanoue T, Nishida E. Docking interactions in the mitogen-activated protein kinase cascades. Pharmacol THERAPEUT (2002) 93(2-3):193–202. doi: 10.1016/S0163-7258(02)00188-2

CrossRef Full Text | Google Scholar

47. Dent P, Yacoub A, Fisher P, Hagan MP, Grant S. MAPK pathways in radiation response. Oncogene (2003) 22:5885–96. doi: 10.1038/sj.onc.1206701

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Harper SJ, Lograsso P. Signalling for survival and death in neurones: the role of stress-activated kinases, JNK and p38. Cell Signal (2001) 13(5):299–310. doi: 10.1016/S0898-6568(01)00148-6

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Kang HM, Jeong CB, Lee YH, Cui YH, Kim DH, Lee MC, et al. Cross-reactivities of mammalian MAPKs antibodies in rotifer and copepod: Application in mechanistic studies in aquatic ecotoxicology. Mar pollut Bull (2016) 124((2):614–23. doi: 10.1016/j.marpolbul.2016.11.049

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Liu Y, Li W, Su J. Progress on the regulatory mechanisms of signals associated with mammalian oocyte maturation. Progress in Veterinary Medicine (2019) 40(9):5. doi: CNKI:SUN:DYJZ.0.2019-09-021

Google Scholar

51. Zhao B, Tumaneng K, Guan KL. The Hippo pathway in organ size control, tissue regeneration and stem cell self-renewal. Nat Cell Biol (2011) 13. doi: 10.1038/ncb2303

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Zhao B, Li L, Guan KL. Hippo signaling at a glance. Company Biologists Ltd (2010) 123((23):4001–6. doi: 10.1242/jcs.069070

CrossRef Full Text | Google Scholar

53. Zhao B, Li L, Lei Q, Guan KL. The Hippo-YAP pathway in organ size control and tumorigenesis: an updated version. Genes Development: J Devoted to Mol Anal Gene Expression Eukaryotes Prokaryotes Viruses (2010) 9):24. doi: 10.1101/gad.1909210

CrossRef Full Text | Google Scholar

54. Tremblay AM, Camargo FD. Hippo signaling in mammalian stem cells. Semin Cell Dev Biol (2012) 23(7):818–26. doi: 10.1016/j.semcdb.2012.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Hatsell S, Rowlands T, Hiremath M, Cowin P. Beta-catenin and Tcfs in mammary development and cancer. J MAMMARY GLAND Biol (2003) 2):8. doi: 10.1023/A:1025944723047

CrossRef Full Text | Google Scholar

56. Yang Y. Wnts and wing: Wnt signaling in vertebrate limb development and musculoskeletal morphogenesis. Birth Defects Res Part C Embryo Today Rev (2003) 69(4):305–17. doi: 10.1002/bdrc.10026

CrossRef Full Text | Google Scholar

57. Kalderon D. Similarities between the Hedgehog and Wnt signaling pathways. Trends Cell Biol (2002) 12(11):523–31. doi: 10.1016/S0962-8924(02)02388-7

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Kawai T, Akira S. TLR signaling. Cell Death Differentiation (2006) 19(1):24–32. doi: 10.1038/sj.cdd.4401850

CrossRef Full Text | Google Scholar

59. Neill LAJ O, Bowie AG. The family of five: TIR-domain-containing adaptors in Toll-like receptor signalling. Nat Rev Immunol (2007) 7(5):353–64. doi: 10.1038/nri2079

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Devi YS, Seibold AM, Shehu A, Maizels E, Halperin J, Le J, et al. Inhibition of MAPK by prolactin signaling through the short form of its receptor in the ovary and decidua: involvement of a novel phosphatase. J Biol Chem (2011) 286(9):7609–18. doi: 10.1074/jbc.M110.166603

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Barclay AN. Membrane proteins with immunoglobulin-like domains–a master superfamily of interaction molecules. Semin Immunol (2003) 15(4):215–23. doi: 10.1016/S1044-5323(03)00047-2

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Simpson ER, Clyne C, Rubin G, Boon WC, Robertson K, Britt K, et al. AROMATASE—A BRIEF OVERVIEW. Annu Rev Physiol (2002) 64(1):93–127. doi: 10.1146/annurev.physiol.64.081601.142703

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Tiwari M, Prasad S, Tripathi A, Pandey AN, Ali I, Singh AK. Apoptosis in mammalian oocytes: a review. Apoptosis (2015) 20(8):1019–25. doi: 10.1007/s10495-015-1136-y

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Huang XY, Xiong XR, Han J, Yang XY, Wang Y, Wang B, et al. Expression of KDM1A during the Follicular Development Process in Yak Oocytes. Chinese Agricultural Science (2019) 52(24):8. doi: 10.3864/j.issn.0578-1752.2019.24.016

CrossRef Full Text | Google Scholar

65. Zhang J, Zhang JY, Yue YQ, Zhao D, Fan YL, Ma Y, et al. Regulatory Role of the LKB1 Gene in Steroid Hormone-Related Gene Expression in Ovarian Granulosa Cells. Chinese Agricultural Science (2022) 55(10):10. doi: 10.3864/j.issn.0578-1752.2022.10.015

CrossRef Full Text | Google Scholar

66. Karnoub AE, Weinberg RA. Ras oncogenes: split personalities. Nat Reviews: Mol Cell Biol (2008) 7):9. doi: 10.1038/nrm2438

CrossRef Full Text | Google Scholar

67. Van Der Heide LP, Hoekman MFM, Smidt MP. The ins and outs of FoxO shuttling: mechanisms of FoxO translocation and transcriptional regulation. Biochem J (2004) 380(2):297–309. doi: 10.1042/bj20040167

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Glaviano A, Foo ASC. PI3K/AKT/mTOR signaling transduction pathway and targeted therapies in cancer. Mol Cancer (2023) 22(1):138. doi: 10.1186/s12943-023-01827-6

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Engelman JA, Luo J, Cantley LC. The evolution of phosphatidylinositol 3-kinases as regulators of growth and metabolism. Nat Rev Genet (2006) 7:606–19. doi: 10.1038/nrg1879

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Ruf F, Sealfon SC. Genomics view of gonadotrope signaling circuits. Trends Endocrin Met (2004) 15(7):331–8. doi: 10.1016/j.tem.2004.07.008

CrossRef Full Text | Google Scholar

71. Harrison GS, Wierman ME, Nett TM, Glode LM. Gonadotropin-releasing hormone and its receptor in normal and Malignant cells. Endocr-Relat Cancer (2004) 4):11. doi: 10.1677/erc.1.00777

CrossRef Full Text | Google Scholar

72. Flier AVD, Sonnenberg A. Function and interactions of integrins. Cell Tissue Res (2001) 305(3):285–98. doi: 10.1007/s004410100417

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Reddy KVR, Mangale SS. Integrin receptors: the dynamic modulators of endometrial function. Tissue Cell (2003) 35(4):260–73. doi: 10.1016/S0040-8166(03)00039-9

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Ven ATVD, Vivante A, Hildebrandt F. Novel insights into the pathogenesis of monogenic congenital anomalies of the kidney and urinary tract. J Am Soc Nephrol (2018) 29((1):36–50. doi: 10.1681/ASN.2017050561

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Muñoz-Espín D, Serrano M. Cellular senescence: from physiology to pathology. Nat Rev. Mol Cell Biol (2014) 15(7):482–96. doi: 10.3410/f.718463914.793498452

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Childs BG, Durik M, Baker DJ, Van Deursen JM. Cellular senescence in aging and age-related disease: from mechanisms to therapy. Nat Med (2015) 21(12):1424–35. doi: 10.1038/nm.4000

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Wiley C, Campisi J. From ancient pathways to aging cells—Connecting metabolism and cellular senescence. Cell Metab (2016) 23(6):1013–21. doi: 10.1016/j.cmet.2016.05.010

PubMed Abstract | CrossRef Full Text | Google Scholar

78. Kong Q-Q, Wang J, Xiao B, Lin F-H, Zhu J. Cumulus cell-released tumor necrosis factor (TNF)-α promotes post-ovulatory aging of mouse oocytes. Aging (2018) 10((7):1745–57. doi: 10.18632/aging.101507

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Thomson DW, Dinger ME. Endogenous microRNA sponges: evidence and controversy. Nat Rev Genet (2016) 17((5):272–83. doi: 10.1038/nrg.2016.20

PubMed Abstract | CrossRef Full Text | Google Scholar

80. Tamer E C, Ozcivit IB. Anti-Müllerian hormone and ovarian aging. Gynecol Endocrinol (2021) 37(10):867–8. doi: 10.1080/09513590.2021.1977276

PubMed Abstract | CrossRef Full Text | Google Scholar

81. Gore AC, Hall JE, Hayes FJ. (2015).

Google Scholar

82. Park SU, Walsh L, Berkowitz KM. Mechanisms of ovarian aging. Reproduction: Off J Soc Study Fertility (2021) 2):162. doi: 10.1530/REP-21-0022

CrossRef Full Text | Google Scholar

83. Uliani RC, Conley AJ, Corbin CJ, Friso AM, Maciel LFS, Alvarenga MA. Anti-Müllerian hormone and ovarian aging in mares. J Endocrinol (2019) 2:147–56. doi: 10.1530/JOE-18-0391

CrossRef Full Text | Google Scholar

84. De Kat AC. Female Aging : Quantification of ovarian reserve and its association with cardiovascular health. Dissertation, The Netherlands, Utrecht University. (2017).

Google Scholar

85. Shimasaki ES. The physiology of folliculogenesis: the role of novel growth factors. Fertil Steril (2001) 76. doi: 10.1016/s0015-0282(01)02859-x

CrossRef Full Text | Google Scholar

86. Wood JR, Strauss JF. Multiple signal transduction pathways regulate ovarian steroidogenesis. Rev Endocr Metab Disord (2002) 3(1):33–46. doi: 10.1023/A:1012748718150

PubMed Abstract | CrossRef Full Text | Google Scholar

87. Halder G, Johnson RL. Hippo signaling: Growth control and beyond. DEVELOPMENT (2011) 138(1):9–22. doi: 10.1242/dev.045500

PubMed Abstract | CrossRef Full Text | Google Scholar