Whole-Transcriptome Analysis of LncRNAs Mediated ceRNA Regulation in Granulosa Cells Isolated From Healthy and Atresia Follicles of Chinese Buffalo

Granulosa cells (GCs) are the main supporting cells in follicles and play an important role in the regulation of oocyte maturation and follicular atresia. Accumulating evidence indicates that non-coding RNAs participate in regulation of the physiological function of GCs. However, whole-transcriptome analysis for GCs of buffalo has yet to be reported. In this study, healthy follicles (HFs) and atretic follicles (AFs) were defined according to the apoptosis rate of GCs and the hormone level in follicular fluid. GCs were collected from HFs and AFs (n = 15, 5 < n < 8 mm) for whole-transcriptome analysis using second-generation high-throughput sequencing. A total of 1,861 and 1,075 mRNAs, 159 and 24 miRNAs, and 123 and 100 lncRNAs, were upregulated and downregulated between HFs and AFs, respectively. Enrichment of functions and signaling pathways of these differentially expressed (DE) genes showed that most of DEmRNAs and targets of DEmiRNAs were annotated to the categories of ECM–receptor interaction and focal adhesion, as well as PI3K-AKT, mTOR, TGF-beta, Rap1, and estrogen signaling pathways. The competing endogenous RNA (CeRNA) network was also constructed based on the ceRNA theory which further revealed regulatory roles of these DERNAs in GCs of buffalo follicles. Finally, we validated that lnc4040 regulated the expression of Hif1a as miR-709 sponge in a ceRNA mechanism, suggesting their critical functions in GCs of buffalo follicles. These results show that lncRNAs are dynamically expressed in GCs of HFs and AFs, and interacting with target genes in a ceRNA manner, suggesting their critical functions in buffalo follicular development and atresia.

Introduction

In the process of mammalian follicle development, more than 99% of follicles become atretic and very few of follicles develop to ovulation (1). Chinese buffalo is an important large domestic animal distributed in the tropical and subtropical regions of China. It is generally believed that the sexual maturity of buffalo is later, and its postpartum estrus period is longer than that of cattle (2, 3). As a mono-ovulating animal, follicular atresia is also an important factor restricting its reproductive performance, and the follicular atresia rate of buffalo is higher than that of cattle (4). During each estrus period, only one follicle matures and ovulates and other follicles undergo atresia during follicular development, and this occurs repeatedly throughout the developmental stages. However, the exact mechanism of cell death during this process is not fully understood. Therefore, elucidating the molecular mechanism of follicle atresia is expected to reverse the fate of follicles in buffalo, thus effectively utilizing more follicle resources and improving the reproductive efficiency of Chinese buffalo.

Previous studies have shown that follicular atresia is closely related to the apoptosis of ovarian granulosa cells (GCs) (5–8). GCs gradually shrink in size, and apoptotic bodies are formed. The chromatin is then condensed, a large amount of which is shed into the follicular fluid of atretic follicles (AFs) and the blood supply around the follicles is greatly reduced (9–12). In addition, interferon and some growth factors can also induce apoptosis. The intrinsic factors related to apoptosis are usually those associated with stress, such as nutritional deficiency, oxidative damage, and genetic damage as well as molecular stress (13). Recent studies have shown that the autophagy of GCs partially dominates the occurrence of follicular atresia (14–16).

However, the molecular mechanism of follicular atresia needs to be further explored. The development of transcriptome techniques provides an effective route to study the mechanism of follicular development and atresia. Early transcriptomics studies of follicles focused on mRNA, and researchers conducted studies at different stages of follicular development in cattle (17–20). Subsequent research showed that miRNA expression was spatiotemporal specific during bovine follicular development. In contrast to small follicles, bta-mir-144, bta-mir-202, bta-mir-451, bta-mir-652, and bta-mir-873 were all upregulated in large healthy follicles (HFs) and they targeted the signaling pathways involved in follicular cell proliferation and steroid generation (21). Long non-coding RNAs (lncRNAs) have received considerable attention for their important roles in epigenetic regulation, chromatin modification, genomic imprinting, transcriptional control, and pre-translational as well as posttranslational mRNA processing (22, 23).

The main mechanism of competing endogenous RNAs (ceRNAs) is to competitively bind miRNAs, remove the inhibitory effect of miRNAs on target genes, and improve the expression levels of target genes (24). The development of this hypothesis has given mRNAs and non-coding RNAs new and broader biological functions (25). The research between compact cumulus cells and expanded GCs found 89 lncRNAs, of which 12 are encoded within introns of genes known to be involved in GC processes. This suggested that unique non-coding RNA transcripts may contribute to the regulation of cumulus expansion and oocyte maturation (26).

In the study of mammalian follicles, researchers have obtained a large amount of transcriptomics data, which has a crucial regulatory role at all stages of follicular development (27). Several miRNAs and lncRNAs were screened for specific signal pathway regulation analysis. However, for large domestic animals, there are few such studies. Therefore, we used high-throughput sequencing to perform a whole-transcriptome analysis of GCs isolated from AFs and HFs of the Chinese buffalo, in order to obtain various types of DE-RNA molecules, and this allowed the construction of a ceRNA network to explore the molecular regulatory pathways involved. The non-coding RNAs obtained in this study can provide a strong basis for the subsequent study of the altered physiological functions of GCs during follicular development and atresia in buffalo. The aim of the study was to provide a theoretical basis and clues for follicular atresia observed in the Chinese buffalo.

Results

Appearance and Classification of HFs and AFs in Buffalo

By observing exposed follicles, the surfaces of HFs were seen to be rich in capillaries and the follicular fluid was luminous yellow, while little or no capillaries were observed on the surface of AFs, and the follicular fluid was grayish white (Figure 1A). Follicular fluid and GCs were isolated from both HFs and AFs. The concentrations of estradiol (E₂) and progesterone (PROG) in isolated follicular fluid were measured by ELISA, and these were 126.21 ± 9.28 and 21.88 ± 1.78 ng·ml⁻¹, respectively, and the ratio was 5.77. The concentrations of E₂ and PROG in AFs were 76.63 ± 1.09 and 53.43 ± 2.68 ng·ml⁻¹, respectively, and the ratio was 1.43 (Table 1). The apoptosis rate of GCs in HFs was 4.93 ± 1.59, while that in AFs was 21.31 ± 1.40 (Figure 1B and Table 2). Transmission electron microscope observations showed that the morphology and submicrostructure of each organelle in HF were more normal than those in AF, there were more mitochondria in HF compared to the AF, chromatin breaks down to form granules, and a large number of vacuoles and apoptotic bodies appear in the cytoplasm of GCs in AFs (Figure 1C).

FIGURE 1

Figure 1. Phenotypic detection of healthy and atretic follicles. (A) The appearance of follicles with scale bars = 1.5 cm: the capillaries were rich on the surface of healthy follicles, and the follicular fluid was yellow, while little or no capillaries were seen on the surface of atretic follicles, and the follicular fluid was grayish white. (B) The apoptosis rate of the granulosa cells (n = 3) of healthy and atretic follicles was assessed by flow cytometry. (C) Transmission electron microscope of granulosa cells in healthy and atretic follicles with scale bars = 5 and 1 μm, respectively. M indicates mitochondria, and V represents vacuoles.

TABLE 1

Table 1. Hormone levels in follicular fluid of HFs and AFs.

TABLE 2

Table 2. Apoptotic rate of GCs in HFs and AFs.

Overview of mRNA-Seq Data

From mRNA-seq data, 40,990 protein-coding genes were identified in the buffalo genome. The Pearson correlation coefficient calculated according to the gene fragments per kilobase of exon model per million reads mapped (FPKM) between the HFs and AFs samples showed that the consistency among the three repeated samples in the same group was more than 0.9, indicating that the repeatability between each two samples was dependable (Figure 2A). Their log2FC values are presented as volcano plot diagrams in Figure 2B. Among all of these genes, 2,936 (1,861 upregulated and 1,075 downregulated) DEmRNAs (fold change ≥ 2 and p-value < 0.05) were identified in GCs between HFs and AFs. A heat map of DEmRNAs showed the general expression profiles of the DEmRNAs in each group (Figure 2C).

FIGURE 2

Figure 2. Identification and analysis of DEmRNAs in granulosa cells between healthy and atretic follicles. (A) The Pearson correlation heat map between each sample. (B) Volcano plot diagrams showing log2FC and p-values of DEmRNAs. (C) Clustering heat map for DEmRNAs in each sample between the two groups. The expression of DEmRNAs was calculated with 2 as the base logarithm, and different samples and transcripts were analyzed by hierarchical cluster analysis; red and blue represent upregulation and downregulation, respectively. (D) Top 20 of the KEGG pathway bubble chart of DEmRNAs. The bubble size represents the number of mRNAs enriched in the pathway, and the bubble color represents the p-value. (E) GO enrichment histogram of DEmRNAs; blue and yellow represent upregulation and downregulation, respectively.

To explore the functions of the DEmRNAs, KEGG was the main public database used regarding the pathway. The results showed that there were 42 and 72 DEmRNAs significantly enriched to ECM–receptor interactions and focal adhesion, respectively. In addition, there were 102, 48, 29, 35, and 33 enriched to the PI3K-AKT, thyroid hormone, TGF-beta, mTOR, and estrogen signaling pathways, respectively (Figure 2D). GO annotation and GO enrichment analysis were also performed. Based on GO annotation, it was found that DEmRNAs were annotated to the cellular process, biological regulation, single-organism processes, metabolic processes, and response to stimulus under biological processes (BP), to the cells, cell parts, organelles, and organelle parts under cellular components (CC) and to the binding, catalytic activities, and molecular transducer activities under molecular functions (MF) (Figure 2E). In organisms, different transcripts coordinated with each other to exercise their biology and pathway-based analysis is helpful to further understand the biological functions of the transcripts.

Overview of miRNA-Seq Data

A total of 704,220 tags were generated. The results were filtered based on length (18–35 nt), and most selected tags were 22 nt in length in both the HF and AF groups (Figure 3A). All the clean tags were aligned with the reference genome. The tags mapped to repeat sequences were removed. Ultimately, 2,446 mature miRNAs (1,173 known and 1,273 novel) were detected, and these were drawn as a scatter plot using the expression between the different groups (Figure 3B). We identified miRNAs with a fold change ≥ 2 and p-values < 0.05 in a comparison of significant DEmiRNAs (159 upregulated and 24 downregulated) (Figure 3C). The heat maps of all DEmiRNAs were drawn to display the miRNA expression levels in different samples and to cluster miRNAs with a similar expression pattern (Figure 3D). Targeted mRNAs of these DEmiRNAs are listed in Supplementary Table 1. Pathway enrichment analysis of the targets of DEmiRNAs showed that there were 930, 884, 702, and 407 mRNAs linked to PI3K-Akt, MAPK, Ras, and mTOR signaling pathways, respectively. In addition, 799 and 400 mRNAs were linked to endocytosis and autophagy, respectively (Figure 3E).

FIGURE 3

Figure 3. Identification and analysis of DEmiRNAs in granulosa cells between healthy and atretic follicles. (A) Length statistic histogram of tags between the two groups. (B) The expression of identified miRNAs between the two groups was upregulated (30) and downregulated (blue). (C) Volcanic diagram of the expression of DEmiRNAs in different groups. (D) Clustering heat map for DEmiRNAs in each sample between the two groups with red for upregulation and blue for downregulation. (E) Top 20 of the KEGG pathway bubble chart of target mRNAs; the bubble size represents the number of mRNAs enriched in the pathway, and the bubble color represents the p-value.

Overview of lncRNA-Seq Data

The expressions of 3,488 known and 1,338 novel transcripts were identified. Two different softwares, namely, CNCI (version 2) (28) and CPC (29) (http://cpc.cbi.pku.edu.cn/), were used to assess the protein-coding potential of novel transcripts by default parameters. Novel transcripts were mapped with the SwissProt database to assess protein annotation. The intersection results of both non-protein-coding potential and non-protein annotation were chosen as the new lncRNAs (1,344) (Figure 4A). The lncRNAs obtained were classified into five classes according to their location relative to protein-coding genes, and these were named intergenic, bidirectional, intronic, antisense, and sense overlapping lncRNAs. The different types of lncRNAs probably have different biological functions (Figure 4B). A comparison of the genomic characterizations of the lncRNAs with mRNAs showed that their transcripts were similar in length distribution, except that lncRNA had relatively higher numbers of long transcripts (>4,500 bp) than mRNAs. For the number of exons, a higher percentage of lncRNAs had two to four exons. In addition, lncRNAs tended to have a shorter ORF length and a lower FPKM value (Figures 4C–E). Transcripts with a fold change ≥ 2 and a p-value < 0.05 were identified in a comparison as significant DElncRNAs (209 in final, 122 upregulated, and 87 downregulated) (Supplementary Table 4) were seen when a volcano plot was drawn (Figure 4F). Based on the expression of the DElncRNAs, the relationship between the sample and the DElncRNAs was clustered hierarchically, and a heat map is used to represent the clustering results (Figure 4G).

FIGURE 4

Figure 4. Identification and analysis of DElncRNAs in granulosa cells between healthy and atretic follicles. (A) Venn diagram of the lncRNA annotation results of CPC, CNCI, and SwissProt. (B) LncRNA classification histogram. (C–E) Comparison of lncRNA with mRNA with respect to the transcript length, exon number, ORF length, and FPKM values. (F) Volcano plot of the expression of DElncRNAs in different groups. (G) Heat map of all DElncRNAs.

CeRNA Regulatory Network in GCS Between HFs vs. AFs

To reveal the global regulatory network of protein-coding RNAs and ncRNAs in GCs between HFs vs. AFs, a ceRNA network was constructed using DEmiRNAs, DEmRNAs, and DElncRNAs based on the ceRNA theory. In total, 1,458 DEmRNAs and 134 DElncRNAs were predicted as targets for 135 miRNAs. The lncRNA–miRNA–mRNA network was constructed by assembling all the co-expression competing triplets, which has been identified above. The connectivity of RNA in the ceRNA network was defined as the number of co-expressed targeted miRNAs. So ceRNAs with the highest connectivity were regarded as hub genes, which are more essential in biological networks.

In this ceRNA network, miR-424-3p, miR-542-5p, miR-735-5p, mi-503-3p, miR-212-3p, miR-216-3p, and miR-2440-5p are involved in morethan 100 nodes, suggesting that they may act as core regulators. In addition, lncRNAs including XR_003111160.1, XR_003106664.1, XR_003103086.1, XR_327919.2, XR_003111331.1, XR_003103012.1, XR_003111622.1, and TCONS_00165726 participated in the network.

Because of the large number of node genes and the tedious gene interaction network, the top 30 DEmiRNAs were selected as the source genes to simplify visualization using Cytoscape software (v3.6.0) (http://www.cytoscape.org/) (Figure 5A) with their target mRNAs. These 30 miRNAs included miR-146-3p, miR-450-3p, miR-424-3p, miR-212-3p, and miR-29-3p. GO enrichment analysis of all target mRNAs showed that most mRNAs were enriched in cell parts, binding, cellular processes, and organelles. KEGG pathways were enriched in ECM–receptor interactions, PI3K-AKT signaling pathway, protein digestion and absorption, focal adhesion, Hippo signaling pathway, mTOR signaling pathway, and various other metabolic pathways (Figures 5B,C).

FIGURE 5

Figure 5. CeRNA regulatory network and functional analysis. (A) CeRNA network built with the top 30 DEmiRNAs. Red and blue represent the upregulated and downregulated levels, respectively. (B) Level 2 GO enrichment of target mRNAs. (C) Top 20 KEGG pathways of target mRNAs.

Verification of the Regulation of Hif1a by ceRNA Circuitry in GCs

To further support the ceRNA hypothesis, a regulatory circuitry containing TCONS_00104040 (lnc4040), miR-709-5p (miR-709), and TCONS_00111150 (Hif1a) was selected for further verification. Hif1a can regulate the proliferation, autophagy, and cell cycle of GCs, and it also plays an important role in ovulation selection (31–33). Our initial experiments confirmed that the expression of lnc4040 was positively correlated with Hif1a and negatively correlated with miR-709; Western blot analysis revealed that the protein expression of Hif1a was significantly higher in the HF group than in the AF group (Figures 6A,B). We hypothesized that lnc4040 function as miRNA (miR-709) sponges to regulate the expression of Hif1a indirectly. To further confirm the binding events among miR-709, lnc4040, and Hif1a, we predicted the potential miR-709-binding site within lnc4040 and Hif1a by RNA22 version 2.0 (http://cm.jefferson.edu/rna22/Interactive/) (Figure 6C). Subsequent experiments validated the direct regulatory effect of miR-709 on Hif1a in 293T cells. Dual luciferase expression vectors carrying the wild-type (WT) miR-709-binding sites or one of the mutated versions of the miR-709-binding sites were constructed. MiR-709 mimics by co-transfection inhibited luciferase activity by the luciferase vectors carrying the mutations in site 1 (mut-1) or both sites 1 and 2 (mut-3), but not by the luciferase vectors carrying the WT Hif1a-binding sites or the mutation in binding site 2 (mut-2) (Figure 6D). This suggests that the binding site 1 of Hif1a might play the major role in miR-709 regulation of Hif1a. We also found that miR-709 can significantly inhibit the luciferase activity (Figure 6E). The HIF1A protein expression level in cultured GCs was decreased after miR-709-mimic transfection and conversely was increased after miR-709-inhibitor transfection (Figure 6F). Moreover, after overexpression of lnc4040, the downregulation of HIF1A protein expression induced by miR-709-mimic was alleviated (Figure 6G). These results evidenced that lnc4040 may function as a sponge of miR-709 to regulate the expression of Hif1a. In summary, our results showed that three-node lncRNA-mediated regulatory circuitry might play an essential role in GCs during follicle development.

FIGURE 6

Figure 6. Experimental validated of ceRNAs regulation with miR-709. (A) Relative expression levels of miR-709, Hif1a, and lnc4040 in GCs between HFs and AFs. (B) Expression of HIF1A protein in GCs between HFs vs. AFs and quantified by ImageJ (v1.45), error bars represent SEM (n = 3) **p < 0.01. (C) Predicted miR-709-binding site on Hif1a and lnc4040, and design of luciferase reporter. (D) 293T cells were co-transfected with wild-type (WT) or mutant (MUT) luciferase reporters of Hif1a with miR-709 mimics, mimics NC, or miRNA inhibitors. The relative levels of firefly luminescence normalized to Renilla luminescence are plotted. Error bars represent SEM (n = 3) *p < 0.05. (E) 293T cells were co-transfected with wild-type (WT) or mutant (MUT) luciferase reporters of lnc4040 with miR-709 mimics, mimics NC, or miRNA inhibitors. The relative levels of firefly luminescence normalized to Renilla luminescence are plotted. Error bars represent SEM (n = 3) *p < 0.05. (F) Expression of HIF1A protein in GCs infected with miR-709 mimics, mimics NC, miR-709 inhibitor, and inhibitor NC and quantified by ImageJ (v1.45); error bars represent SEM (n = 3) **p < 0.01. (G) Expression of HIF1A protein in GCs infected with miRNA mimics, mimics NC, and lnc4040 overexpression vector and quantified by ImageJ (v1.45); error bars represent SEM (n = 3) **p < 0.01.

Materials and Methods

Tissues

All experiments regarding animals were performed in the State Key Laboratory for Conservation and Utilization of Subtropical Agro-bio-resources and were conducted in accordance with its guidelines for the care and use of laboratory animals. Ovaries were collected from a local abattoir in Nanning (from non-pregnant female buffaloes) and transported to the laboratory in sterile saline maintained at 38°C. After rinsing with 75% ethanol, the ovaries were placed on a gauze after high-pressure treatment, then sliced open, and single follicles (5 < n < 8 mm, n = 20) were separated from the ovarian tissue with ophthalmic scissors and repeatedly rinsed with 0.1 M phosphate buffer (pH 7.25).

In this experiment, a total of 15 follicles were successfully collected, and eight HFs and seven AFs were identified. Because a certain amount of follicular fluid is needed to detect hormone levels, we selected six HFs and five AFs for testing. Since the GCs collected need to be detected for apoptosis rate, electron microscopy, and subsequent RNA sequencing, we only collected enough GCs from three HFs and three AFs for experiments.

Classification of HFs and AFs

We judged the status of a follicle by the appearance of the follicle. We classified the follicles with yellowish surfaces, rich capillaries and they bright red color, uniform distribution, meanwhile clear follicular fluid as HFs. When the surfaces appeared gray to white, with fewer capillaries and turbidity in the follicular fluid in conjunction with dark masses inside the follicles or the follicles having severe internal flocculation, these were classified as AFs (34). Then, we detected the hormone content of follicular fluid (estrogen and PROG) by ELISA, the apoptosis of GCs by flow cytometric analysis, and the microstructure of GCs by TEM between HFs and AFs.

Detection of E₂ and PROG in Follicular Fluid

The follicular fluid was collected from HFs (n = 6) and AFs (n = 5), centrifuged for 2,000 g for 10 min, and E₂ and PROG were measured in the supernatants using ELISA kits (Cusabio, CSB-E08893b, Wuhan, China), and an enzyme-labeled instrument (Tecan Trading AG, Mannedorf, Switzerland) was used for measurements. All assays were validated in our laboratory by showing parallelism between serial sample dilutions and the provided assay standard curves (ranges, 40–1,000 pg/ml for E₂ and 0.5–30 ng/ml for PROG). Sensitivities of the assays were 40 pg/ml and 0.2 ng/ml, and the intra-assay coefficients of variations were 15% for each.

Detection of Apoptosis of GCs

The FITC Annexin V Apoptosis Detection Kit I (BD Pharmingen, 556547, Franklin Lakes, NJ, USA) was used to measure apoptosis. GCs were collected and washed twice with cold PBS and then resuspended in 1 × binding buffer at a concentration of 1 × 10⁶ cells/ml. Five microliters of FITC Annexin V and 5 μl PI were added, and the cells were gently vortexed and then incubated for 15 min at room temperature (25°C) in the dark. Four hundred microliters of 1 × binding buffer was then added to each tube per test, and the samples were analyzed by flow cytometry (ACCURI C6, BD, Franklin Lakes, NJ, USA) within 1 h.

TEM Staining of GCs

GCs were collected and washed twice with cold PBS, and the cell suspension was fixed at 4°C by mixing with an equal volume of fixing solution: 2.5% glutaraldehyde and 0.1 M phosphate buffer as solvent (Servicebio, G1102, Wuhan, China) for 2–4 h. The cells were transferred to centrifuge tubes and spun to get the cell pellets, and then washed three times with 0.1 M PBS. Dehydration was performed by placing the samples in the following solvents for 15 min each: 50% ethanol, 70% ethanol, 80% ethanol, 90% ethanol, 95% ethanol, two changes of 100% ethanol, and two changes of acetone. The samples were then infiltrated with 1:1 acetone: EMbed 812 (SPI, West Chester, PA 90529-77-4, USA) for 2–4 h followed by 2:1 acetone: EMbed 812 overnight in a desiccator without covering and then pure EMbed 812 for 5–8 h in a 37°C oven. Then, the samples were embedded by baking in a 60°C oven for 48 h. Ultrathin sections (60–80 nm) were cut with an ultramicrotome (Leica UC7, Leica, Wetzlar, Germany). Sections were stained with uranyl acetate in pure ethanol for 15 min and then rinsed with distilled water. Then, they were stained with lead citrate for 15 min and again rinsed with distilled water. The sections were air-dried overnight and observed on a TEM (HT7700, Hitachi, Tokyo, Japan).

RNA Preparation and Sequencing Analysis

Total RNA was extracted from the GCs of HFs and AFs by using the TRIzol™ Reagent (Life Technologies, Mt Waverley, VIC, Australia) according to the manufacturer's protocol. The protocol was as follows: drastic mechanic vibration after add 1 mL of TRIzol^® reagent (50–100 mg tissue or 5–10 × 10⁶ cells), add 0.2 ml of chloroform and leave for 2–3 min (room temperature), centrifuge (12,000 × g, 4°C for 15 min), add 0.5 ml of isopropanol and leave for 5–10 min (room temperature), centrifuge (12,000 × g, 4°C for 10 min) for RNA deposition, and then add 1 ml 75% ethanol for RNA washing. RNA quality and integrity were estimated with an Agilent 2100 Bioanalyzer and RNA 6000 Nano Kit Reagents (Agilent Technologies, Santa Clara, CA, USA). Only high-quality RNA samples (concentration ≥ 250 ng/μl, RIN ≥ 7.0, total content ≥ 20 ng) were used to construct the sequencing libraries.

RNA-Seq

Details of the mRNA-seq, miRNA-seq, and lncRNA-seq methods are described in detail in Supplementary Material.

CeRNA Regulatory Network

Correlations of expression between mRNAs, miRNAs, and lncRNAs were evaluated using the Spearman rank correlation coefficient (SCC). Pairs with SCC < −0.7 were selected such as those with negatively co-expressed lncRNA–miRNA pairs or mRNA–miRNA pairs, or those where both mRNA and lncRNA were miRNA target genes and also where all the RNAs were differentially expressed. Pairs with Pearson correlation coefficients (PCC) > 0.9 were selected as co-expressed lncRNA–mRNA pairs, where both the mRNAs and lncRNAs in the pair were targeted and co-expressed negatively with a common miRNA. As a result, only the gene pairs with a p-value < 0.05 were selected. Pathway significant enrichment analysis used the KEGG pathway as a unit, and the hypergeometric test was used to determine whether the pathway was significantly enriched in differentially expressed transcripts when compared with the background of the whole transcripts.

Quantitative Real-Time PCR (qRT-PCR) Experimental Validation

qRT-PCR was performed to validate the expression levels of DEmRNAs, DElncRNAs, and DEmiRNAs. Total RNA was extracted by using the TRIzol™ Reagent (Life Technologies, Mt Waverley, VIC, Australia) according to the manufacturer's protocol. First-strand cDNA was synthesized from 1 μg of total RNA with the HiScript^® II QRT SuperMix (Vazyme, Cat#R223, Nanjing, China) for qPCR. In addition, 1 μg of small RNA was used for cDNA synthesis using a miRNA 1st-Strand cDNA Synthesis Kit (reaction contains 10 μl of 2 × RT Mix, 2 μl of random hexamers, 2 μl of HiScript II Enzyme Mix, and 1 pg−1 μg total RNA from the samples, added with ddH₂O to a final volume of 20 μl) (Vazyme, Cat#MR101, China) with the stem-loop primers designed using the stem-loop sequence (GTCGTATCCAGGGTCCGAGGTATTCGCACTGGATACGAC) except for the internal reference, U6. qPCR was performed on the Applied Biosystems™ 7500 Real-Time PCR Systems using 10 μl of ChamQ™ Universal SYBR^® qPCR Master Mix (Vazyme, Cat#Q711), lnRcute lncRNA qPCR Detection Kit (Cat#FP402, Tiangen, Beijing, China), and miRNA Universal SYBR^® qPCR Master Mix (Vazyme, Cat#MQ101) in each 20-μl reaction volume per well by following the manufacturer's instructions. The 2^−ΔΔCT method was used to normalize and determine the RNA level relative to an internal reference gene, beta-actin (Cs1g05000.1) or U6 (35). All primers used are listed in Supplementary Table 3.

Dual-Luciferase Reporter Assay

PsiCHECK 2.0-Hif1a WT and psiCHECK 2.0-Hif1a mut 1, 2, and 3 and psiCHECK 2.0-lnc4040 WT and psiCHECK 2.0-lnc4040 mut were constructed and validated by DNA sequencing. 293T cells (5 × 10⁴ per well) were seeded in 24-well plates and grown to a density of 70–80%. Cotransfection of cells was performed using 50 nM miR-709 mimic, NC, and 25 nM miR-709-inhibitor, NC (RiboBio, Guangzhou, China), or a 0.5 μg reporter vector by X-tremeGENE™ HP DNA Transfection Reagent (6366244001, Roche, Basel, Switzerland). Twenty-four hours after transfection, culture medium was removed, cells were gently washed with PBS once, and 100 μl of passive lysis buffer was added and incubated at 37°C for 10–12 min. Cells were lysed by pipetting up and down several times and centrifuged at 5,000 rpm for 4 min to remove debris, and 10–20 μl was used to assay for luciferase activity using dual luciferase reporter assay (GeneCopoeia, LF004, Rockville, MD USA) in a single injector luminometer. The method refers to the previous study (36).

Vector Construction and Transfection

The lncRNA overexpression vector was constructed using the linear sequence of lnc4040 amplified from GCs by PCR. They were then cloned into the pcDNA3.1-EGFP vector in accordance with the manufacturer's protocol using the NheI and KpnI restriction sites. GCs were collected from follicular fluid and cultured to P2 with 10% FBS (Gibco-BRL, Grand Island, NE, USA) added to normal DMEM (Gibco-BRL, Grand Island, USA). miRNA mimics, mimics NC, and lncRNA overexpression vector were transiently transfected into GCs using Lipofectamine 3000 reagent (Invitrogen, Carlsbad, CA, USA). GCs were collected 48 h after transfection, and proteins were extracted for further experiments.

Western Blot Analysis

Total proteins from GCs were homogenized using RIPA buffer (Servicebio, G2002, Wuhan, China). Protein concentrations were determined using the BCA Protein Assay Kit (Solarbio, PC0020, Wuhan, China). Proteins were separated by 10% sodium dodecyl sulfate polyacrylamide gel electrophoresis, transferred to a polyvinylidene fluoride membrane (Millipore, IEVH00005, Burlington, MA, USA), and then incubated with antibodies (HIF1A, Novus, NB100, Littleton, CO, USA; GAPDH, Proteintech, 60004-1-Ig, Wuhan, China) overnight at 4°C and then with HRP-conjugated secondary antibody for 1 h at room temperature. Pictures were captured by an imaging system (UVP).

Statistical Analysis

The quantitative results are presented as means ± SEM based on at least three independent experiments. Significant variance in multiple comparisons was performed using multiple t-tests and one-way ANOVA with GraphPad Prism version 8.0 (GraphPad Software, La Jolla, CA, USA). Differences were considered significant when p-values ≤ 0.05.

Discussion

At present, there is no exact method to judge the difference between healthy and atresia-associated buffalo follicles. Previous judgment on the appearance of sheep follicles was reported (37), and this was applied to naked follicles collected and observed in these studies. The hormone levels in follicular fluid were also used as an important basis for judging whether the follicle was in the process of atresia. The ratio of E₂ to PROG in follicular fluid has been discussed differently in different species. The PROG/E₂ ratio of ≥10 usually indicated follicle atresia in bovine (38), and the E₂/PROG ratio of healthy camel follicles was 22.8 while that of AFs was 0.08 (39), and in swine, follicles with a PROG/E₂ ratio of <5 were always classified as HF (40). Through a comprehensive comparison of the hormone ratios and appearance, it was found that the ratio of E₂ to PROG in follicular fluid for HFs was 5.77 and that for AFs was 1.43 in our study. Follicular atresia is mainly attributed to GC apoptosis (1, 41–43), and the apoptosis rate of GCs in HFs is basically <10%, while that in AFs can be as high as 30% (44–46). In bovine, the apoptosis rate of follicular GCs with <10% apoptotic cells was designated as HFs (44). In this study, the apoptosis rate of GCs in HFs was 4.93 ± 1.59; that is, the proportion of healthy cells was about 95%, which was consistent with previous studies. The accuracy of sequencing and grouping was guaranteed by the comprehensive determination of three indicators.

In this study, whole-transcriptome RNA sequencing was performed for the first time on buffalo ovarian GCs, and the differences between HFs and AFs at the RNA level were analyzed systematically. Our study found 2,936 DEmRNAs (Supplementary Table 2 for details), among which Serpine2, Inha, Hif1a, Fst, Jak3, Cyp19a1, and Vegfa were the top ones with the highest expression levels, and these were highly expressed in HF. Serpine2 was at the highest level in GCs of growing dominant bovine follicles (47). This is in accordance with the report that the expression of Serpine2 was regulated by follicle-stimulating hormone (FSH) and growth factors in non-luteinizing bovine GCs, and the authors proposed that it could regulate atresia in bovine follicles (48). This is consistent with the high expression of Serpine2 in AFs, which further indicates that Serpine2 is needed in GCs of AFs. In vitro studies on follicle wall explants confirmed the significant differential expression of Inha and Fst during follicle development in the laying hen (49). Moreover, several specific mRNAs expressed in GCs during follicular development detected in previous studies were confirmed in this experiment (50–52).

There were also many reports on the regulation of follicular development by miRNA in GCs. A previous study confirmed that miR-21-3p inhibits bovine GC autophagy by targeting VEGFA so as to regulate follicular atresia (53). In addition, miRNA-10001-3p and miRNA-12030-3p were shown to regulate VEGFA and FST, which suggests that both these factors probably act synergistically during buffalo follicular atresia. Hif1a, as a regulator of VEGF, has been shown to induce autophagy of follicular GCs via FSH (31). In our study, it was found that Hif1a expressed highly in HFs, suggesting that it might promote autophagy of GCs in order to maintain the dynamic balance of follicular development and protect the follicles from entering atresia.

One of the uppermost markers of ovarian differentiation, FOXL2, is necessary for the normal development of GCs (30, 54), and its expression was the highest in these cells, followed by stromal cells, and it was not expressed in oocytes (54, 55). It was predicted that miR-212-3p may target FOXL2 in buffalo GCs, and the results of sequencing and quantitative analysis confirmed the high expression of FOXL2 and the low expression of miR-212-3p in HFs. The results seen for SOX9, which has the opposite effect to FOXL2, showed that its expression in AFs was higher than that in HFs and it is a target by miR-424-3p. This finding is important within the context of the regulation of follicular atresia.

By visualizing the network diagram, several pairs of important ceRNA relationships are seen. These include a key miRNA, miR-424-3p, which is predicted to target four DEmRNAs in GCs between HFs and AFs, including FGFR1 (TCONS_00004077), SGK1 (XM_006059282.2), SOX9 (XM_025279387.1), and THBS (XM_025294260.1). In addition, miR-212-3p is predicted to target FOXL2 (XM_025290957.1), JAK3 (TCONS_00098588), and LDHA (XM_006056061.2). Also, miR-709-5p is seen to target Hif1a (TCONS_00111150) and INHBB (XM_025277483.1). All these target mRNAs regulate cell proliferation, apoptosis, and response to hormones and play regulatory roles at different stages of follicular development.

The ceRNA hypothesis was first proposed by Pandolf in 2011. He presented a unifying hypothesis regarding how mRNAs, transcribed pseudogenes, and lncRNAs “communicate” to each other using microRNA response elements (MREs) (56). MiRNAs can bind to the target mRNAs, inhibiting their translation or leading to their degradation, thus achieving the function of posttranscriptional regulation of gene expression. In addition to pseudogenes, mRNAs and lncRNAs can function as ceRNAs.

In this study, 49.6% of DEmRNAs were predicted as the targets of 135 DEmiRNAs, which constituted 4,309 DEmiRNA–DEmRNA interactions. Half of the DEmRNAs were not used during construction of the ceRNA network because the SCCs of expression between mRNA and miRNA were < -0.7. Also, the PCC of the targeting relationship with miRNA was >0.9 and the p-values were < 0.05. These DEmRNAs competitively bound 134 DElncRNAs as ceRNAs. Because miRNAs play a critical role in the connection and regulation of different ceRNAs, particular attention was paid to the following miRNAs: miR-212-3p, miR-424-3p, and miR-709-5p. MiR-212 had been shown to be highly induced by LH in GCs, thus reducing the expression of CTBP1 protein (57). An increase in follicular levels of miR-212 was also observed following administration of an ovulatory dose of hCG (58). Studies on miR-424 showed that it is associated with therapeutic targets for a variety of diseases, such as muscular dystrophy, liver cancer, HPV, and vascular remodeling as well as regulating the STAT, NOTCH, and HIF1A signaling pathways (59–62).

In our study, XR_003106012.1 and XR_003106717.1 were found to competitively bind to miR-424-3p and TCONS_00104040 was competitively bound to miR-709-5p. It suggests that lncRNA may function as miRNA sponges to participate the regulation of follicular development. However, it should be pointed out that most of DElncRNAs were not included in the ceRNA network. Therefore, it is assumed that these lncRNAs are involved in the regulation of follicular development via other mechanisms. These may include a requirement for transcription or splicing of an RNA, due to DNA elements within the lncRNA promoter or gene body that function independently of the transcription, induction of DNA methylation, and modification of chromatin or serve as transcription enhancers (63).

During follicular development, miRNAs act as key regulators to modulate the upregulation or downregulation of important GC processes such as proliferation, apoptosis, differentiation, transcriptional regulation, and modification. As a consequence, through the regulation of GCs, it may be possible to control the development of dominant follicles to mature follicles and ovulate, while the non-dominant or unselected follicles gradually proceed to atresia, so as to regulate the dynamic development of follicles and maintain the balance of ovarian development. In this process, some lncRNAs can act as ceRNAs to competitively bind the MRE of miRNAs, which may indirectly affect the expression of particular mRNAs. Our experimental validation demonstrated that lnc4040 function as miRNA (miR-709) sponges to regulate the abundance of HIF1α. Thus, we proposed that lncRNA can regulate the abundance of key transcription factors in a ceRNA manner, then fine-tune the expression of cell-related genes in GCs during follicle development. However, further work is needed to elucidate their functions in buffalo follicular atresia.

In general, we have carried out an in-depth molecular analysis of GCs in HFs and AFs. Taking the buffalo as an entry point, we can temporarily fill in the gaps in the non-coding RNA data of this species in the field of reproduction. The buffalo could also serve as a model for mono-ovulating mammals, which could be useful for other mono-ovulating livestock. Next studies will investigate the role of non-coding RNAs in regulating the physiological state and function of GCs, thus affecting follicle development and atresia in Chinese buffalo.

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 study was reviewed and approved by Animal Care & Welfare Committee of Guangxi University.

Author Contributions

YP, SY, DS, and YD designed the research. YP, JC, QL, and QX performed the research. YP, RZ, and JL analyzed the data. SY provided the experimental material. YP, DS, and YD wrote the paper. DS and YD provided the supervision and funding. All authors contributed to the article and approved the submitted version.

Funding

The research was funded by the National Natural Science Foundation of China (grant nos. 31860644 and 31760334), Natural Science Foundation of Guangxi Province (grant nos. 2019GXNSFDA185002, AA17204051, and 2018GXNSFAA281007), and Nanning Scientific Research Technological Development Foundation (grant nos. 20192087 and 20194147).

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.

Acknowledgments

We are also grateful to Dr. Dev Sooranna, Imperial College London, for the English language edits of the manuscript.

Supplementary Material

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

References

1. Tilly JL, Kowalski KI, Johnson AL, Hsueh AJ. Involvement of apoptosis in ovarian follicular atresia and postovulatory regression. Endocrinology. (1991) 129:2799–801. doi: 10.1210/endo-129-5-2799

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Warriach HM, McGill DM, Bush RD, Wynn PC, Chohan KR. A review of recent developments in buffalo reproduction - a review. Asian Austral J Anim. (2015) 28:451–5. doi: 10.5713/ajas.14.0259

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Barile VL. Improving reproductive efficiency in female buffaloes. Livest Prod Sci. (2005) 92:183–94. doi: 10.1016/j.livprodsci.2004.06.014

CrossRef Full Text | Google Scholar

4. Perera B. Reproductive cycles of buffalo. Anim Reprod Sci. (2011) 124:194–9. doi: 10.1016/j.anireprosci.2010.08.022

CrossRef Full Text | Google Scholar

5. 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:493–514. doi: 10.1262/jrd.50.493

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Bhardwaj JK, Sharma RK. Scanning electron microscopic changes in granulosa cells during follicular atresia in Caprine ovary. Scanning. (2011) 33:21–4. doi: 10.1002/sca.20217

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Sharma R, Bhardwaj J. Granulosa cell apoptosis in situ in caprine ovary. Biology. (2007) 7:1111–4.

Google Scholar

8. Sharma R, Bhardwaj J. In situ evaluation of granulosa cells during apoptosis in caprine ovary. Int J Integr Biol. (2009) 5.

Google Scholar

9. Takagi K, Yamada T, Miki Y, Umegaki T, Nishimura M, Sasaki J. Histological observation of the development of follicles and follicular atresia in immature rat ovaries. Acta Med Okayama. (2007) 61:283–98. doi: 10.18926/AMO/32892

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Rodgers RJ, Irving-Rodgers HF. Morphological classification of bovine ovarian follicles. Reproduction. (2010) 139:309–18. doi: 10.1530/REP-09-0177

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Maeda A, Goto Y, Matsuda-Minehata F, Cheng Y, Inoue N, Manabe N. Changes in expression of interleukin-6 receptors in granulosa cells during follicular atresia in pig ovaries. J Reprod Develop. (2007) 53:727–36. doi: 10.1262/jrd.19011

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Eppig JJ, Schroeder AC. Capacity of mouse oocytes from preantral follicles to undergo embryogenesis and development to live young after growth, maturation, and fertilization in vitro. Biol Reprod. (1989) 41:268–76. doi: 10.1095/biolreprod41.2.268

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Quirk SM, Harman RM, Cowan RG. Regulation of Fas antigen (Fas, CD95)-mediated apoptosis of bovine granulosa cells by serum and growth factors. Biol Reprod. (2000) 63:1278–84. doi: 10.1095/biolreprod63.5.1278

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Choi J, Jo M, Lee E, Choi D. Induction of apoptotic cell death via accumulation of autophagosomes in rat granulosa cells. Fertil Steril. (2011) 95:1482–6. doi: 10.1016/j.fertnstert.2010.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Choi JY, Jo MW, Lee EY, Yoon BK, Choi DS. The role of autophagy in follicular development and atresia in rat granulosa cells. Fertil Steril. (2010) 93:2532–7. doi: 10.1016/j.fertnstert.2009.11.021

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Choi J, Jo M, Lee E, Choi D. AKT is involved in granulosa cell autophagy regulation via mTOR signaling during rat follicular development and atresia. Reproduction. (2014) 147:73–80. doi: 10.1530/REP-13-0386

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Evans AC, Ireland JL, Winn ME, Lonergan P, Smith GW, Coussens PM, et al. Identification of genes involved in apoptosis and dominant follicle development during follicular waves in cattle. Biol Reprod. (2004) 70:1475–84. doi: 10.1095/biolreprod.103.025114

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Girard A, Dufort I, Douville G, Sirard M-A. Global gene expression in granulosa cells of growing, plateau and atretic dominant follicles in cattle. Reprod Biol Endocrinol. (2015) 13:17. doi: 10.1186/s12958-015-0010-7

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Sisco B, Hagemann LJ, Shelling AN, Pfeffer PL. Isolation of genes differentially expressed in dominant and subordinate bovine follicles. Endocrinology. (2003) 144:3904–13. doi: 10.1210/en.2003-0485

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Zielak AE, Forde N, Park SD, Doohan F, Coussens PM, Smith GW, et al. Identification of novel genes associated with dominant follicle development in cattle. Reprod Fertil Dev. (2007) 19:967–75. doi: 10.1071/RD07102

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Sontakke SD, Mohammed BT, McNeilly AS, Donadeu FX. Characterization of microRNAs differentially expressed during bovine follicle development. Reproduction. (2014) 148:271–83. doi: 10.1530/REP-14-0140

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Lee JT. Epigenetic regulation by long non-coding RNAs. Science. (2012) 338:1435–9. doi: 10.1126/science.1231776

CrossRef Full Text | Google Scholar

23. Zhang F-L, Li N, Wang H, Ma J-M, Shen W, Li L. Zearalenone exposure induces the apoptosis of porcine granulosa cells and changes long noncoding RNA expression to promote antiapoptosis by activating the JAK2-STAT3 pathway. J Agric Food Chem. (2019) 67:12117–28. doi: 10.1021/acs.jafc.9b05189

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Feng X, Li F, Wang F, Zhang G, Pang J, Ren C, et al. Genome-wide differential expression profiling of mRNAs and lncRNAs associated with prolificacy in Hu sheep. Biosci Rep. (2018) 38:BSR20171350. doi: 10.1042/BSR20171350

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Yerushalmi GM, Salmon-Divon M, Yung Y, Maman E, Kedem A, Ophir L, et al. Characterization of the human cumulus cell transcriptome during final follicular maturation and ovulation. Mol Hum Reprod. (2014) 20:719–35. doi: 10.1093/molehr/gau031

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Zhang Y, Yan Z, Qin Q, Nisenblat V, Chang HM, Yu Y, et al. Transcriptome landscape of human folliculogenesis reveals oocyte and granulosa cell interactions. Mol Cell. (2018) 72:1021–34 e4. doi: 10.1016/j.molcel.2018.10.029

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, et al. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. (2013) 41:e166. doi: 10.1093/nar/gkt646

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Kong L, Zhang Y, Ye Z-Q, Liu X-Q, Zhao S-Q, Wei L, et al. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. (2007) 35:W345–9. doi: 10.1093/nar/gkm391

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Schmidt D, Ovitt CE, Anlag K, Fehsenfeld S, Gredsted L, Treier AC, et al. The murine winged-helix transcription factor Foxl2 is required for granulosa cell differentiation and ovary maintenance. Development. (2004) 131:933–942. doi: 10.1242/dev.00969

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Zhou J, Yao W, Li C, Wu W, Li Q, Liu H. Administration of follicle-stimulating hormone induces autophagy via upregulation of HIF-1alpha in mouse granulosa cells. Cell Death Dis. (2017) 8:e3001. doi: 10.1038/cddis.2017.371

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Szymanska M, Manthe S, Shrestha K, Girsh E, Harlev A, Kisliouk T, et al. Sirtuin-1 inhibits endothelin-2 expression in human granulosa-lutein cells via HIF1A and epigenetic modifications. Biol Reprod. (2020) 2020:ioaa199. doi: 10.1093/biolre/ioaa199

CrossRef Full Text | Google Scholar

33. Yalu R, Oyesiji AE, Eisenberg I, Imbar T, Meidan R. HIF1A-dependent increase in endothelin 2 levels in granulosa cells: role of hypoxia, LH/cAMP, and reactive oxygen species. Reproduction. (2015) 149:11–20. doi: 10.1530/REP-14-0409

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Sharma RK, Singh R, Bhardwaj JK, Saini S. Topographic and ultrastructural variations in isthmus segment of oviduct during oestrous cycle in Caprines. Scanning. (2013) 35:344–8. doi: 10.1002/sca.21073

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2[-Delta Delta C(T)] method. Methods. (2001) 25:402–8. doi: 10.1006/meth.2001.1262

CrossRef Full Text | Google Scholar

36. Hong LJ, Gu T, He YJ, Zhou C, Hu Q, Wang XW, et al. Genome-wide analysis of circular RNAs mediated ceRNA regulation in porcine embryonic muscle development. Front Cell Dev Biol. (2019) 7:289. doi: 10.3389/fcell.2019.00289

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Moor RM, Hay MF, Dott HM, Cran DG. Macroscopic identification and steroidogenic function of atretic follicles in sheep. J Endocrinol. (1978) 77:309–18. doi: 10.1677/joe.0.0770309

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Grimes RW, Matton P, Ireland JJ. A comparison of histological and non-histological indices of atresia and follicular function. Biol Reprod. (1987) 37:82–8. doi: 10.1095/biolreprod37.1.82

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Kafi M, Maleki M, Davoodian N. Functional histology of the ovarian follicles as determined by follicular fluid concentrations of steroids and IGF-1 in Camelus dromedarius. Res Vet Sci. (2015) 99:37–40. doi: 10.1016/j.rvsc.2015.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Zhang J, Liu Y, Yao W, Li Q, Liu H, Pan Z. Initiation of follicular atresia: gene networks during early atresia in pig ovaries. Reproduction. (2018) 156:23–33. doi: 10.1530/REP-18-0058

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Palumbo A, Yeh J. In situ localization of apoptosis in the rat ovary during follicular atresia. Biol Reprod. (1994) 51:888–95. doi: 10.1095/biolreprod51.5.888

PubMed Abstract | CrossRef Full Text | Google Scholar

42. McRae RS, Johnston HM, Mihm M, O'Shaughnessy PJ. Changes in mouse granulosa cell gene expression during early luteinization. Endocrinology. (2005) 146:309–17. doi: 10.1210/en.2004-0999

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Irving-Rodgers HF, van Wezel IL, Mussard ML, Kinder JE, Rodgers RJ. Atresia revisited: two basic patterns of atresia of bovine antral follicles. Reproduction. (2001) 122:761–75. doi: 10.1530/reprod/122.5.761

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Feng WG, Sui HS, Han ZB, Chang ZL, Zhou P, Liu DJ, Tan JH, et al. Effects of follicular atresia and size on the developmental competence of bovine oocytes: a study using the well-in-drop culture system. Theriogenology. (2007) 67:1339–50. doi: 10.1016/j.theriogenology.2007.01.017

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Yu YS, Sui HS, Han ZB, Li W, Luo MJ, Tan JH. Apoptosis in granulosa cells during follicular atresia: relationship with steroids and insulin-like growth factors. Cell Res. (2004) 14:341–6. doi: 10.1038/sj.cr.7290234

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Blondin P, Dufour M, Sirard MA. Analysis of atresia in bovine follicles using different methods: flow cytometry, enzyme-linked immunosorbent assay, and classic histology. Biol Reprod. (1996) 54:631–7. doi: 10.1095/biolreprod54.3.631

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Bédard J, Brûlé S, Price CA, Silversides DW, Lussier JG. Serine protease inhibitor-E2 (SERPINE2) is differentially expressed in granulosa cells of dominant follicle in cattle. Mol Reprod Dev. (2003) 64:152–65. doi: 10.1002/mrd.10239

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Cao M, Nicola E, Portela VM, Price CA. Regulation of serine protease inhibitor-E2 and plasminogen activator expression and secretion by follicle stimulating hormone and growth factors in non-luteinizing bovine granulosa cells in vitro. Matrix Biol. (2006) 25:342–54. doi: 10.1016/j.matbio.2006.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Lovell TM, Gladwell RT, Groome NP, Knight PG. Ovarian follicle development in the laying hen is accompanied by divergent changes in inhibin A, inhibin B, activin A and follistatin production in granulosa and theca layers. J Endocrinol. (2003) 177:45–55. doi: 10.1677/joe.0.1770045

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Richani D, Constance K, Lien S, Agapiou D, Stocker WA, Hedger MP, et al. Cumulin and FSH cooperate to regulate inhibin B and activin B production by human granulosa-lutein cells in vitro. Endocrinology. (2019) 160:853–62. doi: 10.1210/en.2018-01026

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Ndiaye K, Castonguay A, Benoit G, Silversides DW, Lussier JG. Differential regulation of Janus kinase 3 (JAK3) in bovine preovulatory follicles and identification of JAK3 interacting proteins in granulosa cells. J Ovarian Res. (2016) 9:71. doi: 10.1186/s13048-016-0280-5

PubMed Abstract | CrossRef Full Text | Google Scholar

52. McFee RM, Rozell TG, Cupp AS. The balance of proangiogenic and antiangiogenic VEGFA isoforms regulate follicle development. Cell Tissue Res. (2012) 349:635–47. doi: 10.1007/s00441-012-1330-y

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Ma L, Zheng Y, Tang X, Gao H, Liu N, Gao Y, et al. miR-21-3p inhibits autophagy of bovine granulosa cells by targeting VEGFA via PI3K/AKT signaling. Reproduction. (2019) 158:441–52. doi: 10.1530/REP-19-0285

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Cocquet J, Pailhoux E, Jaubert F, Servel N, Xia X, Pannetier M, et al. Evolution and expression of FOXL2. J Med Genet. (2002) 39:916–21. doi: 10.1136/jmg.39.12.916

CrossRef Full Text | Google Scholar

55. Pannetier M, Fabre S, Batista F, Kocer A, Renault L, Jolivet G, et al. FOXL2 activates P450 aromatase gene transcription: towards a better characterization of the early steps of mammalian ovarian development. J Mol Endocrinol. (2006) 36:399–413. doi: 10.1677/jme.1.01947

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language? Cell. (2011) 146:353–8. doi: 10.1016/j.cell.2011.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Fiedler SD, Carletti MZ, Hong X, Christenson LK. Hormonal regulation of MicroRNA expression in periovulatory mouse mural granulosa cells. Biol Reprod. (2008) 79:1030–7. doi: 10.1095/biolreprod.108.069690

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Schauer SN, Sontakke SD, Watson ED, Esteves CL, Donadeu FX. Involvement of miRNAs in equine follicle development. Reproduction. (2013) 146:273–82. doi: 10.1530/REP-13-0107

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Tian Q, Li Y, Wang F, Li Y, Xu J, Shen Y, et al. MicroRNA detection in cervical exfoliated cells as a triage for human papillomavirus-positive women. J Natl Cancer Inst. (2014) 106:dju241. doi: 10.1093/jnci/dju241

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Han H, Du Y, Zhao W, Li S, Chen D, Zhang J, et al. PBX3 is targeted by multiple miRNAs and is essential for liver tumour-initiating cells. Nat Commun. (2015) 6:8271. doi: 10.1038/ncomms9271

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Ghosh G, Subramanian IV, Adhikari N, Zhang X, Joshi HP, Basi D, et al. Hypoxia-induced microRNA-424 expression in human endothelial cells regulates HIF-α isoforms and promotes angiogenesis. J Clin Invest. (2010) 120:4141–54. doi: 10.1172/jci42980

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Connolly M, Paul R, Farre-Garros R, Natanek SA, Bloch S, Lee J, et al. miR-424-5p reduces ribosomal RNA and protein synthesis in muscle wasting. J Cachexia Sarcopenia Muscle. (2018) 9:400–16. doi: 10.1002/jcsm.12266

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Kopp F, Mendell JT. Functional classification and experimental dissection of long noncoding RNAs. Cell. (2018) 172:393–407. doi: 10.1016/j.cell.2018.01.011

PubMed Abstract | CrossRef Full Text | Google Scholar