Deoxynivalenol- and zearalenone-contaminated feeds alter gene expression profiles in the livers of piglets.

OBJECTIVE: The Fusarium mycotoxins of deoxynivalenol (DON) and zerolenone (ZEN) cause health hazards for both humans and farm animals. Therefore, the main intention of this study was to reveal DON and ZEN effects on the mRNA expression of pro-inflammatory cytokines and other immune related genes in the liver of piglets.
METHODS: In the present study, 15 six-week-old piglets were randomly assigned to the following three different dietary treatments for 4 weeks: control diet, diet containing 8 mg DON/kg feed, and diet containing 0.8 mg ZEN/kg feed. After 4 weeks, liver samples were collected and sequenced using RNA-Seq to investigate the effects of the mycotoxins on genes and gene networks associated with the immune systems of the piglets.
RESULTS: Our analysis identified a total of 249 differentially expressed genes (DEGs), which included 99 upregulated and 150 downregulated genes in both the DON and ZEN dietary treatment groups. After biological pathway analysis, the DEGs were determined to be significantly enriched in gene ontology terms associated with many biological pathways, including immune response and cellular and metabolic processes. Consistent with inflammatory stimulation due to the mycotoxin-contaminated diet, the following Kyoto encyclopedia of genes and genomes pathways, which were related to disease and immune responses, were found to be enriched in the DEGs: allograft rejection pathway, cell adhesion molecules, graft-versus-host disease, autoimmune thyroid disease (AITD), type I diabetes mellitus, human T-cell leukemia lymphoma virus infection, and viral carcinogenesis. Genome-wide expression analysis revealed that DON and ZEN treatments downregulated the expression of the majority of the DEGs that were associated with inflammatory cytokines (interleukin 10 receptor, beta, chemokine [C-X-C motif] ligand 9), proliferation (insulin-like growth factor 1, major facilitator superfamily domain containing 2A, insulin-like growth factor binding protein 2, lipase G, and salt inducible kinase 1), and other immune response networks (paired immunoglobulin-like type 2 receptor beta, Src-like-adaptor-1 [SLA1], SLA3, SLA5, SLA7, claudin 4, nicotinamide N-methyltransferase, thyrotropin-releasing hormone degrading enzyme, ubiquitin D, histone H2B type 1, and serum amyloid A).
CONCLUSION: In summary, our results demonstrated that high concentrations DON and ZEN disrupt immune-related processes in the liver.

INTRODUCTION

Fusarium mycotoxins are secondary metabolites produced by various moulds that are commonly found in foods, predominantly cereals and forages [1]. Mycotoxins cause major economic losses due to crop infection and frequently contaminate food and animal feeds, thereby posing health hazards to humans and livestock. The toxic effects of Fusarium mycotoxins in animals include feed rejection, reduced growth, immunosuppression, gastrointestinal lesions, and reproductive and neurological disorders [2].

Deoxynivalenol (DON) and zerolenone (ZEN) are the most frequently detected mycotoxins among the Fusarium mycooxins. Therefore, it is interesting to determine their toxic effects when present in the animal feed. In particular, wheat, maize, and triticale grains are susceptible to Fusarium infection and have also been reported to be more frequently contaminated with DON and ZEN than other cereal grains. Among farm animals, the pig is relatively more sensitive to higher mycotoxin concentrations. In particular, DON causes reduced voluntary feed intake, while ZEN is associated with fertility disorders and hyperestrogenism in pigs. On the other hand, ruminants and poultry were observed to be less sensitive to mycotoxins [3].

DON exerts toxic and immunotoxic effects in various animal species. However, compared to other species, pigs are more sensitive to DON, partly because of differences in DON metabolism and also due to the low de-epoxidation activity of gut microbiota. DON acts as a protein synthesis inhibitor and thus reduces claudin and occluding expression in epithelial cells, which in turn results in the deterioration of the intestinal barrier and increases permeability to bacteria [4]. DON exerts many negative effects, especially on cytokines, and also causes anorexia because of changes in immune defense regulation, oxidative status, barrier functions, and various regulatory mechanisms [5]. Chronic low dietary concentrations of DON do not alter animal performance, hematological and immunological variables, and biochemical parameters. However, higher acute doses can cause vomiting, diarrhea, hemorrhagic, and circulatory shock [6].

ZEN severely disrupts estrogenic mechanisms due to its structural similarities with estrogen. ZEN is known to cause estrogenic effects on several animal species, especially pigs. The most common pathological effects of ZEN include hyperestrogenism, reduced fertility, abortion, and anoestrus, as well as higher rates of embryonic and fetal death [7]. In addition, ZEN toxicity is linked to decreased litter size, smaller pituitary and thyroid glands in offspring, and altered serum levels of progesterone and estradiol. Several previous studies have reported the immunotoxic effects of ZEN in swine, particularly in inflammatory processes. ZEN and its metabolites were confirmed to exert different effects on innate immunity-related processes in the pig and can act as inducers or suppressors of inflammatory cytokine expression in peripheral blood cells [8].

Several groups have demonstrated the adverse effects of DON or ZEN on different animals, which result from complex, multi-pathway mechanisms that are regulated at various developmental stages. Numerous studies have identified the molecules involved in mediating the signal transduction and toxicity pathways of these toxins and have investigated both the complexities and commonalities in their molecular pathways in animal cells. For example, Brzuzan et al [9] demonstrated that exposure of immature gilts to DON and ZEN contaminated diet shows the pathways which are involved in cell proliferation and survival were disordered. Based on previous research reports on pigs; in the present study we are using acute toxic levels of 8 mg/kg and 0.8 mg/kg, DON and ZEN, respectively, for 4 weeks treatment. In the present study, RNA transcriptome analysis was conducted on liver samples of piglets fed with DON- and ZEN-contaminated feeds. RNA-Seq analysis is especially useful as a screening technique for studying biological systems with little known information. Gene expression profiling via RNA-Seq should help gain insights into the molecular effects of mycotoxins and identify metabolism-related genes that can be targeted for DON and ZEN management. Putative marker genes and pathways influenced by DON and ZEN should be identified to obtain a more detailed understanding of the mechanisms by which chronic levels of DON and ZEN in contaminated feeds affect the livers of pigs. In addition, the biologically relevant concentrations at which DON and ZEN influence gene expression in the liver should be evaluated.

MATERIALS AND METHODS

Ethics statement

The protocols used for the animal experimental procedures were reviewed and approved by the Institutional Animal Care and Use Committee of the National Institute of Animal Science (No. 2015-147).

Animal exposure to DON and ZEN and experimental design

The present study was conducted using 15 male castrated 6-week-old piglets (~19 kg), which were purchased from a commercial pig farm. Each piglet was housed in a separate pen (2,100 cm×1,400 cm). Piglets were allowed to acclimate to their new housing conditions for one week at 25°C±1°C and were subsequently assigned to one of three treatment groups (5 piglets/group). Piglets in the three treatment groups had approximately equal body weights. Piglets in the control group were provided with a standard diet to meet the nutritional requirements for piglets [10], and those in the treatment groups were provided with the standard diet added with DON or ZEN. In particular, commercial DON and ZEN were properly mixed into the diet at 8 mg/kg and 0.8 mg/kg, respectively, added as purified toxins (Biomin Pte. Ltd. Singapore). Throughout the 4-week experimental period, the control and DON- and ZEN-contaminated diets and water were provided to the pigs ad libitum. Ethical guidelines for animal protection rights were observed.

Mycotoxin analysis

Quantitative determination of DON and ZEN were analyzed in DON and ZEN mixed corn feeds by using ultra performance liquid chromatography (UPLC). A homogenized DON mixed grain sample (1 g) was extracted with 20 ml of distilled water and shake for 30 min, and for ZEN mixed corn sample (1 g), add 0.5 g of NaCl and 20 mL acetonitrile (ACN) then shake for 1 hour. After filtering the extract through a Whatman paper (No. 1), 5 mL of DON mixed filtrate sample was diluted in 20 mL of phosphate-buffered saline (PBS), and 5 mL of ZEN mixed grain filtrate sample diluted in 20 mL of 1% tween 20 solution. DON and ZEN mixed extract samples were loaded separately onto the appropriate immunoaffinity chromatography columns. The DON mixed sample loaded columns were allowed for completely dried and washed with 10 mL of PBS and distilled water in the order, and then eluted with 0.5 mL of MeOH and 1.5 mL of ACN. In the case of ZEN, the column washed with 10 mL of distilled water, then eluted with 1.5 mL of MeOH. Elutes were dried under N2 gas and dissolved in 1 mL of 50% MeOH then filtered by using a 0.2 μm syringe filter, and to be injected into UPLC (Water Acquity UPLC H Class Milford, MA, USA). The photodiode array and fluorescence detector detector were used for DON and ZEN, respectively. The Waters Acquity UPLC BEH C18 column (2.1×100 mm, 1.7 μm particle size) was used for both analyses. We found that, the amounts of DON and ZEN in mixed corn feed were very close to the original concentrations (8 mg DON/kg feed and 0.8 mg ZEN/kg feed). There are no DON and ZEN contaminations were observed in the control feed sample.

Sample collection and processing

After 4 weeks of treatment, all control and DON- and ZEN-treated pigs were sacrified via anesthetic overdose using a mixture of barbiturates and pentobarbital. Immediately after cardiac attest, liver samples were collected, rapidly frozen in liquid nitrogen, and stored at −80°C for further analysis.

RNA extraction, mRNA library construction, and sequencing analysis

For cDNA library construction from liver samples, RNAs were isolated using the TruSeq RNA library kit and purified with Maxwell 16 LEV simply RNA Purification Kit (Promega, Madison, WI, USA) according to the manufacturer’s instructions. First, mRNAs were purified using poly-T oligo-attached magnetic beads. After purification, mRNAs were fragmented via incubation with divalent cations at room temperature. cDNAs were synthesized from the cleaved RNA fragments using reverse transcriptase and random primers. Next, second strand cDNA synthesis was performed using DNA Polymerase I and RNase H. cDNA fragments were end repaired by addition of a single ‘A’ base and were subseqently ligated with adapters. The resulting products were then purified and enriched by polymerase chain reaction (PCR) to generate the final cDNA library. An RNA-Seq library was constructed for each liver sample and was sequenced following the Illumina HiSeq2000 protocol to generate 90-bp paired-end reads. The quality of the RNA-Seq reads from all liver samples were assessed using FastQC. All reads passed the quality control critera and were mapped to Sus scrofa sequences obtained from UCSC using Tophat2 (v2.0.2). The number of mapped reads was counted using HTseq (v0.5.3p3).

Identification of differentially expressed genes from control and treatment groups

edgeR was used to identify the differentially expressed genes (DEGs) between control and DON and ZEN treatment groups [11]. edgeR was designed to analyze replicated count-based expression data and is based on a negative binomial model. Cox-Reid profile methods are incorporated in the edgeR package to estimate the dispersion for pairwise comparisons between control and DON, control and ZEN, and DON and ZEN treatment groups. The remaining Cox-Reid common volume dispersion and statistical analyses were followed as published in our earlier study [12].

Functional annotation of DON and ZEN DEGs

The Ensembl gene IDs of the DEGs were used to obtain the official gene symbols in the pig. The official gene symbols of the corresponding human homologs of pig DEGs were subjected to functional clustering and enrichment analysis using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) tool [13]. Expression patterns of functional groups in DON- and ZEN-treated samples relative to the whole genome were analyzed using the Expression Analysis Systematic Explorer (EASE) tool within DAVID, which implements a modified Fisher’s exact test to measure the enrichment of gene ontology (GO) terms. DAVID was used to assign Kyoto encyclopedia of genes and genomes (KEGG) pathways to the resulting DEGs. To identify enriched GO terms, functionally clustered genes were filtered using an EASE value cutoff of 0.01.

RNA isolation, reverse transcription, and quantitative real-time polymerase chain reation

To validate the results from differential expression analysis, RNAs were isolated from pig liver samples from the control, DON, and ZEN treatment groups using TRIzol reagent (Ambion, Carlsbad, CA, USA) following the manufacturer’s instructions. Total RNA (2 μg) was used as template for cDNA synthesis using GoTaq 2-step reverse transcription-quanititative polymerase chain reation (RT-qPCR) system (Promega, USA) kit with a final reaction volume of 20 μL. Quantitative RT-PCR (qRT-PCR) was performed to measure the mRNA expression levels of 12 DEGs in either the DON or ZEN treatment group. Next, significant DEGs were selected based on a fold change (FC) cutoff of logFC ≥ 2 in at least one of the DON and ZEN treatment groups. Primers specific to the selected genes of interest were designed based on available pig (Sus scrofa) sequences using Integrated DNA Technologies Oligo Analyzer Tool. Primer sequences, accession numbers, and amplified product lengths are shown in Supplementary Table S1. The qRT-PCR was carried out using a GoTaq qPCR master mix (Promega, USA) on an Applied Biosystems 7500 Fast Real-Time PCR System (Foster, CA, USA). Reactions were run according to the following amplification profile: 1 cycle of 95°C for 5 min, 40 cycles of 95°C for 15 s, 55°C for 30 s, and 72°C for 30 s, followed by a melting curve. Glyceraldehyde-3-phosphate dehydrogenase and β-actin were used as internal controls for normalization. Relative abundance between experimental and control groups was calculated using the 2−ΔΔCt method.

Statistical analysis

All sample data were expressed as mean±standard error of the mean. The analysis of variance, and t-test analyses were carried out to determine the statistical differences between control and mycotoxin-treated groups for all parameters analyzed. Further variations between means were analyzed using Fisher’s least significant difference method. Statistical significance was accepted at p≤0.05.

RESULTS

Identification of differentially expressed genes

To examine the volume dispersion and expression profile differences between pig liver tissues, we identified DEGs among the control and DON- and ZEN-treated samples (Supplementary Figure S1). We identified a total of 101 (37 upregulated and 64 downregulated), 78 (21 upregulated and 57 downregulated), and 70 (41 upregulated and 29 down regulated) DEGs after comparing DON and control, ZEN and control, and DON and ZEN treatment groups, respectively (false discovery rate [FDR] <0.01, FC≥2, and p≤0.05, Figure 1). We also generated a two-way hierarchical clustering dendrogram of gene expression (logFC>2) comparing the DON, ZEN, and control treatments, which clearly demonstrated that some genes were upregulated, but most genes were downregulated in DON and ZEN treatment groups relative to the control group (Supplementary Figure S2). Expression levels of DEGs from each pairwise comparison are shown in Supplementary Figure S3. We identified expression level, fold change, p-value, and annotation of all the 20,596 genes in the DON and ZEN treatment groups (data not shown). Considering the similar effects of DON and ZEN, RNA-Seq data from both DON and ZEN treatment groups and their corresponding DEGs were pooled and examined.

Figure 1

Expression counts of differentially expressed genes (DEGs) in control, deoxynivalenol (DON), and zerolenone (ZEN) treatment groups. The y-axis represents upregulation and downregulation of DEGs between control livers and mycotoxin-treated livers. The x-axis represents the total number of transcripts (DEGs for each pairwise comparison were considered at fold change [FC]>2 and p≤0.05).

Gene ontology and the functional annotation of DEGs

We performed GO classification of the DEGs obtained by comparing expression levels of samples in the DON and ZEN treatment groups relative to those of the control group. GO assignments under the biological process category are shown in Table 1. The most significantly enriched biological process GO terms in the DON- and ZEN-regulated DEGs included immune system process (p = 0.003), response to stimulus (p = 0.013), and cellular catabolic process (p = 0.005). The top three molecular function GO terms were binding (p = 0.050), transition metal ion binding (p = 0.020), and hormone activity (p = 0.007). The top two GO terms under the cellular component category were extracellular space (p = 0.029) and cellular region (p = 0.050).

Table 1

Top biological process GO terms associated with the DEGs identified by comparing the DON and ZEN mycotoxin treatment groups with the control group

Name	GOID1)	Term2)	Count3)	Genes4)	p-value
biological_process	GO:0002376	Immune system process	8	SLC7A2, CD207, UBD, LOC100622791, SLA-7, ORM1, IGFBP2, IGF1	0.003
biological_process	GO:1901575	Organic substance catabolic process	7	SDS, OAT, UBD, CCRN4L, LOC100525856, TRHDE, IGF1	0.005
biological_process	GO:0044248	Cellular catabolic process	7	SDS, OAT, UBD, CCRN4L, LOC100525856, TRHDE, IGF1	0.005
biological_process	GO:0044281	Small molecule metabolic process	7	CA2, SDS, OAT, GPT2, FAXDC2, IGF1, HMGCS2	0.011
biological_process	GO:0019752	Carboxylic acid metabolic process	5	SDS, OAT, GPT2, FAXDC2, IGF1	0.012
biological_process	GO:0009056	Catabolic process	7	SDS, OAT, UBD, CCRN4L, LOC100525856, TRHDE, IGF1.	0.013
biological_process	GO:0050896	Response to stimulus	17	SLC7A2, CD207, CA2, UBD, CCRN4L, LOC100515475, GAS1, LOC100525856, TRHDE, ARHGEF6, LOC100622791, SLA-7, GALP, IGFBP2, MT1A, IGF1, ADORA3	0.013
biological_process	GO:0023051	Regulation of signaling	7	CA2, UBD, GAS1, LOC100525856, ARHGEF6, IGFBP2, IGF1	0.049
biological_process	GO:0043436	Oxoacid metabolic process	5	SDS, OAT, GPT2, FAXDC2, IGF1	0.015
biological_process	GO:0006082	Organic acid metabolic process	5	SDS, OAT, GPT2, FAXDC2, IGF1	0.016
biological_process	GO:1901605	Alpha-amino acid metabolic process	3	SDS, OAT, GPT2	0.023
biological_process	GO:1901565	Organonitrogen compound catabolic process	3	SDS, OAT, TRHDE	0.024
biological_process	GO:0042104	Positive regulation of activated T cell proliferation	2	IGFBP2, IGF1	0.025
biological_process	GO:0045321	Leukocyte activation	4	SLC7A2, UBD, IGFBP2, IGF1	0.026
biological_process	GO:0043567	Regulation of insulin-like growth factor receptor signaling pathway	2	IGFBP2, IGF1	0.028
biological_process	GO:0001775	Cell activation	4	SLC7A2, UBD, IGFBP2, IGF1	0.041
biological_process	GO:0006520	Cellular amino acid metabolic process	3	SDS, OAT, GPT2	0.046

GO, gene ontology; DEGs, differentially expressed genes; DON, deoxynivalenol; ZEN, zerolenone.

Gene ontology ID of the biological process.

Name of the biological term.

Number of significant genes which are involved in biological terms.

Name of the genes which are involved in the particular biological term.

Immune- and disease-related DEGs

Our analysis identified interesting upregulated and downregulated DEGs that were related to inflammatory response and disease in both the DON and ZEN treatment groups. Some of the major upregulated and downregulated genes and their corresponding biological process assignments are listed in Tables 2, respectively. Some of the upregulated genes comprised genes with important functions in the immune response and proliferation, such as claudin 4 (CLDN4), orosomucoid 1 (ORM1), insulin-like growth factor binding protein 2 (IGFBP2), cytochrome P450, family 26, subfamily A, polypeptide 1 (CYP26A1), and growth arrest-specific 1 (GAS1). However, the majority of the DEGs that were associated with immune response, proliferation, and cytokine networks were downregulated, including Src-like-adaptor-1 (SLA-1), SLA-3, SLA-5, SLA-7, interleukin 10 receptor, beta (IL10RB), chemokine (C-X-C motif) ligand 9 (CXCL9), ubiquitin D (UBD), paired immunoglobulin-like type 2 receptor beta (PILRB), histone H2B type 1, sphingomyelin phosphodiesterase 3 (SMPD3), solute carrier family 7 member 2 (SLC7A2), thyrotropin-releasing hormone degrading enzyme (TRHDE), and serum amyloid A. We also identified certain diseases associated with the following downregulated DEGs: ornithine aminotransferase, nicotinamide N-methyltransferase (NNMT), insulin-like growth factor 1 (IGF-1), hydroxylmethylglutaryl-CoA synthase 2 (HMGCS2), GAS1, and major facilitator superfamily domain containing 2A.

Table 2

Upregulated genes in the livers of piglets exposed to both deoxynivalenol (DON) (8 mg/kg) and zerolenone (ZEN) (0.8 mg/kg) dietary treatments

Gene symbol	Gene description	Log2 FC	p-value	Biological process
CLDN4	Claudin 4	2.02	0.044	Regulation of immune system process, acute-phase response
CYP26A1	Cytochrome P450, family 26, subfamily A, polypeptide 1	2.13	0.012	Transport and catabolism, Retinol metabolism
IGFBP2	Insulin-like growth factor binding protein 2	2.75	0.028	Growth regulation, response to insulin
ORM1	Orosomucoid 1	2.25	0.028	Regulation of immune system process
FAXDC2	Fatty acidd hydroxylase domain containing 2	2.07	0.028	Fatty acid biosynthetic process
SMPD3	Sphingomyelin phosphodiesterase 3	2.21	0.033	Hematopoietic progenitor cell differentiation
OAZ3	Ornithine decarboxylase antizyme 3	3.46	0.017	Polyamine metabolic process
SCT	secretin	2.49	0.037	Hormone activity
UBD	Ubiquitin D	3.59	0.028	Positive regulation of apoptotic process
CA2	Carbonic anhydrase 2	2.10	0.048	Morphogenesis of an epithelium
GAS1	Growth arrest specific	2.33	0.047	Regulation of apoptotic process
GPT2	Glutamic-pyruvic transaminase 2	5.90	0.002	2-oxoglutarate metabolic process
PCTP	Phosphatidylcholine transfer protein	2.18	0.001	Intracellular vesicular traffic
IL10RB	Interleukin 10 receptor, beta	−2.10	0.031	Cytokine-cytokine receptor interaction
IGF1	Insulin like growth factor 1	−3.90	0.033	Positive regulation of cell proliferation, activation of MAPK activity
OAT	Ornithine aminotransferase	−2.02	0.047	Arginine catabolic process to proline
NOCT	Nocturnin	−2.41	0.043	Deadenylation-dependent decapping of nuclear-transcribed mRNA
NNMT	Nicotinamide N-methyl transferase	−4.09	0.032	Promotes epigenetic remodeling in cancer
LIPG	Lipase G, endothelial type	−2.27	0.043	Positive regulation of high-density lipoprotein particle clearance
SIK1	Salt inducible kinase 1	−2.24	0.011	Glucagon signaling pathway
TRHDE	Thyrotropin releasing hormone degrading enzyme	−2.17	0.036	Regulation of blood pressure, peptide catabolic process
SDS	Serine dehydratase	−2.52	0.012	Cellular amino acid metabolic process
SLC7A2	Solute carrier family 7 member 2	−2.10	0.017	Nitric oxide production involved in inflammatory response
SAA	Serum amyloid A	−3.44	0.007	Response to inflammatory stimuli
MT1A	Metallothionein 1A	−2.10	0.038	Negative regulation of growth
VSIG4	V-set and immunoglobulin domain containing 4	−2.10	0.006	Negative regulation of interleukin-2 production
ARHGEF6	Rac/Cdc42 guanine nucleotide exchange factor 6	−2.75	0.001	Regulation of Rho protein signal transduction
ADORA3	Adenosine A3 receptor	−2.14	0.009	Neuroactive ligand-receptor interaction
PKD1L3	Polycystic kidney disease 1-like 3	−2.45	0.018	Transient Receptor Potential Interacting
CD207	CD207 molecule, langerin	−5.60	0.034	Defense response to virus
PILRB	Paired immunoglobulin-like type 2 receptor beta	−2.16	0.029	Regulation of the immune system
SLA-1	Src-like-adaptor 1	−3.15	0.049	Immune response
SLA-3	Src-like-adaptor 3	−2.00	0.020	Immune response
SLA-5	Src-like-adaptor 5	−2.46	0.013	Immune response
SLA-7	Src-like-adaptor 7	−3.27	0.012	Antigen processing and presentation
HMGCS2	3-hydroxy-3-methylglutaryl-CoA synthase 2	−10.94	0.036	Isoprenoid biosynthetic process
GALP	Galanin like peptide	−14.31	0.012	Neuropeptide signaling pathway
UBD	Ubiquitin D	−3.53	0.048	Ubiquitin-dependent protein catabolic process
MFSD2A	Major facilitator superfamily domain containing 2A	−2.80	0.009	Establishment of blood-brain barrier
SMPD3	Sphingomyelin phosphodiesterase 3	−2.01	0.036	Hematopoietic progenitor cell differentiation
CGREF1	Cell growth regulator with EF-hand domain 1	−2.47	0.027	Calcium ion binding
CXCL9	Chemokine (C-X-C motif) ligand 9	−3.28	0.017	Inflammatory response
HIST1H2B	Histone H2B type 1	−2.10	0.012	Innate immune response in mucosa

KEGG pathway analysis of DEGs in the DON and ZEN treatment groups

GO enrichment analysis was carried out using the DAVID software and knowledgebase to capture enriched biological terms associated with the DEGs from the DON and ZEN treatment groups. Based on KEGG pathway analysis, DEGs from both mycotoxin treatment groups were mainly clustered under immune response, disease, and metabolism. As shown in Table 3, the majority of DEGs obtained from the DON and ZEN treatment groups relative to the control group were associated with immune, disease, and viral carcinogenesis pathways. The most significantly enriched pathways included metabolic pathways (p = 4.45E-07), immune associated pathways (allograft rejection pathway [2.34925E-06], graft-versus-host disease [GVHD] [p = 1.84987E-06], AITD [p = 5.53981E-06], cell adhesion molecules [CAMs] [p = 2.5123E-06], type I diabetes mellitus [3.14384E-06], and viral myocarditis [p = 1.02775E-05]), infectious disease pathways (herpes simplex infection [p = 6.0546E-06], Epstein-Barr virus infection [p = 0.00025], and human T-cell leukemia type 1 [HTLV-I] [p = 0.000059]), and the cancer pathway viral carcinogenesis (p = 0.00025) (Supplementary Figure S4).

Table 3

KEGG pathways associated with DEGs from the DON and ZEN dietary treatment groups

KEGG ID1)	KEGG map name2)	Count3)	Gene names4)	p-value	Bonferroni	FDR
01100	Metabolic pathways	9	CYP26A1, NNMT, SDS, LIPG, OAT, SMPD3, GPT2, FAXDC2, HMGCS2	4.44551E-07	0.00003.2896	3.28967E-05
05332	Graft-versus-host disease	4	SLA-3, SLA-1, SLA-5, SLA-7	1.84987E-06	0.000136891	4.64776E-05
05330	Allograft rejection	4	SLA-1, SLA-3, SLA-7, SLA-5	2.34925E-06	0.000173845	4.64776E-05
04514	Cell adhesion molecules (CAMs)	5	CLDN4, SLA-3, SLA-1, SLA-5, SLA-7	2.51230E-06	0.000185910	4.64776E-05
04940	Type I diabetes mellitus	4	SLA-5, SLA-7, SLA-3, SLA-1	3.14384E-06	0.000232644	4.65289E-05
05320	Autoimmune thyroid disease	4	SLA-1, SLA-3, SLA-7, SLA-5	5.53981E-06	0.000409946	6.40058E-05
05168	Herpes simplex infection	5	PILRA, SLA-3, SLA-1, SLA-5, SLA-7	6.05460E-06	0.000448041	6.40058E-05
05416	Viral myocarditis	4	SLA-3, SLA-7, SLA-5, SLA-1	1.02775E-05	0.000760534	9.50668E-05
04612	Antigen processing and presentation	4	SLA-7, SLA-5, SLA-3, SLA-1	1.28527E-05	0.000951102	1.05678E-04
04145	Phagosome	4	SLA-1, SLA-3, SLA-7, SLA-5	1.27643E-04	0.009445557	9.44556E-04
05203	Viral carcinogenesis	4	SLA-3, SLA-1, SLA-5, SLA-7	2.59025E-04	0.019167814	1.59731E-03
05169	Epstein-Barr virus infection	4	SLA-7, SLA-1, SLA-3, SLA-5	2.59025E-04	0.019167814	1.59731E-03
05166	HTLV-I infection	4	SLA-5, SLA-1, SLA-3, SLA-7	5.93257E-05	0.043900983	2960192E-03

KEGG, Kyoto encyclopedia of genes and genomes; DEGs, differentially expressed genes; DON, deoxynivalenol; ZEN, zerolenone; FDR, false discovery rate; HTLV-I, human T-cell leukemia lymphoma virus I.

KEGG pathway IDs in DON and ZEN treatment groups.

Name of the KEGG pathway which is identified in DON and ZEN treatment groups.

Number of significant genes which is involved in KEGG pathways.

Name of the genes which are involved in the particular KEGG pathways.

qRT-PCR validation of DEGs

The qRT-PCR was performed to validate the DEGs identified in the DON and ZEN treatment groups, and results are shown in Table 4. We randomly selected the following six upregulated and six downregulated genes from the putative DEGs identified from DON and ZEN treatment groups that were involved in disease or immunity: IGF-1, ORM-1, UBD, SLA-1, HMGCS2, metallothionein 1A (MT1A), nocturnin (NOCT), TRHDE, NNMT, Rho guanine nucleotide exchange factor 6 (ARHGEF6), claudin 4 (CALDN4), and paraoxonase 3 (PON3). Out of these 12 genes, ORM1, HMGCS2, NOCT, NNMT, CLDN4, and PON3 were confirmed to be differentially expressed in both DON- and ZEN-treated liver samples. The IGF-1 gene was only significantly differentially expressed in the DON-treated group, while the UBD and SLA-1 genes were validated only in the ZEN-treated group. There were no significant differences in the expression levels of MT1A, TRHDE, and ARHGEF6 in both the DON and ZEN treatment groups when compared relative to the control group. In summary, most of the DEGs identified based on RNA-Seq analysis were validated via qRT-PCR in both the DON and ZEN treatment groups (Table 4).

Table 4

Comparison of FCs estimated from RNA-Seq data from the DON and ZEN treatment groups and expression validation via quantitative real-time PCR

Gene name1)	RNA-sequencing FC2)				qRT-PCR FC2)
	DON/CON	p value	ZEN/CON	p value	DON/CON	p value	ZEN/CON	p value
IGF-1	−3.90	*	−2.43	NS	−4.32	*	−2.21	NS
ORM1	+2.24	*	+1.62	NS	2.91	*	2.02	*
UBD	+3.59	*	+1.01	*	3.83	NS	1. 32	*
SLA-1	+1.32	NS	−2.38	*	1.16	NS	−2.67	*
HMGCS2	+5.08	*	−2.15	*	3.94	*	−1.93	*
MT1A	−1.46	NS	−2.10	*	−1.73	NS	−2.79	NS
NOCT	−2.41	*	−2.29	*	−3.59	*	−1.93	*
TRHDE	−2.16	*	+1.04	NS	−3.02	NS	0.67	NS
NNMT	−4.08	*	−3.07	*	−3.41	*	−2.16	*
ARHGEF6	−1.30	NS	−2.7	*	−1.03	NS	−2.08	NS
CLDN4	+2.02	*	+1.22	NS	1.92	*	1.19	*
PON3	+3.00	*	+2.99	*	2.01	*	2.31	*

FCs, fold changes; DON, deoxynivalenol; ZEN, zerolenone; PCR, polymerase chain reaction; qRT-PCR, quantitative real-time polymerase chain reaction.

IGF-I, insulin-like growth factor I; ORM1, orosomucoid 1; UBD, ubiquitin D; SLA-1, Src-like-adaptor-1; HMGCS2, hydroxy-3-methylglutaryl-CoA synthase 2; MT1A, metallothionein 1A; NOCT, nocturnin; TRHDE, thyrotropin-releasing hormone degrading enzyme; NNMT, nicotinamdie N-methyltransferase; ARHGEF6, rho guanine nucleotide exchange factor 6; CLDN4, claudin 4; PON3, paraoxonase 3; GAPDH, glyceraldehydes 3-phosphate.

Fold change values expressed as a ratio between DON and CON groups, and ZEN and CON groups, were FC>1.00 indicates an increase in gene expression and FC<–1 indicates a decrease in gene expression n = 5 for each of the CON (0 mg/kg of DON), DON (8 mg/kg of DON) and ZEN (0.8 mg/kg of ZEN) treatments.

p≤0.05;

NS, no significant.

DISCUSSION

Among common farm animals, pigs are the most sensitive to the effects of the mycotoxins DON and ZEN [14]. In any RNA-Seq experiment, the appropriate bioinformatics pipeline used for differential gene expression analysis obviously depends on the final objective of the experiment. To date, little is known about the effects of DON and ZEN on the transcription profiles of pigs. So, in this study, we employed an RNA-Seq approach to investigate the effects of DON and ZEN on the global transcriptome expression profiles in the livers of pigs fed with DON- and ZEN-contaminated diets for 28 days. RNA-seq allowed reliable examination of a broad range of expression levels in DON- and ZEN-treated liver samples. In any RNA-Seq experiment, the appropriate bioinformatics pipeline used for differential gene expression analysis obviously depends on the final objective of the experiment.

RNA-Seq data from both DON- and ZEN-treated samples allowed the identification of genes with altered expression patterns through comparison of expression levels with the control group. Our analysis returned 58 upregulated and 121 downregulated genes. The majority of the DEGs were downregulated because of the predominantly toxic effects of DON and ZEN; these downregulated genes were primarily involved in cell signaling. Table 2 list some of the most significant upregulated and downregulated genes and their associated biological process GO terms. The effects of DON and ZEN treatment on the immune system, proliferation, and cytokine networks and the most significantly altered genes are discussed herein. RNA-Seq analysis revealed significant upregulation of CALDN4 (FC: 2.02), UBD (FC: 3.59), IFGBP2 (FC: 2.02), GAS1 (FC: 2.32), ORM1 (FC: 2.02), and SMPD3 (FC: 2.20). CLDN4 is known to promote pancreatic cancer and gastric carcinoma and has also been shown to cause impaired function of tight junctions, which are in turn related to tumor differentiation [15]. Treatment of pigs with 30 μmol/L DON resulted in the loss of epithelial barrier integrity of the intestine and decreased expression of the CLDN4 gene. Similar to our results, Alizadeh et al [16] data also found upregulation of CLDN4 in the cecum of DON fed (0.9 mg/kg) piglets. The ORM1 gene is responsible for the multifactorial and immunogenetic etiology of sarcoidosis. ORM exhibits anti-inflammatory and immunomodulatory properties, including anti-complement and anti-neutrophil functions [17]. UBD is involved in the immune response and plays a crucial role in the cell cycle by exerting pro-apoptotic (by interacting with HIV) and pro-proliferative (in liver and colon cancers) effects. UBD is overexpressed in various human gastrointestinal and gynecological cancers, as well as in more than 60% of hepatocellular carcinoma and 75% of liver cancer stem cells [18]. IGFBP2 is a member of the insulin-like growth factor-binding family of proteins and is secreted by most sarcoma cell lines. IGFBP2 also appears to be involved in the regulation of growth and transformation. The IGFBP2 levels measured at the end of cancer treatment are correlated with residual tumor load, suggesting that IGFBP2 can serve as an early indicator of therapeutic response [19]. GAS1 plays a role in growth repression and acts by blocking entry into the S phase and suppressing the cell cycle in both normal and transformed cells. GAS1 is a pleiotropic gene that influences apoptosis and cell arrest in various tumors. GAS1 has been associated with peripheral neuropathy, semilobar holoprosencephaly, and alobar holoprosencephaly diseases in humans [20]. In mice, GAS1 overexpression has been demonstrated to reduce the gene size, proliferation activity, and malignancy of liver tumors [21]. The SMPD3 is a transcription factor that mediates transforming growth factor beta 1 signaling. SMAD3 has been recognized to play crucial roles in various cellular process that control inflammation and fibrosis and also serves as a primary regulator of carcinogenesis [22].

We also identified several downregulated genes that play significant roles in immune response, proliferation, and cytokine networks. Our results showed that genes encoding the Src-like adaptor proteins (SLAPs), such as SLA-1, SLA-3, SLA-5, and SLA-7, were significantly downregulated (FC: −3.15, −2.00, −2.46, and −3.27 respectively) in the mycotoxin-treated livers. These SLA proteins act as negative regulators of receptor signaling. In addition to their well-characterized immunoregulatory functions, SLA proteins also have important roles in the pathogenesis of type I hypersensitivity, osteoporosis, and numerous other malignant diseases. SLA proteins are expressed in most tissues, wherein they exert inhibitory effects on multiple intracellular signaling pathways [23]. Similar to our results, SLA mRNA expression levels were previously reported to be downregulated in human B-cells and were strongly expressed in naïve, pre-germanal center, and germinal-center B-cells based on gene expression analysis [24]. The most significantly downregulated gene, galanin-like peptide (GALP, FC: −14.31), is a protein-coding gene involved in the regulation of appetite that also exerts other functions associated with inflammation, sexual behavior, and stress in humans and various animals [25]. Ganglioneuroblastoma and ganglioneuroma are diseases that are linked to GALP. The observed significant downregulation of the GALP gene could have been caused by the high toxicity of Fusarium mycotoxins. IL10RB (FC: −2.10), an important anti-inflammatory marker gene, was also found to be downregulated in mycotoxin-exposed livers. IL10RB is a cytokine signaling protein belonging to the cytokine receptor family. IL10RB serves as an accessory chain necessary for the formation of the active interleukin 10 receptor complex. Coexpression of the L10RB gene was demonstrated to be essential for IL10-induced signal transduction. Some diseases linked to IL10RB include inflammatory bowel disease 25, autosomal recessive, early onset, and colitis. IL10RB is also associated with immune system pathways and the JAK/STAT signaling pathway [26]. Our results were also consistent with a previous study showing that IL10 was significantly downregulated in the livers of pigs fed with naturally contaminated ZEN [1]. The NNMT gene, which was found to be downregulated (FC: −4.68) in our study, is primarily expressed in the liver, where it acts as a significant cytosolic methyltransferase and is categorized as a phase II metabolic enzyme. NNMT is essential for the biotransformation and detoxification of several xenobiotic compounds and acts by catalyzing the N-methylation of nicotinamide, pyridines, and other structural analogs. Abnormal NNMT expression has been recognized in several kinds of diseases and pathophysiological processes, such as cancer and chemotherapy resistance [27]. IGF1 (FC: −3.90) is a protein that acts as an endocrine hormone and is primarily produced by the liver. IGF1 is associated with growth retardation with deafness and mental retardation due to IGF1 deficiency [28]. IGF1 is a key mediator of the effects of the growth hormone (GH). GH is synthesized in the anterior pituitary gland and subsequently released into the bloodstream, wherein it stimulates the liver to produce IGF1. In turn, IGF1 stimulates systemic body growth and exerts growth-promoting effects in many cell types. PILRBs (FC: −2.15) consist of activating and inhibitory receptors that are highly involved in immune regulation. PILRBs have been recognized as cell signaling activators that contain immunoreceptor tyrosine-based inhibitory motifs (ITIM) that harbor adapter molecules and are found on the cell surface. According to Mousseau et al [29], PILRB mRNA levels in various human tissues were expressed or paired. The group concluded that PILRA and PILRB genes represent a novel ITIM-bearing and non-ITIM-bearing receptor pair. The downregulated CXCL9 (FC: −3.28) is a small cytokine related to the CXC chemokine family. CXCL9 is a T-cell chemoattractant that is stimulated by interferon-γ. CXCL9 is associated with endotheliitis and Sydenham chorea, which are regulated by connected pathways via G protein-coupled receptors and peptide ligand-binding receptors. A study by Altara et al [30] demonstrated that the CXCL9 gene is a reliable biomarker for heart failure and left ventricular dysfunction in humans.

In addition to the abovementioned genes, most of the other DEGs are listed in the Table 2. Some of the non-listed significant DEGs also found to be associated with the immune system and proliferation signaling networks. We hypothesized that toxic levels of DON and ZEN lead to significant reductions in the expression of inflammatory cytokines and decreased levels of other inflammatory mediators. Our analysis revealed that most of the DEGs in the DON and ZEN treatments were downregulated. This observation can be explained by the acute toxicity levels of DON and ZEN present in the diet, which in turn inhibited the expression of inflammation- and proliferation-associated genes. Moreover, some of the upregulated genes were involved in various diseases. Functional annotation analysis also revealed that the DEGs in both the DON and ZEN treatment groups were enriched in GO biological process terms and KEGG pathways that were related to immune system, metabolic process, and diseases (Tables 1, 2, 3).

KEGG enrichment analysis identified significant disease- and immune-related pathways (Table 3, Supplementary Figure S4). Here, we predicted the impact of DON and ZEN toxicity on biological processes in humans and identified some human disease-related pathways associated with DON and ZEN. Almost all the genes were downregulated, which are found in the present study pathways. We hypothesized that severe toxic levels of DON and ZEN could modulate the gene expression, and dowdnregulated the various genes might inhibit the activation, apoptosis and the death of immune cells. Alterations in metabolic pathways can give rise to various diseases. In this study, 13 DEGs were found to be involved in metabolic pathways that are related to certain diseases, cancer, and immune suppression. Diseases arise as a result of dysfunction in certain metabolic processes; however, metabolic systems can tolerate a certain degree of imbalance. A disease will develop only when the activity of a significant enzyme is disrupted or when a control mechanism for a metabolic pathway is altered. The molecular mechanisms underlying the metabolism of most mycotoxins in humans and animals remain unclear. Metabolic pathways of trichothecenes, which mostly include DON, T-2 toxin, 3-aDON, nivalenol, fusarenon-X and diacetoxy scirpenol, in swine, humans, poultry, ruminants, and rodents [31]. The allograft rejection pathway was significantly enriched in the DEGs identified in this study and thus provided insights into the molecular interactions involved in the basic adaptive immune response to allograft destruction. The innate immune system is the critical component involved in the early stage of the allogeneic response, in which chemokines and CAMs play crucial roles not only in mediating leukocyte exodus into the graft but also in expediting dendritic cell and T-cell trafficking between the lymph nodes and the transplant, this outcome in a particular and obtained alloimmune response mediated by T-cells [32]. The strong toxicities of DON and ZEN could be responsible for the observed downregulation of the genes SLA-1, SLA-3, SLA-5, and SLA-7, which are involved in the alloimmune response. The GVHD pathway is responsible for lethal complications that arise as a result of allogeneic hematopoietic stem cell transplantation, wherein immunocompetent donor T-cells attack the genetically different host cells. Liver GVHD is difficult to diagnose and is primarily characterized by jaundice or an increase in alkaline phosphatase levels; increased bilirubin levels leading to hepatomegaly can also be noted. However, GVHD is often difficult to differentiate from other liver dysfunctions that develop following allogeneic hematopoietic cell transplantation, such as mycotoxin infection, veno-occlusive disease, and drug toxicity [33]. Based on the abovementioned results, we hypothesized that acute toxicity caused by DON and ZEN treatment played a role in the development of GVHD in the livers of the piglets. The DEGs were also found to be associated with the AITD, which is characterized by severe chronic inflammation of thyroid gland caused by white blood cells and abnormal blood antibodies that erroneously attack and injure healthy thyroid cells. Untreated AITD can lead to the destruction of the thyroid gland, which eventually leads to hypothyroidism. AITD also participates in crosstalk with other autoimmune endocrine disorders. Animals with experimental autoimmune thyroiditis provide the means to investigate disease-associated pathways. Similarly, farm animals with spontaneous thyroid autoimmune diseases can be used to study underlying molecular pathways of these diseases [34].

The DEGs identified in this study were also associated with the CAMs pathway. CAMs are expressed on the cell surface and play crucial roles in a wide range of biological processes, including the immune response, hemostasis, inflammation, embryogenesis, and improvement of neuronal tissue. The CAMs pathway produces specific signals in response to infection and wounding and control the expression and activation of specific adhesion molecules [35]. The CAMs pathway regulates the expression of the immune response-related genes, namely, CLDN4, SLA-1, SLA-3, SLA-5, and SLA-7, which were downregulated mycotoxin-treated livers. The DEGs were also significantly enriched in the viral carcinogenesis pathway. Viruses that can induce carcinogenesis are major pathogens of humans and farm animals. These viruses promote cell transformation and rapid and uncontrollable cell proliferation, which lead to the progression of malignant tumors. Various additional host factors and exogenous cofactors influence the course of viral carcinogenesis and play roles in various processes ranging from manipulating primary infection to immunological responses to chronic infection. The endogenous viral genomes are inherited as host genetic elements and can be considered as host factors that interact with exogenous carcinogens, such as mycotoxins, chemicals, and ionizing radiation [36]. Viral myocarditis pathway, which is related to cardiac disease, is also associated with inflammation of the myocardium. In addition to viral infection, viral myocarditis can also be caused by bacterial infections, chemical toxins, mycotoxins, certain medications, and autoimmune disorders. Viral myocarditis has also been detected in many farm animals often as a result of viral infection or exposure to toxins, such as chemicals and mycotoxins [37]. Consistent with the results of previous studies, we have also identified important pathogen-related pathways in the DEGs, including type 1 diabetes mellitus, Herpes simplex infection, Epstein-Barr virus infection, and HTLV-I [38].

The liver plays a critical role in the inflammatory response, as well as in the development of acquired immunity [1]. However, little is known about the in vivo effects of DON and ZEN on the expression of liver inflammatory markers in pigs. Exposure of piglets to DON- and ZEN-contaminated diets altered the expression of genes involved in the inflammatory response. Our inflammatory gene expression results were slightly similar with those obtained via microarray analysis in the spleen samples of ZEN-fed pigs [39]. The genes GAS1 (FC: 2.33), CLDN4 (FC: 2.02), ORM1 (FC: 2.34), and UBD (FC: 2.34), which are recognized as some of most potent inflammatory markers, were found to be significantly upregulated in the livers of mycotoxin-treated groups. In turn, overexpression of the above genes can potentially worsen the inflammatory response in lesions formed at the tissue stage [40]. In this study, we observed that consumption of DON- and ZEN-contaminated diets severely downregulated the expression of some of the most critical inflammatory cytokines. Similar to our results, Pistol and his team [1] also demonstrated significant downregulation of pro-and anti-inflammatory cytokine genes in the livers of pigs exposed to ZEN contamination. Based on these results, we hypothesized that DON and ZEN can exert in vivo biphasic effects through the suppression of the inflammatory response in piglet livers, thereby altering immune homeostasis.

Previous studies have reported that the genotoxic effects of ZEN are mediated by impaired protein synthesis as a result of DNA fragmentation. DNA fragmentation has been observed in human hepatocarcinoma cells exposed to ZEN [41]. Moreover, DON treatment reduced the expression of interleukin-6 (IL-6), Tumor necrosis factor- α, and IL-1β in the liver of mice spleen [41]. DON also increased the expression of the Th2-mediated cytokines IL-6 and IL-4 and downregulated those of the cytokines Th1 and interferon-γ in Peyer’s patches of mice infected with retrovirus previously exposed to DON contamination [42]. Differences in the observed effects could have been caused by the varying toxicities and use of different animal models. DON and ZEN caused severe hepatic immunosuppression (significant downregulation of some pro- and anti-inflammatory markers), which could have influenced the infection process. DON and ZEN exert their effects by altering the expression of several genes through interactions with transcription factors and signaling factors. When ZEN binds to estrogen receptors, the ZEA-estrogen receptors complex translocates to the nucleus, wherein it binds to steroid-responsive elements and controls the expression and activity of essential molecules involved in cell signaling and transcription. However, DON binds to ribosomal peptidyl-transferase and particularly suppresses protein and DNA synthesis, leading to reduced cell proliferation [43].

In summary, our results demonstrated that diet containing the Fusarium mycotoxins DON and ZEN induced significant changes in the global transcriptome profiles of piglet livers. Based on RNA-Seq analysis, a total of 249 genes were differentially expressed in the DON and ZEN treatment groups relative to the control group. Of the 249 identified DEGs, 39% were upregulated and 61% were downregulated. Pathway analysis showed that DON and ZEN influenced genes involved in diseases and infections. Moreover, the majority of the DEGs that were downregulated upon DON and ZEN treatment were found to be involved in the inflammatory response, cytokine networks, and cell signaling. The expression patterns of the liver inflammatory cytokines IL10RB and CXCL9 observed in this study were similar to those reported by Pistol et al [1], in which the expression of pro-and anti-inflammatory cytokines were determined to be dramatically downregulated upon exposure to ZEN-contaminated diet. Based on these results, we hypothesized that DON and ZEN exerted in vivo biphasic effects and are involved in the suppression of the inflammatory response in pig livers, which in turn altered immune homeostasis. Thus, the pig can serve as a valuable model for investigating early stages of many infections and elucidating the mechanisms underlying immune system dysregulation in various infections. Further research should focus on the identification of genes show expression patterns are regulated by the unusually expressed RNA-Seq, as they might become excellent biomarkers of DON and ZEN toxicity.