Abstract
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 (). 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 (, ). 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 (). 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) (–). 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 (–). 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 (). Recent studies have shown that the autophagy of GCs partially dominates the occurrence of follicular atresia (–).
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 (–). 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 (). 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 (, ).
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 (). The development of this hypothesis has given mRNAs and non-coding RNAs new and broader biological functions (). 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 ().
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 (). 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 (E2) 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−1, respectively, and the ratio was 5.77. The concentrations of E2 and PROG in AFs were 76.63 ± 1.09 and 53.43 ± 2.68 ng·ml−1, 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
Table 1
| Follicle classification | Estrogen content (ng·ml−1) | Progesterone content (ng·ml−1) | E2/PROG (per) | E2/PROG |
|---|---|---|---|---|
| Healthy follicles (n = 6) | 133.83 | 25.74 | 5.20 | 5.90 ±a |
| 124.43 | 23.99 | 5.19 | ||
| 107.76 | 18.29 | 5.89 | ||
| 167.78 | 27.42 | 6.12 | ||
| 115.70 | 17.55 | 6.59 | ||
| 125.74 | 19.66 | 6.40 | ||
| Atretic follicles (n = 5) | 75.24 | 55.79 | 1.35 | 1.35 ±b |
| 72.83 | 48.93 | 1.49 | ||
| 76.59 | 49.07 | 1.56 | ||
| 79.26 | 62.93 | 1.26 | ||
| 82.20 | 74.36 | 1.11 |
Hormone levels in follicular fluid of HFs and AFs.
The data are the means ± SEM.
Indicate p < 0.01 for HFs and AFs, respectively.
Table 2
| GCs of HFs | GCs of AFs | |||||
|---|---|---|---|---|---|---|
| Q1 (%) | 0.18 | 0.07 | 0.13 | 0.91 | 3.60 | 2.30 |
| Q2 (%) | 5.22 | 5.77 | 0.63 | 8.93 | 5.84 | 21.50 |
| Q3 (%) | 0.92 | 1.11 | 1.14 | 13.00 | 12.80 | 1.86 |
| Q4 (%) | 93.70 | 93.00 | 98.10 | 77.10 | 77.70 | 74.30 |
| Q2+Q3 | 6.14 | 6.88 | 1.77 | 21.93 | 18.64 | 23.36 |
| Apoptosis rate (%) | 4.93 ± 1.59a | 21.31 ± 1.40b | ||||
Apoptotic rate of GCs in HFs and AFs.
The data are the means ± SEM, n = 3.
Indicate p < 0.01 for HFs and AFs, respectively.
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
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
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) (
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

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 (
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 (
Detection of E2 and PROG in Follicular Fluid
The follicular fluid was collected from HFs (n = 6) and AFs (n = 5), centrifuged for 2,000 gfor 10 min, and E2 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 E2 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 × 106 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 × 106 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 ddH2O 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 (
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 × 104 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 (
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 (
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 (
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 (
One of the uppermost markers of ovarian differentiation, FOXL2, is necessary for the normal development of GCs (
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) (
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 (
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 (
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.
Statements
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).
Acknowledgments
We are also grateful to Dr. Dev Sooranna, Imperial College London, for the English language edits of 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.
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.
TillyJLKowalskiKIJohnsonALHsuehAJ. Involvement of apoptosis in ovarian follicular atresia and postovulatory regression. Endocrinology. (1991) 129:2799–801. 10.1210/endo-129-5-2799
2.
WarriachHMMcGillDMBushRDWynnPCChohanKR. A review of recent developments in buffalo reproduction - a review. Asian Austral J Anim. (2015) 28:451–5. 10.5713/ajas.14.0259
3.
BarileVL. Improving reproductive efficiency in female buffaloes. Livest Prod Sci. (2005) 92:183–94. 10.1016/j.livprodsci.2004.06.014
4.
PereraB. Reproductive cycles of buffalo. Anim Reprod Sci. (2011) 124:194–9. 10.1016/j.anireprosci.2010.08.022
5.
ManabeNGotoYMatsuda-MinehataFInoueNMaedaASakamakiKet al. Regulation mechanism of selective atresia in porcine follicles: regulation of granulosa cell apoptosis during atresia. J Reprod Dev. (2004) 50:493–514. 10.1262/jrd.50.493
6.
BhardwajJKSharmaRK. Scanning electron microscopic changes in granulosa cells during follicular atresia in Caprine ovary. Scanning. (2011) 33:21–4. 10.1002/sca.20217
7.
SharmaRBhardwajJ. Granulosa cell apoptosis in situ in caprine ovary. Biology. (2007) 7:1111–4.
8.
SharmaRBhardwajJ. In situ evaluation of granulosa cells during apoptosis in caprine ovary. Int J Integr Biol. (2009) 5.
9.
TakagiKYamadaTMikiYUmegakiTNishimuraMSasakiJ. Histological observation of the development of follicles and follicular atresia in immature rat ovaries. Acta Med Okayama. (2007) 61:283–98. 10.18926/AMO/32892
10.
RodgersRJIrving-RodgersHF. Morphological classification of bovine ovarian follicles. Reproduction. (2010) 139:309–18. 10.1530/REP-09-0177
11.
MaedaAGotoYMatsuda-MinehataFChengYInoueNManabeN. Changes in expression of interleukin-6 receptors in granulosa cells during follicular atresia in pig ovaries. J Reprod Develop. (2007) 53:727–36. 10.1262/jrd.19011
12.
EppigJJSchroederAC. 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. 10.1095/biolreprod41.2.268
13.
QuirkSMHarmanRMCowanRG. Regulation of Fas antigen (Fas, CD95)-mediated apoptosis of bovine granulosa cells by serum and growth factors. Biol Reprod. (2000) 63:1278–84. 10.1095/biolreprod63.5.1278
14.
ChoiJJoMLeeEChoiD. Induction of apoptotic cell death via accumulation of autophagosomes in rat granulosa cells. Fertil Steril. (2011) 95:1482–6. 10.1016/j.fertnstert.2010.06.006
15.
ChoiJYJoMWLeeEYYoonBKChoiDS. The role of autophagy in follicular development and atresia in rat granulosa cells. Fertil Steril. (2010) 93:2532–7. 10.1016/j.fertnstert.2009.11.021
16.
ChoiJJoMLeeEChoiD. AKT is involved in granulosa cell autophagy regulation via mTOR signaling during rat follicular development and atresia. Reproduction. (2014) 147:73–80. 10.1530/REP-13-0386
17.
EvansACIrelandJLWinnMELonerganPSmithGWCoussensPMet al. Identification of genes involved in apoptosis and dominant follicle development during follicular waves in cattle. Biol Reprod. (2004) 70:1475–84. 10.1095/biolreprod.103.025114
18.
GirardADufortIDouvilleGSirardM-A. Global gene expression in granulosa cells of growing, plateau and atretic dominant follicles in cattle. Reprod Biol Endocrinol. (2015) 13:17. 10.1186/s12958-015-0010-7
19.
SiscoBHagemannLJShellingANPfefferPL. Isolation of genes differentially expressed in dominant and subordinate bovine follicles. Endocrinology. (2003) 144:3904–13. 10.1210/en.2003-0485
20.
ZielakAEFordeNParkSDDoohanFCoussensPMSmithGWet al. Identification of novel genes associated with dominant follicle development in cattle. Reprod Fertil Dev. (2007) 19:967–75. 10.1071/RD07102
21.
SontakkeSDMohammedBTMcNeillyASDonadeuFX. Characterization of microRNAs differentially expressed during bovine follicle development. Reproduction. (2014) 148:271–83. 10.1530/REP-14-0140
22.
LeeJT. Epigenetic regulation by long non-coding RNAs. Science. (2012) 338:1435–9. 10.1126/science.1231776
23.
ZhangF-LLiNWangHMaJ-MShenWLiL. 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. 10.1021/acs.jafc.9b05189
24.
ThomsonDWDingerME. Endogenous microRNA sponges: evidence and controversy. Nat Rev Genet. (2016) 17:272–83. 10.1038/nrg.2016.20
25.
FengXLiFWangFZhangGPangJRenCet al. Genome-wide differential expression profiling of mRNAs and lncRNAs associated with prolificacy in Hu sheep. Biosci Rep. (2018) 38:BSR20171350. 10.1042/BSR20171350
26.
YerushalmiGMSalmon-DivonMYungYMamanEKedemAOphirLet al. Characterization of the human cumulus cell transcriptome during final follicular maturation and ovulation. Mol Hum Reprod. (2014) 20:719–35. 10.1093/molehr/gau031
27.
ZhangYYanZQinQNisenblatVChangHMYuYet al. Transcriptome landscape of human folliculogenesis reveals oocyte and granulosa cell interactions. Mol Cell. (2018) 72:1021–34 e4.10.1016/j.molcel.2018.10.029
28.
SunLLuoHBuDZhaoGYuKZhangCet al. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. (2013) 41:e166. 10.1093/nar/gkt646
29.
KongLZhangYYeZ-QLiuX-QZhaoS-QWeiLet al. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. (2007) 35:W345–9. 10.1093/nar/gkm391
30.
SchmidtDOvittCEAnlagKFehsenfeldSGredstedLTreierACet al. The murine winged-helix transcription factor Foxl2 is required for granulosa cell differentiation and ovary maintenance. Development. (2004) 131:933–942. 10.1242/dev.00969
31.
ZhouJYaoWLiCWuWLiQLiuH. Administration of follicle-stimulating hormone induces autophagy via upregulation of HIF-1alpha in mouse granulosa cells. Cell Death Dis. (2017) 8:e3001. 10.1038/cddis.2017.371
32.
SzymanskaMMantheSShresthaKGirshEHarlevAKislioukTet al. Sirtuin-1 inhibits endothelin-2 expression in human granulosa-lutein cells via HIF1A and epigenetic modifications. Biol Reprod. (2020) 2020:ioaa199. 10.1093/biolre/ioaa199
33.
YaluROyesijiAEEisenbergIImbarTMeidanR. HIF1A-dependent increase in endothelin 2 levels in granulosa cells: role of hypoxia, LH/cAMP, and reactive oxygen species. Reproduction. (2015) 149:11–20. 10.1530/REP-14-0409
34.
SharmaRKSinghRBhardwajJKSainiS. Topographic and ultrastructural variations in isthmus segment of oviduct during oestrous cycle in Caprines. Scanning. (2013) 35:344–8. 10.1002/sca.21073
35.
LivakKJSchmittgenTD. Analysis of relative gene expression data using real-time quantitative PCR and the 2[-Delta Delta C(T)] method. Methods. (2001) 25:402–8. 10.1006/meth.2001.1262
36.
HongLJGuTHeYJZhouCHuQWangXWet al. Genome-wide analysis of circular RNAs mediated ceRNA regulation in porcine embryonic muscle development. Front Cell Dev Biol. (2019) 7:289. 10.3389/fcell.2019.00289
37.
MoorRMHayMFDottHMCranDG. Macroscopic identification and steroidogenic function of atretic follicles in sheep. J Endocrinol. (1978) 77:309–18. 10.1677/joe.0.0770309
38.
GrimesRWMattonPIrelandJJ. A comparison of histological and non-histological indices of atresia and follicular function. Biol Reprod. (1987) 37:82–8. 10.1095/biolreprod37.1.82
39.
KafiMMalekiMDavoodianN. 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. 10.1016/j.rvsc.2015.01.001
40.
ZhangJLiuYYaoWLiQLiuHPanZ. Initiation of follicular atresia: gene networks during early atresia in pig ovaries. Reproduction. (2018) 156:23–33. 10.1530/REP-18-0058
41.
PalumboAYehJ. In situ localization of apoptosis in the rat ovary during follicular atresia. Biol Reprod. (1994) 51:888–95. 10.1095/biolreprod51.5.888
42.
McRaeRSJohnstonHMMihmMO'ShaughnessyPJ. Changes in mouse granulosa cell gene expression during early luteinization. Endocrinology. (2005) 146:309–17. 10.1210/en.2004-0999
43.
Irving-RodgersHFvan WezelILMussardMLKinderJERodgersRJ. Atresia revisited: two basic patterns of atresia of bovine antral follicles. Reproduction. (2001) 122:761–75. 10.1530/reprod/122.5.761
44.
FengWGSuiHSHanZBChangZLZhouPLiuDJTanJHet 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. 10.1016/j.theriogenology.2007.01.017
45.
YuYSSuiHSHanZBLiWLuoMJTanJH. Apoptosis in granulosa cells during follicular atresia: relationship with steroids and insulin-like growth factors. Cell Res. (2004) 14:341–6. 10.1038/sj.cr.7290234
46.
BlondinPDufourMSirardMA. Analysis of atresia in bovine follicles using different methods: flow cytometry, enzyme-linked immunosorbent assay, and classic histology. Biol Reprod. (1996) 54:631–7. 10.1095/biolreprod54.3.631
47.
BédardJBrûléSPriceCASilversidesDWLussierJG. Serine protease inhibitor-E2 (SERPINE2) is differentially expressed in granulosa cells of dominant follicle in cattle. Mol Reprod Dev. (2003) 64:152–65. 10.1002/mrd.10239
48.
CaoMNicolaEPortelaVMPriceCA. 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. 10.1016/j.matbio.2006.05.005
49.
LovellTMGladwellRTGroomeNPKnightPG. 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. 10.1677/joe.0.1770045
50.
RichaniDConstanceKLienSAgapiouDStockerWAHedgerMPet 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. 10.1210/en.2018-01026
51.
NdiayeKCastonguayABenoitGSilversidesDWLussierJG. 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. 10.1186/s13048-016-0280-5
52.
McFeeRMRozellTGCuppAS. The balance of proangiogenic and antiangiogenic VEGFA isoforms regulate follicle development. Cell Tissue Res. (2012) 349:635–47. 10.1007/s00441-012-1330-y
53.
MaLZhengYTangXGaoHLiuNGaoYet al. miR-21-3p inhibits autophagy of bovine granulosa cells by targeting VEGFA via PI3K/AKT signaling. Reproduction. (2019) 158:441–52. 10.1530/REP-19-0285
54.
CocquetJPailhouxEJaubertFServelNXiaXPannetierMet al. Evolution and expression of FOXL2. J Med Genet. (2002) 39:916–21. 10.1136/jmg.39.12.916
55.
PannetierMFabreSBatistaFKocerARenaultLJolivetGet 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. 10.1677/jme.1.01947
56.
SalmenaLPolisenoLTayYKatsLPandolfiPP. A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language?Cell. (2011) 146:353–8. 10.1016/j.cell.2011.07.014
57.
FiedlerSDCarlettiMZHongXChristensonLK. Hormonal regulation of MicroRNA expression in periovulatory mouse mural granulosa cells. Biol Reprod. (2008) 79:1030–7. 10.1095/biolreprod.108.069690
58.
SchauerSNSontakkeSDWatsonEDEstevesCLDonadeuFX. Involvement of miRNAs in equine follicle development. Reproduction. (2013) 146:273–82. 10.1530/REP-13-0107
59.
TianQLiYWangFLiYXuJShenYet al. MicroRNA detection in cervical exfoliated cells as a triage for human papillomavirus-positive women. J Natl Cancer Inst. (2014) 106:dju241. 10.1093/jnci/dju241
60.
HanHDuYZhaoWLiSChenDZhangJet al. PBX3 is targeted by multiple miRNAs and is essential for liver tumour-initiating cells. Nat Commun. (2015) 6:8271. 10.1038/ncomms9271
61.
GhoshGSubramanianIVAdhikariNZhangXJoshiHPBasiDet al. Hypoxia-induced microRNA-424 expression in human endothelial cells regulates HIF-α isoforms and promotes angiogenesis. J Clin Invest. (2010) 120:4141–54. 10.1172/jci42980
62.
ConnollyMPaulRFarre-GarrosRNatanekSABlochSLeeJet al. miR-424-5p reduces ribosomal RNA and protein synthesis in muscle wasting. J Cachexia Sarcopenia Muscle. (2018) 9:400–16. 10.1002/jcsm.12266
63.
KoppFMendellJT. Functional classification and experimental dissection of long noncoding RNAs. Cell. (2018) 172:393–407. 10.1016/j.cell.2018.01.011
Summary
Keywords
buffalo (Bubalus bubalis), follicles, granulosa cells, whole-transcriptome, ceRNA
Citation
Pan Y, Yang S, Cheng J, Lv Q, Xing Q, Zhang R, Liang J, Shi D and Deng Y (2021) Whole-Transcriptome Analysis of LncRNAs Mediated ceRNA Regulation in Granulosa Cells Isolated From Healthy and Atresia Follicles of Chinese Buffalo. Front. Vet. Sci. 8:680182. doi: 10.3389/fvets.2021.680182
Received
13 March 2021
Accepted
09 June 2021
Published
14 July 2021
Volume
8 - 2021
Edited by
Jordi Roca, University of Murcia, Spain
Reviewed by
Jitender Kumar Bhardwaj, Kurukshetra University, India; Abouzar Najafi, University of Tehran, Iran
Updates

Check for updates
Copyright
© 2021 Pan, Yang, Cheng, Lv, Xing, Zhang, Liang, Shi and Deng.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Deshun Shi ardsshi@gxu.edu.cnYanfei Deng yanfei-dun@163.com
This article was submitted to Animal Reproduction - Theriogenology, a section of the journal Frontiers in Veterinary Science
†These authors have contributed equally to this work
Disclaimer
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.