Integrated Analysis of miRNA-mRNA Interaction Network in Porcine Granulosa Cells Undergoing Oxidative Stress

Oxidative stress (OS), a common intracellular phenomenon induced by excess reactive oxygen species (ROS) generation, has been shown to be associated with mammalian ovarian follicular development blockage and granulosa cell (GC) impairment. However, the mechanism involved in these effects remains unknown, and the effect of OS on the transcriptome profiles in porcine GCs has not been fully characterized. In this study, we found that hydrogen peroxide-mediated oxidative stress induced porcine GC apoptosis and impaired cell viability. Moreover, RNA-seq analysis showed that oxidative stress induced dramatic changes in gene expression in porcine GCs. A total of 2025 differentially expressed genes (DEGs) were identified, including 1940 DEmRNAs and 55 DEmiRNAs. Functional annotation showed that the DEGs were mainly associated with cell states and function regulation. In addition, multiple hub genes (FOXO1, SOD2, BMP2, DICER1, BCL2L11, FZD4, ssc-miR-424, and ssc-miR-27b) were identified by constructing protein-protein interaction and DEmiRNA-DEmRNA regulatory networks. Furthermore, a gene-pathway-function coregulatory network was established and demonstrated that these hub genes were enriched in FoxO, TGF-β, Wnt, PIK3-Akt, MAPK, and cAMP signaling pathways, which play important roles in regulating cell apoptosis, cell proliferation, stress responses, and hormone secretion. The current research provides a comprehensive perspective of the effects of oxidative stress on porcine GCs and also identifies potential therapeutic targets for oxidative stress-induced female infertility.


Introduction
In mammalian ovaries, less than 1 percent of the follicles are mature and capable of ovulation, whereas the majority of the follicles undergo atresia and degeneration during folliculogenesis and follicular development [1]. It is generally believed that the fate of follicles is determined by the state of the follicular granulosa cells (GCs) [2,3]. Recent reports have suggested that follicular atresia is mainly attributed to the apoptosis of GCs [4] and nonapoptotic forms of programmed cell death [5], which are mediated by a complex regulatory network that consists of multiple factors, including environmental factors [6], homeostasis [7], steroid hormones [8], cytokine [9], and epigenetic regulators [10]. Besides, accumulating evidence shows that the reactive oxygen species-(ROS-) induced oxidative stress also plays an important role in regulating the state and function of granulosa cells and even causes several anovulatory disorders [11,12].
Oxidative stress, a common phenomenon in mammalian cells, is mainly caused by excessive ROS accumulation due to redox unbalance and involved in multiple critical biological processes [13,14]. ROS are the natural byproducts of intracellular aerobic metabolism occurring into the mitochondria. Basel ROS concentration in normal cells can be beneficial for maintaining physiological functions, but excessive ROS accumulation disrupts cellular homeostasis and leads to oxidative stress-induced cellular damage and mitochondrial dysfunction [15][16][17]. It has been reported that excessive ROS levels are generated and accumulated in cells undergoing dramatic changes or processes that have a high aerobic energy metabolism requirement, such as autophagy, endoplasmic reticulum stress, carcinogenesis, and reproduction impairment [18][19][20]. During follicular development, metabolic rates accelerate to meet the increased demand for nutrients and energy, which inevitably generates excessive ROS and further induces oxidative stress in follicles [21].
Previous studies using hydrogen peroxide-(H 2 O 2 -) treated mouse models have shown that oxidative damage can block GC development and trigger follicular atresia [22]. However, the underlying mechanism of oxidative stress-induced GC injury and follicular atresia remains largely unknown.
In the current study, we attempted to identify the core RNA molecules and crucial pathways involved in the response of porcine GCs to oxidative stress, by constructing mRNA and miRNA expression profiles through high-throughput sequencing technology. The differentially expressed (DE) miRNA-mRNA regulatory axis and gene-pathway-function interaction network associated with H 2 O 2 -induced oxidative stress were also established. These data will lay a preliminary foundation for further investigation of the biological mechanisms of oxidative stress-induced porcine GC apoptosis and provide opportunities for the future development of molecular-targeted therapy for oxidative stress-induced female infertility.

Materials and Methods
2.1. Cell Culture and H 2 O 2 Treatment. Porcine granulosa cells were derived from healthy ovarian follicles (diameter 3-6 mm) by using syringe with a 22-gauge needle and cultured into a DMEM/F12 medium with 10% fetal bovine serum (FBS) at 37°C in a 5% CO 2 incubator as previously described [10,23]. For H 2 O 2 treatment, the medium was first changed by DMEM/F12 without FBS for 12 h and then H 2 O 2 were added into the medium with the final concentration at 150 μM for 3 h. The morphological features of porcine granulosa cells after H 2 O 2 treatment were observed and recorded in Figure 1S. This study was reviewed and approved by the Animal Ethics Committee of Nanjing Agricultural University, China (SYXK (Su) 2015-0656).

ROS Measurement.
A Reactive Oxygen Species Assay Kit (#S0033, Beyotime, Haimen, China) was used to measure ROS levels in porcine GCs according to the manufacturer's instructions. Briefly, dichloro-dihydro-fluorescein diacetate (DCFH-DA) was diluted 1 : 1000 with a DMEM/F12 medium without FBS, to a concentration of 10 μM. Porcine GCs were then submerged in 1 mL of DCFH-DA (10 μM) for 20 min at 37°C. After incubation, cells were washed three times with a non-FBS medium and H 2 O 2 was added to the medium for 2 h at a final concentration of 150 μM. The entire procedure was performed in a darkroom. ROS levels in porcine GCs after treatment with H 2 O 2 were detected by fluorescence microscopy and flow cytometric analysis.

Cell Apoptosis and Viability Analysis.
To detect the apoptosis rate of porcine GC, Annexin V-FITC and propidium iodide (PI) were used according to the protocol (Vazyme Biotech Co., Ltd). Briefly, 2 × 10 5 cells were collected and dyed with Annexin V-FITC and PI for 10 min in a darkroom, which were then sorted by flow cytometry with a cell counting machine (Becton Dickinson). FlowJo software (TreeStar) was used for data analysis. The apoptosis rate was calculated by the ratio of (cell numbers in Q2 and Q3)/total cells. For cell viability detection, Cell Counting Kit-8 (CCK-8, #K1018, APExBIO, USA) was used following the manufacturer's instructions. Briefly, porcine GCs were seeded into 96-well cell plates, and after treatment with PBS or H 2 O 2 , 10 μL CCK-8 was added into the medium and incubated at 37°C for 2 h. Then, the absorbance (optical density, OD value) of porcine GCs was detected by using a microplate reader.
2.4. RNA Extraction, Library Preparation, and Sequencing. Porcine granulosa cells after H 2 O 2 treatment were collected, and total RNA were extracted using the High purity RNA extraction kit (#RP1202, BioTeke Corporation, Beijing, China). The extracted RNA was run on 1.5% agarose gels to detect degradation and contamination; their quantity and quality were also estimated by the NanoDrop 3000 system (Agilent Technologies, USA). The cDNA libraries for sequencing were prepared according to the modified method [24] and subsequently sent for sequencing by Biomarker Technologies Co. Ltd., Beijing, China. The proportion of each category in relation to total clean tags was determined, and sequences obtained from Sus Scrofa RefSeq (Sscrofa 11.1) databases were used for reads mapping. The raw transcriptome sequencing data have been submitted to Sequence Read Archive (SRA) database of NCBI (accession number SUB6086396).

Bioinformatics Analysis
2.5.1. Differentially Expressed Gene Analysis. Raw data were extracted and low-quality reads were removed through Perl scripts developed by Biomarker Technologies Co. Ltd. (Beijing, China). After quantile normalization, sequencing data were log 2 transformed and expression levels of each gene in all samples were normalized as fragments per kilobase of transcript sequence per million mapped reads (FPKM). Differentially expressed genes (DEGs) between different samples were detected by using the DESeq R package (1.10.1), and P values were adjusted to control for the false discovery rate (FDR). Significant DEGs (DEmRNAs and DEmiRNAs) were identified using | log 2 ðfold changeÞ| ≥ 1 and adjusted FDR < 0:05 as cut-off criteria.

Functional Annotation of Differentially Expressed
Genes. To assess the functions, roles, and biological processes of the DEGs and their enrichment in different biological pathways, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed by using online database, DAVID (the Database for Annotation, Visualization and Integrated Discovery, version 6.7, https://david.ncifcrf.gov/). A significance level of P < 0:05 and an enrichment score > 2 were set as the thresholds. For the functional annotation and miRNA/pathway clustering of DEmiRNAs, the DIANA-miRPath v3.0 database (http://www.microrna.gr/miRPathv3/) was used according to the published instructions [25].

Protein-Protein Interaction Network Construction.
DEG-encoded proteins were chosen for construction of a protein-protein interaction (PPI) network. Briefly, potential or confirmed protein interactions were generated and analyzed using STRING online software (http://string-db .org/), with a minimum confident interaction score > 0:9 (0-1) and the interacted protein amount ≥ 1. The PPI network was then visualized using Cytoscape v3.7.1 software. Hub genes were identified as the nodes with higher degrees (top 5%) using CytoHubba functions. The Cytoscape software MCODE package was performed to analyze the modules in the PPI network.
2.5.4. Construction of DEmiRNA-DEmRNA Regulatory Network and Functional Assessment. Target genes of DEmiR-NAs were first predicted using the miRWalk v3.0 database (http://mirwalk.umm.uni-heidelberg.de/search_mirnas/), microRNA.org (http://microrna.org/), TargetScan (http:// www.targetscan.org/), and miRDB (http://www.mirdb.org/). The common genes both belong to DEmiRNA targets and DEmRNAs which have inverse expression relationship with DEmiRNAs were chosen to analyze miRNA-mRNA pairs and retained for DEmiRNA-DEmRNA regulatory network construction through using Cytoscape software. To assess the functions of these miRNA-mRNA regulatory networks, DAVID was used for GO annotation and KEGG pathway analysis of the differentially expressed target genes. Venn diagram indicating the intersected genes was generated by a Draw Venn Diagram online tool (http://bioinformatics.psb .ugent.be/webtools/Venn/). Hub miRNAs were defined as the miRNA nodes with higher degree (top 5%) in the DEmiRNA-DEmRNA regulatory network.
2.6. Quantitative Real-Time PCR Validation. Total RNA from porcine granulosa cells after H 2 O 2 treatment was reverse-transcribed into cDNA with three biological repeats by using HiScript® III RT SuperMix for qPCR (+gDNA wiper, #R323-01, Vazyme Biotech Co., Ltd.) according to the manufacturer's instructions. Several significant DEGs were chosen for sequencing accuracy detection, and qRT-PCR were performed by using AceQ qPCR SYBR Green Master Mix (#Q111-03, Vazyme Biotech Co., Ltd, Nanjing, China) on a StepOne Plus System (Applied Biosystems) with three biological repeats. The experimental data were analyzed using the 2 -ΔΔCT method. The expression levels of coding genes were normalized to that of GAPDH. U6 was chosen as an internal control of miRNAs' expression levels. The primers used here were designed using the primer 5.0 software and listed in Supplementary Table S1.

Statistical and Data
Analysis. Statistical analyses were performed using GraphPad Prism 7.0 (GraphPad Software, CA, USA) and SPSS 20.0 (SPSS, IL, USA). The comparisons were conducted by a two-tailed Student's t-test. P value < 0.05 was considered as statistically significant.

H 2 O 2 Induced
Oxidative Stress in Porcine GCs. In this study, 150 μM H 2 O 2 was used to establish oxidative stress in a porcine GC model. To confirm that the model was constructed successfully, we first analyzed ROS levels in porcine GCs under different treatment conditions ( Figure 1S(a)). ROS levels in porcine GCs treated with 150 μM H 2 O 2 were significantly upregulated and higher than that in the control group (PBS), indicating that excessive ROS were generated and accumulated after 150 μM H 2 O 2 treatment. Besides, we observed that H 2 O 2treated porcine GCs had shrunken appearance with jagged edges, suggesting the loss of membrane integrity and low cell viability of porcine GCs ( Figure 1S(b)). In addition, H 2 O 2 -mediated oxidative stress significantly upregulated porcine GC apoptosis ( Figure 1S(c)) and dramatically inhibited cell viability ( Figure 1S(d)). These observations suggested that 150 μM H 2 O 2 could induce oxidative stress in porcine GCs.

Identification of Differentially Expressed RNAs in Porcine
GCs Treated with H 2 O 2 . To investigate the crucial RNA molecules and pathways involved in the responses of porcine GCs to oxidative stress, a high-throughput sequencing strategy was employed, as shown in Figure 1 Table S3). Besides, 284 novel and 600 function-unknown genes were identified in the sequencing data. In addition, heat maps of these DEGs were generated with hierarchy cluster analysis to show their expression patterns (Figures 1(c) and 1(d)). The top 10 most up-and downregulated DEmRNAs and DEmiRNAs according to fold change are presented in Tables 1 and 2, respectively. Among these, SYVN1 and COX-3 were the most up-and downregulated mRNAs, whereas novel-miR-336 and novel-miR-418 were the most up-and downregulated miRNAs, respectively. To validate the accuracy of the sequencing data, top 10 significantly changed DEmRNAs and DEmiRNAs were selected for qRT-PCR detection and as shown in Figure 1(e), the changes of their expression level after H 2 O 2 treatment were generally similar in qRT-PCR and sequencing data, indicating the high accuracy of our sequencing analysis.

Functional Analysis of Differentially Expressed mRNAs.
To further investigate the role of these DEGs in porcine GCs under oxidative stress, Gene Ontology (GO) analysis was performed to assess the function of 1270 functionknown DEmRNAs using DAVID. This analysis identified 135 significantly altered GO terms (P < 0:05, Supplementary Table S4). As shown in Figure 2(a), three GO categories, including biological process (BP), cell component (CC), and molecular function (MF), were analyzed. The three most enriched GO terms in the BP category were "negative regulation of transcription from RNA polymerase II promoter," "positive regulation of transcription from RNA polymerase II promoter," and "heart development." In CC category, "nucleoplasm," "nucleus," and "cytoplasm" were the three most enriched GO terms, whereas the top three terms in MF were "zinc ion binding," "transcription coactivator activity," and "GTPase activator activity." Furthermore, KEGG pathway enrichment analyses showed   that these DEmRNAs were mainly enriched in pathways associated with regulation of cell state and functions, such as PI3K-Akt, AMPK, cAMP, TGF-β, and FoxO signaling pathways (Supplementary Table S5). The top 20 of the 38 significantly enriched pathways (P < 0:05) is shown in Figure 2(b).

Functional Annotation of Differentially Expressed miRNAs.
To investigate the functions of these 55 DEmiRNAs in porcine GCs treated with H 2 O 2 , their targets were first predicted and GO analysis was performed. As shown in Figure 2(c), three GO categories including biological process, cell component, and molecular function were analyzed. In BP,   the most enriched GO term was "small molecule metabolic process". In CC, "protein complex" was the most enriched GO term and "ion binding" was the most enriched GO term in MF (Supplementary Table S6). In addition, KEGG pathway enrichment analysis was performed and showed that the targets of DEmiRNAs were mainly enriched in estrogen generation, Hippo, TGF-β, FoxO, Ras, and mTOR signaling pathways (Figure 2(d); Supplementary Table S7). Moreover, DEmiRNA annotation and miRNA/pathway clustering were performed and showed that these DEmiRNAs mainly participated in regulating cell pluripotency, cell growth and death, fatty acid metabolism, estrogen biosynthesis, and diseases ( Figure 2S).

DEG Protein-Protein Interaction Analysis and Hub Gene
Identification. After removing the novel and functionunknown DEGs, a PPI network was built and visualized by using STRING online database and the Cytoscape v3.7.1 visualization tool. As shown in Figure 3, there were 483 nodes (380 up-and 103 downregulated genes) and 2499 edges in the PPI network. Among them, 30 nodes with higher degree (top 5%) were considered as hub genes with CREBBP, HIST1H2BD, CDK1, CDC20, PIK3R1, FOXO1, DYNC1H1, SUMO1, CBL, and FN1 being the most significant 10 node degree genes (Supplementary Table S8). Module analysis was conducted and showed that there were four significant modules existing in the PPI network ( Figure 3). KEGG pathway analyses of the genes in these modules showed that they were mainly enriched in pathways involved in regulating the cell cycle, growth, apoptosis, autophagy, oxidative stress, and ubiquitin modification. The top 6 hub genes were selected for qRT-PCR validation. The qRT-PCR and sequencing results were highly consistent, except for the SUMO1 gene that was not significantly affected by H 2 O 2 treatment ( Figure 3S).

DEmiRNA-DEmRNA Regulatory Network and Functional
Assessment. To establish the DEmiRNA-DEmRNA regulatory network, the potential targets of DEmiRNAs were first analyzed and 5439 genes were identified, including 586 common genes with DEmRNAs in the sequencing data (Figure 4(a)).
With the DEmiRNAs (14 up-and 9 downregulated) and common DEmRNAs (46 up-and 24 downregulated) mentioned above, a regulatory network was established including 93 nodes and 124 edges (Figure 4(b), Supplementary  Table S9). In this work, ssc-miR-424 and ssc-miR-27b with higher degrees (top 5%) were considered as hub miRNAs. Functional assessment showed that their targets were mainly enriched in TGF-β, FoxO, Hippo, Wnt, cAMP, PI3K-Akt, and MAPK signaling pathways (Figure 4(c) upper lane). To further assess the effects of the miRNA-mRNA regulatory network on porcine GCs, gene-pathway-function coexpression patterns were analyzed and indicated that the oxidative stress-induced miRNA-mRNA regulatory axes exerted important roles in regulating states (proliferation, survive, and apoptosis) and functions (stress responses and hormone secretion) of porcine GCs (Figure 4(c) lower lane).

Discussion
Ovarian follicle development is a complex process that has been proved to be regulated by multiple follicular fluid factors, such as FSH, IGF-1, and ROS [26][27][28]. Previous studies have shown that ROS levels in follicular fluid are closely related to follicular development, atresia, and ovarian  [29]. Besides, excess ROS induces autophagy in oocytes, which eventually leads to follicular atresia and premature ovarian failure [30,31]. 3-NP has been used as an oxidant in several studies and proves that it causes excess ROS generation and induces mouse follicular atresia and mouse GC (mGC) apoptosis in vivo [32,33]. Meanwhile, proanthocyanidins and gallic acid have been shown to be excellent antioxidants, which further decrease the follicular ROS levels in mice and suppress atresia by inhibiting mGC apoptosis [34]. In vitro, H 2 O 2 -mediated oxidative stress has been proved to significantly increase mGC apoptosis and impair relative functions [26,35]. In the present study, we found that 150 μM H 2 O 2 induced excess ROS generation and oxidative stress in porcine GCs which further induced porcine GC apoptosis and inhibited cell viability in vitro. Accumulating evidence obtained through sequencing technology has shown that oxidative stress significantly changes gene expression in multiple cell types of different species, such as bovine, mouse, and human [36,37]. In this    Supplementary Table S4-S7. study, DEmRNAs were screened and a total of 1970 DEmRNAs were identified in porcine GCs treated with H 2 O 2 . Moreover, multiple hub genes including FOXO1, SOD2, COX-3, BMP2, and FZD4 were identified. Among which, FOXO1, a sensor of ROS, has been proven to be regulated by oxidative stress at different levels [38,39]. Members of TGF-β (BMP2) and Wnt signaling pathway (TCF7 and FZD4) have also been proven to be regulated by oxidative stress in human bone marrow stromal cells (hBMSCs) [40]. In addition, the expression and activity of SOD2 have been reported to be modulated by oxidative stress during mitophagy and vascular hypertension [41,42]. GO and KEGG pathway enrichment analyses demonstrated that these DEmRNAs mainly serve as regulators of porcine GC states and functions. For instance, CDK6, CDK18, CDKN2B, CDKN2C, CCNB1, CCNJ, and CCNK are associated with cell cycle [43][44][45]. BMP2, GDF15, TGF-β1, SMAD5, and PCNA are involved in cell proliferation [46,47]. BCL2L11, FOXO1, FOXO3, MMP2, and SOS1 are apoptotic factors [48,49]. ACVR2B, BMP1, BMP4, and NR5A2 are associated with hormone secretion and cytokine responses [50,51]. These identified DEmRNAs may partially explain the effects of oxidative stress on porcine GC states and functions. Apart from these well-known genes, many function-unknown DEmRNAs were also identified and their functions will be the focus of future investigations. MicroRNAs (miRNAs), one of the most important endogenous epigenetic factors, are a class of 18~24 nt short noncoding RNA that inhibit gene expression at the posttranscriptional level by complementary binding to the 3 ′ -untranslated regions (3 ′ -UTR) of specific target genes [52][53][54]. Thus, miRNAs are involved in regulating numerous biological processes such as cell death, proliferation, autophagy, and various diseases [55][56][57]. Importantly, miRNAs have been shown to be the major mediators affecting cell function under conditions of oxidative stress [58,59]. For example, oxidative stress significantly upregulated various miRNAs, including miR-30b, miR-194, miR-125, and miR-128 which further suppressed proliferation of human fibroblasts and HN cells [59,60]. Besides, Lee et al. showed that oxidative stress-induced hnRNPA2B1 increases miR-17/93 expression in epithelial cells and elicits an innate immune response [61]. In addition, oxidative stress-induced APE1 is required for miR-221/222 processing and regulating the tumor suppressor, PTEN [62]. A previous SOS1 IRS1  7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 B B B B B B B B B B N N N N N N N N N N N N N N N P214 1 P2 2 21 1 21 2 2 21 2 2 2 2 2 NUP NUP P P P P2 P P P P P N N NU NU UP2 P2 P P2 2 2 2 2 P2 2 2 N N N N N N N N N N P214 2 P2 2 2 2 21 21 14 1 2 2 2 2 2 2 NUP N P P NU N P P P P GATAD2B

MYC T T T T T T T T T T T T T T T T T T T T T T T T T T T T T T T T TN N TN T T T T T T T TN N T T T T T TN T T T T T T T T TN N N T TN N N TN T T TN N N N N N N T T T T T TN N N N N N N T T TN AGO2
DPY30 ETS2 PHC3 P P P P P P P P 1 1 1 1 CDKN2D RXRA SUN1 CDKN2C CD40 PIK3CD PI P P P P P P P P PI PI PI P P P P P PI P PI P PIK K IK PI P P P P P P PIK K PLXNA4 P P PL PL L L L P P P P P P P P PL  I I I I I I I I I I I  1 1 1 N N N N N N N N N N

GAN H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H HE H H H H HE HE H H
GLI3 P P P P P P P P P P P P P P P P P P CB CB C 2C 2C C 2 2CB C C P2 P2 2 2 P2 P2 P P P P P P PP P P PP P P P P P PP P P P P P P P P P P P P B B CB C C C CB B B B B B B CB CB CB CB B B B B B 2C C 2C C C 2C C 2C C C C P2 2 P2 P P2 2 2 2 P P PP P PP P P P PP PP P P P P P P P P P P P P P P 2 2 2 2 2 2 2 P P P P P PP PP PP PP B B B B B B B B B B B B B B B B B B B B B B B C C C C C C C C C C C C C C C C C C C C 2C C C C C C C C C C CB B B B P2 2 P2 P P PPP P PP PP P P PP P P P P P P P 2 2 P P P P P B CB CB CB C C C CDC20   study using high-throughput sequencing technology showed that oxidative stress modulates the expression of miRNAs in bovine GCs [63]. In the present study, we identified 55 DEmiRNAs (38 up-and 17 downregulated) in oxidative stress porcine GCs through RNA-seq with the criteria | log 2 ðfold changeÞ| ≥ 1 and FDR < 0:05. Interestingly, we observed that the miRNAs mentioned above also existed in our DEmiRNA library. Apart from these well-known miRNAs, multiple newly identified OS-induced DEmiRNAs may also play vital roles in ovaries. miR-126 has been shown to induce GC apoptosis by directly targeting FSHR in pigs [64]. miR-142 regulates the proliferation and apoptosis of human GC by inhibiting TGFBR1 [65]. The miR-182/183/96 cluster actively participates in sexual differentiation in primordial germ cells [66]. Additionally, in another ongoing study, we showed that miR-130a significantly induces porcine GC apoptosis and follicular atresia by targeting the TGF-β signaling pathway (data not shown). Moreover, several DEmiRNAs, such as miR-141 and miR-210, have been proven to control oxidative stress responses in ovarian cancer cells and endometriotic cells by targeting p38α and BARD1, respectively [67,68], suggesting that several identified DEmiRNAs may function as modulators in response to oxidative stress. We also identified 19 novel, functionunknown DEmiRNAs which require further investigation. It is worth noting that DROSHA and DICRE1 (miRNA-  processing enzyme encoding genes) were significantly upregulated during this process which may partially explain the multiple miRNAs differentially expressed in porcine GCs undergoing oxidative stress. Previous studies have demonstrated that oxidative stress often induces multiple signaling pathways during the regulation of different cellular biological processes. These mainly include the FOXO [26], TGF-β [13], Wnt [69], Hippo [70], and PI3K-Akt signaling pathways [71], which were also enriched in the miRNA-mRNA pathway-function regulatory network. Furthermore, cAMP signaling, MAPK signaling, regulation of pluripotency, and miRNAs in cancer were also enriched. Therefore, we hypothesized that oxidative stress might participate in regulating the morphology, pluripotency, energy metabolism, and small noncoding RNA processing in porcine GCs, which requires confirmation in future research. Although the expression profiles of RNAs (mRNAs and miR-NAs) were discussed and multiple differentially expressed RNAs were verified by qRT-PCR in the present study, further studies are required to support the results in vivo. Moreover, the mechanism involved in the differential expression of RNAs in responses to oxidative stress and how their interaction network affects the state and function of porcine GCs needs to be more precisely described in the future.

Conclusion
In summary, we constructed differentially expressed RNA profiles using RNA-seq technology in this study and demonstrated that dramatic changes in gene expression occurred in porcine GCs under oxidative stress. Functional annotation analysis showed that these DEGs were mainly involved in regulating the state and function of porcine GCs, which was highly consistent with our observations that oxidative stress significantly induced porcine GC apoptosis and dramatically impaired cell viability. A DEmiRNA-DEmRNA pathwayfunction coregulatory network and a PPI network were established and indicated that multiple physiological processes and signaling pathways were involved in the response of porcine GCs to oxidative stress ( Figure 5). The integrated analysis of miRNA-mRNA interaction networks also provides a series of potential therapeutic targets for oxidative stress-induced female infertility.

Data Availability
The data used to support the findings of this study are included within the article and supplementary information file. DEmiRNA/KEGG pathway clustering was analyzed in porcine GCs undergoing oxidative stress. 11 significant enrichment signaling pathways existed in the heat map. Figure  3S: hub gene expression validation, related to Figure 3. The expression levels of 6 hub genes were verified by qRT-PCR.
Black columns indicate data from RNA-seq, and red columns indicate qRT-PCR data. Supplementary