Analysis of the miRNA Expression Profiles in the Zearalenone-Exposed TM3 Leydig Cell Line

Zearalenone (ZEN), an important environmental pollutant, can cause serious harm to human and animal health. The aim of our study was to examine the effect of zearalenone (ZEN) on miRNA expression profiles in the mouse Leydig cell line (TM3 Leydig cell line) by miRNA sequencing. The effect of ZEN on the viability of TM3 Leydig cells was verified by Cell Counting Kit-8 (CCK-8). MiRNA sequencing was performed 24 h after the exposure of TM3 Leydig cells with 50 μmol/L of ZEN. Bioinformatics predicted the miRNA target genes, performed Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses, and conducted miRNA-gene-pathway mapping to show the relationship between miRNA, the target gene, and the signalling pathway. The expression levels of miRNA and the miRNA target genes associated with ZEN toxicology were verified by quantitative real-time polymerase chain reaction. The miRNA sequencing revealed a significant change (p < 0.05) in the 197 miRNAs in the ZEN-treated and control groups, among which 86 were up-regulated and 111 were down-regulated. GO analysis of the target genes of these miRNAs indicated various biological functions. KEGG analysis showed that the predicted miRNA target genes were involved in signalling pathways, such as cancer, apoptosis, and oxidation, namely, the Ras signalling pathway, Rap1 signalling pathway, PI3K-AKT signalling pathway, Foxo signalling pathway, and AMPK signalling pathway. These results suggest that ZEN, as an estrogen-like toxin, is regulated by microRNAs. Our results can help to examine the toxicological effects of ZEN-regulated miRNAs on germ cells.


Introduction
Zearalenone (ZEN) is a mycotoxin produced by Fusarium fungi [1,2]. The structure of ZEN is similar to that of 17β-oestradiol: ZEN competitively binds to estrogen receptors and activates the transcription of estrogen-responsive genes [3,4]. Therefore, ZEN plays a role by interfering with the physiological estrogen signalling pathway. ZEN may cause reproductive problems, such as ovarian dysfunction, decreased fertility, early abortion, reduced litter size, lower testicular weight, decreased motility of spermatozoa, and a lower total motile sperm count [5][6][7]. These reproductive toxicities

Differential Expression of miRNAs in the ZEN Exposure Groups
The miRNA deep sequencing analysis was used to identify the differential expression of miRNAs in the TM3 Leydig cells between the ZEN-treated and control groups. In the miRNA sequencing results, 197 miRNAs changed significantly (fold change > 1.5, p < 0.05), with 86 miRNAs up-regulated and 111 miRNAs down-regulated ( Figure 1B).

Prediction and Functional Classification of Target Genes
For the target animal species, Target Gene and MiRanda software were used to perform target gene predictions for miRNAs with significant differences. The predicted genes were classified into 3 Gene ontology (GO) categories and 50 terms were enriched from gene ontology analysis ( Figure 2). Twenty-five terms were enriched to the biological process category, including "biological process" and "regulation of transcription". Fifteen terms were enriched to the cellular component category, including "cytoplasm", "nucleus", and "integral components of membrane". Ten terms were enriched to the molecular function category, including "protein binding", "molecular function", and "metal ion binding".

Prediction and Functional Classification of Target Genes
For the target animal species, Target Gene and MiRanda software were used to perform target gene predictions for miRNAs with significant differences. The predicted genes were classified into 3 Gene ontology (GO) categories and 50 terms were enriched from gene ontology analysis ( Figure 2). Twenty-five terms were enriched to the biological process category, including "biological process" and "regulation of transcription". Fifteen terms were enriched to the cellular component category, including "cytoplasm", "nucleus", and "integral components of membrane". Ten terms were enriched to the molecular function category, including "protein binding", "molecular function", and "metal ion binding". Pathway significance enrichment analysis revealed that the ZEN-affected miRNA target genes have a close relationship with signal pathways, such as the Ras signalling pathway, the pathway in cancer, Foxo signalling pathway, and endocytosis. The number of miRNA target genes in these pathways is significant. These pathways have a regulatory effect on germ cell apoptosis, autophagy, oxidative stress, cancer development, invasion, and differentiation. The above signal pathways were the focus of the research on the effects of ZEN on germ cells (Figure 3). Some miRNAs likely to be associated with ZEN toxicology were screened by GO and Kyoto Encyclopedia of Genes and . Gene ontology (GO) term enrichment of the predicted targets of differentially expressed miRNAs. GO term enrichment showed that the putative targets of the differentially expressed miRNAs were associated with diverse functional terms, including binding, membrane, and reproduction.
Pathway significance enrichment analysis revealed that the ZEN-affected miRNA target genes have a close relationship with signal pathways, such as the Ras signalling pathway, the pathway in cancer, Foxo signalling pathway, and endocytosis. The number of miRNA target genes in these pathways is significant. These pathways have a regulatory effect on germ cell apoptosis, autophagy, oxidative stress, cancer development, invasion, and differentiation. The above signal pathways were the focus of the research on the effects of ZEN on germ cells ( Figure 3). Some miRNAs likely to be associated with ZEN toxicology were screened by GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. The results are shown in Table 1. Note: a log2 (fold change) refers to the relative expression levels between ZEN-exposed to control groups generated by sequencing data. Positive Expressions Up and Negative Expressions Down.
We mapped the miRNA-gene-pathway interaction map based on the GO and KEGG analyses ( Figure 4). One miRNA can regulate multiple target genes, and one target gene is regulated by multiple miRNAs. The target genes regulated by miRNAs regulate the corresponding pathways. Through the miRNA-gene-pathway interaction map, the relationship between the miRNA, target genes, and signalling pathways can be clearly seen. However, these signalling pathways are closely related to the toxicological mechanisms of ZEN. KEGG pathway analysis showed that the predicted targets were involved in cancer, apoptosis, and signalling pathways, such as the rat sarcoma signalling pathway (Ras signalling pathway), Ras-related protein1 signalling pathway (Rap1 signalling pathway), Phosphatidylinositide 3-kinases/AKT signalling pathway (PI3K-AKT signalling pathway), Forkhead boxo signalling pathway (Foxo signalling pathway), Adenosine 5'-monophosphate (AMP)-activated protein kinase (AMPK signalling pathway), and mitogen-activated protein kinase signalling pathway (MAPK signalling pathway). . Pathway enrichment of the predicted targets of the differentially expressed miRNAs. KEGG pathway analysis showed that the predicted targets were involved in cancer, apoptosis, and signalling pathways, such as the rat sarcoma signalling pathway (Ras signalling pathway), Ras-related protein1 signalling pathway (Rap1 signalling pathway), Phosphatidylinositide 3-kinases/AKT signalling pathway (PI3K-AKT signalling pathway), Forkhead boxo signalling pathway (Foxo signalling pathway), Adenosine 5'-monophosphate (AMP)-activated protein kinase (AMPK signalling pathway), and mitogen-activated protein kinase signalling pathway (MAPK signalling pathway).  Note: a log2 (fold change) refers to the relative expression levels between ZEN-exposed to control groups generated by sequencing data. Positive Expressions Up and Negative Expressions Down.
We mapped the miRNA-gene-pathway interaction map based on the GO and KEGG analyses ( Figure 4). One miRNA can regulate multiple target genes, and one target gene is regulated by multiple miRNAs. The target genes regulated by miRNAs regulate the corresponding pathways. Through the miRNA-gene-pathway interaction map, the relationship between the miRNA, target genes, and signalling pathways can be clearly seen. However, these signalling pathways are closely related to the toxicological mechanisms of ZEN.

Verification of the miRNA Target Genes
The gene expression levels represent the mRNA expression levels relative to the control levels (the values represent mean ± SD). The asterisks are used to indicate a statistically significant difference: * p < 0.05; ** p < 0.01.

Verification of the miRNA Target Genes
As shown in Figure 5B, we verified the expression levels of the miRNA target genes in the Ras signalling pathway, Rap1 signalling pathway, PI3K-AKT signalling pathway, Foxo signalling pathway, AMPK signalling pathway, and MAPK signalling pathway by qRT-PCR. The findings are consistent with the sequencing results.

Discussion
The mechanism of ZEN toxicity to humans and animals is complicated. Elucidating the mechanism of ZEN toxicity to humans and animals is of great clinical significance to prevent and treat ZEN. Recent studies have shown that the addition of reproductive hormones to TM3 cells affects the expression of miRNAs in TM3 Leydig cells [17]. Therefore, we speculate that miRNAs are involved in the regulation of ZEN toxicology after ZEN is exposed to TM3 cells.
GO enrichment analyses were performed for the miRNA target genes. The results showed that a lot of GO terms were enriched. Specifically, the GO term "cytoplasm", "nucleus", "integral component of membrane", and "protein binding" were the most significantly enriched. Coincidently, ZEN can affects the cellular structure and cellular connections of germ cells [18]. Therefore, we hypothesized a correlation between the toxicological effects of ZEN and these GO terms.
We selected some signalling pathway related to ZEN toxicology through GO and KEGG analyses. The Ras signalling pathway, Rap1 signalling pathway, PI3K-AKT signalling pathway, Foxo signalling pathway, AMPK signalling pathway, MAPK signalling pathway, and other signalling pathways were used to map the miRNA-gene-pathway interaction. The miRNA-gene-pathway interaction map clearly showed the mutual regulation among these six signalling pathways. As indicated by the KEGG analysis, miRNA could be involved in the regulation of target genes in these six signalling pathways by ZEN and affect the cell phenotype of germ cells.
Ras oncoprotein plays a key role in the development and maintenance of many tumor types [19,20]. Rap1 and Ras have a high degree of sequence similarity, have overlapping binding partners and have been shown to antagonize and mimic the Ras-driven cancer phenotype [21]. Studies have shown that ZEN, as a natural female-like toxin, mainly produces carcinogenic effects by interfering with endocrine balance [22]. In our results, the results of the KEGG analysis were significant in the Ras and Rap1 signalling pathways; the expressions of H-ras, K-ras, Rap1a, and Rap1b genes in the Ras and Rap1 proteins were significantly regulated, and H-ras, K-ras, Rap1a, and Rap1b were regulated. The expression of the isogenic miRNA (Figures 4 and 5) was also significant. This finding suggests that miRNAs may be involved in the toxicological process of ZEN carcinogenesis through the Ras and Rap1 signalling pathways. In addition, some scholars reported that ZEN influences the tumorigenic mechanism of TM3 Leydig cells by affecting proto-oncogenes [23]. Our research provides a new idea for the tumorigenic mechanism of ZEN.
The PI3K-AKT signalling pathway is an anti-apoptotic pathway that regulates the expression of downstream target proteins, such as Bax and Bcl-2, and thus participates in the regulation of cell growth, apoptosis, differentiation, migration, invasion, and angiogenesis [24,25]. The PTEN protein upstream of the PI3K-AKT signalling pathway can dephosphorylate phosphatidylinositol triphosphate to phosphatidylinositol diphosphate, thereby negatively regulating the PI3K-AKT signalling pathway and playing an important role in apoptosis [26,27]. The Foxo signalling pathway is located downstream of the PI3K-AKT signalling pathway and is regulated by the PI3K-AKT signalling pathway. The Foxo signalling pathway is closely related to apoptosis and oxidative stress [28,29]. A large number of studies have shown that ZEN can produce toxicological effects, such as oxidative stress and apoptosis, on germ cells. Long et al. investigated the oxidative damage of testis induced by ZEN in male mice and the apoptosis of ZEN induced by the Nrf2/ARE signalling pathway [11,30]. ZEN is known to cause changes in cell phenotypes, such as oxidative stress and apoptosis, in germ cells through the signalling pathways. In our experiments, genes such as PTEN, AKT2, AKT3, Foxo1, Foxo3, and Foxo6 were significantly expressed in the PI3K-AKT signalling pathway and the Foxo signalling pathway. The expression of miRNAs (Figures 4 and 5) regulating PTEN, AKT2, AKT3, Foxo1, Foxo3, and other genes was significant. This finding suggests that ZEN may cause germ cell apoptosis through the PI3K-AKT signalling pathway and the Foxo signalling pathway, and that miRNA may participate in the regulation of germ cell apoptosis through ZEN by affecting the target genes in the PI3K-AKT signalling pathway and the Foxo signalling pathway. Although no studies have shown that ZEN induces germ cell apoptosis through the PI3K-AKT signalling pathway and the Foxo signalling pathway, our study provides strong evidence of ZEN being able to induce germ cell apoptosis through the PI3K-AKT signalling pathway and the Foxo signalling pathway.
In the AMPK signalling pathway, the expression of the Prkaa2 genes that regulate the AMPK protein is significant, and the differential expression of the miRNAs that regulate the Prkaa2 genes (Figures 4 and 5) is also significant. Currently, studies have shown that in addition to the classical pathway of endoplasmic reticulum-induced apoptosis, the ATP/AMPK signalling pathway is regulated by endoplasmic reticulum stress during apoptosis induced by ZEN [31]. Therefore, through our experimental results, we boldly speculate that miRNA may be involved in the regulation. In addition, Pistol et al. [32] verified by transcriptome sequencing that the ZEN inflammatory stimuli and immunotoxicity in pig spleen cells could be the result of the JNK pathway activation but not that of the p-38/MAPK and NF-kB genes and proteins. On the basis of their results, we validated the JNK and JUN genes Mapk-8 and Mapk-10 as well as the JUN genes and their miRNAs (Figures 4 and 5) in the MAPK signalling pathway, and the results were all significant. This outcome suggests that the inflammatory stimuli and immunotoxicity produced by the ZEN activation of the JNK pathway may involve miRNAs. Our research provides a theoretical basis for this assumption.
Our results of this study reveal a complex network of miRNA associated with diverse molecular functions, which may be engaged in cellular response to ZEA exposure. However, more detailed aspects of these findings should be elucidated in future research, such as luciferase assay or gene knockout experiment, for better characterization of genes selected in this study.

Cell Culture and ZEN Exposure
Mouse TM3 Leydig cells were obtained from the American Type Culture Collection (Beijing, China). The cells were seeded as monolayer cultures in Dulbecco's modified Eagle 's medium/F-12 (HyClone, Logan, UT, USA) with 10% (v/v) inactivated fetal bovine serum and 1% (v/v) penicillin/streptomycin incubated at 37 • C with 5% CO 2 . When the volume of adherent cells occupied 80% of the dish, the cells were treated with 50 µmol/L of ZEN for 24 h. Each treatment was replicated thrice.

Cell Viability Assay
The cytotoxic effects of ZEN on TM3 cells were determined using a Cell Counting Kit-8 (CCK-8) (Solarbio, Beijing, China) assay. Cells were plated at a density of 2 × 10 5 per well in 96-well plates. After treatment with 0-90 µmol/L ZEN for 24 h, 10 µL of CCK-8 was added to each well, and the cells were incubated for 2 h at 37 • C. Non-treated cells served as a negative control. The absorbance was determined at a wavelength of 450 nm using a microplate reader (Infinite 200 PRO, ABI, New York, NY, USA). The results are presented as percentage of the values measured for untreated control cells. It has been reported that the IC 50 of ZEN-treated TM3 cells is 50 µmol/L [33].

Small RNA Sequencing and Bioinformatics Analysis
Total RNA was extracted using Trizol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's procedure. The total RNA quantity and purity were analyzed by Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, Palo Alto, CA, USA) with RIN number greater than 7.0. Approximately 1 µg of the total RNA was used to prepare a small RNA library according to the protocol of TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, CA, USA). We performed single-end sequencing (36 bp) on an Illumina Hiseq2500 in the Lianchuan Biotechnology (Hangzhou, China) following the vendor's recommended protocol.
Raw reads were subjected to an in-house program, ACGT101-miR (LC Sciences, Houston, TX, USA) to remove adapter dimers, junk, low complexity, common RNA families (rRNA, tRNA, snRNA, snoRNA), and repeats. Subsequently, unique sequences with length in 18~26 nucleotide were mapped to specific species precursors in miRBase 21.0 (http://www.mirbase.org/) by BLAST search to identify known miRNAs and novel 3p-and 5p-derived miRNAs. Length variation at both 3 and 5 ends and one mismatch inside of the sequence were allowed in the alignment. The unique sequences mapping to specific species mature miRNAs in hairpin arms were identified as known miRNAs. The unique sequences mapping to the other arm of the known specific species precursor hairpin opposite to the annotated mature miRNA-containing arm were considered to be novel 5p-or 3p-derived miRNA candidates. The remaining sequences were mapped to other selected species precursors (with the exclusion of specific species) in miRBase 21.0 by BLAST search, and the mapped pre-miRNAs were further BLASTed against the specific species genomes to determine their genomic locations. The above two were defined as known miRNAs. The unmapped sequences were BLASTed against the specific genomes, and the hairpin RNA structures containing sequences were predicated from the flank 80 nt sequences using RNAfold software (http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi).
To predict the genes targeted by the differentially expressed miRNAs, two computational target prediction algorithms (TargetScan 50 and miRanda 3.3a) were used to identify the miRNA binding sites. The data predicted by both algorithms were then combined, and the overlaps were calculated. The GO terms and KEGG pathway of these differentially expressed miRNA targets were also annotated. Then, we mapped the predicted target genes to map the miRNA-gene-pathway interactions.

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)
Total RNA was extracted using Trizol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's procedure. The total RNAs from each sample were reverse-transcribed to cDNA using the miRNA first-strand cDNA synthesis kit (by stem-loop) (Vazyme, Nanjing, China). The qRT-PCR was performed using the miRNA Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China). Reverse transcription reactions were performed at the following parameters: 25 • C (mRNA) for 5 min, 50 • C for 15 min, and 85 • C for 5 min. PCR reactions were performed at the following parameters: 95 • C for 5 min followed by 40 cycles of 95 • C for 10 s and 60 • C for 30 s. A U6 small nuclear RNA was used as the endogenous control for the data normalization of miRNAs, and β-actin was then adopted as the internal control of mRNA expression. The relative expression was calculated by the comparative threshold cycle method. The sequences of the primers used for reverse transcription and qRT-PCR were purchased from Sangon Biotech Company, Shanghai, China (Tables 2 and 3).

Conclusions
In summary, differentially expressed miRNAs in TM3 Leydig cells after ZEN exposure were identified. Sixteen miRNAs and target genes were identified to exhibit differential expression. GO enrichment analysis and pathway interaction analysis showed that miRNA and target genes are closely related to the toxicological mechanism of ZEN. QRT-PCR analysis of representative genes indicated that Ras signalling pathway, Rap1 signalling pathway, PI3K-AKT signalling pathway, Foxo signalling pathway, AMPK signalling pathway, and MAPK signalling pathway lead to oxidative stress, apoptosis, and carcinogenesis of germ cells. This study provides a molecular basis and new insights into the toxicological mechanisms of ZEN.