Zearalenone was purchased from Sigma Company (St. Louis, MO, United States). Stock solutions of ZEA were fixed by dissolving ZEA in dimethyl sulfoxide (DMSO). DMSO (D12345), fetal bovine serum (FBS; 10100147), M-199 medium (11150-059), penicillin and streptomycin were procured from Gibco Company (Carlsbad, CA, United States).

The mature donkeys’ ovaries used in the experiment were obtained from the Dong E Donkey Production Company (Qingdao, Shandong, China). The ovaries were collected from the slaughterhouse of the company and maintained at 32–35°C for the isolation of GCs. All procedures of animal handling in this study were reviewed and approved by the Ethical Committee (Agreement No. 2017-18) of Qingdao Agricultural University.

Donkey GCs (dGCs) were aspirated from the antral follicles using a 10 ml syringe (Zhu et al., 2012). Standing for 15–18 min, the dGCs were centrifuged at 300 g for 5 min in accordance with previous report (Qin et al., 2015). Then the dGCs were cultured in DMEM medium (HyClone, SH30022.01, Beijing, China) supplement with 10% FBS (10099-141, Gibco, Australia) and 1% penicillin-streptomycin (Hyclone, SV30010) in incubator with 5% CO₂ at 37°C (Duda et al., 2011).

The primary GCs were passaged after culture 48–36 h. To avoid the stress of passage response, drug exposure was performed until 12 h later. The GCs were inoculation in a 6 cm petri dish (Corning, 430166, United States) at a density of 1 × 10⁶ cells. ZEA was added to the cultured medium at final concentrations of the 10 or 30 μM, then the cells were incubated with ZEA for 72 h. The control and 10 μM ZEA group added the same dose of DMSO to 30 μM ZEA group for accuracy.

The GCs were collected and fixed in 4% paraformaldehyde for 2 h, then heated at 42°C for 2 h, finally attached onto a poly-lysine coated slide. For immunofluorescence, the sections were blocked with the BDT (10% goat serum in TBS, 3% BSA) for 35 min, and then incubated overnight with primary antibodies at 4°C (Supplementary Table S1). The sections were then incubated with Cy3/FITC-conjugated goat anti-rabbit secondary antibody (Beyotime, A0208, Nantong, China) at a dilution of 1:50 at 37°C for 1.5 h. Finally, the sections were incubated with Hoechst33342 (Beyotime, C1022) to visualize nuclei of GCs for 3 min at room temperature. The immunosignal was detected using a fluorescence microscope (Olympus, XB51, Japan), the images were captured and analyzed in accordance with cellSens Standard.

TUNEL BrightRed Apoptosis Detection Kit was utilized to evaluate the GCs apoptosis (Vazyme, A11302, Nanjing, China). Briefly, the GCs were fixed with 4% paraformaldehyde for 2 h after ZEA treatment. Observation under fluorescence microscopy was carried out after TUNEL reaction following manufacturer’s recommendations [60 min at 37°C without light, and DNA staining with Hoechst 33342 (Beyotime, C1022) incubation at room temperature for 3 min]. The TUNEL positive cells were detected under the fluorescence microscope (Olympus, XB51). To analyze TUNEL positive cell ratio, more than 2,000 cells were counted in each group, and each group included three biological replicates at least.

Total RNA was extracted from the cultured GCs by the use of an RNAprep pure MicroKit in line with the manufacturer’s instruction (Aidlab, RN07, Beijing, China). And the first-strand cDNA acquisition was utilized a cDNA Synthesis Kit (TransGen, AT311-03, Beijing, China) with reference to previous study (Pan et al., 2011). The reaction program was setup at: 42°C for 15 min, 85°C for 30 s, and 4°C for cooling. Then, RNA sequencing was performed with Hiseq 4000 platform by Novogene (Beijing, China). Three biological replicates were manipulated in each group. Primary RNA-seq data were uploaded to the SRA repository. The SRA index number was SRP139976. The RNA-seq data matrix used for DESeq2 analysis were uploaded to the GEO repository and the GEO accession number was GSE116950¹.

R Bioconductor/DESeq2 package was applied to analyze the dGCs groups (control, 10 and 30 μM ZEA) to identify the difference of gene expression. DESeq2, the negative binomial was applied as the reference distribution and taken own normalization approach for raw counts in differential expression analysis. In other methods to avoid possible biases, the data of differential expression analysis has to be previously normalized (Benjamini et al., 2001; Luo and Brouwer, 2013). Because the design of sequencing contains biological replicates in each group, log2|fold change| is not setup as the filter condition. So p_adj < 0.01 was really considered as statistical significance.

R Bioconductor/clusterProfiler package was applied to analyze functional profiles (KEGG) of differentially expressed genes (DEGs) (Yu et al., 2012). Furthermore, R Bioconductor/Pathview package was applied to visualize KEGG enrichment results. According to the DEGs log2|fold change| value shows enrichment signaling pathway active status. The p-value adjusted used the Benjamini & Hochberg (BH) method (Luo and Brouwer, 2013). Still, the p_adj < 0.05 was considered as statistical significance.

Gene set enrichment analysis (GSEA) does not need to specify a clearly differential expression of gene threshold, the algorithm performs GSEA based on overall trend of the gene raw read count. The file has to be made suitable for GSEA software analysis with R script (Subramanian et al., 2005). The data for GSEA analysis contained nine samples within three groups. The GSEA software option “Collapse dataset to gene symbols” parameters was “false” and option “Permutation type” parameters were “gene set.” Other parameters were set as default as software manual references. The GSEA gene set with FDR q-value <0.05 were defined as significant difference.

Search Tool for the Retrieval of Interacting Genes/Proteins² (STRING) is a database of protein–protein interaction (PPI). This database contains the direct and physically related interactions between known and predicted protein and genes. The R Bioconductor/STRINGdb was applied for PPI of interested DEGs (Franceschini et al., 2013). The Cytoscape software was applied to visualize PPI results (Shannon et al., 2003).

Total RNA and cDNA reverse transcription was extracted as described above. Quantitative real-time PCR (RT-qPCR) was performed using Light-Cycler^® 480 SYBR Green I Master Kit (Roche, Germany) on LightCycler 480 real-time PCR instrument. The RT-qPCR reaction was set as: 10 min at 95°C, followed by 45 cycles of 95°C for 10 s, 60°C for 30 s, 72°C for 30 s, and cooling step at 4°C. Quantitative RT-PCR primers are listed in Supplementary Table S2 and the primer efficiency in the Supplementary Figure S1. The standard curve gene GAPDH in dGCs was used as the reference to normalize the related genes’ mRNA expression. The difference multiples = 2^-ΔΔCT method was employed for relative quantification of PCR. Each gene was expressed as the mean ± standard deviation (SD), which was calculated from independent biological replicates at least three times.

Protein lysates isolated from dGCs were used for western blotting according to the standard protocol (Zhang et al., 2010; Chao et al., 2012). Proteins from each ZEA treatment were separated by 10% SDS-PAGE, then transferred onto the PVDF membranes. After, the membranes were incubated with primary antibodies (Supplementary Table S1) at 4°C overnight. Then the membranes were incubated with secondary antibodies (Beyotime, A0208) at 37°C for 2 h in TBST after rinsing three times with TBST. The final detection of related genes was carried out by AlphaImager^® HP (ProteinSimple, 92-13824-00, United States). The band intensity was quantified using GAPDH as internal control and measured with IPWIN software.

Data are presented by mean ± SD. Different effects between the control and ZEA treatment groups of donkey was statistically determined by One-way ANOVA for multiple comparisons. All analyses were conducted using Graph-Pad Prism analysis software (San Diego, CA, United States). All experiments were repeated at least three times unless otherwise noted. Results were considered statistically significant at p-value <0.05.

Donkey granulosa cells were cultured in vitro and exposed to 10 or 30 μM ZEA for 72 h (Figure 1A). The percentages of TUNEL positive dGCs significantly increased as a result of exposure of ZEA (10 μM: 36.04 ± 1.52%; 30 μM: 44.30 ± 1.33%) compared to that of the control (0 μM: 9.83 ± 0.21%; P < 0.01; Figure 1B). As shown in Figure 1C, the ratios of BAX/BCL2 mRNA expression significantly increased in 10 and 30 μM ZEA exposed dGCs.

We performed RNA-seq to verify the effect of ZEA exposure on dGCs. Based on the criterion FDR <0.05, 14,506 DEGs was observed between the control and the ZEA-treated dGCs (Figure 2A). A total of 7,253 and 6,984 DEGs were noted between the control and 10 and 30 μM ZEA treatment groups, respectively. A total of 269 DEGs were observed between the 10 and 30 μM ZEA treatment groups. In this study, the DEGs between the control and ZEA treatment groups and among the treatment groups were obtained from: 0 μM vs. 10 or 30 μM and 10 vs. 30 μM ZEA group, respectively (Figure 3A).

We obtained three RNA-seq replicates for the 0, 10, and 30 μM ZEA-treated dGCs. The variation of the replicates from the control and ZEA treatment groups are shown in Figure 2B. A heat map was drawn from the DEG results (Figure 2C). In this study, we chose the DEGs from the 0 μM vs. 10 and 30 μM groups.

To explore the potential mechanism of ZEA exposure in dGCs, the STRING database was applied in annotating functional interactions DEGs between the control and ZEA-treatment group. Also PPI networks was visualized by Cytoscape (Supplementary Figures S2A). The center nodes of the networks, some interesting genes, were observed, such as PTEN, AKT, ATM, TGFβ, PI3K, CCND2, HEY2, CDK2, and FDX1. In addition, the cBioPortal was used to provide visualization, analysis of DEGs related to the ovarian cancer from large-scale cancer genomics data sets (Supplementary Figures S2C, S3).

In order to attain functional insights of DEGs, the R package of clusterProfiler was carried out to establish extremely affected KEGG pathways. The Venn diagram were constructed by the results of DEGs (Figure 3A). A total of 5,863 DEGs containing in control and 30 μM ZEA groups were quantified (Figure 3B). The DEGs were significantly enriched in the PI3K-AKT signaling pathway (Count = 125, p_adj = 0.0049), endocytosis (Count = 91, p_adj = 0.049), regulation of actin cytoskeleton (Count = 90, p_adj = 0.0001), and proteoglycans in cancer (Count = 83, p_adj = 0.0009) in ZEA treated dGCs (Figure 3C).

The GSEA analysis showed that OVARIAN_CANCER_LMP and TNF_SIGNALING_VIA_NFKB gene set was enriched by 10 μM ZEA vs. the control of dGCs (Figures 4A,B), while CANCER_HEAD_AND_NECK_VS_CERVICAL and BREAST_CANCER_16Q24_AMPLICON gene set was enriched by the 30 μM ZEA vs. the control in dGCs (Figures 4C,D), respectively. The GSEA enrichment plot revealed a regulated tendency for most of enriched genes in ZEA-treatment compared with the control. The GSEA analysis also verified the KEGG analysis.

In addition, six genes were identified by performing the KEGG and GSEA pathway analysis in ZEA treatment, the PI3K and AKT genes were involved in PI3K-AKT signaling pathway of dGCs upregulating. DGCs in ZEA treatment, another four genes involved in anti-oncogene process downregulating were also identified, such as PTEN, TGFβ, CDK2, and ATM. Finally, CCND2, CDK6, TNF, and TP53 genes in treatment dGCs, were identified involving in cancer process either upregulation or downregulation.

It is hypothesis that 10 and 30 μM ZEA exposure may affect the proliferation and carcinogenesis of dGCs (Figure 4). RT-qPCR was conducted to verify the hypothesis by evaluating the expression of different transcript together with enzymes in the pathway of dGCs between the control and ZEA treatments. RT-qPCR analyses indicated that 10 and 30 μM ZEA exposure significantly downregulated mRNA abundance of PTEN while upregulated PI3K and AKT genes in dGCs (P < 0.05 or P < 0.01; Figure 5A). As illustrated in Figure 5B, ZEA-treated dGCs exhibited lower protein levels of PTEN but higher protein levels of PI3K and AKT compared with that of the control dGCs (P < 0.05 or P < 0.01). Moreover, 10 and 30 μM ZEA exposure significantly down-regulated the mRNA abundance and protein levels of CDK2, TGFβ, and ATM genes than the control dGCs (P < 0.05 or P < 0.01; Figures 6A,B). Interestingly, the ZEA treatments significantly decreased the number of PTEN and TGFβ but significantly increased the number of PI3K immunofluorescence positive dGCs, in comparison to the control (Figures 7A,B).