The study was approved by the Ethics Committee of the Shengjing Hospital of China Medical University (Shenyang, China; approval number 2022PS854K). This study conformed to the relevant ethical standards. The rats used in this study (Beijing HFK Bioscience Co., Ltd.; Beijing, China) had a Wistar genetic background. The rats were adolescent, male, and pathogen-free; weighed approximately 160-180 g; and were housed in cages at 24°C with a 12h alternating light/dark cycle and free access to water and food. An adolescent rat model of endotoxic shock was generated by intraperitoneal (i.p.) injection of LPS, as we described previously (6). Briefly, rats were anesthetized with 20% urethane (1g/kg, i.p.). The left femoral artery was cannulated (Biopac MP150; Biopac Systems, Goleta, CA, USA) to continuously measure the mean arterial pressure (MAP). After the MAP stabilized, LPS (L-2880; Sigma-Aldrich, St. Louis, MO, USA) was administered i.p. (20 mg/kg, 10 mg LPS dissolved in 1 ml of 0.9% saline). 2 hours after LPS injection, MAP began to decline, and the model was considered successful when MAP decreased by 20-25 percent. The control group was injected with an equal volume of 0.9% normal saline solution. The preliminary experiment found that the death rate of rats was relatively high 24 hours after injection of LPS, but could basically be stabilized to 12 hours. In order to obtain the pathological changes in the late stage of LPS-induced myocardial injury, we selected 12 hours after injection of LPS to conduct this study. Twelve hours post-LPS or normal saline injection, euthanasia was performed, and the left ventricle was excised and placed in an RNase-free centrifuge tube, quickly placed in liquid nitrogen, and stored at -80°C. Blood samples were collected from the abdominal aorta at 12h post-LPS or saline administration and centrifuged (3,000 rpm, 10 min) after incubation at room temperature (15°C-25°C) for 15 min. Supernatants were collected and stored at -80°C for enzyme-linked immunosorbent assay (ELISA) analysis. Six pairs of left ventricles with endotoxaemia or controls were selected for MeRIP-seq and RNA-seq, and the remaining samples were saved for validation experiments.

Rat heart samples were fixed in 4% paraformaldehyde for 48 h, dehydrated, permeabilized, and embedded in paraffin. Paraffin-embedded tissues were sectioned (0.35 µm) and stained with H&E. Pathological changes were observed microscopically.

Twelve hours after the LPS injection, echocardiography was performed. A special ultrasound instrument (Vinno Technology, China) for small animals was used to obtain echocardiographic images. Two-dimensional guided M-mode measurements of the left ventricular (LV) internal diameter were obtained from the long-axis view at the level of the papillary muscles. The interventricular septal thickness at diastole (IVSd), left ventricular intra-diameter at diastole (LVIDd), left ventricular end-diastolic posterior wall thickness (LVPWd), interventricular septal thickness at systole (IVSs), left ventricular intra-diameter at systole (LVIDs), left ventricular end-systolic posterior wall thickness (LVPWs), left ventricular end-systolic volume (LVESV), left ventricular end-diastolic volume (LVEDV), left ventricular ejection fraction (LVEF), stroke volume (SV), left ventricular fractional shortening (LVFS), heart rate (HR) and cardiac output (CO) were recorded.

Cardiac troponin I (c-TNI), creatine kinase myocardial band (CK-MB) and interleukin (IL)-6 levels in blood supernatants from both groups were determined via ELISA analysis. IL-6 ELISA kits were purchased from Wanlei Biological Company (Shenyang, China), whereas CK-MB and c-TNI ELISA kits were purchased from Elabscience Biotechnology Co., Ltd. (Wuhan, China). The biomarker concentrations in each sample were estimated based on optical density values at 450 nm and a standard curve.

Total RNA was isolated from heart tissue using the TRIzol Reagent (Invitrogen, Life Technologies, Carlsbad, CA, USA). The RNA quantity and purity of each sample were determined using Nanodrop 2000 (Thermo Scientific, USA). m⁶A RNA Methylation Quantification Kit (ab185912, Abcam, Shanghai, China) was used to detect the total m⁶A levels in 200ng RNA by measuring absorbance at 450nm.

RNA extraction and purity measurement are the same as 2.5. RNA-integrity numbers were assessed using a Bioanalyzer 2100 (Agilent, CA, USA) and exceeded 7.0, as confirmed by denaturing agarose-gel electrophoresis. Poly (A) RNA was purified from 50μg total RNA using Dynabeads Oligo (dT)25-61005 (Thermo Fisher, CA, USA) and two rounds of purification. Then, the poly(A) RNA was fragmented using a Magnesium RNA Fragmentation Module (NEB, USA) at 86°C for 7min. The RNA fragments were incubated for 2h at 4°C with an m⁶A-specific antibody (Synaptic Systems, Germany) in IP buffer (50mM Tris-HCl, 750mM NaCl, and 0.5% Igepal CA-630). Both the input samples (without IP) and the m⁶A-IP samples were used to construct RNA-seq libraries with the NEBNext Ultra RNA Library Prep Kit for Illumina. Finally, we performed 2×150bp paired-end sequencing (PE150) on an Illumina Novaseq 6000 instrument (LC-Bio Technology Co., Ltd., Hangzhou, China).

Fastp software (12) (default parameters) was used to remove reads containing adaptors, low-quality bases, and undetermined bases. The quality of the IP and input sequences was verified using FastQC and RseQC. Filtered high-quality data were mapped to the Rattus norvegicus genome (Ensemble release 101) using HISAT2 (13). Peak calling and differential-peak analysis were performed using the exomePeak2 package (14), and peaks were annotated using ANNOVAR. The peaks in each group with a fold change (FC) of ≥2 was considered significant. HOMER (15) was used for de novo motif finding, followed by motif localization with respect to the peak summit. Then, StringTie (16) was used to determine the expression levels of all mRNAs from the input libraries using the following relationship: (total exon fragments/mapped reads [millions] × exon length [kB]).

Genes with differential expression or m⁶A modification sites were imported into the STRING database (17), and the interaction relationships between the imported genes were determined. A PPI network was constructed by importing the data into Cytoscape 3.8.2 (18)and analyzing the data using a network analyzer. The genes imported with combined interaction scores of >0.4 were selected to construct the PPI network.

MeRIP assays were performed using the m⁶A RNA Enrichment Kit (Epigentek, USA). Target m⁶A-containing fragments were immunoprecipitated using a bead-bound m⁶A capture antibody, and RNA sequences at both ends of the m⁶A-containing regions were cleaved using Cleavage Enzyme Mix. The enriched RNA was released, purified, and eluted. Quantitative reverse transcriptase-PCR (qRT-PCR) analysis was performed following MeRIP to quantify changes in target-gene m⁶A modification. Primers ( Table S1 ) were synthesized by Sangon Biotech Co., Ltd. (Shanghai, China).

RNA was reverse transcribed to complementary DNA using the PrimeScript™ RT Reagent Kit with gDNA Eraser (Takara Bio, Being, China), with a reaction time of 15min at 37°C, followed by denaturation for 5s at 85°C. A 7500 Real-Time PCR System (Applied Biosystems, USA) was used for PCR. The primers used are listed in Table S2 .

Frozen ventricular tissue samples were homogenized in RIPA buffer (Beyotime Biotechnology, China) containing protein enzyme inhibitors, and protein concentration was quantified using a BCA Protein Assay Kit (Beyotime Biotechnology, China). The lysates were separated on 8–15% SDS-PAGE gels and then electro-transferred to 0.45μm polyvinylidene difluoride membranes (Millipore, USA). After blocking with 5% nonfat milk in Tris-buffered saline at room temperature for 2h, the membranes were incubated for overnight at 4°C with primary antibodies. The primary antibody included IL-6 (Wanlei, China), IL-1β (Affinity, China), IL-18 (Affinity, China) and Tubulin (Affinity, China). After washing three times for 10min per wash, the membranes were incubated with an HRP-conjugated secondary antibody for 2h at room temperature. Subsequently, the blots were imaged using a Bio-Rad imaging system (Bio-Rad, USA) and protein expression levels were determined using ImageJ software.

All data were obtained from three or more independent experiments and are presented as mean ± standard deviation. Statistical analyses were conducted using GraphPad Prism (version 8.0). Student’s t-test for unpaired data was used to compare two groups. Statistical significance was set at P < 0.05. Peaks and transcripts with |FC| ≥1.5 and P < 0.05 between control group and LPS group were considered differential peaks and transcripts.

In this study, we conducted a rat model of endotoxic shock by intraperitoneal injection of LPS. The MAP in rats decreased by 25-30% 2h after LPS administration ( Supplementary Figure 1 ). Firstly, we assessed cardiac function using echocardiograph. As shown in Table 1 , the IVSd, LVIDd, LVIDs, LVEDV, LVESV, SV, HR and CO were deceased significantly in the LPS group while the LVPWd and LVSs were increased significantly in the LPS group. However, no difference in LVPWs, EF and FS were observed between two groups. Then, we detected the pathological changes in the cardiac tissue founding that myocardial cells showed greater disorder and inflammatory cell infiltration in the LPS group than in the control group (H&E staining; Figure 1A ). Strikingly, serum c-TNI and CK-MB (markers of injury to cardiomyocytes) were significantly elevated in the LPS group ( Figure 1B, C ), as was the serum inflammatory biomarker, IL-6 ( Figure 1D ). We also analyzed the protein expression levels of inflammatory cytokines (IL-6, IL-18, IL-1β) in the heart using western blot. The protein expression levels of IL-6, IL-18, and IL-1β were upregulated in the LPS group compared to the control group ( Figure 1E ), suggesting inflammatory injury of cardiac tissue in the rat model induced by LPS.

In addition, higher total m⁶A modification levels were found in the LPS group than that in the control group ( Figure 1F ). Thus, we performed transcriptome-wide MeRIP-seq and RNA-Seq to generate an m⁶A modification map in cardiac tissue of rats with LPS-induced myocardial injury.

Next, we compared the m⁶A peaks between LPS group and control group. We identified 14,360 common m⁶A peaks, representing 9,840 genes in both groups using exomePeak2 ( Figure 2A ). In addition, 4,497 unique peaks and 3,589 unique genes appeared in the LPS group, whereas 5,162 unique peaks and 4,102 unique genes appeared in the control group ( Figure 2B ). Chromosome 1 showed the largest number of m⁶A modification sites with 2,272 and 2,355 m⁶A peaks in the LPS group and the control group, respectively, followed by chromosome 10 with 1,613 and 1,662 m⁶A peaks in the LPS group and the control group, respectively ( Supplementary Figure 2A ). Most genes had one to three m⁶A modification sites ( Supplementary Figure 2B ).

In addition, we identified the top five consensus motifs for the m⁶A peaks in the LPS group and the control group ( Figure 2C ). RRACH motif was identified in both groups, where R = G or A; A = m⁶A, and H = U, A, or C, as described (19). The m⁶A peaks were predominantly distributed in coding sequences and the 3’ untranslated region (UTR) ( Supplementary Figure 2C ). Finally, m⁶A peaks were analyzed in whole transcriptome data and divided into exon, intergenic, intron, 3’UTR, 5’UTR regions, based on their location in RNA transcripts. The peaks were mainly distributed in 3’UTR, exon and 5’UTR and different distribution patterns were observed between the two groups. The LPS group showed slightly more m⁶A peaks in the exon (27.97% versus. 28.15%) and 3′UTR (56.24% versus. 56.62%) and slightly less m⁶A peaks in the 5′UTR (15.56% versus. 14.97%) than the control group ( Figure 2D ).

Furthermore, we found 14,360 m⁶A peaks in the LPS group and the control group, 5,740 of which were differentially modification (P < 0.05). Among them, 512 hypermethylated and 658 hypomethylated m⁶A modification sites (|FC| ≥ 1.5, P < 0.05) were found in the LPS group versus. the control group. The 10 genes with the most significant differences were identified ( Figure 2E ). Table 2 showed the top 20 differentially upregulated and downregulated m⁶A peaks.

To determine the possible functions of the differentially m⁶A-methylated genes involved in regulation, we performed Gene Ontology (GO) (20), Kyoto Encyclopedia of Genes and Genomes (KEGG) (21), and PPI network analyses.

Most genes with differential m⁶A modification were related to apoptosis and inflammatory responses by GO analysis ( Figures 3A, D ). Furthermore, KEGG pathway analysis of genes with differentially m⁶A modification revealed enrichment for terms related to (i) inflammatory responses, such as Staphylococcus aureus infection, Fc gamma R-mediated phagocytosis, complement and coagulation cascades, Th17 cell differentiation, chemokine signaling pathway, IL17 signaling pathway and antigen processing and presentation; (ii) important mediators of apoptosis, such as Hippo signaling pathway, AGE-RAGE signaling pathway in diabetic complications, MAPK signaling pathway and PI3K/Akt signaling pathway ( Figures 3B, E ).

We then performed PPI-network analysis of genes with differential m⁶A modification using Cytoscape (18). Stat3, Hif1a, and Creb1 were the most central genes with upregulated m⁶A modification, whereas Tp53, Mapk11, and Hsp90aa1 were the most central genes with downregulated m⁶A modification ( Figures 3C, F ). These central proteins are particularly associated with inflammatory responses and apoptosis.

RNA-seq data showed that 1,836 genes were significantly upregulated and that 2,467 genes were significantly downregulated in the LPS group compared to the control group (|FC| ≥1.5, P < 0.05). Figure 4A shows the 10 most upregulated and downregulated genes. In addition, principal component analysis showed that samples from the LPS group and control group clustered separately within each group ( Supplementary Figure 3A ). The heatmap of Pearson correlation coefficient matrix also revealed good intragroup correlations in both groups ( Supplementary Figure 3B ). Moreover, GO analysis suggested that differentially expressed genes were associated with inflammatory responses and apoptosis ( Figure 4B ). KEGG pathway analysis revealed enrichment for terms related to (i) inflammatory responses terms, such as the chemokine signaling pathway, Toll-like receptor signaling pathway, Th17 cell differentiation, natural killer cell-mediated cytotoxicity, NOD-like receptor signaling pathway, cytokine-cytokine receptor interaction and TNF signaling pathway; (ii) apoptosis terms such as PI3K-Akt signaling pathway and MAPK signaling pathway ( Figure 4C ). Furthermore, PPI-network analysis demonstrated that the central proteins were Stat1, Tnf, and IL-6 among the upregulated genes (which are critical for inflammatory responses), and Erbb3, Igf1, and Pik3r3 among the downregulated genes (which are vital for involvement of the PI3k-Akt signaling pathway) ( Figures 4D, E ).

Seventy-two genes with m⁶A hypermodification and 103 genes with m⁶A hypomodification were upregulated. In addition, 137 genes with m⁶A hypermodification and 133 genes with m⁶A hypomodification were downregulated. Notably, 4 upregulated genes had both m⁶A upregulated and downregulated modification at different sites, and 9 downregulated genes had both m⁶A upregulated and downregulated modification sites ( Figure 5A ).

GO analysis suggested that genes with both differential m⁶A modification and differential expression were related to inflammatory responses and apoptosis ( Figure 5B ). KEGG terms related to (i) certain immune pathways, such as Staphylococcus aureus sepsis, Th17 cell differentiation, natural killer cell-mediated cytotoxicity, cell adhesion molecules and antigen processing and presentation and (ii) several pathways involved in apoptosis, such as PI3K-Akt signaling pathway and MAPK signaling pathway, were significantly enriched ( Figure 5C ). In addition, Tp53, stat3, Mst1r, and Hsp90aa1 were the most central proteins in the PPI network associated with immune response and apoptosis ( Figure 5D ).

Lasso screening identified Utrn, Stat3, Stard9, Pcdh12, LOC108348047, and Dnmt3a as key m⁶A-hypermodified genes ( Figure 6A ), whereas Tnfrsf9, Sele, Pcdh12, Myrip, Gng2, Gng12, Bdkrb2, Batf, and Agrn were identified as key m⁶A-hypomodified genes in the LPS group compared to the control group ( Figure 6B ). Furthermore, we analyzed the m⁶A modification sites with the Integrative Genomics Viewer (IGV) (22) ( Figure 6C ), verified the m⁶A modification and gene expression levels of leucine-rich repeat kinase 2 (Lrrk2), selectin E (Sele), TNF receptor superfamily member 9 (Tnfrsf9), bradykinin receptor B2 (Bdkrb2), and heat shock protein 90 alpha family class A member 1 (Hsp90aa1) in the LPS group and the control group ( Figures 6D, E ), which showed that they had the same expression tendencies, consistent with our sequencing data.

We examined the expression of m⁶A methyltransferase (“writers”), including Mettl3, Mettl14, Wtap, Mettl16, Rbm15, Rbm15b, Zc3h13 and Virma and found that Mettl16 was upregulated whereas Rbm15 was significantly downregulated in the LPS group versus. the control group ( Figure 7A ). In addition, among the m⁶A demethylases (“erasers”) studied (Fto and Alkbh5), Fto was significantly downregulated in the LPS group ( Figure 7B ). When studying m⁶A-recognition factors (“readers”), we discovered that the YTH domain family protein, Ythdc2, was significantly upregulated in the LPS group versus. the control group and that heterogeneous nuclear ribonucleoprotein G (Hnrnpg) was downregulated ( Figure 7C ).