Genome-Wide DNA Methylome and Transcriptome Analysis of

Subscriber access provided by OCCIDENTAL COLL

Omics Technologies Applied to Agriculture and Food

Genome-wide DNA methylome and transcriptome analysis of porcine intestinal epithelial cells upon deoxynivalenol exposure Haifei Wang, Qiufang Zong, Shiqin Wang, Chengxiang Zhao, Shenglong Wu, and wenbin bao J. Agric. Food Chem., Just Accepted Manuscript • DOI: 10.1021/acs.jafc.9b00613 • Publication Date (Web): 23 Apr 2019 Downloaded from http://pubs.acs.org on April 24, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 30

Journal of Agricultural and Food Chemistry

1

Genome-wide DNA methylome and transcriptome analysis of porcine intestinal epithelial cells

2

upon deoxynivalenol exposure

3 4

Haifei Wang†, Qiufang Zong†, Shiqin Wang†, Chengxiang Zhao†, Shenglong Wu†‡, Wenbin Bao†‡*

5 6

†Key

7

Animal Science and Technology, Yangzhou University, No.48 Wenhui East Road, Yangzhou

8

225009, China

9

‡Joint

Laboratory for Animal Genetics, Breeding, Reproduction and Molecular Design, College of

International Research Laboratory of Agriculture & Agri-Product Safety, Yangzhou

10

University, No.48 Wenhui East Road, Yangzhou 225009, China

11

*Email:

[email protected]

12 13 14 15 16 17 18 19 20 21 22 1

ACS Paragon Plus Environment


Page 2 of 30

23

ABSTRACT: Deoxynivalenol (DON) is a type of mycotoxin disruptive to intestinal and immune

24

systems. To better understand the molecular effects of DON exposure, we performed genome-wide

25

comparisons of DNA methylation and gene expression from porcine intestinal epithelial cell IPEC-J2

26

upon DON exposure using reduced representation bisulfite sequencing and RNA-seq technologies.

27

We characterized the methylation pattern changes and found 3030 differentially methylated regions.

28

Moreover, 3226 genes showing differential expression were enriched in pathways of protein and

29

nucleic acid synthesis, and ribosome biogenesis. Integrative analysis identified 29 genes showing

30

inverse correlations between promoter methylation and expression. Altered DNA methylation and

31

expression of various genes suggested their roles and potential functional interactions upon DON

32

exposure. Our data provided new insights into epigenetic and transcriptomic alterations of intestinal

33

epithelial cells upon DON exposure, and may advance the identification of biomarkers and drug

34

targets for predicting and controlling the toxic effects of this common mycotoxin.

35

KEYWORDS: mycotoxin, cell viability, DNA methylation, gene expression, cytotoxic effects

36 37 38 39 40 41 42 43 44 2


Page 3 of 30


45

INTRODUCTION

46

Deoxynivalenol (DON) is a type of B trichothecene mycotoxin that primarily produced by

47

plant pathogens of Fusarium graminearum and Fusarium culmorum, and occurs predominantly in

48

grains such as wheat, corn, oats, and barley.1 DON is absorbed from gastrointestinal tract through

49

digestion of contaminated food and then reaches the blood compartment and whole body through

50

blood circulation. DON can restrain protein synthesis and disturb immune and metabolic functions,

51

leading to toxic symptoms including vomiting, anorexia, malnutrition, leukocytosis, and diarrhea.

52

Presence of DON was widely detected in food materials, processed food products, animal products,

53

and animal feed,2-3 which represents a risk to food safety and the health of humans and farm animals.

54

Moreover, it is difficult to degrade DON during processing and cooking due to its resistance to high

55

temperature and acidic conditions. Cereal and cereal products contaminated with DON usually have

56

to be used as fuels and even directly abandoned, resulting in huge economic losses for agricultural

57

industry. To control the potential harm caused by DON contamination food, many counties have

58

established DON maximum level in various cereal and cereal products.4 In addition, a maximum

59

tolerable daily intake for DON of 1 ug/kg body weight per day has been proposed by the Joint

60

FAO/WHO Expert Committee on Food Additives (WHO, 2017). Therefore, identifying control

61

strategies for DON contamination is highly important for enhancing food safety in the food chain

62

and attenuating economic losses.

63

DNA methylation is one of the epigenetic modifications that mainly occurs at cytosine at position

64

C5 in CpG dinucleotides in mammal genome. It is reported to be involved in a number of processes

65

such as genomic imprinting, gene transcriptional regulation, development, and tumorigenesis.5

66

Among different kinds of epigenetic marks, DNA methylation is characterized as the most stable and 3



Page 4 of 30

67

easily accessible biomarker candidate. Proper gene expression pattern will be established under the

68

influence of environment stimuli by changing DNA methylation state of the genome, and thereby

69

leading to phenotypic consequences. For instance, pathogenic stimulation, drug treatment and food

70

supply could induce changes in gene-specific or global DNA methylation that further regulate

71

expression of the responsive genes.6-8 Moreover, DON can trigger expression of proinflammatory

72

genes (e.g., TNFα, IL2, IL6), expression and activation of transcription factors (e.g., JunB, EGR1,

73

NF-κB) involved in regulating expression of inflammation and immune related genes.9 Expression of

74

the IL6 and IL8 genes are epigenetically modulated by promoter DNA methylation in Caco-2 cells.10

75

However, systematic investigations on global DNA methylation changes induced by DON exposure

76

and the methylation pattern of responsive genes are still scant until now.

77

In this study, to characterize genome-wide DNA methylation and gene expression in mammal

78

cells exposed to DON, we comprehensively analyzed the methylomic and transcriptomic changes

79

induced by DON exposure in the porcine small intestine epithelial cell IPEC-J2 using reduced

80

representation bisulfite sequencing (RRBS) and RNA-seq technologies. A subset of genes involved

81

in biological processes and pathways for responses to DON exposure were identified and their

82

diverse DNA methylation and expression patterns were observed. The results enriched our

83

understanding of epigenetic and transcriptomic changes of intestinal cells in response to DON

84

exposure and may contribute to the identification of biomarkers and drug targets for predicting and

85

controlling DON contamination.

86

MATERIALS AND METHODS

87

Cell Viability Measurement. IPEC-J2 cells were seeded in a 96-well plate at a density of

88

5×104 cells/ml in 100 μL of culture medium and incubated in a CO2 incubator at 37°C for 24 h. Cell 4


Page 5 of 30


89

were then exposed to different concentrations (0, 300, 500, 1000, 2000, 2500, and 3000 ng/ml) of

90

DON (Sigma, Germany) and cultured for 24 h, 48 h, and 72 h. Cell viability was determined using

91

Cell Counting Kit-8 following the manufacturer’s protocols (Dojindo, Japan). The absorbance at 450

92

nm was measured on a Tecan Infinite 200 microplate reader (Tecan, Switzerland).

93

Sample Preparation and Nucleic Acid Isolation. Cells were seeded in 6-well plate at a density

94

of 5×104 cells/ml and cultured overnight. DON was then added to the medium of experimental wells

95

at a concentration of 1000 ng/ml, and an equal amount of phosphate buffer saline (PBS) was added

96

to the medium of control wells. DON-treated and control cells were cultured for 48 h and collected

97

for nucleic acid isolation. Total RNA and genomic DNA were respectively isolated using the TRizol

98

reagent (Thermo Fisher Scientific, Waltham, USA) and QIAamp DNA extraction kit (Qiagen,

99

Hilden, Germany) according to the manufacturer’s guidelines. Concentration of RNA and DNA

100

samples were quantified with Qubit Fluorometer (Thermo Scientific, USA). Integrity of nucleic acid

101

samples were further checked by electrophoresis in 1% agarose gel and Agilent 2100 Bioanalyzer

102

(Agilent Technologies, Palo Alto, USA). For RNA samples, only those with a RNA integrity number

103

greater than 7.0 were retained in subsequent analysis.

104

Methylation Library Construction, Sequencing and Data Analysis. Genomic DNA of the cell

105

samples were subjected to library construction and reduced representation bisulfite sequencing

106

(RRBS). The MspI enzyme (recognition site CC^GG) was used to digest genomic DNA, and the

107

products were processed with end repairing, adenylation, and 5-methylcytosine-modified adapter

108

ligation. DNA fragments in the length of 40-220 bp were selected with gel extraction. Bisulfite

109

treatment was then performed using the EZ DNA Methylation Gold Kit (Zymo Research, CA, USA)

110

according to the manufacturer’s protocols. DNA fragments were enriched by PCR amplification. For 5



Page 6 of 30

111

library quality control, the insert size of libraries was examined using Agilent 2100, and the library

112

concentration was quantified by quantitative PCR. The libraries were then subjected to Illumina

113

sequencing.

114

Raw reads were filtered by trimming to remove sequencing adaptors and low quality reads.

115

Sequence

116

(https://www.ncbi.nlm.nih.gov/genome/?term=pig) was conducted using Bismark.11 Methylation

117

level of individual cytosine was calculated as the ratio of the number of methylated reads to the total

118

reads number covering the tested cytosine. Differentially methylated loci and differentially

119

methylated region (DMR) between the two groups were called using the DSS software.12 Regions

120

with the minimal length of 50 bp, more than three CG sites, and more than 50% of differentially

121

methylated loci in the regions with P < 1E-05 were considered as DMRs. Adjacent DMRs were

122

combined when the distance between two DMRs was less than 100 bp.

alignment

to

the

genome

assembly

Sscrofa11.1

123

RNA-seq Library Construction and Sequencing. For each sample, 3 μg total RNA was used to

124

enrich mRNA with the poly-T oligo-attached magnetic beads. The purified mRNA was then

125

randomly cleaved into small fragments using NEBNEXT RNA fragmentation buffer (NEB, Beijing,

126

China) according to the manufacturer's instructions. The first strand cDNA was synthesized using

127

random hexamer primer and M-MulV Reverse Transcriptase, and the second strand was synthesized

128

using DNA Polymerase I and RNase H (NEB, Beijing, China). Double-strand cDNA was purified,

129

adenylated at 3’ ends, and ligated with sequencing adaptor. The cDNA fragments with the length of

130

250-300 bp were selected using with AMPure XP system (Beckman Coulter, Beverly, USA) and

131

PCR amplified using Phusion High-Fidelity DNA polymerase (NEB, Beijing, China). PCR products

132

were then purified using AMPure XP system (Beckman Coulter, USA) and library quality was 6


Page 7 of 30


133

assessed on the Agilent Bioanalyzer 2100 system (Agilent Technologies, USA). After cluster

134

generation of the index-coded samples, library preparations were sequenced utilizing the Illumina

135

HiSeq-PE150 high-throughput sequencing platform.

136

RNA-seq Data Quality Control and Gene Expression Quantification. Raw reads were firstly

137

processed through in-house perl scripts to remove reads containing adapter or ploy-N or with base

138

quality score lower than 20. The clean reads were aligned to the genome assembly Sscrofa11.1

139

(https://www.ncbi.nlm.nih.gov/genome/?term=pig) using TopHat2,13 and reads numbers mapped to

140

each gene were calculated using the HTSeq program.14 The fragments per kilobase of transcript

141

sequence per million base pairs of each gene was determined by the length of the gene and reads

142

count mapped to this gene. Differential expression analyses of DON-treated and control groups were

143

performed using the DESeq of R package.15 Benjamini and Hochberg’s method was applied to adjust

144

the resulting P-values for controlling the false discovery rate. Genes with an adjusted P-value 1.5 were defined as differentially expressed.

146

Gene Set Enrichment Analysis (GSEA) and Functional Annotation. Biological processes and

147

pathways potentially perturbed by DON treatment were identified using GSEA software.16 Pathways

148

with FDR q value < 0.2 were considered to be statistically significant. Gene Ontology enrichment

149

analyses were implemented with GOseq, which is based on the Wallenius non-central

150

hypergeometric distribution in R package.17 The GO terms with corrected P-values less than 0.05

151

were regarded as significantly enriched by the tested genes. We employed the KOBAS software to

152

examine the statistical enrichment of genes in KEGG database.18 Pathways with corrected P-values

153

less than 0.05 were considered to be statistically significant.

154

Validation of RNA-seq Data by Quantitative Real-time PCR (qRT-PCR). Total RNA of the 7



Page 8 of 30

155

samples was purified and reversely transcribed into cDNA by using PrimerScript RT Reagent Kit

156

with gDNA Eraser (Takara Biotechnology (Dalian) Co., Ltd, Dalian, China) following the

157

manufacturer’s instructions. The qRT-PCR assays were performed with a volume of 20 μL

158

containing 10 μL SYBR Green Mixture, 1 μL of each primer, 0.4 μl of 50× ROX Reference Dye II, 1

159

μL template of cDNA, and 6.6 μL deionized water. The thermal conditions were as follows: 95°C for

160

15 s, 40 cycles of 95°C for 5 s, 60°C for 30 s. The GAPDH gene was used as the internal control.

161

Primer sequences used for qRT-PCR assays are listed in Table S1. Each qRT-PCR assay was carried

162

out in triplicate, and relative gene expression level was calculated by using the 2-ΔΔCt method.19

163

RESULTS

164

DON Exposure Reduced the Viability of IPEC-J2 Cells. To investigate the cytotoxic effects of

165

DON on cell viability, the IPEC-J2 cells were treated with a serial concentrations of DON (0, 300,

166

500, 1000, 2000, 2500, and 3000 ng/ml) and cultured for different periods of time (24, 48, and 72 h).

167

We measured the cell viability (Figure 1A) and performed least square fitting to predict responses to

168

DON treatment (Figure 1B), which indicated that the cell viability tended to decrease over time with

169

increasing DON dose. As shown in Figure 1A, at 48 h of incubation and under the DON

170

concentration of 1000 ng/ml, cell viability was sharply decreased to 59.1% of control. The cytotoxic

171

effects of DON on IPEC-J2 cell growth were further evidenced by observations of cell morphology

172

and cell confluency (Figure S1). Therefore, cell samples treated with 1000 ng/ml DON and cultured

173

for 48 h were used for subsequent experiments.

174

Profiling of DNA Methylome of IPEC-J2 Cells in Response to DON Treatment. To detect

175

DNA methylation changes induced by DON treatment, four samples under DON treatment and four

176

controls were examined by RRBS. An average of 38.5 million clean reads per sample were yielded, 8


Page 9 of 30


177

and the average bisulfite conversion rate of C to T was greater than 99.5% (Table S2). Among the

178

clean reads, 67.4% on average could be mapped to the pig reference genome (Table S2). The regions

179

with at least 5×sequencing depth were analyzed for methylation level measurement in each cytosine,

180

which covered an average of 19.67% CpG, 8.31% CHG, and 5.44% CHH cytosines in the pig

181

genome (Table S3). As DNA methylation predominantly occurs at CpG cytosines (average

182

methylation level of 41.9% for CpG, of 0.71% for CHG, and of 0.50% for CHH) (Table S3), we only

183

focused on the methylation level of CpG cytosines in subsequent analyses. The distribution of CpG

184

cytosine methylation level in the tested samples appeared to be the bimodal distribution (Figure S2),

185

consistent with the previous findings in porcine tissues and human cells.20,21 To reflect the potential

186

effects of DON exposure on DNA methylome, we performed principal component analysis. The

187

results showed that the samples exposed to DON form a cluster that clearly separates from the

188

control samples (Figure S3). To examine CpG methylation changes at distinct genomic contexts, we

189

computed the mean methylation level at eight genomic contexts. The promoter and UTR5 regions

190

exhibited lower methylation levels than other genomic contexts (Figure 2A) and the methylation

191

level tended to sharply decrease toward to transcription start sites and increase toward to gene bodies

192

(Figure 2B). The pattern of changes in DNA methylation at transcription start sites was consistent

193

with previous findings on other kinds of cells and tissues.22-24

194

To detect changes in DNA methylome induced by DON treatment, we performed differential

195

methylation analysis with smoothing-based methods. A total of 3030 DMRs were identified between

196

the DON-treated and control groups, of which 2093 DMRs showed differential higher methylation

197

levels and 937 differential lower methylation levels in the DON-treated group (Table S4). DNA

198

methylation level in the DMRs and the differences between the two groups were displayed in Figure 9



Page 10 of 30

199

S4. The DMR length was mainly distributed in the range from 50 to 200 bp (Figure S5). DMR

200

distribution analysis demonstrated that the majority of DMRs were located at the CGI, CGI shore,

201

exon, intron, and repeat regions (Figure S6). Based on the value of areaStat, we presented circular

202

visualization that reflects the significance of DMRs and its distribution on each chromosome (Figure

203

3). The identified DMRs were mapped to promoters associated with 328 genes and gene bodies of

204

1717 genes based on their genomic locations (Table S5). A collection of genes including STAT6,

205

SOCS2, and MAPK1 involved in the mitogen-activated protein kinase pathway,25,26 through which

206

DON induces stimulated cellular processes, demonstrated differential methylation changes post

207

DON treatment. Functional annotation of genes associated with DMRs showed that these genes were

208

significantly involved in biological process including “regulation of small GTPase mediated signal

209

transduction (GO:0051056, 28 genes)”, “cellular metabolic process (GO:0044237, 607 genes)”, and

210

“Ras protein signal transduction (GO:0007265, 24 genes)”; cellular component including “cell

211

periphery (GO:0071944, 85 genes)”, “host cell part (GO:0033643, 27 genes)”, and “nucleus

212

(GO:0005634, 198 genes)”; molecular function including “binding (GO:0005488, 1016 genes)”, “ion

213

binding (GO:0043167, 535 genes)”, and “protein binding (GO:0005515, 561 genes)” (Table S6).

214

Moreover, two pathways of “axon guidance (ss04360)”, “Rap1 signaling pathway (ss04015)” were

215

significantly enriched for genes associated with DMRs (Table S7).

216

Characterization of Transcriptomic Changes Induced by DON Treatment. Four IPEC-J2 cell

217

samples treated with DON and four control samples were used for transcriptomic analyses by

218

RNA-seq (Table S8). A total of approximately 609.3 million raw reads were yielded using

219

next-generation RNA sequencing of all the samples, from which 597.2 million clean reads were

220

obtained after quality control analysis, with an average of 74.6 million clean reads per sample (Table 10


Page 11 of 30


221

S8). Alignment analysis showed that about 552 million reads (92.4%) were mapped to the pig

222

genome, of which 528.6 million reads (88.5%) were uniquely mapped (Table S9). Analysis of reads

223

distribution across the genomic regions demonstrated that most of the reads (> 86%) were mapped to

224

exons, and a minority of reads were mapped to introns and intergenic regions (Table S10), indicating

225

our data can efficiently reflect genome-wide gene expression profiling in the tested samples.

226

To unravel the transcriptomic differences between DON-treated and control samples, differential

227

gene expression analysis was conducted. In total, 3226 differentially expressed annotated genes

228

(|log2 fold change| > 1.5, corrected P < 0.05) were identified, comprising 1004 upregulated and 2222

229

downregulated genes (Figure 4). Genes showing differential expression between the two groups are

230

listed in Table S11. To validate the expression pattern of differentially expressed genes from

231

RNA-seq data, nine genes (APOC3, ANXA4, ENTPD5, AIG1, GDPD2, TLR3, ISG15, SDC2, and

232

KLF4) were randomly selected and quantified by qRT-PCR. The results indicated concordant

233

expression patterns of these genes between RNA-seq data and qRT-PCR analyses (Figure 5), with

234

the Pearson’s correlation coefficient of 0.823 (P = 0.006), indicating the high reliability and accuracy

235

of transcriptomic data analyses.

236

GSEA was used to examine biological processes and pathways associated with DON exposure.

237

The results showed that genes significantly associated with pathways such as valine, leucine and

238

isoleucine degradation, drug metabolism by cytochrome P450, glutathione metabolism were

239

downregulated in DON-treated samples (Figure 6, Table S12). The upregulated gene sets were

240

significantly enriched in pathways including ribosome biogenesis in eukaryotes, mRNA surveillance

241

pathway, RNA transport, and TNF signaling pathways (Figure 6, Table S12). DON inhibits protein

242

and nucleic acid synthesis through interacting with ribosome. A subset of genes exhibiting distinct 11



Page 12 of 30

243

expression changes and significantly enriched in protein degradation and ribosome biogenesis were

244

identified herein (Figure 7A and 7B), which indicates their roles in response to the cytotoxic effects

245

of DON.

246

Integrative Analysis of Transcriptome and Methylation Data. To explore the effects of DNA

247

methylation on gene expression, we conducted integrative analyses between DMRs and gene

248

expression profiles. DNA methylation occurring at the gene promoters are usually involved in

249

inhibiting the expression of the corresponding genes.27 We explored the relationship between

250

differential methylation changes at the promoter regions and differential gene expression changes.

251

We found 29 DMRs at the promoter regions displayed inverse correlations with differential

252

expression of the corresponding genes (Table S13). These DMRs associated genes including ASAP3,

253

ST3GAL5, and SLC4A11 which are essential for cell growth and cell proliferation28-30 were

254

hypermethylated and downregulated post DON exposure. As the CGI is an intriguing genomic

255

element that crucial for gene expression regulation,31 we examined whether these DMRs are enriched

256

in CGIs. The results showed that most of these promoter DMRs (27 out of 29) were located in CGIs

257

or CGI shores (Table S13), suggesting the potential roles of promoter CGI methylation changes in

258

controlling the expression of the corresponding genes in response to DON exposure. Moreover, we

259

identified a total of 258 differentially expressed genes were differentially methylated at their gene

260

bodies (Table S14). We classified the methylation level of these genes into three groups, low (0-0.4),

261

medium (0.4-0.7) and high (0.7-1) groups (Figure S8A), and found that the expression level changes

262

of genes in the three groups showing consistent trend with the methylation changes at gene bodies.

263

This result was concordant with the idea that DNA methylation at gene bodies is usually positively

264

correlated with gene expression.27 12


Page 13 of 30


265

DISCUSSION

266

The contamination of food and feeds with DON poses a significant threat to human and animal

267

health. Alterations in epigenetic modifications may be important mechanisms involved in the

268

development of diseases induced by toxicants exposure. 32,33 As DNA methylation encodes important

269

information for normal biology and disease, characterizing DNA methylation pattern across the

270

genome is crucial for understanding the influence of epigenetics.34 Therefore, we herein explored

271

genome-wide changes in DNA methylation from porcine intestinal epithelial cell post DON

272

exposure. We identified more than 3000 DMRs and profiled their distribution across the genome. A

273

large proportion of the DMRs (2093 out of 3030) were hypermethylated, which indicated DON can

274

lead to increased global DNA methylation of the intestinal epithelial cells. The induced epigenetic

275

effects were consistent with the observations of aberrant epigenetic modifications in porcine oocytes

276

induced by DON exposure.35 The DMRs associated genes were highly enriched in physiological

277

processes such as cellular metabolic process, nucleic acid biosynthesis, and protein signaling

278

transduction, through which DON exerts toxic effects at the cellular level. Recent reports mainly

279

focused on the expression of genes involved in these processes and the phenotypic consequences of

280

cell viability and intestinal immunity and functions induced by DON exposure.36,37 Our data set

281

offers a layer of epigenetic insights for comprehensively deciphering the molecular mechanisms

282

underlying DON toxic effects on the intestinal cells.

283

We identified 27 CGI-associated promoters showing inverse correlations between DNA

284

methylation and the corresponding gene expression. Among these genes, VAV3 is an exchange

285

factor for GTPase pathway and stimulates the transcription of IL-2 via activation of the transcription

286

factor NF-AT.38 SLC4A11 and ST3GAL5 have been found to be involved in cell growth and cell 13



Page 14 of 30

287

proliferation.15,31 MGAT3 encodes a glycosyltransferase that catalyzes the addition of

288

β1,4-bisecting-N-acetylglucosamine on N-glycans.39 DUSP9 encodes the dual specificity protein

289

phosphatase that negatively regulate members of the MAPK superfamily, which is associated with

290

cellular proliferation and differentiation.40 The aforementioned genes were previously found to

291

display inverse associations between promoter DNA methylation and their expression in human

292

tissues and cells.41-45 After DON exposure, these genes were all hypermethylated at promoter regions

293

and downregulated (Table S13). It can be speculated that the genes may be functionally linked and

294

regulated by promoter methylation in responses to DON-induced cytotoxic effects.

295

We observed significant reductions in IPEC-J2 cell growth after DON exposure. It has been

296

shown that DON exposure affects cell growth by triggering a significant increase of caspase 3

297

activity which is one marker for the induction of apoptosis and inducing a cell cycle arrest in G2/M

298

phase in IPEC-J2 cells.46 We also identified the significant increased expression of the BCL10 and

299

AEN genes that function in the induction and enhancement of apoptosis47,48. DON is known to induce

300

apoptosis by increasing caspase 3 activity in human colon carcinoma cells and porcine hepatocyte,

301

repress cell growth by retarding cell cycle in the G2/M phase in human epithelial cells, and trigger

302

autophagy and early apoptosis by affecting expression of the genes including LAMP2, LC3, and

303

mTOR in porcine oocyte.35,49,50 These indicated the toxicological effects of DON on cell growth of

304

different cell types by activating the pathways of cell cycle arrest, autophagy and apoptosis.

305

Moreover, recent studies showed that DON exposure could potentiate pro-inflammatory cytokine

306

expression and further stimulate IgA production in vivo and vitro.51,52 IL6 can drive the

307

differentiation of IgA-committed B-cells and further stimulate IgA production.9 We herein found the

308

significant upregulation of IL6 and IL8 after DON exposure (Table S11), consistent with previous 14


Page 15 of 30


309

investigations on human and porcine intestinal cell lines treated with DON in different dose and

310

culture time.44,45,51-53 Furthermore, the expression of IL11 and IL15 were also found to be markedly

311

upregulated after DON exposure (Table S11). These indicated the roles of these cytokines in

312

mucosal immunity against DON exposure. However, the expression of several mediators including

313

TNFSF13B and PIGR that critical to IgA immunity were significantly decreased after DON exposure

314

(Table S11). Expression of TNFSF13B and PIGR were also shown to be suppressed by

315

DON-induced ribosomal stress in human and murine enterocytes.37,54 Given the crucial role of PIGR

316

in transporting IgA into mucosal secretions and acting as the precursor of secretory component,55

317

these findings revealed that DON may influence IgA secretion and further disrupt the mucosal IgA

318

responses and homeostasis by suppressing the expression of PIGR.

319

It has been demonstrated that absorption of DON causes disorders of the digestive system,

320

nervous and endocrine systems,56 which poses the risk of animal production and human health.

321

Biological control is proposed as a promising strategy to control DON contamination. Therefore,

322

identification of biomarker for DON exposure and genes involved in DON detoxification will be of

323

great significance to DON control. DNA methylation has proven to be effective biomarkers for

324

disease diagnosis and prognosis, and response to chemotherapeutic drugs.57 Some cytochrome P450

325

genes involved in the xenobiotic detoxification were not induced in their expression post DON

326

exposure, which was inconsistent with their expressions induced by xenobiotics.58,59 The varied roles

327

of cytochrome P450 genes in different species and under different xenobiotics dose as well as

328

genetically and epigenetically multiplex modulation of their expression may account for the

329

discrepancies of cytochrome P450 genes expression.60,61 We identified increased DNA methylation

330

at gene bodies of the ESR1 gene (Table S5), which plays key roles in reproductive 15



Page 16 of 30

331

and endocrine systems.62 DNA methylation of the ESR1 gene is associated with ESR1 binding and

332

capable of predicting drug response in breast cancer.63,64 DON exposure has been shown to reduce

333

the oocytes maturation capability and embryo development in different animal models.65 These

334

reflected the potential applications of the ESR1 gene methylation as biomarker candidates for

335

DON-induced reproductive and endocrine disorders. In addition, a cluster of genes encoding the

336

DNA repair factor MGMT, nuclear receptor NSD1, apoptotic activator APAF1, and transcriptional

337

regulator PITX2 showed markedly altered methylation after DON exposure (Table S5). DNA

338

methylation alterations of these genes have been shown the clinical value as biomarkers in human

339

diseases.57 The findings indicated the effects of DON on DNA methylation patterns of these genes

340

and may shed lights on their potential use as biomarkers for DON-induced toxicity.

341

In conclusion, we profiled the overview landscape of DNA methylation and gene expression in

342

responses of porcine intestinal cells to DON exposure. Integrative analysis of methylation and

343

transcriptome data allowed the detection of a subset of genes implicated in response to DON

344

cytotoxicity, and the expression changes of these genes may be driven by DNA methylation status.

345

Our findings provided new insights into the molecular effects of DON and contributed to further

346

studying the role of epigenetic modifications in the responses of intestinal cells to DON exposure.

347

Supporting Information

348

DON-treated and control IPEC-J2 cells cultured for 48 h; distribution of CpG cytosine

349

methylation percentage per base of DON-treated and control groups; principal component analysis

350

based

351

methylation between DON-treated and control groups; distribution of DMRs length; number of

352

DMRs mapped to different genomic contexts; pearson correlation and principal component analysis

on

DNA

methylation

level

of

the

16


samples;

heatmap of differential

Page 17 of 30


353

of DON-treated and control samples; analysis of changes in gene body methylation and the

354

corresponding gene expression (PDF)

355

Primer sequences used for qRT-PCR; statistics for RRBS data of all samples; sequencing depth of

356

cytosine methylation; list of DMRs; DMRs mapped to promoters and gene bodies; functional

357

annotation of DMR associated genes; KEGG pathways enriched for DMR associated genes; statistics

358

for RNA-seq data; summary of RNA-seq mapping results; distribution of reads across different

359

genomic regions; differentially expressed genes between DON-treated and control groups; GSEA for

360

transcriptome of DON-treated cells; genes showing inverse correlations between promoter

361

methylation and expression; differentially expressed and differentially methylated genes at gene

362

bodies (XLSX)

363

Funding Sources

364

This work was supported by the National Natural Science Foundation of China (31702082,

365

31772560), the Independent Innovation Fund Project of Agricultural Science and Technology in

366

Jiangsu [CX (16)1003], the China Postdoctoral Science Foundation (2017M621842, 2018T110564),

367

Qing Lan Project of Yangzhou University, and the Priority Academic Program Development

368

(PAPD) of Jiangsu Higher Education Institutions.

369

Notes

370

The authors declare no competing financial interest.

371

REFERENCE

372

(1) Creppy, E. E. Update of survey, regulation and toxic effects of mycotoxins in Europe. Toxicol.

373

Lett. 2002, 127, 19–28.

374

(2) Turner, P. C.; Burley, V. J.; Rothwell, J. A.; White, K. L. M.; Cade, J. E.; Wild, C. P. Dietary 17



Page 18 of 30

375

wheat reduction decreases the level of urinary deoxynivalenol in UK adults. J. Expo. Sci. Environ.

376

Epidemiol. 2008, 18, 392–399.

377

(3) Zain, M. E. Impact of mycotoxins on humans and animals. J. Saudi Chem. Soc. 2011, 15, 129–

378

144.

379

(4) Urbanek, K. A.; Habrowska ‐ Górczyńska, D. E.; Kowalska, K.; Stańczyk, A.; Domińska, K.;

380

Piastowska ‐ Ciesielska, A. W. Deoxynivalenol as potential modulator of human steroidogenesis. J.

381

Appl. Toxicol. 2018, 38, 1450-1459.

382

(5) Schübeler, D. Function and information content of DNA methylation. Nature 2015, 517,

383

321-326.

384

(6) Kiga, K.; Mimuro, H.; Suzuki, M.; Shinozaki-Ushiku, A.; Kobayashi, T.; Sanada, T.; Kim,

385

M.; Ogawa, M.; Iwasaki, Y. W.; Kayo, H.; Fukuda-Yuzawa, Y.; Yashiro, M.; Fukayama, M.; Fukao

386

T.; Sasakawa, C. Epigenetic silencing of miR-210 increases the proliferation of gastric epithelium

387

during chronic Helicobacter pylori infection. Nat. Commun. 2014, 5, 4497.

388

(7) Jiang, S.; Yan, K.; Sun, B.; Gao, S.; Yang, X.; Ni, Y.; Ma, W.; Zhao, R. Long-term high-fat diet

389

decreases hepatic iron storage associated with suppressing TFR2 and ZIP14 expression in rats. J.

390

Agric. Food Chem. 2018, 66, 11612-11621.

391

(8) Swathy, B.; Saradalekshmi, K. R.; Nair, I. V.; Nair, C.; Banerjee, M. Understanding the influence

392

of antipsychotic drugs on global methylation events and its relevance in treatment

393

response. Epigenomics 2018, 10, 233-247.

394

(9) Pestka, J. J. Deoxynivalenol-induced IgA production and IgA nephropathy-aberrant mucosal

395

immune response with systemic repercussions. Toxicol. Lett. 2003, 140, 287–295.

396

(10) Caradonna, F.; Cruciata, I.; Schifano, I.; La Rosa, C.; Naselli, F.; Chiarelli, R.; Perrone, A., 18


Page 19 of 30


397

Gentile, C. Methylation of cytokines gene promoters in IL-1β-treated human intestinal epithelial

398

cells. Inflamm. Res. 2018, 67, 327-337.

399

(11) Krueger, F.; Andrews, S. R. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq

400

applications. Bioinformatics 2011, 27, 1571-1572.

401

(12) Park, Y.; Wu, H. Differential methylation analysis for BS-seq data under general experimental

402

design. Bioinformatics 2016, 32, 1446-1453.

403

(13) Kim, D.; Pertea, G.; Trapnell, C.; Pimentel, H.; Kelley, R.; Salzberg, S. L. TopHat2: accurate

404

alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol.

405

2013, 14, R36.

406

(14) Anders, S.; Pyl, P. T.; Huber, W. HTSeq-a Python framework to work with high-throughput

407

sequencing data. Bioinformatics 2015, 31, 166-169.

408

(15) Anders, S.; Huber, W. Differential expression of RNA-Seq data at the gene level–the DESeq

409

package. Euro. Mol. Biol. Lab. 2012.

410

(16) Subramanian, A.; Tamayo, P.; Mootha, V. K.; Mukherjee, S.; Ebert, B. L.; Gillette, M.

411

A.; Paulovich, A.; Pomeroy, S. L.; Golub, T. R.; Lander, E. S.; Mesirov, J. P. Proc. Natl. Acad. Sci.

412

2005,102:15545–50.

413

(17) Young, M. D.; Wakefield, M. J.; Smyth, G. K.; Oshlack, A. Gene ontology analysis for

414

RNA-seq: accounting for selection bias. Genome Biol. 2010, 11, R14.

415

(18) Mao, X.; Cai, T.; Olyarchuk, J. G.; Wei, L. Automated genome annotation and pathway

416

identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics 2005, 21,

417

3787-3793.

418

(19) Livak, K. J.; Schmittgen, T. D. Analysis of relative gene expression data using real-time 19



Page 20 of 30

419

quantitative PCR and the 2− ΔΔCT method. Methods 2001, 25, 402-408.

420

(20) Meissner, A.; Mikkelsen, T. S.; Gu, H.; Wernig, M.; Hanna, J.; Sivachenko, A.; Zhang,

421

X.; Bernstein, B. E.; Nusbaum, C.; Jaffe, D. B.; Gnirke, A.; Jaenisch, R.; Gnirke, A. Genome-scale

422

DNA methylation maps of pluripotent and differentiated cells. Nature 2008, 454, 766.

423

(21) Choi, M.; Lee, J.; Le, M. T.; Nguyen, D. T.; Park, S.; Soundrarajan, N.; Schachtschneider, K.

424

M.; Kim, J.; Park, J. K.; Kim, J. H.; Park, C. Genome-wide analysis of DNA methylation in pigs

425

using reduced representation bisulfite sequencing. DNA Res. 2015, 22, 343-355.

426

(22) Li, Y.; Zhu, J.;Tian, G.; Li, N.; Li, Q.; Ye, M.; Zheng H.; Yu, J.; Wu, H.; Sun, J.; Zhang,

427

H.; Chen, Q.; Luo, R.; Chen, M.; He, Y.; Jin, X.; Zhang, Q.; Yu, C.; Zhou, G.; Sun, J.; Huang,

428

Y.; Zheng, H.; Cao, H.; Zhou, X.; Guo, S.; Hu, X.; Li, X.; Kristiansen, K.; Bolund, L.; Xu, J.; Wang,

429

W.; Yang, H.; Wang, J.; Li, R.; Beck, S.; Wang, J.; Zhang, X.; Zhang, H. The DNA methylome of

430

human peripheral blood mononuclear cells. PLoS Biol. 2010, 8, e1000533.

431

(23) Huang, Y. Z.; Sun, J. J.; Zhang, L. Z.; Li, C. J.; Womack, J. E.; Li, Z. J.; Lan X. Y.; Lei, C.

432

Z.; Zhang, C. L.; Zhao, X.; Chen, H. Genome-wide DNA methylation profiles and their relationships

433

with mRNA and the microRNA transcriptome in bovine muscle tissue (Bos taurine). Sci. Rep.

434

2014, 4, 6546.

435

(24) Wang, H.; Wang, J.; Ning, C.; Zheng, X.; Fu, J.; Wang, A., Zhang Q.; Liu, J. F. Genome-wide

436

DNA methylation and transcriptome analyses reveal genes involved in immune responses of pig

437

peripheral blood mononuclear cells to poly I: C. Sci. Rep. 2017, 7, 9709.

438

(25) Shuai, K.; Liu, B. Regulation of JAK–STAT signalling in the immune system. Nat. Rev.

439

Immunol. 2003, 3, 900.

440

(26) Zhang, W.; Liu, H. T. MAPK signal pathways in the regulation of cell proliferation in 20


Page 21 of 30


441

mammalian cells. Cell Res. 2002, 12, 9.

442

(27) Jones, P. A. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat.

443

Rev. Genet. 2012, 13, 484-492.

444

(28) Gracheva, E. V.; Samovilova, N. N.; Golovanova, N. K.; Kashirina, S. V.; Shevelev, A.;

445

Rybalkin, I.; Gurskaya, T.; Vlasik, T. N.; Andreeva, E. R.; Prokazova, N. V. Enhancing of GM3

446

synthase expression during differentiation of human blood monocytes into macrophages as in vitro

447

model of GM3 accumulation in atherosclerotic lesion. Mol. Cell. Biochem. 2009, 330, 121-129.

448

(29) Park, M.; Li, Q.; Shcheynikov, N.; Zeng, W.; Muallem, S. NaBC1 is a ubiquitous electrogenic

449

Na+-coupled borate transporter essential for cellular boron homeostasis and cell growth and

450

proliferation. Mol. Cell 2004, 16, 331-341.

451

(30) Ha, V. L.; Bharti, S.; Inoue, H.; Vass, W. C.; Campa, F.; Nie, Z.; de Gramont, A.; Ward, Y.;

452

Randazzo, P. A. ASAP3 is a focal adhesion-associated Arf GAP that functions in cell migration and

453

invasion. J. Biol. Chem. 2008, 283, 14915-14926.

454

(31) Deaton, A. M.; Bird, A. CpG islands and the regulation of transcription. Genes Dev. 2011, 25,

455

1010-1022.

456

(32) Bezek, Š.; Ujházy, E.; Mach, M.; Navarová, J.; Dubovický, M. Developmental origin of chronic

457

diseases: toxicological implication. Interdiscipl. Toxicol. 2008, 1, 29-31.

458

(33) Dai, Y.; Huang, K.; Zhang, B.; Zhu, L.; Xu, W. Aflatoxin B1-induced epigenetic alterations: an

459

overview. Food Chem. Toxicol. 2017, 109, 683-689.

460

(34) Laird, P. W. Principles and challenges of genome-wide DNA methylation analysis. Nat. Rev.

461

Genet. 2010, 11, 191-203.

462

(35) Han, J.; Wang, Q. C.; Zhu, C. C.; Liu, J.; Zhang, Y.; Cui, X. S.; Kim, N. H.; Sun, S. C. 21



Page 22 of 30

463

Deoxynivalenol exposure induces autophagy/apoptosis and epigenetic modification changes during

464

porcine oocyte maturation. Toxicol. Appl. Pharmacol. 2016, 300, 70-76.

465

(36) Pinton, P.; Oswald, I. Effect of deoxynivalenol and other Type B trichothecenes on the intestine:

466

a review. Toxins 2014, 6, 1615-1643.

467

(37) Do, K. H.; Park, S. H.; Kim, J.; Yu, M.; Moon, Y. Ribosome inactivation leads to attenuation of

468

intestinal polymeric Ig receptor expression via differential regulation of human antigen R. J.

469

Immunol. 2016, 197, 847-858.

470

(38) Bustelo, X. R. Vav proteins, adaptors and cell signaling. Oncogene 2001, 20, 6372–6381.

471

(39) Zhao, Y.; Nakagawa, T.; Itoh, S.; Inamori, K. I.; Isaji, T.; Kariya, Y.; Kondo, A.; Miyoshi,.

472

E.; Miyazaki, K.; Kawasaki, N.; Taniguchi, N.; Gu, J. N-acetylglucosaminyltransferase III

473

antagonizes the effect of N-acetylglucosaminyltransferase V on α3β1 integrin-mediated cell

474

migration. J. Biol. Chem. 2006, 281, 32122-32130.

475

(40) Dickinson, R. J.; Delavaine, L.; Cejudo-Marín, R.; Stewart, G.; Staples, C. J.; Didmon, M. P.;

476

Trinidad, A. G.; Alonso, A.; Pulido, R.; Keyse, S. M. Phosphorylation of the kinase interaction motif

477

in mitogen-activated protein (MAP) kinase phosphatase-4 mediates cross-talk between protein kinase

478

A and MAP kinase signaling pathways. J. Biol. Chem. 2011, 286, 38018-38026.

479

(41) Qin, L.; Li, T.; Liu, Y. High SLC4A11 expression is an independent predictor for poor overall

480

survival in grade 3/4 serous ovarian cancer. PLoS One 2017, 12, e0187385.

481

(42) Hui, L.; Zhang, J.; Ding, X.; Guo, X.; Jang, X. Identification of potentially critical differentially

482

methylated genes in nasopharyngeal carcinoma: A comprehensive analysis of methylation profiling

483

and gene expression profiling. Oncol. Lett. 2017, 14, 7171-7178.

484

(43) Pinho, S. S.; Oliveira, P.; Cabral, J.; Carvalho, S.; Huntsman, D.; Gärtner, F.; Seruca, R.; Reis, 22


Page 23 of 30


485

C. A.; Oliveira, C. Loss and recovery of Mgat3 and GnT-III Mediated E-cadherin N-glycosylation is

486

a mechanism involved in epithelial-mesenchymal-epithelial transitions. PLoS one 2012, 7, e33191.

487

(44) Peralta-Arrieta, I.; Hernández-Sotelo, D.; Castro-Coronel, Y.; Leyva-Vázquez, M. A.;

488

Illades-Aguiar, B. DNMT3B modulates the expression of cancer-related genes and downregulates

489

the expression of the gene VAV3 via methylation. Am. J. Cancer Res. 2017, 7, 77-87.

490

(45) Wu, F.; Lv, T.; Chen, G.; Ye, H.; Wu, W.; Li, G.; Zhi, F. C. Epigenetic silencing of DUSP9

491

induces the proliferation of human gastric cancer by activating JNK signaling. Oncol. Rep. 2015, 34,

492

121-128.

493

(46) Diesing, A. K.; Nossol, C.; Panther, P.; Walk, N.; Post, A.; Kluess, J.; Kreutzmann, P.; Danicke,

494

S.; Rothkotter, H. J.; Kahlert, S. Mycotoxin deoxynivalenol (DON) mediates biphasic cellular

495

response in intestinal porcine epithelial cell lines IPEC-1 and IPEC-J2. Toxicol. Lett. 2011, 200,

496

8-18.

497

(47) Kawase, T.; Ichikawa, H.; Ohta, T.; Nozaki, N.; Tashiro, F.; Ohki, R.; Taya, Y. p53 target gene

498

AEN is a nuclear exonuclease required for p53-dependent apoptosis. Oncogene 2008, 27, 3797-3810.

499

(48) Yui, D.; Yoneda, T.; Oono, K.; Katayama, T.; Imaizumi, K.; Tohyama, M. Interchangeable

500

binding of Bcl10 to TRAF2 and cIAPs regulates apoptosis signaling. Oncogene 2001, 20,

501

4317-4323.

502

(49) Bensassi, F.; Golli-Bennour, E.; Abid-Essefi, S.; Bouaziz, C.; Hajlaoui, M.R.; Bacha, H.

503

Pathway of deoxynivalenol-induced apoptosis in human colon carcinoma cells. Toxicology 2009,

504

264, 104–109.

505

(50) Mikami, O.; Yamamoto, S.; Yamanaka, N.; Nakajima, Y. Porcine hepatocyte apoptosis and

506

reduction of albumin secretion induced by deoxynivalenol. Toxicology 2004, 204, 241–249. 23



Page 24 of 30

507

(51) Van De Walle, J.; Romier, B.; Larondelle, Y.; Schneider, Y. J. Influence of deoxynivalenol on

508

NF-kappaB activation and IL-8 secretion in human intestinal Caco-2 cells. Toxicol. Lett. 2008, 177,

509

205–214.

510

(52) Wan, L. Y.; Turner, P. C.; El-Nezami, H. Individual and combined cytotoxic effects of

511

Fusarium toxins (deoxynivalenol, nivalenol, zearalenone and fumonisins B1) on swine jejunal

512

epithelial cells. Food Chem. Toxicol. 2013, 57, 276–283.

513

(53) Li, M.; Cuff, C. F.; Pestka, J. J. Modulation of murine host response to enteric reovirus infection

514

by the trichothecene deoxynivalenol. Toxicol. Sci. 2005, 87, 134–145.

515

(54) Do, K. H.; Choi, H. J.; Kim, J.; Park, S. H.; Kim, K. H.; Moon, Y. SOCS3 regulates BAFF in

516

human enterocytes under ribosomal stress. J. Immunol. 2013, 190, 6501-6510.

517

(55) Johansen, F. E.; Kaetzel, C. S. Regulation of the polymeric immunoglobulin receptor and IgA

518

transport: new advances in environmental factors that stimulate pIgR expression and its role in

519

mucosal immunity. Mucosal Immunol. 2011, 4, 598-602.

520

(56) Maresca, M. From the gut to the brain: Journey and pathophysiological effects of the

521

food-associated trichothecene mycotoxin deoxynivalenol. Toxins 2013, 5, 784-820.

522

(57) Heyn, H.; Esteller, M. DNA methylation profiling in the clinic: applications and

523

challenges. Nat. Rev. Genet. 2012, 13, 679-692.

524

(58) Pizzo, F.; Caloni, F.; Schreiber, N. B.; Cortinovis, C.; Spicer, L. J. In vitro effects of

525

deoxynivalenol and zearalenone major metabolites alone and combined, on cell proliferation, steroid

526

production and gene expression in bovine small-follicle granulosa cells. Toxicon 2016, 109, 70-83.

527

(59) Smith, M. C.; Madec, S.; Pawtowski, A.; Coton, E.; Hymery, N. Individual and combined

528

toxicological effects of deoxynivalenol and zearalenone on human hepatocytes in in vitro chronic 24


Page 25 of 30


529

exposure conditions. Toxicol. Lett. 2017, 280, 238-246.

530

(60) Naselli, F.; Catanzaro, I.; Bellavia, D.; Perez, A.; Sposito, L.; Caradonna, F. Role and

531

importance of polymorphisms with respect to DNA methylation for the expression of CYP2E1

532

enzyme. Gene 2014, 536, 29-39.

533

(61) Wu, Q. H.; Wang, X.; Yang, W.; Nüssler, A. K.; Xiong, L. Y.; Kuča, K.; Dohnal, V.; Zhang,

534

X.J.; Yuan, Z. H. Oxidative stress-mediated cytotoxicity and metabolism of T-2 toxin and

535

deoxynivalenol in animals and humans: an update. Arch. Toxicol. 2014, 88, 1309-1326.

536

(62) Arevalo, M. A.; Azcoitia, I.; Garcia-Segura, L. M. The neuroprotective actions of oestradiol and

537

oestrogen receptors. Nat. Rev. Neurosci. 2015, 16, 17-29.

538

(63) Stone, A.; Zotenko, E.; Locke, W. J.; Korbie, D.; Millar, E. K.; Pidsley, R.; Stirzaker,

539

C.; Graham, P.; Trau, M.; Musgrove, E. A.; Nicholson, R. I.; Gee, J. M.; Clark, S. J. DNA

540

methylation of oestrogen-regulated enhancers defines endocrine sensitivity in breast cancer. Nature

541

Commun. 2015, 6, 7758.

542

(64) Fan, M.; Yan, P. S.; Hartman-Frey, C.; Chen, L.; Paik, H.; Oyer, S. L.; Salisbury JD.; Cheng, A.

543

S.; Li, L.; Abbosh, P. H.; Huang, T. H.; Nephew, K. P. Diverse gene expression and DNA

544

methylation profiles correlate with differential adaptation of breast cancer cells to the antiestrogens

545

tamoxifen and fulvestrant. Cancer Res. 2006, 66, 11954-11966.

546

(65) Yu, M.; Chen, L.; Peng, Z.; Nuessler, A. K.; Wu, Q.; Liu, L.; Yang, W. Mechanism of

547

deoxynivalenol effects on the reproductive system and fetus malformation: current status and future

548

challenges. Toxicol. In Vitro 2017, 41, 150-158.

549 550 25



Page 26 of 30

551 552 553

Figure Captions

554

Figure 1. Cytotoxic effects of DON on IPEC-J2 cells measured by cell viability assay. (A) Cell

555

viability of cells treated with a serial concentrations of DON (0, 300, 500, 1000, 2000, 2500, and

556

3000 ng/ml) and cultured for different periods of time (24, 48, and 72 h). Cell viability is calculated

557

as percentage of controls and expressed as mean ± standard derivation of three independent

558

experiments. (B) The least square fitting model to predict responses of cells to DON treatment. The x

559

axis represents culture time; the y axis, concentration of DON; the z axis, cell viability.

560

Figure 2. Distribution of DNA methylation in different genomic elements (A) and in upstream and

561

downstream of gene bodies (B). The genomic elements (CGI, CGI shore, promoter, UTR5, exon,

562

intron, UTR3, and repeat) of each gene were divided into 20 equal bins, and average methylation

563

level was calculated from the corresponding genomic elements of all genes. The upstream 2 K, gene

564

bodies, and downstream 2 K regions were divided into 50 bins, and the methylation levels of all

565

regions were respectively averaged. CGI: CpG island; UTR5: 5’-untranslated region; UTR3:

566

3’-untranslated region.

567

Figure 3. Circular display of DMR distribution across the pig genome. The outermost circle

568

represents the autosomes and X chromosome. The red and blue dots denote hypermethylated and

569

hypomethylated DMRs, respectively. The larger the dots, the greater the methylation differences

570

between the two groups. TE: transposable element.

571

Figure 4. Volcano plot of differentially expressed genes between DON-treated and control group.

572

The red dots represent significantly upregulated genes; the green dots, significantly downregulated 26


Page 27 of 30


573

genes; the black dots, no significant differential expression.

574

Figure 5. Expression changes of nine genes determined by RNA-seq and qRT-PCR. Fold changes

575

are expressed as ratio of gene expression in DON-treated samples to control samples. The orange and

576

blue bars represent the results of RNA-seq and qRT-PCR, respectively.

577

Figure 6. Gene set enrichment analysis of transcriptome of DON-treated cells. The red means

578

positive correlation; and the blue, negative correlation. NES: normalized enrichment score.

579

Figure 7. Expression of genes quantified by RNA-seq. Gene expression is denoted by color: red,

580

high relative expression; blue, low relative expression. Treated: DON-treated samples; Control:

581

untreated samples.

582 583

584 585

Figure 1

586

27



587 588

Figure 2

589

590 591

Figure 3 28


Page 28 of 30

Page 29 of 30


592 593

Figure 4

594 595

596 597

Figure 5

598 599 600

29



601 602

Figure 6

603

604 605

Figure 7

606 30


Page 30 of 30

Genome-Wide DNA Methylome and Transcriptome Analysis of

Recommend Documents