MAC-T cells (an immortalized bovine mammary epithelial cell line) were donated by Professor Mark Hanigan of Virginia Tech University. MAC-T cells were cultured by adapting the procedure used by Li et al. (40).

S. aureus (ATCC 29213) and E. coli (ATCC 25922) were donated by Professor Zhou Rui and Professor Wang Xiangru of Huazhong Agricultural University. S. aureus and E. coli were resuscitated, single colony purified, cultured, and heat-inactivated by adapting the procedure used by Li et al. (40).

MAC-T cells were seeded in 6-well plates (Corning, US) at a density of 2×10⁵ cells/well. After 12 h of culture, inactivated bacteria were added at an MOI of 10:1. Inactivated bacteria and cells were cocultured for 24 h. Cells in each well were washed 3 times with cold PBS (HyClone, China) before TRIzol (Invitrogen, US) was added. There were 3 replicate wells for each treatment and 3 replicate wells for the control, for a total of 9 samples.

Total RNA of the above 9 samples was extracted according to the instructions of the commercial reagent manufacturer and the RNA isolation steps. A NanoDrop ND-1000 instrument (Agilent Inc. USA) was used to measure RNA concentration and purity, and the quality control index of RNA purity was based on the OD260/OD280 value between 1.8 and 2.0. Denaturing agarose gel electrophoresis was used to measure RNA integrity and gDNA contamination. The sample was stored at −80°C for later use.

The RNA samples obtained were reverse transcribed into cDNA using Vazyme HiScript II Reverse Transcriptase (+gDNA wiper) (Vazyme, Nanjing, China), and the expression of related RNA was detected by Vazyme AceQ^® SYBR qPCR Master Mix (Vazyme) in a ViiA7 Real-time PCR System (Applied Biosystems Inc., Foster City, CA, United States). The 2^−△△ct method was used to analyze the fluorescence quantitative data, and GraphPad Prism 7.0 was used to process the data. Related primer sequences are shown at Table S1 .

After MAC-T cells were digested, collected, and washed three times with PBS, the cells were resuspended in 100 µl of Binding Buffer. Then, 5 µl each of FITC and PI dyes (Vazyme, China) was 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).

MAC-T cells were treated with S. aureus and E. coli and fixed with 2% glutaraldehyde fixation buffer. The cells were used for TEM slide preparation. Sections were observed via a 100 KV H7650 transmission electron microscope (HITACHI, Japan).

Nine samples under different treatments were sent to Cloud-Seq Biotech (Shanghai, China) for MeRIP-CircRNA sequencing. The operation was briefly described as follows: (1) To enrich circRNA, a circRNA Enrichment Kit (Cloud-seq, USA) was used according to the supplier’s instructions. (2) An NEBNext^® Ultra™ II Directional RNA Library Prep Kit (New England Biolabs, Inc., Massachusetts, USA) was used to pretreat RNA, and a sequencing library was constructed by immunoprecipitation. (3) Quality control and quantification of the library were performed with the BioAnalyzer 2100 system (Agilent Technologies, USA), and the Illumina HiSeq instrument was used to perform 150 bp paired-end sequencing.

After paired-end sequencing, raw circRNA data were obtained. According to Q30 for quality control, cutadapt software (v1.9.3) was used to obtain high-quality reads. STAR software (v2.5.1b) was used to compare high-quality reads to the reference genome/transcriptome (bosTau9), and DCC software (v0.4.4) was used for circRNA detection and identification. The circBase database and Circ2Traits were used to annotate the identified circRNA. Then, edgeR software (v3.16.5) was used for data standardization and differential expression circRNA screening (fold change ≥ 1.5, p-value ≤ 0.05). The standardized number of reads was used to calculate the differential expression of circRNA between the three groups of samples. A fold change ≥ 2.0 and p-value ≤ 0.05 were considered the thresholds for differential circRNA.

According to the circRNA verified by sequencing, the differentially expressed circRNA was selected. Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/keg, accessed on 27 April 2020) and Gene Ontology (GO; http://www.geneontology.org, accessed on 27 April 2020) enrichment analyses for the circRNA source gene were performed using the DAVID biometric analysis tool. The downstream mRNAs were predicted using miRBase (http://mirbase.org/index.shtml).

Statistical analysis of the data was performed using GraphPad Prism 7.0 (GraphPad Software, La Jolla, CA, United States). Data are represented as the mean with SD, and the significant differences were analyzed by Student’s t-test; p < 0.05 was considered statistically significant.

Additionally, RNAhybrid (2.1.2) was used to predict circRNA and miRNA interaction (41). Mirbase 22 (https://mirbase.org/, accessed on March 3, 2018) was used to predict the downstream mRNA of the predicted miRNA.

To determine whether the stimulation of inactivated bacteria on MAC-T cells can cause cell damage, we added bacteria to stimulate the cells for 24 h. The results showed that IL-1β, IL-6, and TNF-α were significantly upregulated in the S. aureus group and E. coli group, and the expression of these three genes in the E. coli group was significantly higher than that in the S. aureus group ( Figures 1A–C ). Additionally, cell apoptosis was significantly increased ( Figures 1D–F ). It showed almost 5 times the apoptosis rate in the S. aureus group and E. coli group compared to the control group. In contrast to the trend in mRNA expression associated with inflammation, flow cytometry illustrated a similar apoptosis rate when cells were treated with S. aureus and E. coli ( Figure 1G ).

To more intuitively observe the cell injury in S. aureus treatment and E. coli treatment, transmission electron microscopy was used ( Figures 1H–M ). Interestingly, we found the nuclear membrane disappearance, chromatin condensation ( Figure 1K ), and autophagosomes ( Figure 1L ) in MAC-T cells injured by S. aureus ( Figure 1I ), which indicated minor damage in these cells. However, cells with acute cytosolic voids ( Figure 1J ) and apoptotic bodies ( Figure 1M ) were found in MAC-T cells injured by E. coli.

To investigate the circRNA profile in MAC-T cells injured by heat-inactivated S. aureus or E. coli, purified cellular RNA was subjected to circRNA-seq and MeRIP-seq. The sequencing raw reads were generated from the control group, S. aureus group, and E. coli group. Bioinformatic approaches, including Cutadapt and Bowtie2 software, were performed for data filtering and quality control. We obtained clean reads and the proportions of net reads were between 99.122% and 99.995% with stringent quality control ( Table S2 ).

To clarify the relationship between m⁶A modification and normal circRNA, we first analyzed the circRNA profile from circRNA-seq and found that most circRNAs were less than 2,000 bp and derived from exons ( Figures S1A, B ). Moreover, genes were mainly distributed on chromosome 3 ( Figure S1E ). When MAC-T cells were injured by S. aureus and E. coli, the source and length of the differential circRNA showed similar patterns ( Figures S1C, D, F ). Most circRNAs were located on chromosome 3 in the S. aureus group, while most circRNAs were located on chromosome 8 in the E. coli group ( Figures S1F, G ).

Further analysis of the MeRIP-seq data was performed, and the m⁶A methylation peaks were compared and analyzed based on the sequencing data. There were 1,599 m⁶A methylation peaks within 1,035 circRNAs in the control group, 35 peaks within 32 circRNAs in the S. aureus group, and 1,016 within 728 circRNAs in the E. coli group ( Figures 2A, B ).

A total of 17 m⁶A peaks and 18 circRNA transcripts were observed in all three groups, whereby 933 methylation peaks and 687 circRNAs occurred in both the E. coli group and the control group. The number of m⁶A methylation peaks was obviously different among the three groups. There were 18 m⁶A peaks in the S. aureus group instead of the control group and 83 m⁶A peaks in the E. coli group instead of the control group.

Further analysis was performed to assess the features of m⁶A-modified circRNAs. The number of m⁶A methylation peaks in each circRNA was highly similar in the control group, S. aureus group, and E. coli group ( Figure 2C ). Similar to recent studies, we have found that most circRNAs have one m⁶A methylation peak, but there are also a large number of circRNAs with more than one m⁶A methylation peak, which indicates that m⁶A modification sites are not unique in circRNAs. Moreover, the lengths of m⁶A-modified circRNAs in each group were analyzed. The length of most m⁶A-modified circRNAs was less than 5,000 bp, and similar characteristics of length were shown in all three groups ( Figure 2D ). In our studies, the sources of the circRNA in the three groups were subsequently analyzed. CircRNAs in the control group and E. coli group derived from exons were the most abundant. However, in the S. aureus group, there was a certain amount of m⁶A-modified circRNA derived from antisense and sense overlapping ( Figure 2E ). Finally, chromosome distribution also revealed that m⁶A-methylated circRNA is more likely to be present on chromosome 3 ( Figure 2F ).

To identify the role of m⁶A-modified circRNA in regulating MAC-T cell injury induced by S. aureus and E. coli, the various m⁶A modification levels between S. aureus/E. coli and the control group were analyzed (fold change > 2, p-value < 0.05). Compared to the control group, a total of 259 differential m⁶A modification peaks within 200 circRNAs were found. Further analysis confirmed that 67 m⁶A methylation peaks within 63 circRNAs were significantly different between the S. aureus group and the control group, of which 26 hypermethylated peaks were within 25 circRNAs (e.g., circRNA_OPA3), and 41 hypomethylated peaks were within 38 circRNAs (e.g., circRNA_AKT3). Meanwhile, 192 m⁶A methylation peaks within 137 circRNAs were significantly different between the E. coli group and the control group, of which 75 hypermethylated peaks were within 47 circRNAs (e.g., circRNA_GRB10), and 117 hypomethylated peaks were within 90 circRNAs (e.g., circRNA_SMAD1). Data visualization analysis was performed by IGV to show the differential m⁶A peaks between the control group and treatment group ( Figure 3A ). The top 20 differentially methylated circRNAs with hypermethylation or hypomethylation in the S. aureus and E. coli group compared to the control group are shown ( Tables 1 , 2 ).

Additionally, the feature of the differentially m⁶A-methylated circRNA in the S. aureus/E. coli group compared to the control group is also shown. Instead of being a major part of short m⁶A-modified circRNA derived from exon ( Figure 2D ), it is longer and derived more from sense overlapping when MAC-T cells are injured by S. aureus and E. coli ( Figures 3B, C ). Additionally, the differential circRNA with methylation showed totally reverse patterns in the S. aureus group compared to the E. coli group, which were mainly hypomethylation in the S. aureus group and hypermethylation in the E. coli group. The most differentially methylated circRNA in the S. aureus group was found on chromosome 10, while in the E. coli group, it was on chromosome 19 ( Figure 3D ).

Recent studies have shown that circRNAs can affect the expression of its cis genes (42). Therefore, to clarify the essential role of m⁶A-modified circRNAs in MAC-T cell injury induced by S. aureus and E. coli, GO and KEGG pathway enrichment analyses for the m⁶A-modified circRNA source genes were performed (p-value < 0.05) ( Tables S3, S4 ). In the S. aureus group, the GO annotation of circRNA with hypermethylation illustrated that it was mainly enriched in transmembrane, cell part morphogenesis, plasma membrane-bounded cell projection organization, P-body, and protein tyrosine kinase activator activity ( Figure 4A ). Furthermore, the KEGG analyses indicated that they were enriched in endocytosis, antigen processing and presentation, and ubiquitin-mediated proteolysis ( Figure 5A ). Likewise, the GO analysis of circRNA with hypomethylation showed that it was closely related to clathrin, SMAD binding, integrin binding, and PI3K ( Figure 4B ). KEGG analysis of hypomethylated circRNAs revealed enrichment in focal adhesion, the RAP1 signaling pathway, and regulation of the actin cytoskeleton ( Figure 5B ).

In the E. coli group, the GO analysis predicted that the functions of circRNAs with hypermethylation were significantly enriched in cytoskeleton, and actin filament ( Figure 4C ). The KEGG pathway analysis indicated that circRNAs with hypermethylation were significantly related to the Ras signaling pathway, endocytosis, and PI3K-Akt signaling pathway ( Figure 5C ). Correspondingly, the GO analysis of differential m⁶A-modified circRNA with hypomethylation illustrated that the meaningful terms may be related to regulation of the retinoic acid receptor signaling pathway, SMAD binding, Ras GTPase binding, and small GTPase binding ( Figure 4D ). The KEGG analysis also revealed enrichment in focal adhesion, regulation of actin cytoskeleton, bacterial invasion of epithelial cells, and Rap1 signaling pathway ( Figure 5D ).

Also, some representative genes closely related to apoptosis and inflammation in the KEGG pathways analysis were shown, such as BOLA, AKT3, MAP2K4, and JAK2 ( Figures 5A–D ). Given the above, we found that the functions of differentially expressed m⁶A-modified circRNAs were highly similar when MAC-T cell injury was induced by either S. aureus or E. coli. These findings may suggest that these m⁶A-modified circRNAs regulated immunity, cell junctions, growth metabolism, and resistance to bacterial invasion through a similar pathway in both induced groups.

CircRNAs are able to regulate the expression of target genes as sponges for miRNAs based on complementary base pairing. By predicting the target miRNA between circRNA and mRNA, we constructed a circRNA–miRNA–mRNA interaction network ( Table S5 ). In this study, 10 circRNAs (5 hypermethylation and 5 hypomethylation), 10 miRNAs, and 22 mRNAs were included in the S. aureus group ( Figure 6A ), and 10 circRNAs (5 hypermethylation and 5 hypomethylation), 8 miRNAs, and 19 mRNAs were included in the E. coli group ( Figure 6B ).

In the interactive network of the S. aureus group, many inflammation-related miRNAs were predicted. For example, miR-2305 and miR-1777a, which have been considered to have an effect on immunity (43), were also shown to enrich for transcripts related to inflammation and apoptosis in our study (including key genes such as MAPK8IP2 and WNT10A). Additionally, we predicted that some miRNAs may be regulated by multiple circRNAs at the same time. For instance, miR-11988 was predicted to be regulated by circ_FAM214A, circ_NEDD4 L, and circ_SETD2 simultaneously.

Similarly, in the E. coli group, a large amount of mRNA was also predicted to be associated with inflammation, mainly in the MAPK signaling pathway and the WNT signaling pathway. In addition, BCL2, which is the key molecule in apoptosis, was predicted to be regulated by miR-2305. Some natural immune-related mRNAs such as TNFAIP8L2, FCER1G, and MSR1B were present in the network. It is predicted that the miRNAs may regulate multiple mRNAs. For example, miR-2305 could be able to interact with BCL2, FCER1G, and MAX. Interestingly, in both groups, we found a large number of similar miRNAs, suggesting that circRNAs modified with different m⁶A modifications may play a similar role in S. aureus- and E. coli-induced cell injury.

To further explore the potential function of circRNAs with m⁶A modification in S. aureus- and E. coli-induced cell injury, a conjoint analysis of circRNA-seq and MeRIP-seq was performed.

In addition, the differential expression of circRNAs with significant changes in the S. aureus and the E. coli groups compared to the control group was analyzed ( Figures 7A–D ). We found 3,914 novel circRNAs, of which 1,532 were upregulated and 11 were downregulated in the S. aureus group, and 536 were upregulated and 11 were downregulated in the E. coli group. The top 20 most significantly differentially expressed circRNAs in the S. aureus group and E. coli group are listed separately ( Tables S6, S7 ).

To predict the especially important circRNA in MAC-T injury and reveal the relationship between m⁶A modification and expression, a conjoint analysis was performed. Finally, 3 significantly upregulated circRNAs with hypermethylation and 4 upregulated circRNAs with hypomethylation were found in the S. aureus group ( Figure 7E ). Two upregulated circRNA with hypomethylation were found in the E. coli group ( Figure 7F ).