Integrated Profiles of Transcriptome and mRNA m6A Modification Reveal the Intestinal Cytotoxicity of Aflatoxin B1 on HCT116 Cells

Aflatoxin B1 (AFB1) is widely prevalent in foods and animal feeds and is one of the most toxic and carcinogenic aflatoxin subtypes. Existing studies have proved that the intestine is targeted by AFB1, and adverse organic effects have been observed. This study aimed to investigate the relationship between AFB1-induced intestinal toxicity and N6-methyladenosine (m6A) RNA methylation, which involves the post-transcriptional regulation of mRNA expression. The transcriptome-wide m6A methylome and transcriptome profiles in human intestinal cells treated with AFB1 are presented. Methylated RNA immunoprecipitation sequencing and mRNA sequencing were carried out to determine the distinctions in m6A methylation and different genes expressed in AFB1-induced intestinal toxicity. The results showed that there were 2289 overlapping genes of the differentially expressed mRNAs and differentially m6A-methylation-modified mRNAs. After enrichment of the signaling pathways and biological processes, these genes participated in the terms of the cell cycle, endoplasmic reticulum, tight junction, and mitophagy. In conclusion, the study demonstrated that AFB1-induced HCT116 injury was related to the disruptions to the levels of m6A methylation modifications of target genes and the abnormal expression of m6A regulators.


Introduction
Aflatoxins (AFTs) are natural toxins produced by fungi such as Aspergillus flavus and Aspergillus parasiticus and they are widely present in the production and preservation of tree nuts, cereal crops, spices, and animal feed [1]. The concentration of AFTs in cereal products intended for human consumption should not exceed 20 µg/kg according to U.S. food supervision system, while the maximum limited dosage was 4 µg/kg for cereals foods in the European Union [2]. According to global mycotoxin occurrence data, in the 74,821 feed samples collected from 100 countries, 23% of them were positive for AFTs, with the positive samples showing a median concentration of 4 µg/kg. Among the regions where aflatoxin B1 (AFB1) concentrations exceed 20 µg/kg, south Asia, sub-Saharan Africa, and southeast Asia have reported the highest positivity ratios of 41.1%, 38.5%, and 20.9%, respectively [3]. The worldwide AFB1 occurrence in highly used food commodities during 2008 to 2018 was reported, and sorghum, spices, and rice had the highest frequencies of 67.3%, 64.4%, and 54.5%, respectively; among the positive samples, the concentration ranges were 0.2-83.6, 0. 2-25.4, and 0.76-73.2 µg/kg, respectively [4]. AFB1 can be transferred (C) representative fluorescence images of ROS accumulation and quantification analysis of fluorescence intensity. * p < 0.05, ** p < 0.01 or *** p < 0.001 indicates statistically significant difference.

Cell Viability
The HCT116 cells were treated with 0, 25, 50, 100, and 200 μM AFB1 for 48 or 72 h. The cell viability was measured using the CCK-8 assay. In brief, after treatment of AFB1, the CCK-8 solution was added for incubation at 37 °C for 2 h and followed by the measurement of absorbance at 450 nm of a microplate reader (BioTek, Santa Clara, CA, USA).

RNA Preparation and m6A MeRIP-Seq of HCT116
After treatment with 50 μM AFB1 for 72 h, the HCT116 cells were lysed by TRIzol

Cell Viability
The HCT116 cells were treated with 0, 25, 50, 100, and 200 µM AFB1 for 48 or 72 h. The cell viability was measured using the CCK-8 assay. In brief, after treatment of AFB1, the CCK-8 solution was added for incubation at 37 • C for 2 h and followed by the measurement of absorbance at 450 nm of a microplate reader (BioTek, Santa Clara, CA, USA).

RNA Preparation and m6A MeRIP-Seq of HCT116
After treatment with 50 µM AFB1 for 72 h, the HCT116 cells were lysed by TRIzol (Invitrogen, Carlsbad, CA, USA). The total RNA was extracted, and the concentration and integrity were quantified. The enrichment of RNAs with m6A methylation modifications and sequencing were performed by Seqhealth (Wuhan, China). Briefly, polyadenylated RNA enrichment was conducted by beads from 10 µg of total RNA. Then, mRNA was cleaved into fragments of 100-200 nt. The 0.5 mg/mL m6A antibody (Synaptic Systems, Goettingen, Germany) was used for m6A IP. Consequently, the RNA library was constructed for sequencing in DNBSEQ-T7 (MGI Technology Co., Ltd., Shenzhen, China).

Bioinformatic Analysis
The data were analyzed in Trim Galore and mapped to hg19 genome [19] with default parameters. The RNA expression level and the differential expression were analyzed by StringTie (Baltimore, MD, USA) [20] and DEseq [21], respectively. The m6A peak calling were identified using exomepeak2 [22,23]. STREME [24] was applied to identify the m6A motif sequences. The position weight matrix Markov model was integrated in the STREME algorithm. MetaTX [25] (Suzhou, China) was utilized to visualize the distribution of the epitranscriptome profiles. The m6A genes of the healthy human colon were obtained from m6A-TSHub [26]. The GO annotations and KEGG analysis were performed using the DAVID database (Frederick, MD, USA) [27]. The result of m6A conservation and disease association were obtained from ConsRM and RMDisease (Suzhou, China) [28,29]. The substrates of m6A regulators were obtained using starBase v2.0 (Guangzhou, China) [30].

Molecular Docking
Molecular docking was performed in SYBYL (Tripos, St Louis, MO, USA) to explore whether differentially expressed m6A regulators could interact with AFB1. The structures of human proteins were obtained from AlphaFold (DeepMind, London, UK) [31] and PDB and were preprocessed by applying SYBYL to extract ligand substructures, add hydrogen atoms, remove heteroatoms and water molecules, and carry out termini treatment [32]. The AFB1 structure was retrieved from PubChem, then converted to 3D structures using Chem3D (Waltham, MA, USA) and saved as a mol2 file. The AFB1 structure was set in an energy-minimized status [33]. The proteins were docked using the semi-flexible docking method in Surflex-Dock Geom mode. When the total score was higher than 5, the interaction between proteins and molecules was considered a stable model [34].

Intracellular ROS Accumulation
The intracellular ROS was evaluated after cells were treated with 50 µM AFB1 for 72 h. Briefly, the cells were incubated with DCFH-DA probe then washed with PBS to remove probes which did not enter the cells, followed by Hoechst 33342 for nuclear staining. The fluorescence images were analyzed in ImageJ to quantify the level of intracellular ROS.

Quantitative Real-Time PCR
The total RNA was extracted and reverse transcribed into cDNA. The AriaMx (Agilent, Palo Alto, CA, USA) provided quantitative detection and data analysis of mRNA expression levels. The human IGF2BP3 primers sequences used in RT-qPCR were as follows (5 -3 ): forward GAGGCGCTTTCAGGTAAAATAG, reverse AATGAGGCGGGATATTTCGTAT.

Western Blotting
Cells were lysed with RIPA and the protein concentrations were measured by bicinchoninic acid. Then the proteins were separated and transferred onto the membranes. After blocking with skimmed milk, the membranes were incubated with primary and secondary antibodies (Proteintech and Abclonal, Wuhan, China), and then exposed with ECL reagents (Advansta, San Jose, CA, USA).

Statistical Analysis
The data analysis was conducted by SPSS (IBM, New York, NY, USA). One-way analysis of variance was performed. Difference of p < 0.05 was considered statistically significant. At least three independent experiments were performed in each assay.

Cell Viability and ROS Accumulation
The cytotoxic effects of AFB1 on HCT116 cells were measured. As shown in Figure 1B, the AFB1-induced inhibition of cell viability was shown to be dose-and time-dependent. At 72 h, the cell viability was 62.4% in the 50 µM AFB1 group. The effect of ROS accumulation on HCT116 and SW480 cells induced by AFB1 was measured by using a DCFH-DA fluorescent probe. ROS accumulations were significantly increased in the AFB1 groups ( Figure 1C).

Transcriptome-Wide MeRIP-Seq Reveals m6A Modification Pattern after AFB1 Treatment of HCT116 Cells
The summary of reads in MeRIP-seq were shown in Table S1. A total of 104,794 peaks were found in 12,848 genes of the control group, and 94,622 m6A peaks were found in 12,214 m6A genes of the AFB1 group. In the analysis of m6A-TSHub, 6222 m6A genes were reported in the healthy human colon (Figure 2A-C). The STREME was used to define the conserved sequence of RRACH, and the results in both the control and AFB1 groups showed the classical motifs ( Figure 2D). To identify the distribution of m6A modifications throughout the transcriptome, the m6A-methylated transcripts were quantified in each gene, and it was found that more than 4000 genes contained 1-2 m6A peaks in each group ( Figure 2E).

Analysis of m6A Modification Distribution
The distribution pattern of m6A methylation was analyzed in the control and AFB1 groups. As a result, the m6A peaks density were increased between 5′ untranslated regions (5′-UTRs) and the start codon region and remained high in the coding sequences

Analysis of m6A Modification Distribution
The distribution pattern of m6A methylation was analyzed in the control and AFB1 groups. As a result, the m6A peaks density were increased between 5 untranslated regions (5 -UTRs) and the start codon region and remained high in the coding sequences (CDS) region ( Figure 2F).

Differentially Expressed Genes (DEGs) and Differentially Methylated Genes
A total of 5487 DEGs were identified through the RNA-seq data and 5312 differentially m6A-modified mRNAs were identified between the control and AFB1 groups after the analysis of the MeRIP-seq data (p < 0.05). A total of 2289 genes were repeatedly detected as overlapping ( Figure 3A). When the statistical threshold was adjusted to fold change > 2 and p < 0.05, a total of 248 mRNAs were up-regulated and 432 mRNAs were down-regulated ( Figure 3B). genes from the HCT116 control group, AFB1 group, and the healthy human colon data obtained by m6A-TSHub; (D) the motifs enriched from identified m6A peaks based on STREME; (E) the m6A-modified peaks number in genes; (F) the density of m6A-modified peaks in mRNA transcripts.

Analysis of m6A Modification Distribution
The distribution pattern of m6A methylation was analyzed in the control and AFB1 groups. As a result, the m6A peaks density were increased between 5′ untranslated regions (5′-UTRs) and the start codon region and remained high in the coding sequences (CDS) region ( Figure 2F).

Differentially Expressed Genes (DEGs) and Differentially Methylated Genes
A total of 5487 DEGs were identified through the RNA-seq data and 5312 differentially m6A-modified mRNAs were identified between the control and AFB1 groups after the analysis of the MeRIP-seq data (p < 0.05). A total of 2289 genes were repeatedly detected as overlapping ( Figure 3A). When the statistical threshold was adjusted to fold change > 2 and p < 0.05, a total of 248 mRNAs were up-regulated and 432 mRNAs were down-regulated ( Figure 3B).

Biological Pathways of DEGs and Differentially m6A-Modified mRNAs
KEGG and GO enrichment was performed using the DAVID bioinformatic database. The GO enrichment analysis contained subtypes of biological process (BP), cellular component (CC), and molecular function (MF). KEGG of 2289 overlapping genes was involved in the cell cycle, endoplasmic reticulum (ER), autophagy, tight junction, mitophagy, and apoptosis ( Figure 4A). The BP terms were cell cycle, autophagy, and response to ER stress. The CC terms included focal adhesion, cell junction, mitochondrion, and ER lumen. The MF terms contained cadherin binding ( Figure 4B).

Conservation and Disease Association of m6A-Modified Genes
The conservation of m6A sites and disease association were explored in the 2289 overlapping genes. The result showed that 92.4% of the m6A-methylation-modified sites were non-conserved, and 43.5% genes and 37.5% peaks were diseased-associated ( Figure 5).

Conservation and Disease Association of m6A-Modified Genes
The conservation of m6A sites and disease association were explored in the 228 overlapping genes. The result showed that 92.4% of the m6A-methylation-modified sit were non-conserved, and 43.5% genes and 37.5% peaks were diseased-associated (Figu 5).

Proteins Expression Levels
In the Western blotting assay, tight junction and mitophagy biomarkers were detected. As shown in Figure 6, after treatment of AFB1, the CLDN3 was decreased in HCT116 cells, and ZO-1, BECN1, LC3-II and SQSTM1 were increased. The CLDN3, ZO-1 and SQSTM1 were decreased, and BECN1 and LC3-II were increased in SW480 cells. The OCLN was not altered by AFB1 in HCT116 and SW480 cells. Y = Yes, and N = No. Whether the mRNA is the substrate of the m6A regulator.

Proteins Expression Levels
In the Western blotting assay, tight junction and mitophagy biomarkers were d tected. As shown in Figure 6, after treatment of AFB1, the CLDN3 was decreased HCT116 cells, and ZO-1, BECN1, LC3-II and SQSTM1 were increased. The CLDN3, ZO and SQSTM1 were decreased, and BECN1 and LC3-II were increased in SW480 cells. T OCLN was not altered by AFB1 in HCT116 and SW480 cells.

Potential RNA m6A Regulators of Differentially Methylated Genes
Through the analysis of mRNA-seq data, a total of 23 mRNA of m6A regulato were differentially expressed (Table S2 and Figure 7). The writers VIRMA, METT RBM15B, and ZC3H13 were upregulated, and the reader IGF2BP3 was downregulat (Figure 7). The mRNA expression of IGF2BP3 has been verified by RT-qPCR AFB1-treated HCT116 and SW480 cells ( Figure S1).

Potential RNA m6A Regulators of Differentially Methylated Genes
Through the analysis of mRNA-seq data, a total of 23 mRNA of m6A regulators were differentially expressed (Table S2 and Figure 7). The writers VIRMA, METTL5, RBM15B, and ZC3H13 were upregulated, and the reader IGF2BP3 was downregulated (Figure 7). The mRNA expression of IGF2BP3 has been verified by RT-qPCR in AFB1-treated HCT116 and SW480 cells ( Figure S1). nes 2023, 14, x FOR PEER REVIEW Figure 7. The significantly differential expressed m6A regulators in HCT116 treatment. Data were expressed as means ± SD. * p < 0.05, ** p < 0.01, *** p < sus the control.
The GO and KEGG terms were enriched in cell cycle, ER, tight junc We collected the genes that were enriched in the mentioned terms and tory relationship between differentially expressed m6A regulators and Figure 7. The significantly differential expressed m6A regulators in HCT116 with or without AFB1 treatment. Data were expressed as means ± SD. * p < 0.05, ** p < 0.01, *** p < 0.001 or p < 0.0001 versus the control.
The GO and KEGG terms were enriched in cell cycle, ER, tight junction, and mitophagy. We collected the genes that were enriched in the mentioned terms and explored the regulatory relationship between differentially expressed m6A regulators and these target genes by STRING to perform protein-protein interactions network analysis. As shown in Figure 8, the top five genes in cell cycle terms were CCNB1, MYC, CDC25C, ATM, and CHEK2, with the degrees of 33, 31, 29, 28, and 28 ( Figure 8A); in terms of protein processing in ER, the top five genes were HSPA5, CANX, P4HB, HSP90B1, and PDIA3, with the degrees of 30, 25, 23, 21, and 21 ( Figure 8B); in the tight junction terms, the top five genes were ACTB, SRC, ACTG1, MYH9, and ACTR23, with the degrees of 27,21,20,19, and 18 ( Figure 8C); in the mitophagy terms, the top five genes were BECN1, SQSTM1, GABARAPL1, MAP1LC3B, and MFN2, with the degrees of 16, 16, 13, 13, and 12 ( Figure 8D). The gene set enrichment analysis of these four terms in HCT116 cells treated with AFB1 were shown in Figure S2. Figure 7. The significantly differential expressed m6A regulators in HCT116 with or with treatment. Data were expressed as means ± SD. * p < 0.05, ** p < 0.01, *** p < 0.001 or p < 0 sus the control.
The GO and KEGG terms were enriched in cell cycle, ER, tight junction, and mi We collected the genes that were enriched in the mentioned terms and explored th tory relationship between differentially expressed m6A regulators and these target STRING to perform protein-protein interactions network analysis. As shown in Figu top five genes in cell cycle terms were CCNB1, MYC, CDC25C, ATM, and CHEK2, degrees of 33, 31, 29, 28, and 28 ( Figure 8A); in terms of protein processing in ER, the genes were HSPA5, CANX, P4HB, HSP90B1, and PDIA3, with the degrees of 30, 2 and 21 ( Figure 8B); in the tight junction terms, the top five genes were ACTB, SRC, MYH9, and ACTR23, with the degrees of 27,21,20,19, and 18 ( Figure 8C); in the m terms, the top five genes were BECN1, SQSTM1, GABARAPL1, MAP1LC3B, and MF the degrees of 16, 16, 13, 13, and 12 ( Figure 8D). The gene set enrichment analysis four terms in HCT116 cells treated with AFB1 were shown in Figure S2.  To detect the substrates of m6A methylation regulators, the writers RBM15B, VIRMA, METTL5, and ZC3H13 and the readers YTHDC2, YTHDC1, YTHDF3, and IGF2BP3 were selected as well as the cell cycle genes CCNB1, MYC, CDC25C, ATM, and CHEK2 ( Figure 9A); the genes of protein processing in ER HSPA5, CANX, P4HB, HSP90B1, and PDIA3 ( Figure 9B); the tight junction genes ACTB, SRC, ACTG1, MYH9, and ACTR2 ( Figure 9C); and the mitophagy genes BECN1, SQSTM1, GABARAPL1, MAP1LC3B, and MFN2 ( Figure 9D). The m6A methylation levels of transcripts were observed by IGV and shown in Figure 9E.

Interactions between AFB1 and m6A Regulators
The interactions between AFB1 and m6A regulators were predicted and shown in Table 2 and Figure 10. Theoretically, AFB1 bound to regulators by hydrogen bonds and hydrophobic contacts. All the total scores except ZC3H13 were higher than 5, indicating that most of these proteins could form stable and direct interactions with AFB1.

Interactions between AFB1 and m6A Regulators
The interactions between AFB1 and m6A regulators were predicted and sh Table 2 and Figure 10. Theoretically, AFB1 bound to regulators by hydrogen bo hydrophobic contacts. All the total scores except ZC3H13 were higher than 5, in that most of these proteins could form stable and direct interactions with AFB1.

Discussion
Generally, the digestive system represents the first barrier against exposure to foods contaminated with AFB1, so the biological interaction between the intestine and AFB1 is an important issue to be clearly elucidated. To date, most studies related to the toxicity of AFB1 have focused on the metabolic process and hepatotoxicity. The liver and hepatocytes are the main experimental models in toxicological research. However, the liver is not the sole toxicity target of AFB1, and the intestine injury induced by AFB1 merits more attention. In this study, human intestinal cells were exposed to AFB1, and the alterations of m6A modifications and the mRNA expression levels were analyzed through MeRIP-seq.
According to the present sequencing data, the m6A-modified transcripts were mainly enriched near the stop codons, which were consistent with the topology of human RNA m6A methylomes [35]. Additionally, by exploring the distribution of m6A-modified peaks, it was found that most genes contain one or two m6A peaks, which was consistent with a previous description [15].
The KEGG enrichments reported terms of cell cycle, focal adhesion, protein processing in ER, autophagy, mitophagy, hepatitis B, etc. Many vital terms were enriched based on GO analysis. The cell cycle, autophagy, and ER stress were enriched in BP; the actin cytoskeleton, mitochondrion, autophagosome, and ER lumen were enriched in CC; and ATP binding, cadherin binding, and ubiquitin protein ligase binding were enriched in MF. The intestine is particularly susceptible to xenobiotic compounds. When the intestine is exposed to mycotoxins, the intestinal epithelial barrier, which consists of various cellular junctions, is regarded as a mechanical barrier to the toxins [36]. The tight junction, focal adhesion, and adherens junction are important components of cellular junctions. The mycotoxins can affect tight junction proteins, thus compromising the integrity of the intestinal barrier [37]. In this study, the cellular junctions and actin cytoskeleton were damaged in HCT116 cells, which may have resulted in abnormal intestinal permeability, absorption, and efflux.
The toxicity of AFB1 is involved in ROS production, which is one of the main factors in the toxicity toward organs [38]. It has been proven that AFB1-induced oxidative stress can activate ER stress (ERS) [39]. The ER is the major membrane-bound organelle for intracellular calcium storage and protein synthesis, folding, and transport. ROS can initiate ERS and mediate the unfolded protein response, resulting in calcium imbalance [40]. Existing studies have reported that ERS could block the cell cycle to prevent daughter cell generation [41]. Previous studies also demonstrated that exposure to AFB1 could induce cell cycle arrest at different phases among diverse cells, such as arrest at the G0/G1 or S phases in human liver cells [42][43][44]. After the GO and KEGG enrichment in the present study, AFB1 was found to arrest cell cycle at G1/S transition, as shown in Figure 4B. In addition, the ER and mitochondria are highly interconnected organelles, since the ER membrane and the outer mitochondrial membrane are tightly connected to each other, and this region is called the mitochondria-associated ER membrane (MAM) [45]. ERS is associated with mitochondrial dysfunction. When the ER homeostasis is disrupted, calcium ions are released from the ER into the mitochondria, leading to the production of large amounts of ROS in the mitochondria, and mitochondrial bioenergetics are regulated by effecting the MAM. In turn, ER failure is accelerated, and mitochondrial damage is aggravated [46]. Then, autophagy is triggered to selectively identify and eliminate damaged mitochondria [47].
RNA m6A methylation involves in the post-transcriptional modification and affects diverse, fundamental biological processes and functions of mRNA, such as stability, translation, and decay. RNA m6A methylation is regulated by methyltransferases, demethylases, and methylation recognition proteins [48][49][50]. In this study, eight regulators were abnormally expressed. Among them, IGF2BP3 was the m6A-methylation-recognition protein that enhanced mRNA stability. As determined using computational prediction technology, AFB1 was able to bind stably to m6A regulators excluding ZC3H13. Then, this study explored the endpoint effects of the differentially expressed m6A regulator genes on the enriched pathways of cell cycle, mitophagy, and protein processing in the ER and tight junction. For example, in the term of cell cycle, cyclin B1 (CCNB1) and MYC promote the progression of the cell cycle [51,52]. For the term of mitophagy, BECN1 and GABA type A receptor-associated protein such as 1 (GABARAPL1) are involved in the formation of mature autophagosome membranes [53,54]. Once autophagy occurs, light chain 3 (LC3) acts as a biomarker, and the lipidation of LC3-I to LC3-II is observed during autophagosome formation [55]. For the term of protein processing in the ER, the heat shock proteins are essential molecular chaperones for regulating protein homeostasis both under normal physiological conditions and in ERS [56,57]. In this study, the m6A methylation levels of CCNB1 and MYC were reduced in the AFB1 group, along with decreased binding to IGF2BP3; conversely, the m6A methylation levels of BECN1, GABARAPL1, and MAP1LC3B were elevated, along with increased binding to IGF2BP3. As a result, the stability of CCNB1 and MYC mRNA was destroyed, leading to decreased mRNA expression and translation. On the contrary, the mRNAs of BECN1, GABARAPL1, and MAP1LC3B were highly expressed and translated.

Conclusions
In this study, the cytotoxicity of AFB1 toward HCT116 cells was observed. The potential intestinal toxicological effects of AFB1, and the GO and KEGG enriched terms of cell cycle, mitophagy, protein processing in ER, and tight junction, are reported. These injuries result in structural abnormalities and dysfunction of the intestine. This study suggests that the underlying molecular toxicological mechanisms are associated with mRNA m6A methylation; namely, AFB1 might decrease CCNB1 and MYC mRNA expression by reducing the levels of mRNA m6A methylation, resulting in cell cycle arrest. On the other hand, AFB1 induced mitophagy, probably by elevating the level of m6A methylation modifications, mRNA expression, and translation of BECN1, GABARAPL1, and MAP1LC3B. This study provides evidence of AFB1-induced intestinal cytotoxicity and a potential m6A epigenetic regulation mechanism for mRNA.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes14010079/s1, Figure S1: The mRNA expression levels of IGF2BP3 in AFB1-treated HCT116 and SW480 cells; Figure S2: The GSEA of (A) cell cycle, (B) protein processing in endoplasmic reticulum, (C) tight junction, and (D) mitophagy-animal in HCT116 cells treated with AFB1; Table S1: Summary of reads of MeRIP-seq and mRNA-seq in HCT116 cells with or without AFB1 treatment; Table S2: The mRNA expression levels of m6A regulators in AFB1-treated HCT116 cells.

Conflicts of Interest:
The authors declare no conflict of interest.