Transcriptome Profiling of m6A mRNA Modification in Bovine Mammary Epithelial Cells Treated with Escherichia coli

Mastitis is a common disease in dairy cows that is mostly caused by E. coli, and it brings massive losses to the dairy industry. N6-Methyladenosine (m6A), a methylation at the N6 position of RNA adenine, is a type of modification strongly associated with many diseases. However, the role of m6A in mastitis has not been investigated. In this study, we used MeRIP-seq to sequence the RNA of bovine mammary epithelial cells treated with inactivated E. coli for 24 h. In this in vitro infection model, there were 16,691 m6A peaks within 7066 mRNA transcripts in the Con group and 10,029 peaks within 4891 transcripts in the E. coli group. Compared with the Con group, 474 mRNAs were hypermethylated and 2101 mRNAs were hypomethylated in the E. coli group. Biological function analyses revealed differential m6A-modified genes mainly enriched in the MAPK, NF-κB, and TGF-β signaling pathways. In order to explore the relationship between m6A and mRNA expression, combined MeRIP-seq and mRNA-seq analyses revealed 212 genes with concomitant changes in the mRNA expression and m6A modification. This study is the first to present a map of RNA m6A modification in mastitis treated with E. coli, providing a basis for future research.


Introduction
Mastitis is one of the most severe diseases in dairy farming, and its diagnosis and treatment present many challenges. Common measures to prevent mastitis include antibiotic treatments, the discarding of milk, and the elimination of sick cows, all of which cause serious losses [1]. Over the past decades, studies of mastitis have mainly focused on the reaction between the pathogen and the host [2]. The innate immune system of the mammary gland has been shown to minimize mastitis caused by pathogenic microorganisms, thereby cooperating with the acquired immune system. Mammary epithelial cells (MECs), as the first barrier of nonspecific immunity of mastitis, play a crucial role in mastitis following infection by pathogenic microorganisms. These cells secrete a variety of inflammatory factors and chemokines [3] that regulate cell apoptosis, the Toll-like receptor (TLR) pathway, the mitogen-activated protein kinase (MAPK) pathway, and other signaling pathways [4]. Mastitis-associated Escherichia coli (E. coli) [5] could result in an acute course of mastitis with severe clinical symptoms [6,7]. Lipopolysaccharides (LPS), a component of E. coli, can activate TLR signal transduction and cause a series of reaction cascades [8][9][10][11]. It has also been reported that E. coli can induce apoptosis in mastitis [12]. These factors eventually lead to inflammation and tissue damage. Although many studies have been conducted on mastitis, knowledge of its diagnosis and treatment remains insufficient. Therefore, new research on pathogen-specific mastitis may provide new directions for diagnosis and treatment.
N 6 -Methyladenosine (m 6 A) is the most common RNA modification in eukaryotes, and it involves methylation at the N 6 position of RNA adenine [13,14]. In the 1970s, m 6 A modification was detected in eukaryotic mRNA and lncRNA [15], and subsequent studies found that about 1-2% of mRNAs contain m 6 A modification [16]. Studies have demonstrated that m 6 A modification is vital to splicing and editing of mRNA, degradation of polyadenylation, and other RNA processing events [17]. Furthermore, m 6 A modification may promote mRNA nuclear export and translation initiation, and it maintains the structural stability of mRNA with poly A binding proteins. m 6 A modification is mainly involved in three enzymes: "writers", methylation transferases responsible for the methylation modification of RNA, such as the METTLE family and WTAP; "erasers", demethylation enzymes, including FTO and ALKBH5, which clear the methylation modification of RNA; "readers", methylation-reading proteins, such as YTH family protein and elF3 [18][19][20][21][22]. m 6 A modification is involved in many pathophysiological processes in mammals [23]. Chen et al. found that m 6 A modification of circular RNA inhibited innate immunity [24]. Yu et al. reported that m 6 A reader protein YTHDF was closely correlated with inflammation induced by LPS via the expression of inflammatory genes such as nuclear factor kappa B (NF-κB) and MAPK [25]. Zhu and Lu reported that the "eraser" ALKBH5 could regulate the cell apoptosis induced by LPS [26]. However, the nature of the relationship between mastitis and m 6 A modification is still unclear.
In this study, we determined potential m 6 A modification in the inflammation of bovine mammary epithelial cells treated with inactivated E. coli using high-throughput sequencing (MeRIP-seq). Differential m 6 A-modified transcripts were analyzed in E. coliinduced MAC-T. Furthermore, we performed biological function analysis of the differential m 6 A modification and clustered signaling pathways of differentially methylated mRNAs. We also investigated the relationship between m 6 A modification and mRNA transcription. These findings lay a foundation for further exploration of m 6 A RNA modification in mastitis treated with E. coli.

Establishment of Mastitis Model In Vitro
In this study, MAC-T cells were treated with heat-inactivated E. coli with an MOI of 10:1 for 24 h. We detected the mRNA and protein expression of IL-1β, IL-6, and TNF-α by RT-qPCR (Figure 1a-c) and ELISA (Figure 1d-f). The results show that the expression levels of proinflammatory factors IL-1β, IL-6, and TNF-α were significantly increased in the MAC-T cells treated with inactivated E. coli (p < 0.05). Meanwhile, we used flow cytometry to detect the apoptosis of MAC-T cells induced by inactivated E. coli. The results show that the percentage of apoptotic cells in the control group was 1.92% and that in the E. coli group was 25.21% (Figure 1g), suggesting that inactivated E. coli could cause MAC-T cell apoptosis. (a-c) Relative mRNA levels of IL-1β, IL-6, and TNF-α measured using real-time PCR in the Con and E. coli groups. (d-f) Concentrations of IL-1β, IL-6, and TNF-α measured using ELISA in the Con and E. coli groups. (g) The apoptotic cell rate in the Con and E. coli groups according to flow cytometry. Data are expressed as the mean ± standard deviation (SD) and analyzed using one-way analysis of variance; *, **, and *** represent p < 0.05, p < 0.01, and p < 0.001, respectively.

Profile of the m 6 A Modification in MAC-T Cells Treated with E. coli
In order to obtain a map of m 6 A modification in MAC-T cells treated with E. coli, we performed a transcriptome profiling of m 6 A modification analysis using meRIP-seq. Comparing the clean reads to the reference genome, there were 16,691 m 6 A peaks within 7066 mRNA transcripts in the Con group and 10,029 peaks within 4891 mRNA transcripts in the E. coli group ( Figure S1a,b, Supplementary Materials). A total of 9005 m 6 A modification peaks and 4675 mRNA transcripts existed in the two groups, whereby 7692 peaks and 2391 mRNA transcripts existed in the Con group and 1030 peaks and 215 mRNA transcripts existed in the E. coli group (Figure 2a,b). The most enriched motif sequence [27,28] m 6 A peak in the Con and E. coli groups was GGACU (Figure 2c).
An analysis of the specific positions of peaks on modified genes revealed that methylation in the CDS region was significantly greater than in other regions in the Con and E. coli groups (Figure 2d), consistent with the results reported in previous studies [29]. Furthermore, the distribution density of methylation peaks in the Con and E. coli groups was highly similar (Figure 2e). Studies have confirmed that there is sometimes more than (a-c) Relative mRNA levels of IL-1β, IL-6, and TNF-α measured using real-time PCR in the Con and E. coli groups. (d-f) Concentrations of IL-1β, IL-6, and TNF-α measured using ELISA in the Con and E. coli groups. (g) The apoptotic cell rate in the Con and E. coli groups according to flow cytometry. Data are expressed as the mean ± standard deviation (SD) and analyzed using one-way analysis of variance; *, **, and *** represent p < 0.05, p < 0.01, and p < 0.001, respectively.

Profile of the m 6 A Modification in MAC-T Cells Treated with E. coli
In order to obtain a map of m 6 A modification in MAC-T cells treated with E. coli, we performed a transcriptome profiling of m 6 A modification analysis using meRIP-seq. Comparing the clean reads to the reference genome, there were 16,691 m 6 A peaks within 7066 mRNA transcripts in the Con group and 10,029 peaks within 4891 mRNA transcripts in the E. coli group ( Figure S1a,b, Supplementary Materials). A total of 9005 m 6 A modification peaks and 4675 mRNA transcripts existed in the two groups, whereby 7692 peaks and 2391 mRNA transcripts existed in the Con group and 1030 peaks and 215 mRNA transcripts existed in the E. coli group (Figure 2a,b). The most enriched motif sequence [27,28] m 6 A peak in the Con and E. coli groups was GGACU (Figure 2c). one methylation peak of mRNA, but one peak is usually the most common. In this study, a count of the methylation peaks of mRNAs revealed that the number of mRNAs with one methylation peak was the largest, although genes with multiple methylation peaks were not uncommon (Figure 2f).

Differential m 6 A Modification between the E. coli Group and Control Group
In order to further study the role of m 6 A modification in dairy bovine mastitis, we analyzed the differences between the m 6 A peaks of RNA in the E. coli and Con groups. Compared with the Con group, there were 2904 significantly differential m 6 A peaks within 2101 mRNAs in the E. coli group (Table S2 Table 1 shows the top 10 hypermethylated or hypomethylated An analysis of the specific positions of peaks on modified genes revealed that methylation in the CDS region was significantly greater than in other regions in the Con and E. coli groups (Figure 2d), consistent with the results reported in previous studies [29]. Furthermore, the distribution density of methylation peaks in the Con and E. coli groups was highly similar (Figure 2e). Studies have confirmed that there is sometimes more than one methylation peak of mRNA, but one peak is usually the most common. In this study, a count of the methylation peaks of mRNAs revealed that the number of mRNAs with one methylation peak was the largest, although genes with multiple methylation peaks were not uncommon (Figure 2f).

Differential m 6 A Modification between the E. coli Group and Control Group
In order to further study the role of m 6 A modification in dairy bovine mastitis, we analyzed the differences between the m 6 A peaks of RNA in the E. coli and Con groups.
Compared with the Con group, there were 2904 significantly differential m 6 A peaks within 2101 mRNAs in the E. coli group (Table S2, Supplementary Materials), of which 644 hypermethylated peaks were within 474 mRNAs such as BCL2 ( Figure 3a) and 2260 hypomethylated peaks were within 1627 mRNAs, such as TIRAP and TLR4 (Figure 3a) (p < 0.00001, fold change >2.0). Table 1 shows the top 10 hypermethylated or hypomethylated genes in E. coli groups compared with the Con group. This study also demonstrated that BCL2 was m 6 A hypermethylated (Figure 3b), while TLR4 was m 6 A hypomethylated (Figure 3c) in E. coli-induced MAC-T cells. Additionally, compared with the Con group, several differential methylation sites in the E. coli group were located on chromosomes 7, 18, 19, and 25 ( Figure 3d). genes in E. coli groups compared with the Con group. This study also demonstrated that BCL2 was m 6 A hypermethylated (Figure 3b), while TLR4 was m 6 A hypomethylated (Figure 3c) in E. coli-induced MAC-T cells. Additionally, compared with the Con group, several differential methylation sites in the E. coli group were located on chromosomes 7, 18, 19, and 25 ( Figure 3d). Figure 3. Distribution of significantly differential peaks between the Con and E. coli groups. (a) Data visualization analysis of differential m 6 A peaks in selected mRNAs (BCL2, TIRAP, and TLR4) in the E. coli group compared with the Con group. MeRIP-qPCR verified the differential m 6 A modification in BCL2 (b) and TLR4 (c). (d) Distribution of differential m 6 A sites on bovine chromosomes. Figure 3. Distribution of significantly differential peaks between the Con and E. coli groups. (a) Data visualization analysis of differential m 6 A peaks in selected mRNAs (BCL2, TIRAP, and TLR4) in the E. coli group compared with the Con group. MeRIP-qPCR verified the differential m 6 A modification in BCL2 (b) and TLR4 (c). (d) Distribution of differential m 6 A sites on bovine chromosomes.

Differential m 6 A Modifications Participate in Important Biological Pathways
In order to demonstrate the important function of m 6 A modification in MAC-T cells induced with E. coli, GO and KEGG enrichment analyses were performed on the above molecules with differential m 6 A methylation. Regarding the biological process function (BP) of GO enrichment analysis, we observed that genes with m 6 A hypermethylation in the E. coli group were mainly enriched in processes related to the regulation of signal transduction, protein modification, cell differentiation, etc. (Figure 4a), whereas genes with m 6 A hypomethylation were mainly enriched in RNA biosynthetic and cellular macromolecule biosynthetic processes (Figure 4b).
Furthermore, KEGG analysis annotated m 6 A hypermethylated genes in the E. coli group as participating in the MAPK signaling pathway, circadian rhythm, Ras signaling pathway, nitrogen metabolism, sphingolipid signaling pathway, etc. (Figure 4c), whereas m 6 A hypomethylated genes were significantly associated with the sulfur relay system, ribosome biogenesis in eukaryotes, the Wnt signaling pathway, the NF-κB signaling pathway, the Hippo signaling pathway, etc. (Figure 4d).

Conjoint Analysis of mRNA-Seq and MeRIP-Seq
To further explore the relationship between m 6 A modification and mRNA expression, we used RNA sequencing data from the input samples to analyze the mRNA expression of the control and E. coli groups. Volcano plots showed significantly different mRNA expression distribution in the E. coli group (p < 0.05, fold change >2.00) (Figure 5a). Hierarchical clustering showed the similarity of relative expression levels between the Con and E. coli groups (Figure 5b). On the basis of this analysis, we identified 190 upregulated and 595 downregulated genes in the E. coli group. Table 2 shows the top 20 most significantly different genes in the Con and E. coli groups.

Conjoint Analysis of mRNA-Seq and MeRIP-Seq
To further explore the relationship between m 6 A modification and mRNA expression, we used RNA sequencing data from the input samples to analyze the mRNA expression of the control and E. coli groups. Volcano plots showed significantly different mRNA expression distribution in the E. coli group (p < 0.05, fold change >2.00) (Figure 5a). Hierarchical clustering showed the similarity of relative expression levels between the Con and E. coli groups (Figure 5b). On the basis of this analysis, we identified 190 upregulated and 595 downregulated genes in the E. coli group. Table 2 shows the top 20 most significantly different genes in the Con and E. coli groups.
A combined analysis showed differential mRNA expression levels (p <0.05, fold change >1.5) and m 6 A methylation peaks (p < 0.00001, fold change >1.5). The correlation coefficient between m 6 A modification and mRNA expression in the E. coli group compared with the Con group was 0.14 ( Figure 5c), indicating a weak connection between m 6 A and the overall mRNA expression level. In addition, these genes were divided into four parts, including 25 hypermethylated and upregulated genes (hyper-up), 154 hypomethylated and downregulated genes (hypo-down), 15 hypermethylated but downregulated genes (hyper-down), and 18 hypomethylated but upregulated genes (hypo-up) (Figure 5d). Table 3 shows the 20 genes with significantly differential m 6 A modification (p <    A combined analysis showed differential mRNA expression levels (p <0.05, fold change >1.5) and m 6 A methylation peaks (p < 0.00001, fold change >1.5). The correlation coefficient between m 6 A modification and mRNA expression in the E. coli group compared with the Con group was 0.14 ( Figure 5c), indicating a weak connection between m 6 A and the overall mRNA expression level. In addition, these genes were divided into four parts, including 25 hypermethylated and upregulated genes (hyper-up), 154 hypomethylated and downregulated genes (hypo-down), 15 hypermethylated but downregulated genes (hyper-down), and 18 hypomethylated but upregulated genes (hypo-up) (Figure 5d). Table 3 shows the 20 genes with significantly differential m 6 A modification (p < 0.00001, fold change >1.5) and mRNA expression (p < 0.05, fold change >1.5) in the Con and E. coli groups.

Discussion
Bovine MECs are the first barrier against invading microorganisms in the mammary gland. The pathogen recognition receptors (PPRs) on the cell surface activate pathogenassociated molecular patterns (PAMPs) after recognizing related pathogens. This major interaction initiates a series of downstream regulatory expression, leading to the expression and release of antimicrobial molecules, chemokines, and cytokines by host cells, even causing apoptosis [30]. Despite the volume of studies that exist on mastitis, there are still uncertainties about the pathogenesis of mastitis caused by E. coli. Therefore, the prevention, diagnosis, and treatment of mastitis remain challenging. The abnormality of m 6 A-modified enzymes can cause a series of diseases [31]; however, the mechanism of m 6 A modification in bovine mastitis remains unclear. This study is the first to analyze the relationship between the m 6 A modification profile and E. coli-induced mastitis.
Using MeRIP-seq to assess m 6 A modification in a mastitis model, we obtained an overview of m 6 A modification in mastitis treated with inactivated E. coli infection. The total m 6 A modification peak numbers revealed significant differences in m 6 A modification between the control and E. coli groups. Therefore, we assume that m 6 A modification may be related to E. coli-induced mastitis.
It is known that m 6 A modification of mRNA often affects changes in and the development of diseases. In this study, we identified a total of 2904 distinct m 6 A methylation peaks in 2101 mRNAs; 474 mRNAs had 622 hypermethylated peaks and 1627 mRNAs had 2260 hypomethylated peaks. Our analysis of the differential methylation patterns revealed some key mRNAs that may be closely related to mastitis (Figure 3a). BCL2 was found to be hypomethylated. As one of the key factors in the study of apoptosis, BCL2 can inhibit cell apoptosis and protect tissue [32,33]. Therefore, we suspect that m 6 A modification of BCL2 is involved in the process of apoptosis in mastitis. TLR4 is an important membrane protein receptor that can recognize lipopolysaccharides on the surface of Gram-positive bacteria such as E. coli [10,34,35], and TIRAP, a Toll/IL-1 receptor (TIR) domain of cohesion in cells, upon identifying PAMPs, can activate a series of downstream immune responses [36,37]. We speculate that the m 6 A hypomethylation in these two key molecules (TLR4 and TIRAP) by MeRIP-seq and IGV visualization reduces the rapid inflammatory reaction of MAC-T treated with E. coli. The findings regarding m 6 A modification of these genes suggest that differential m 6 A modification somewhat enriches the pathogenic mechanism of E. coli infection of bovine mammary glands. Although previous studies have confirmed the important role of some genes in the pathophysiology of mastitis, few studies have been conducted on m 6 A modification. Whether the changes resulting from the m 6 A modification of genes play an important role in mastitis caused by E. coli requires further verification.
E. coli can not only activate the natural immunity of dairy cow mammary gland tissues through the nucleotide binding oligomerization domain (NOD)-like receptor signaling pathway, the TLR signaling pathway, and the NF-κB signaling pathway, but it can also cause cell apoptosis [3,4,8]. Studies have confirmed that the Ras-MAPK signaling pathway affects proliferation, differentiation, and apoptosis by affecting gene transcription and regulation [38,39]. The NF-κB signaling pathway is a key regulatory pathway for immune response, apoptosis, and differentiation [40]. In our study, pathway enrichment analysis showed that the m 6 A differential methylation peaks were mainly enriched in the MAPK signaling pathway, NF-κB signaling pathway, Ras signaling pathway, and Hippo signaling pathway. This evidence also indicates that m 6 A modification is probably associated with mastitis.
Nevertheless, in the combined analysis of mRNA-seq and MeRIP-seq, we found a weak correlation [41] of m 6 A modification and mRNA expression between the E. coli group and the control group. Meanwhile, there were 212 genes in the E. coli group that exhibited differential methylation with significant differential expression of mRNA compared with the Con group (Figure 5a,b). The change in expression of these genes may be closely related to m 6 A modification [42]. These 212 genes (Table 3) included BAD, a protein in the BCL-2 family which plays a key role in mitochondrial-dependent apoptosis and participates in the development of many diseases by regulating cell death [43][44][45]. TMEM214 was also included, which acts as an anchor for the recruitment of procaspase 4 to the endoplasmic reticulum and its subsequent activation, and which is essential for endoplasmic reticulum stress-induced apoptosis [46]. MAP3K2, MAP2K1, and MAPK12 were also identified. Belonging to the MAPK cascade family of molecules, these proteins participate in important biological processes such as cell proliferation, differentiation, and immune response [34,35]. Whether the difference in m 6 A modification of these molecules affected the changes in their mRNA levels and, in turn, affected the release of inflammatory factors and the occurrence of cell apoptosis requires further experimental confirmation. In addition to affecting mRNA expression, m 6 A modification has many other functions, such as affecting the splicing of mRNA precursors, regulating the nuclear export of RNA, regulating mRNA translation, and affecting the stability of mRNAs [47,48]. The genes that did not exhibit expression changes may have one or more of these functions, which should be confirmed by further experiments.
The model established using inactivated bacteria retained the main infectious components of E. coli, while avoiding bacterial growth and apoptosis. Inactivated bacterial infection model experiments have been widely recognized and applied in many mastitis studies [49][50][51]. Although MAC-T cells cannot fully represent the microenvironment during the natural infection process of mastitis, MECs are the first barrier in mastitis, not macrophages or neutrophils [30,51,52]. Therefore, MECs are suitable in vitro models for mastitis. We suspect that m 6 A modification affects the transcription and translation of mRNA in MAC-T cells treated with E. coli and, in turn, affects physiological and pathological processes such as inflammatory release and cell apoptosis ( Figure 6). However, the mechanism of m 6 A-regulated mastitis in dairy cows is not clearly understood. Through the use of advanced technologies, we provide the first m 6 A transcriptome profile of mastitis and an initial map revealing the function of m 6 A modification in mastitis, thereby contributing critical insights for further research on the role of m 6 A in mastitis.

Bacteria Strains and Cell Line
E. coli (ATCC 25922) was donated by associated Professor Wang Xiangru of Huazhong Agricultural University, respectively. E. coli was resuscitated in LB (Luria-Bertani) solid medium. Single colonies were grown at 37 °C overnight. A single bacteria colony was used to inoculate culture bottle containing 10 mL LB broth and incubated in a shaker for 220 r/min. After counting the bacteria, the bacteria were heat-inactivated at 63 °C for 30 min.

Bacteria Strains and Cell Line
E. coli (ATCC 25922) was donated by associated Professor Wang Xiangru of Huazhong Agricultural University, respectively. E. coli was resuscitated in LB (Luria-Bertani) solid medium. Single colonies were grown at 37 • C overnight. A single bacteria colony was used to inoculate culture bottle containing 10 mL LB broth and incubated in a shaker for 220 r/min. After counting the bacteria, the bacteria were heat-inactivated at 63 • C for 30 min.
MAC-T cells (an immortalized bovine mammary epithelial cell line) were donated by Professor Mark Hanigan of Virginia Tech University. The culture medium formulations of MAC-T cell differed from those reported in the literature [10]. MAC-T cells were cultured in DME/F12 medium (Hyclone, Tauranga, New Zealand) supplemented with 10% fetal bovine serum (Gibco, New York, NY, USA), and incubated at 37 • C in a 5% CO 2 humidified incubator. The cells were digested and passaged using 0.25% trypsin and 0.02% EDTA.

Sample Collection and RNA Extraction
Trypsinized MAC-T cells were counted and seeded in a cell culture dish (Corning, New York, NY, USA) with 10 6 cells per well. After 12 h, the DME/F12 maintenance medium containing 2% FBS was replaced. Three groups of cells were set up in triplicate. To the control group (Con) was added only the LB broth. The E. coli group was induced with 10 7 inactivated E. coli at a multiplicity of infection (MOI) of 10:1 and incubated for 24 h. The medium was discarded, and the cells were gently washed thrice using cold PBS (Hyclone, Tauranga, New Zealand), before adding 1 mL of TRIzol reagent (Invitrogen, Carlsbad, CA, USA) to each well for RNA cell lysate collection.
Total RNA was extracted according to the commercial reagent manufacturer's instructions (Invitrogen) and RNA isolation procedures. The NanoDrop 2000 instrument (Thermo Fischer Scientific, Waltham, MA, USA) was used to measure the RNA concentrations, and the RNA purity index was denoted by an OD 260 /OD 280 value between 1.80 and 2.10. The integrity of the sample RNA and potential gDNA contamination were assessed using agarose gel electrophoresis.

Enzyme-Linked Immunosorbent Assay (ELISA)
Cell-free supernatants were collected from the groups of bacteria treatments and assayed for proinflammatory cytokines (IL-1β, IL-6, TNF-α) using an ELISA kit (Cusabio, Wuhan, China). Based on the manufacturer's instructions, the OD of each sample was assessed with a microplate reader at 450 nm wavelength. A standard curve was constructed to calculate the concentration of the samples.

Flow Cytometry
MAC-T cells were digested with trypsin without EDTA for 3 min. After stopping the digestion with DMEM/F12 containing 10% FBS, the cells were collected by centrifugation. After washing three times with PBS, the cells were re-suspended with 100 µL of Binding Buffer. Then, 5 µL each of FITC and PI dyes (Vazyme) were added and left for 10 min at room temperature in the dark. After adding 400 µL of Binding Buffer, the cells were detected using cytoflex-LX (Beckman Coulter, Indianapolis, IN, USA).

Methylated RNA Immunoprecipitation (MeRIP)
Firstly, the collected RNA was fragmented. Secondly, protein A magnetic beads (Thermo Fisher Scientific, Waltham, MA, USA) were incubated with anti-m 6 A antibody (abcam, Cambridge, UK) at room temperature for 1 h. Thirdly, the recovered RNA fragments were co-incubated with the mixture of magnetic beads and antibodies at 4 • C for 3 h. Fourthly, the above mixture was eluted with elution buffer, and the collected RNA was reverse-transcribed by HiScript reverse transcriptase. Lastly, the corresponding mixture of cDNA was configured according to AceQ SYBR ® qPCR Master Mix (Vazyme), and the expression levels of BCL2 and TLR4 were detected using a ViiA7 Real-Time PCR System instrument. The percentage input method was used for data analysis, where % input = 2 −(Average CTRIP − AverageCTinput − log2(input dilution factor)) . Related primer sequences were as follows: TLR4 F, CCGGCTGGTTTTGGGAGAAT and R, ATGGTCAGGTTGCACAGTCC; BCL2 F, CAGTTGCTCTGCTGTTTGAGG and R, CATTACTCTAGTGCTCCCCGC.

MeRIP-seq and mRNA-seq
The extracted RNA samples were sent to Cloud-Seq Biotech (Shanghai, China) for MeRIP-seq and mRNA-seq analyses (GSE 161050). Firstly, the fragmented RNA and m 6 A antibodies (Millipore, Burlington, MA, USA) were incubated in IP buffer at 4 • C for 2 h. Secondly, the mixture was immunoprecipitated with protein A magnetic beads (Thermo Fisher Scientific, Waltham, MA, USA) at 4 • C for 2 h. Thirdly, the RNA on the magnetic beads was eluted with a free m 6 A adenosine analogue and further extracted with TRIzol reagent. Lastly, the m 6 A IP samples and the input samples without immunoprecipitation were used for the NEBNext ® Ultra II Directional RNA Library Prep Kit (New England Biolabs, Inc., Ipswich, MA, USA), and double-ended sequencing was performed on an Illumina Hiseq (Illumina, Inc., San Diego, CA, USA).

Bioinformatics Analysis
Image analysis, base recognition, and quality control were performed after sequencing to produce raw data of original reads. Next, Q30 quality control and Cutadapt software (V1.9.3) (http://code.google.com/p/cutadapt/, accessed on 25 March 2020) were used to remove the connector and low-quality reads, thereby obtaining high-quality clean reads (Table S1, Supplementary Materials). The reads were compared on the basis of the genome/transcriptome (bosTau9) using Hisat2 software (v2.0.4) [53], and the MACS software (v1.4.2) [54] was used to identify methylated mRNA peaks in each sample. IGV software (v2.4.10) [55] was used to visualize the matching of MeRIP and input on the genome. The Gene Ontology (GO) (http://www.geneontology.org, accessed on 27 April 2020) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/keg, accessed on 27 April 2020) databases were used to analyze differentially methylated genes.

Statistical Analysis
Data analysis was performed using GraphPad Prism v7.0 software, and values are reported as the mean ± standard deviation (SD) of at least three replicates in each group. A p-value of <0.05 was considered statistically significant (*, **, and *** represent p < 0.05, p < 0.01, and p < 0.001, respectively).

Conclusions
Our findings clearly profile the m 6 A methylation of E. coli-induced MAC-T cells and identify many differentially methylated genes, which may contribute to the diagnosis and treatment of mastitis. These results are the first to elucidate a potentially meaningful relationship between m 6 A modification and dairy cow mastitis, laying the foundation for future studies of m 6 A modification in mastitis.