Transcriptome Changes in Porcine Granulosa Cells in Response to Gossypol Cytotoxicity


 Background: Gossypol (GP) is a polyphenolic compound in cottonseed. In porcine, GP affects female reproduction and the respiratory system. The aim of this study was to identify differentially expressed genes in porcine granulosa cells (GCs) treated with GP (6.25 and 12.5 μM) for 72 h, in vitro . Results: In GP-treated groups (6.25 and 12.5 μM), the expressions of a number of genes were found to be significantly changed. Gene ontology analysis indicated that the identified DEGs were primarily related to the mitotic cell cycle, chromosome, centromeric region and protein binding. Kyoto Encyclopedia of Genes and Genomes pathway analysis of the GP6.25 group revealed that pathways related to the cell cycle, oocyte meiosis, progesterone-mediated oocyte maturation and p53 signaling were the most changed, while that of the GP12.5 group revealed that pathways related to the PI3K-Akt signaling, focal adhesion, HIF-1 signaling, cell cycle and ECM-receptor interaction were the most changed. Genes associated with female reproductive function ( CDK1 , CCNB1 , CPEB1 and MMP3 ), cellular component organization ( BIRC5 , CYP1A1 , TGFβ3 and COL1A2 ) and oxidation-reduction processes ( PRDX6 , MGST1 and SOD3 ) were confirmed to be differently expressed in GP-treated groups. Conclusion: Our findings provide insight regarding changes in GC gene expression in porcine exposed to GP.


Introduction
Cottonseed and cottonseed meal (CM) are used as protein sources for livestock. While CM is most commonly used as feed for ruminants, its use has been increasing in poultry and monogastric animal species. However, both cottonseed and CM are known to contain a toxic substance, namely, gossypol (GP), which limits the use of CM as animal feed. GP (2,2-bi(8-formyl-1,6,7-trihydroxy-5-isopropyl-3-methylnaphthalene)) is a yellowish polyphenolic compound found in the roots, leaves, stems and seeds of the cotton plant genus Gossypium subspecies. It is known to affect male reproduction, inhibiting spermatogenesis and reducing spermatozoid motility and viability [1,2], and the use of GP as a contraceptive agent in humans has been documented [3]. GP is absorbed into the body, and it does not decompose well and continues to accumulate. Absorbed GP appears to have a long half-life and has been detected in various organs and tissues, including plasma, heart, liver, kidney, muscle, testes and ovaries. In general, acute toxicity symptoms of GP in livestock include respiratory distress, impaired body weight gain, anorexia, weakness, apathy and death after several days [4][5][6]. In addition, certain clinical signs of GP poisoning have been attributed to reduced levels of antioxidants in tissues and increased reactive oxygen species (ROS) formation, resulting in lipid peroxidation [7,8]. At high concentrations, GP has also been reported to impair energy production from oxidative metabolism by interfering with the mitochondrial electron transport chain and oxidative phosphorylation [9].
Ovarian follicles are composed of granulosa and theca cells, and oocytes surrounded by them. The growth of granulosa cells (GCs) and oocytes supports early follicle growth; in this process, GCs surround the oocytes and produce ovarian steroid hormones such as estrogen and progesterone [10,11]. Ovarian steroids play important roles for the normal development and function of several organs, including the brain, uterus and mammary glands [12]. However, few studies have investigated the toxicity of GP using GCs, and among them, only changes in reproductive hormones or specific genes have been investigated [13][14][15].
In recent years, with development of the high-throughput technology capable of processing large amounts of data, and cost reduction of next-generation sequencing (NGS), sequence-based transcriptome analysis such as RNA sequencing (RNA-seq) is used not only in human but also in various animals and plants. RNA-seq has been shown to have considerable advantages such as examining transcriptome profile and differentiate between different transcriptional and splicing isoforms [16,17]. Moreover, RNA-seq technology is known to be useful for analyzing expression patterns of gene populations and understanding the regulatory networks of various metabolic pathways of the analyzed genes [18].
To our knowledge, no studies have used RNA-seq to identify how GP toxicity affects GCs in livestock. Especially, previous studies investigating the effects of GP toxicity, with the goal of identifying differences in gene expression, have not been conducted in porcine. Therefore, this study aims to systematically examine changes in gene expression as a result of GP treatment in porcine GCs using RNA-seq.

Methods
All experimental procedures were conducted in accordance with guidelines approved by the Institutional Animal Care and Use Committee of Kangwon National University, Chuncheon, Korea (KIACUC-18-061).
Isolation and GCs primary culture. Prepubertal gilt ovaries were obtained from a local slaughterhouse and transported to the laboratory within 2 h of isolation in a vacuum thermos flask in 0.9% sterile physiological saline at 37°C. Follicular fluid and GCs were aspirated from follicles with a diameter of 3-6 mm that contained clear follicular fluid using 10 ml syringe with 18-gauge needle. The GCs were then transferred to a 15-ml centrifuge tube and centrifuged at 500 × g for 5 min for precipitation. Precipitated cells were resuspended with ammonium chloride (NH4Cl) at 37°C for 1 min to remove red blood cells. The cells were then washed twice with phosphatebuffered saline (PBS) containing 1% penicillin-streptomycin solution. GCs were cultured in DMEM/F12 (Gibco, Carlsbad, CA, USA) medium supplemented with 10% FBS (Gibco) and 1% penicillin-streptomycin solution in a humidified incubator with 5% CO2 at 37°C for 3 days. When a complete monolayer had formed in the primary culture, the GCs were washed with PBS, trypsinised and harvested for additional experiments.
Subculture and GP treatments. Stock solutions of GP were prepared by first dissolving GP in DMSO at a 10mM concentration and it was stored in aliquots at −20°C until further use. The stock solution was then diluted to concentrations of 0.625, 1.25 and 2.5 mM in DMSO before start of experiments. For experiments, the diluted GP was further diluted (1 in 100) when added to the cell culture medium. The obtained cells from primary culture were subcultured into wells or dishes in suitable amount for each experiment with 1 × 10 6 cells/dish in 100 mm culture dishes (SPL Life science, Pocheon, Korea) for RNA-seq, 1 × 10 5 cells/well in 6-well culture plates (SPL Life science) for qRT-PCR and Western blot analysis and 0.5 × 10 4 cells/well in 96-well culture plates (SPL Life science) for cellular proliferation assays. All culture plates and dishes were incubated at 5% CO2 at 37°C. After seeding the cells for 24 h, to study the cytotoxic effects of GP in GCs, GP was added to the medium at final concentrations of 6.25, 12.5 and 25 µM and the cells were incubated for 72 h. The GP untreated group (Ctrl) was incubated with only DMSO at the same concentrations as the GP treatment groups to ensure the accuracy of the experimental design. The experiment was repeated three times with identical conditions to provide three biological replicates for each time point and treatment. RNA extraction, library construction and sequencing analysis. After GP treatment, total RNA was extracted from the GCs using Trizol reagents (Invitrogen, Carlsbad, CA, USA). RNA extraction was performed in duplicate for each sample. The obtained total RNA was assessed using the Agilent 2100 Bioanalyzer (Agilent technologies, CA, USA). The Illumina TruSeq Stranded mRNA Sample Preparation kit (Illumina, San Diego, USA) was used for preparing mRNA sequencing library from the isolated total RNA. All libraries were quantified by qPCR using a CFX96 Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) and sequenced with a paired-end 75 bp plus single 8 bp index reads using a NextSeq 500 sequencers (Illumina). Potentially existing sequencing adapters and low-quality bases (5' and 3' ends) in the raw reads were trimmed by Skewer software [19]. The cleaned high-quality reads after trimming the low-quality bases and sequencing adapters were mapped to the reference genome by STAR software [20]. The calculation of Q30, GC content and sequence duplication levels were based on the clean reads using by FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/).

Identification of DEGs.
The mapped reads for gene expression values were assembled with Cufflinks following the porcine reference genome (susScr11) annotation [21]. The gene annotation of the reference genome from the UCSC genome (https://genome.ucsc.edu) in GTF format was used as a gene model and the expression values were calculated in the fragments per kilobase of transcripts per million mapped reads (FPKM) unit. The differentially expressed genes (DEGs) between the two selected biological conditions were analyzed by Codify software in Cufflinks package [22]. To compare the expression profiles among the samples, the normalized expression values of the selected few hundred DEGs were unsupervised clustered by in-house R scripts. (https://biit.cs.ut.ee/gprofiler), was performed [23]. A corrected p-value of 0.05 was set as the threshold to identify significantly different pathways [24] and p-values were adjusted using the Benjamini-Hochberg method [25].   Table   1). The reads were obtained for about 88.33% of the porcine reference genome (sus scrofa) in NCBI. In addition, the Q30 range was determined to be from 89.5% to 92.1%. These finding suggests that the library quality for each group was good and suitable for analysis. The results of the hierarchical clustering of DEGs between the three groups showed a clear discrimination which is presented as a heat map in Figure 2a. The distribution of FPKM for all the samples were indicated as similar expression level, and most of the mRNAs detected in the RNA-seq analysis were also found to be protein-coding genes. (Figure 2b).  Table 2 shows the most upregulated and downregulated genes overlapping with the GP6.25 and GP12.5 groups.

Discussion
In recent years, RNA-seq, as a method for the sequence-based analysis of transcriptomes, has become a useful tool for analyzing gene expression patterns with high-throughput NGS technology. In addition, RNA-seq allows for the analysis of complex whole eukaryotic transcriptomes with higher reproducibility, wider dynamic range, less bias and lower frequency of false positive signals compared with traditional cDNA microarray technologies [27]. An earlier study has reported that the Pearson correlation between RNA-seq and qRT-PCR could reach > 0.9 [28], suggesting that RNA-seq is an advanced technology that is reliable and reproducible. In addition, several studies using RNA-seq have been conducted on the porcine transcriptome such as the effects of ageing or various antioxidants in GCs [29][30][31]. However, to date, there exist limited studies that have investigated the toxic effects of GP in porcine, and the majority of these studies have focused on differences in internal organ size and changes in reproductive hormones [6,32]. As such, the studies of gene effects on GCs by GP toxicity are currently almost unknown. Therefore, in this study, we investigated the transcriptome profile of the porcine GCs using RNA-seq, with the aim of identifying candidate genes for major functional changes caused by GP cytotoxicity.
Previous studies have shown that cellular proliferation is affected by GP toxicity. Indeed, GP has been shown to inhibit cell growth in a dose-dependent manner in various species, including human. The majority of the studies evaluating the toxic effects of GP at the cellular level have been performed on cancer cells, including those of leukaemia [33], prostate cancer [34], breast cancer [35], ovarian cancer [36], carcinoma [37] and other malignancies [38], with the focus of preventing cancer metastasis. As a result, it has been reported that the proliferation of cancer cells is inhibited by GP toxicity. However, unlike cancer cells, some studies using primary cells have shown a different pattern of cell growth by GP toxicity. Lin et al. [14] reported that GP at 10, 50 and 150 µg/ml significantly decreased the number of bovine oocytes undergoing cumulus expansion cultured at 24 h, by 15%, 50% and 80%, respectively. In another study, however, cellular proliferation was shown to be significantly increased by GP at 5 and 25 µg/ml in porcine GCs cultured for 44 h [39]. In this study, we found that the number of viable cells was dramatically decreased a dose-dependent manner in GP-treated porcine GCs.
A possible reason for the discrepancy between the results of the present study and those from previous studies using primary cells may be related to differences in the culture period. However, Hsiao et al. [40] reported that the population of apoptotic cells increased significantly with increasing GP dose, and the cells were more likely to enter late apoptosis and undergo cell arrest compared with the control group. Collectively, our results suggested that the inhibition of proliferation by GP toxicity may occur after exposure to GP for a sufficient period of time.
In females, granulosa and theca cells, and oocytes are the constituents of mammalian ovarian follicles. GCs play a complex and central role in the development of the ovarian follicle During this process, GCs surround and nurse the oocyte, supporting its maturation. In this study, we identified DEGs from an RNA-seq library of porcine GCs treated with two GP concentrations (6.25 and 12.5 µM). These DEGs were then divided into three groups (GP6.25, GP12.5 and IPC), and we investigated the differences in gene expression caused by GP toxicity through GO and KEGG pathway analysis. We have found that identified DEGs are associated with processes related to cellular component organisation, oxidation-reduction and female reproductive functions. From these results, we selected 11 genes to investigate differences in mRNA and protein expression as a result of GP exposure.
Meiotic progression and oocyte maturation are regulated by activity of the universal cell cycle regulator called the maturation-promoting factor (MPF). In addition, MPF is important for the precise regulation of mitosis and enables chromatin condensation and alignment through the histone phosphorylation regulation. As such, MPF, which plays an important role in the cell cycle, is a heterodimer mainly composed of a regulatory subunit, cyclin B (CCNB) and a catalytic subunit, cyclin-dependent kinase 1 (CDK1) [41,42]. Adhikari et al. [43] reported that the CDK1 gene is necessary to promote the resumption of meiosis. Jiang et al. [44] reported that cyclin and CDK accelerate cell cycle progression, whereas CDK inhibitors slow the progression. CCNB1 is known as a major protein involved in the regulation of oocyte maturation in mammals. The synthesis of CCNB1 is necessary for germinal vesicle breakdown (GVBD) induction in the second meiosis [45,46]. Stanley et al. [47] reported that the levels of cyclin A, B1 and CDK1 were reduced by hexavalent chromium in rat GCs, delaying cell entry and progression through the G2-M phase. In the present study, the mRNA expression of CDK1 and CCNB1 was decreased in a dose-dependent manner in GP-treated GCs, as observed from both RNA-seq and qRT-PCR. We found that protein expression of CCNB1 also decreased in dose-dependently. These results suggest that GP indirectly inhibits the maturation and development of oocytes by delaying cell cycle progression by reducing the levels of CDK1 and CCNB1 in porcine GCs.
Translational activation and polyadenylation of cell cycle regulators including cyclins and CDK, are found primarily in the 3'-UTR of the mRNAs; the cytoplasmic polyadenylation element (CPE) is located in these sequences [48]. The CPE binds to the CPE-binding protein 1 (CPEB1), the CPEB1 is involved in translational activity or inhibition, according to the phosphorylation state [49]. In oocytes, the non-phosphorylated CPEB1 protein is known to interact with Maskin and Pumilio proteins, which are involved in translational repression [50,51]. In previous studies of various species, including human, CPEB has been reported to be involved in both the inhibition and stimulation of CCNB1 mRNA. In this regard, Nakahata et al. [51] reported that the overexpression of Pumilio in Xenopus oocytes leads to inhibition of CCNB1 levels. Taken together, the phosphorylation of CPEB1 suggests that it may be involved in the regulation of CCNB1 activity. In our results, the mRNA expression of CPEB1 was upregulated and that of CCNB1 was downregulated. We have not been able to identify exactly how GP acts on the phosphorylation of CPEB1. However, our results suggested that Pumilio bound to the CPEB1 gene likely regulated the levels of CCNB1. As a result, GP is involved in the expression of CPEB1 and CCNB1 genes in GCs; thus, it seems to be affecting the growth and maturation of oocytes.
In general, the extracellular matrix (ECM) is a natural substrate that supports the function of cells such as intercellular communication, adhesion, migration and proliferation. In particular, it affects cell morphology, follicle development, aggregation, communication and steroidogenesis in the ovary. ECM metabolism, which causes periodic changes throughout the reproductive cycle, is caused by the action of certain types of proteins known as the matrix metalloproteinases (MMPs). MMPs are enzymes that degrade ECM proteins throughout the body to facilitate tissue remodelling [52]. Many studies have shown that MMPs play an important role in follicle development, ovulation, luteinization and luteolysis associated with follicular ECM alterations [53,54]. Zhu et al. [55] reported that the gene expression of MMP1, MMP3, and MMP9 were significantly increased in the ovary during sexual maturation in chicken. In addition, Hui et al. [29] reported that MMP3 expression in the GCs decreased significantly during ageing in porcine. In this study, we identified that the expression of MMP3 protein and mRNA were significantly increased in GP-treated GCs in a dose-dependent manner. Kuittinen et al. [56] reported that high MMP expression was found in lymphoid and malignant transformed cells. In summary, our results suggested that the upregulation of MMP3 in GCs is the result of GP cytotoxicity rather than the maturation of the cells.
Survivin (BIRC5) is an inhibitor of apoptosis and has been shown to play an important role in cellular apoptosis by blocking the downstream step of mitochondrial cytochrome c release [57]. Several studies have shown that the BIRC5 protein is an inhibitor of apoptosis and is overexpressed in human cancer cells [58,59]. Jiang et al. [44] reported that the toxicity of GP was mediated by the expression of genes such as BIRC5 involved in cell cycle and survival, thereby inhibiting the proliferation of prostate cancer cells. Our mRNA expression results showed that BIRC5 was downregulated by GP cytotoxicity, indicating that GP induces apoptosis by inhibiting cell proliferation of porcine GCs.
The catalytic reaction of cytochrome P450 proteins (CYPs) is not essential for the maintenance of cellular activity, the CYP family are called phase I enzymes; they monooxygenase, reduce, and hydrolyze various substances such as lipids, steroidal hormones, and xenobiotics [61]. Phase I enzymes are expressed primarily in the liver, but they are also present in other tissues such as uterus, placenta, brain, kidney, and testis [62]. Previous studies have reported that certain CYP isoforms (CYP1A1, CYP1A2 and CYP2B) are present in porcine prepubertal ovary cells [63,64]. The Ah receptor (AhR) is a transcription factor that mediates the toxic effects of chemicals. Pocar et al. [65] reported that treatment with TCDD (2,3,7,8 tetrachlorodibenzo-p-dioxin) and PCB126 (polychlorinated biphenyl 126) induced the AhR signal transduction response via the expression of CYP1A1, an AhR target gene in porcine thyrocyte. In addition, Barć et al. [60] reported that the activity and expression of CYP1A1 dosedependently increased in polychlorinated naphthalene-treated porcine GCs. According to the qRT-PCR results in the present study, CYP1A1 mRNA was significantly increased in GP-treated GCs, which is likely due to GP cytotoxicity. However, RNA-seq results did not show a significant difference in CYP1A1 mRNA expression with GP treatment, although an upregulation pattern was observed.
The transforming growth factor β (TGFβ) gene is involved in several important cellular functions, including proliferation and differentiation. It is also a molecular marker related to embryonic growth and development of mammalian oocytes [66]. Dragovic et al. [67] reported that TGFβ gene is an important marker of cumulus expansion, and this process has been shown to be determined by the TGFβ protein signaling pathway. COL1A2 is involved in TGFβ signaling, which plays an important role in mediating ovarian age and oocyte quality [68,69]. The COL1A2 gene is involved in cellular proliferation and regulation of the cell cycle, and decreased levels of COL1A2 in old mares indicate the potential deregulation of TGFβ signaling as an underlying factor in the agerelated decline in oocyte quality [70]. In the present study, we found that the mRNA expression of the COL1A2 and TGFβ3 genes were significantly downregulated. In view of the function of these two genes, these results suggest that GP inhibits the proliferation and differentiation of porcine GCs and thus could affect the growth and development of oocyte.
To evaluate the effect of GP toxicity on intracellular oxidative damage in porcine GCs, the relative mRNA and protein abundance of three genes related to ROS production was analyzed. Peroxiredoxins (PRDXs) are peroxidases that work together to detoxify ROS and provide cytoprotection from internal and external environmental stress [71,72]. Kubo et al. [73] showed that PRDXs are down-regulators of various apoptotic pathways. In particular, since cumulus-oocyte complexes form a functional entity, the upregulation of PRDX6 in oocytes at the protein level and the accumulation of deadenylated transcripts have been shown to play an important role not only in embryo development but also oocyte maturation [74]. Microsomal glutathione transferase 1 (MGST1) is a member of the membrane-associated proteins in eicosanoid and glutathione metabolism (MAPEG) superfamily [75]. MGST1 is a membrane protein located in the endoplasmic reticulum and the outer mitochondrial membrane that protects cells from oxidative stress [76]. Several studies have reported that MGST1 is activated by a variety of factors, including reactive oxygen (O2¯, H2O2) and nitrogen (ONOO¯, NO) species, sulfhydryl reagents, proteolysis, heating, radiation and dimerization of the homotrimers [77,78]. In addition, MGST1 has been shown to play an important role in the detoxication of drugs and xenobiotics and may protect membranes against lipid peroxidation [79,80]. SOD3 is an important antioxidant enzyme located in the ECM of numerous tissues and the glycocalyx of cell surfaces. SOD3 removes superoxide and ROS to prevent cell death and protect normal tissues [81,82]. In addition, SOD3 may play a critical role in the protection against external environmental factors such as microbiological pathogens [83]. Basina et al. [39] reported that GP stimulates O2 generation by SOD activity inhibition in porcine GCs. The authors proposed that the increased O2 levels play a role in reducing steroid production in GCs. In this study, we found that SOD3 gene was downregulated according to GP concentration. This result is consistent with the SOD activity reported by Basini et al. [39]. Furthermore, in our study, the MGST1 and PRDX6 genes were found to be upregulated at both the mRNA and protein levels.
Collectively, these results indicate that these genes may be part of the cellular defense against oxidative stress caused by GP toxicity in the differentiation process of GCs.
In conclusion, we used RNA-seq to examine the gene expression profiles of GP-treated porcine GCs. To our knowledge, this is the first time that the effects of GP cytotoxicity in porcine GCs were studied using RNA-seq.
We selected a total of 11 genes and classified their functions using GO and KEGG analyses of the DEGs obtained with GP treatment. qRT-PCR and Western blot analysis showed significant differences in the expression of selected genes. However, further studies are warranted to identify the mechanisms underlying the functions of each gene. The results of this study provide a basis for examining the effects of GP toxicity in porcine. In addition, our findings provide greater insight into the various aspects of GP cytotoxicity in GCs.

Supplementary information
Additional file 1 Supplemental Table 1 Primer sequence for qRT-PCR Supplemental Figure 1 Uncropped Western blot images of Figure 6.