Transcriptome-wide assessment of the m6A methylome of intestinal porcine epithelial cells treated with deoxynivalenol


 Background: Deoxynivalenol (DON) is a cytotoxic compound found in various food and feed products. N6-Methyladenosine (m6A) is a highly abundant epitranscriptomic marker that modifies a wide range of mRNA molecules in mammalian cells. However, the role of the m6A methylome in DON-induced damage remains poorly understood.Results: In this study, we assessed the transcriptome-wide m6A profile of intestinal porcine epithelial cells (IPEC-J2) treated with 1000 ng/mL DON by m6A sequencing and RNA sequencing. Overall, 5406 new m6A peaks appeared with the disappearance of 2615 peaks in DON-treated IPEC-J2 cells. Genes that were uniquely m6A-modified following DON treatment were found to be associated with the tumor necrosis factor (TNF) signaling pathway. On comparing DON-treated and control cells, we identified 733 differentially expressed mRNAs bearing hyper- or hypomethylated m6A peaks. Further experimental data suggested that METTL3-dependent m6A methylation might also play a role in DON-induced inflammatory response, and CSF2 marker is key functional relevance in the context of DON-induced toxicity. Conclusions: This is the first study to perform a transcriptome-wide assessment of the m6A methylome of IPEC-J2 cells treated with DON. We believe that our findings should be useful for identifying mechanisms whereby m6A modifications influence the outcomes of DON exposure.


Background
RNA modi cations can modulate the stability, function, and metabolic processing of modi ed RNAs. To date, ve primary types of RNA methylation have been described: 5-methylcytosine (m 5 C), N 7 -methylguanine (m 7 G at the 5′-cap), N 6 -methyladenosine (m 6 A), N 1 -methyl adenosine (m 1 A), and pseudouridine (ψ) [1]. Among these, m 6 A modi cations are the most common and functionally important internal RNA modi cations within eukaryotic cells, wherein they serve as key regulators of gene expression [2]. In general, the distribution of m 6 A modi cations is most commonly associated with the stop codon (stopC), transcriptional start site (TSS) sequences, coding sequence (CDS), 3′-untranslated region (3′-UTR) sequences, and 5′-UTR sequences in modi ed RNAs [3]. These m 6 A modi cations can profoundly impact RNA splicing, nuclear export, translation, and degradation [4]. As such, they are closely linked to normal physiological processes, including gametogenesis and cell division, and to pathological conditions, such as cancer and infertility [5]. There is also evidence that m 6 A modi cations are involved in viral replication [6] and the modulation of host-pathogen interactions [7,8]. While many studies on m 6 A modi cations have been conducted in humans and mice, few have assessed the relationship among these modi cations and vital traits in poultry and livestock. Analyses of m 6 A methylation in pigs have largely focused on the impact of these modi cations on fat deposition [9], lipid metabolism [10], stem cell differentiation [11], and liver development [12]. No comprehensive studies have as yet characterized the transcriptomewide distribution of m 6 A modi cations in the context of most porcine diseases, underscoring a novel direction for the ongoing epigenetic studies on economically important livestock.
Mycotoxins are secondary metabolites derived from fungi that often contaminate human and animal food supplies and induce signi cant cytotoxic damage following consumption [13]. Deoxynivalenol (DON) is among the most important mycotoxins associated with impaired intestinal integrity; moreover, it evidently exhibits proin ammatory and immunomodulatory properties [14]. In some previous studies involving an intestinal porcine epithelial cell line (IPEC-J2), exposure to DON was found to induce oxidative stress and in ammation, eventually causing accelerated apoptotic cell death and impaired functionality [15][16][17][18]. The m 6 A status of DON-treated IPEC-J2 cells has not been assessed so far. Thus, we herein conducted an m 6 A-speci c RNA immunoprecipitation assay using the methylated RNA immunoprecipitation sequencing (MeRIP-seq) approach [19] to examine the transcriptome-wide m 6 A modi cation pro le of control and DON-treated cells. Accordingly, we identi ed distinct m 6 A modi cation patterns associated with DON treatment; we further used bioinformatics and qPCR to identify key signaling pathways and genes related to dysregulated m 6 A methylation in IPEC-J2 cells. We believe that our ndings should serve as a foundation for future studies on the roles of m 6 A modi cation in the context of DON-induced cell damage.

Methods
Cell Culture IPEC-J2 cells were obtained from the University of Pennsylvania (PA, USA) and cultured in DMEM/F12 containing 10% fetal bovine serum at 37℃ in a 5% CO 2 incubator. For viability analyses, the cells were added to 96-well plates (2 × 10 4 /well) for 12 h, and then treated with 1000 ng/mL DON (Sigma, Germany) for 24, 48, or 72 h. The Cell Counting Kit-8 (Dojindo, Japan) assay was performed to examine cell viability, as per manufacturer instructions. To assess apoptosis, the cells were treated with DON (1000 ng/mL) for 48 h, followed by staining using the Annexin V-FITC/PI Apoptosis Detection Kit (Solarbio, Beijing, China). For RNA-based experiments, IPEC-J2 cells were added to 6-well plates (1 × 10 6 /well) and incubated until 80% con uence; the cells were subsequently treated with DON (1000 ng/mL) for 48 h. They were then collected for downstream analyses. m 6 A MeRIP-seq and RNA-seq For RNA-seq, three pairs of control and DON-treated cell samples were harvested. TRIzol (Invitrogen, CA, USA) was used for total RNA extraction, and RNA quality was assessed using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scienti c, MA, USA). MeRIP-seq and RNA-seq were performed by CloudSeq Biotech Inc. (Shanghai, China) using previously reported methods [19]. Brie y, the GenSeq™ m 6 A-MeRIP Kit (GenSeq Inc., China) was used, as per manufacturer instructions, to conduct m 6 A RNA immunoprecipitation. The immunoprecipitated (IP) and control input samples were then used for RNA-seq library preparation with the NEBNext® Ultra II Directional RNA Library Prep Kit (New England Biolabs, Inc., USA). An Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., USA) was used to con rm the quality of the prepared libraries, and an Illumina HiSeq 4000 sequencer (Illumina, Inc.) was used for library sequencing with 150 bp paired-end reads. The raw high-throughput m 6 A and RNA-seq data have been uploaded to the Gene Expression Omnibus database (accession no.: GSE156078).

Sequencing data analysis
After sequencing, the samples were subjected to Q30-based quality control. Following 3′-adaptor-trimming, low-quality reads were removed using Cutadapt v1.9.3 [20]. Clean library reads were aligned to a reference genome (Sscrofa11.1, https://www.ncbi.nlm.nih.gov/genome/?term=pig) with Hisat2 v2.0.4 [21]. MACS2 v.2.1.2 was then used to identify m 6 A peaks in MeRIP-seq data [22]. m 6 A peaks that overlapped with transcript exons were selected for further evaluation. Differences in m 6 A methylation rates between DON-treated and control cells were compared using diffReps [23]. Peaks that were signi cantly different between the groups met the following criteria: log2|fold change| ≥ 1 and FDR ≤ 0.05. These peaks were then subjected to gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses.

Statistical analyses
The 2 −ΔΔCt method [25] was used to assess relative gene expression. SPSS 18.0 (SPSS, Inc., IL, USA) was used for statistical analyses. Data are expressed as mean ± standard deviation (SD) and were compared using Student's t-test. P < 0.05 indicated statistical signi cance.

Results
Effects of DON on IPEC-J2 cell viability, in ammation, apoptosis In this study, DON treatment was associated with a marked reduction in the viability of IPEC-J2 cells (Fig. 1a). In comparison with the control group, cell viability signi cantly decreased upon DON treatment for 48 h (P < 0.01, Fig. 1a); moreover, cell con uency of the DON-treated group was lower than that of the control group, and cell morphology showed an obvious alteration (Fig. 1b). We also found that DON-treated cells exhibited signi cantly higher expression levels of IL-6, IL-12, IL-1β, and TNF-α relative to control cells (P < 0.05; Fig. 1c). Consistent with these ndings, DON-treated cells showed signi cantly higher rates of apoptotic cell death relative to control cells, as assessed via ow cytometry (P < 0.01; Fig. 1d-e).

Transcriptome-wide assessment of DON-related m 6 A methylation pro les
Transcriptome-wide m 6 A-sequencing (m 6 A-seq) and RNA-seq of DON-treated and control cells were performed. These two datasets yielded 33.62-43.03 million and 42.17-43.31 million clean reads, respectively, with the RNA-seq dataset also serving as the m 6 A-seq input library in this study (Table S2). Of these clean reads, 80.22% and 85.58% mapped to the reference genome (Sscrofa11.1), respectively.
Model-based ChIP-seq (MACS) analyses led to the identi cation of 15,879 and 13,088 m 6 A peaks in DON-treated and control samples, corresponding to the transcripts of 7265 and 6354 genes, respectively ( Fig. 2a-b). Ben diagram analyses revealed that 10,473 m 6 A peaks and 5649 m 6 A-modi ed genes were detectable in both the groups. In comparison with the control group, the DON-treated group exhibited 2615 new peaks and 5406 peaks that were no longer detectable, con rming that there were substantial differences in global m 6 A modi cation patterns between the groups (Fig. 2c). In addition, m 6 A peaks were found to be characterized by the canonical GGACU motif in the groups (Fig. 2d).
We next sought to understand m 6 A-modi ed peak distributions on a per-gene basis and determined that approximately 30% of all modi ed genes (2321/7970) exhibited a unique m 6 A modi cation peak. Most modi ed genes (5077/7970) exhibited between one and three m 6 A-modi ed sites (Fig. 2e), while genes unique to the DON-treated group typically exhibited only one m 6 Amodi ed site relative to m 6 A-modi ed genes that were unique to the control group (Fig. 2f).
In addition, for both the groups, we assessed m 6 A peak distributions throughout the entire transcriptome to determine the presence of any treatment-speci c sites of preferential modi cation. The m 6 A peaks were thus separated into the start codon (startC), 5′-UTR, CDS, stopC, and 3′-UTR peaks based on their transcript locations. A Metagene pro le of peak density revealed these m 6 A modi cations to be speci cally enriched in startC and stopC regions (Fig. 2g). Total and unique m 6 A peaks were particularly enriched in the vicinity of CDS and start/stopC regions (Fig. 2h-j).

Enrichment of m 6 A-modi ed genes in signaling pathways
We next compared the differences in m 6 A peak abundance between the DON-treated and control groups. Of the 10,473 m 6 A peaks that were evident in both the groups, we selected 3871 differentially methylated sites for further analyses. Relative to the control group, 2156 signi cantly hypermethylated m 6 A peaks (Table S3) and 1715 signi cantly hypomethylated m 6 A peaks (Table S4) were detected in the DON-treated group (fold change ≥ 2 and P < 0.05; Fig. 3a). The Integrative Genomics Viewer tool was used to visualize these methylated sites, which revealed that the GGACU motif surrounded the corresponding m 6 A peaks in the DON-treated and control groups (Fig. 3b).
To assess the potential biological impact of m 6 A modi cations on DON-treated cells, we performed GO and KEGG pathway enrichment analyses of differentially methylated genes. We found that the genes harboring hypermethylated m 6 A sites in DONtreated samples were enriched in many signaling pathways, including the TNF, TGF-β, and NF-κB signaling pathways ( Fig. 3c-d, Table S5-S6). Moreover, the transcripts with DON-speci c m 6 A peaks were found to be involved in several biological processes, including the TNF signaling pathway (Fig. 3e-f).

Identi cation of differentially expressed genes (DEGs) in DON-treated cells
The aforementioned m 6 A-seq input library was utilized as an RNA-seq dataset for elucidating global mRNA expression patterns in DON-treated and control cells. Scatter plots were constructed and hierarchical clustering analyses were performed using these RNA-seq data ( Fig. 4a-b). In total, 1776 DEGs were identi ed upon comparing DON-treated and control samples; 728 genes were upregulated (Table S7) and 1048 were downregulated (Table S8) (fold change > 2 and P < 0.05) in DON-treated samples. GO enrichment analyses indicated that upregulated genes in DON-treated cells were involved in processes such as gene expression and RNA metabolism (Fig. 4c), while KEGG pathway enrichment analyses revealed these genes to be signi cantly enriched in the TNF, NOD-like receptor, and MAPK signaling pathways (Fig. 4d-e, Table S9-S10).
Conjoint analysis of m 6 A-RIP-seq and RNA-seq data Conjoint analyses of our m 6 A-RIP-seq and RNA-seq datasets revealed a strong positive correlation (R = 0.93, P < 0.0001) between differential m 6 A methylation and mRNA expression (Fig. 5a). Of the 2156 hypermethylated m 6 A sites detected via m 6 A-seq, 364 genes were also upregulated ("hyper-up") and 10 were downregulated ("hyper-down") at the mRNA level. In contrast, only three genes exhibiting hypomethylated m 6 A sites were upregulated ("hypo-up"), while 356 were downregulated ("hypo-down") at the mRNA level (Fig. 5b). Notably, 364/367 (99%) upregulated mRNA transcripts in DON-treated samples were associated with hypermethylated m 6 A peaks. In addition, substantially more "hyper-up" and "hypo-down" genes were detected relative to "hyperdown" and "hypo-up" genes ( Fig. 5b). Together, these data suggested that m 6 A hypermethylation is typically associated with the upregulation of gene expression in DON-treated cells.
We then selected 41 candidate genes associated with DON-related pathways (TNF, NF-κB, NOD-like receptor, and MAPK signaling pathways) that exhibited signi cant changes in both m 6 A abundance and mRNA transcript expression in DON-treated cells relative to control cells (Table 1, Fig. 5c). Next, m 6 A peak locations in RNA transcripts were assessed by analyzing the overall expression levels of m 6 A-methylated and non-m 6 A-methylated transcripts in DON-treated and control samples (Fig. 5d). We further analyzed the transcripts that had been classi ed into subgroups (startC, 5′-UTR, CDS, stopC, and 3′-UTR) based on their m 6 A modi cation sites. All m 6 A peaks were analyzed, as were those that were unique to DON-treated and control samples. We found that the transcripts with m 6 A modi cations proximal to CDS were typically expressed at signi cantly lower levels relative to those harboring stopC modi cations (P < 0.01, Fig. 5e). Different genes exhibit variable numbers of m 6 A-modi ed sites, as mentioned earlier. We found that a higher number of m 6 A-modi ed sites per gene was associated with upregulated gene expression (Fig. 5f). While gene expression is controlled by many complex mechanisms, differential m 6 A modi cations seem to clearly impact gene expression; further studies are thus warranted.  Table S11), and it was visualized using Cytoscape [26]. CSF2 was the most central protein in this network (Fig. 6), and CSF2 was a DEM and DMG in RNA-seq and m 6 A-seq analyses, respectively, exhibiting a signi cant difference in its expression level in DON-treated cells (Table 1). To further validate our sequencing results, we performed qPCR to assess the expression level of several hyper-and hypomethylated genes (CSF2, KLF6, CCR7, CCL5, KLF10, PIGR, CD86, DUSP6, TLR3, and STAT2), and we found the expression levels to be similar to those detected using RNA-seq (Fig. 7). Notably, the expression of CSF2 was markedly upregulated in DON-treated cells (P < 0.01).

Expression change of m 6 A methyl marks in DON-treated IPEC-J2 cells
To gain insights into how the m 6 A methyl marks regulated the DON-induced damage, we analyzed the mRNA and protein expression pro les of three m 6 A writers (METTL3, METTL14, and WTAP) and two m 6 A erasers (ALKBH5 and FTO) in DON-treated IPEC-J2 cells. qPCR detection showed (Fig. 8a) methylation patterns in these species [3]. Some other studies were able to characterize m 6 A methylation pro les in Arabidopsis thaliana [27], Oryza sativa [28], and Saccharomyces cerevisiae [29]. However, only few studies have assessed porcine m 6 A methylation pro les to date. In this study, we characterized transcriptome-wide m 6 A modi cation patterns in IPEC-J2 cells using the MeRIP-seq approach. We found that the m 6 A modi cation sites were primarily enriched in CDS and start/stopC regions, with relatively few sites in 3′-UTR and 5′-UTR, and this result was consistent with that of an earlier study involving porcine liver tissues [12]. Similarly, m 6 A peaks have been reported to be primarily concentrated around TSS and stopC regions in human cells. In mice, these peaks were more abundant throughout the internal region of modi ed mRNAs, with the following distribution frequencies: 37% CDS, 30% stopC, 20% 3′-UTR, and 13% TSS [30]. Overall, these ndings demonstrate that at a high level, m 6 A modi cation pro les in mammalian transcriptomes are similar across species. In contrast, in case of A. thaliana, the m 6 A modi cation sites were primarily enriched in 3′-UTR and stop/startC regions, and distinct motifs in m 6 A-modi ed startC regions were noted in A.
thaliana relative to those observed in mammals [27]. Further, in case of O. sativa, the majority of the m 6 A modi cation sites were found to be distributed in intergenic regions, with 70% of such sites being found in CDS and 3′-UTR, and 20% in intronic regions and 5′-UTR [28]. These ndings suggest that m 6 A methylation patterns are species speci c. In general, most m 6 A methylation sites in mammals and plants were signi cantly enriched in stopC regions and 3′-UTR, suggesting that these modi cation patterns are consistent with the typical topological m 6 A patterning of mature mRNAs [27,28,30]. The m 6 A enrichment in these terminal regions seems to be associated with RNA stability and signal transduction, and may further regulate the recruitment of speci c factors necessary for mRNA transport or translation (Wang et al., 2014;Li et al., 2014). The RRACH consensus sequence is characteristic of the canonical conserved domains at which m 6 A methylation occurs [31][32][33], and this was con rmed via our m 6 A high-throughput sequencing results [34][35][36]. No such RRACH sequence, however, was detected in O. sativa mRNAs [28].
Further studies should be conducted to explore differences in species-speci c m 6 A consensus sequences.
At a functional level, m 6 A modi cations modulate diverse processes, including metabolism, differentiation, and disease pathology [11,37]. including heat shock, interferons, and ultraviolet radiation [3]. Environment-and stimuli-speci c m 6 A methylation patterns have also been studied using murine embryonic stem cells, embryonic broblasts, and preadipocytes [38,39]. In this study, on comparing differences in m 6 A methylation patterns between DON-treated and control IPEC-J2 cells, genes that were both upregulated and more frequently m 6 A modi ed were found to be associated with various immunological signaling pathways, such as the TNF and MAPK pathways. Wang et al. (2018) revealed that DON treatment induces toxicity and apoptotic cell death in piglet hippocampal nerve cells via the MAPK signaling pathway [40]. Further, Zhang et al. (2020) reported that in IPEC-J2 cells, DON treatment induced the production of proin ammatory factors by activating p38 and ERK1/2 [18], while  showed that treating IPEC-J2 cells with DON caused in ammatory injury via activation of the NF-κB signaling pathway [41]. With regard to the m 6 A-modi ed genes identi ed in this study, our RNA-seq and qPCR analyses indicated that CSF2 was differentially expressed the most signi cantly in DON-treated cells, and CSF2 was the most central protein in the constructed PPI network (Fig. 6). CSF2 is involved in the TNF signaling pathway and thus regulates human disease pathology [42][43][44], which suggests that CSF2 is of key functional relevance in the context of DON-induced toxicity. Few studies to date, however, have analyzed porcine CSF2, indicating that it is pivotal to systematically assess the association between CSF2 expression and DON treatment via Cas9-mediated knockout and overexpression experiments.
To explore the role of m 6  attenuate LPS induced in ammatory response in macrophages through NF-κB signaling pathway [48]. METTL3 knockdown inhibits osteoblast differentiation and activates the in ammatory response by regulating MAPK signaling in LPS-induced in ammation [49]. Moreover, the m 6 A methylome analysis showed that the upregulated m 6 A-modi ed genes in DON-treated IPEC-J2 cells mainly participate in the TNF and MAPK signaling pathways. Accordingly, the similar expression pattern of the m 6 A methyltransferase METTL3 and m 6 A-modi ed genes suggested that METTL3 might play a functional role in DON-induced in ammatory response, but further veri cation will be required in the future.

Conclusions
To summarize, this study characterized m 6        PPI network constructed with dysregulated m6A-modi ed genes. Cytoscape was used to visualize a network incorporating m6Amodi ed genes/proteins based on the conjoint analysis of m6A-RIP-seq and RNA-seq data. Figure 7 qPCR validation of dysregulated m6A-modi ed genes. Relative mRNA expression levels of 10 representative genes were assessed using qPCR. Bars correspond to mean values of four technical replicates. **P < 0.01.