All procedures in this study were carried out in accordance with the Chinese Guidelines for Animal Welfare and Experimental Protocol and were approved by the Institutional Animal Care and Use Committee of Henan Institute of Science and Technology. Sixty fertilized Ross eggs were incubated in an automatic incubator with a temperature of 37°C ± 0.5°C and a relative humidity of 65%. Two weeks after hatching, chicks were randomly divided into Control and LPS group. Chicks in the LPS group were injected intraperitoneally with LPS from Escherichia coli 055:B5 (L4524, Sigma-Aldrich) at a dose of 1 mg/kg body weight. Chicks in the Control group were injected with the same volume of PBS. Two and 24 h after the treatment, blood samples were obtained and the plasma were separated, and livers were obtained and frozen in liquid nitrogen and stored at −80°C. Triglycerides and cholesterol content in plasma and liver were measured using the GPO-PAP method.

Total RNA was isolated by using TRIzol reagent (15596026, Invitrogen) according to the manufacturer’s introduction. By treating with DNase I (D2215, Takara), the potential genomic DNA contamination was eliminated. Two micrograms of DNase I treated RNA was used to obtained cDNA by using the GoScript Reverse Transcription System (A5001, Promega). cDNA was diluted (1:15) and proceed to real-time PCR analysis, and the TB Green Premix Ex Ta Ⅱ (RR820A, TaRaKa) was used for PCR amplification. The following standard PCR conditions were carried out for all transcripts: 95°C for 2 min; 35 cycles of 95°C for 15 s, 60°C for 30 s, and 72°C for 30 s; one cycle of 95°C for 15 s, 60°C for 15 s, and 95°C for 15 s. The last cycle was carried out to provide the melt curve for assessment of the specificity of amplification. Primers used in this study were listed in Supplementary Table S1, and chicken β-actin was selected as a reference gene. The method of 2^−ΔΔCt was used to analyze the data. Sequences of the primers used in real-time PCR were listed in Supplementary Table S1.

Two hundred micrograms of total RNA from the liver of both the Control and LPS groups were subjected to m⁶A quantification using the EpiQuik m⁶A RNA Methylation Quantification Kit (Colorimetric) (Epigentek, P-9005) following the manufacturer’s protocol.

Total protein was extracted from frozen liver with RIPA lysate. The protein concentration was then detected by using the Pierce BCA Protein Assay Kit (23227, Thermo Scientific). Equal amount of protein extracts was separated by electrophoresis on a 9% SDS-PAGE gel, and proteins were then transferred onto nitrocellulose membranes and immunoblotted with primary antibodies at 4°C for overnight. Primary antibodies used were as follows: anti-METTL3 (ProteinTech, 15073-1-AP); anti-METTL14 (ABclonal, A8530); anti-FTO (ProteinTech, 27226-1-AP); anti-YTHDF3 (AVIVA, ARP55530); and anti-β actin (Bioworld, LM1713). The membranes were then immunoblotted with the second antibodies at room temperature for 2 h. Reactive signals were detected by Pierce ECL Western Blotting Substrate (32109, Thermo Scientific). ECL signals were visualized with an imaging system (Bio-Rad, USA) and normalized to β actin by using ImageJ software.

For genome-wide m⁶A profiling, approximately 400 µg of total RNA was extracted from liver using TRIzol in accordance with the manufacturer’s instructions. RNA quality was controlled by using Nano Drop ND-2000 and agarose gel electrophoresis. High-throughput m⁶A and mRNA sequencing were performed by LC Sciences, Inc. (Hangzhou, China) following the published procedure. Briefly, mRNA was firstly purified by using oligo (dT)-coupled magnetic beads. Following purification, the poly(A) mRNA fractions were then fragmented into ∼100-nt-long oligonucleotides using Magnesium RNA Fragmentation Module (NEB, cat. e6150, USA) under 86°C for 7 min. Then, equal amount of cleaved RNA fragments (in biological duplicates) from Control and LPS-treated group were subjected to m⁶A immunoprecipitation by incubating with anti-m⁶A antibody (No. 202003, Synaptic Systems, Germany) for 2 h at 4°C in IP buffer (50 mM Tris-HCl, 750 mM NaCl, and 0.5% IGEPAL CA-630). Rabbit normal IgG was used as the negative control. An equal amount of fragmented mRNA from each biological sample was reserved as the input control (Input). The eluted m⁶A-containing fragments (IP) and the input were converted by SuperScript^TM Ⅱ Reverse Transcriptase (Invitrogen, cat. 1896649, United States) to obtain cDNA, which was then used to synthesize U-labeled second-stranded DNAs with E. coli DNA polymerase Ⅰ (NEB, cat. m0209, United States), RNaseH (NEB, cat. m0297, United States) and dUTP solution (Thermo Fisher, cat. R0133, United States). A-bases were then added to the blunt ends of each strand for the subsequent ligation to the indexed adapters. Each adapter contains a T-base overhang for ligating the adapter to the A-tailed fragmented DNA. Single- or dual-index adapters were ligated to the fragments, and size selection was performed with AMPureXP beads. After the heat-labile UDG enzyme (NEB, cat. m0280, United States) treatment of the U-labeled second-stranded DNAs, the ligated products are amplified with PCR by the following conditions: initial denaturation at 95°C for 3 min; eight cycles of denaturation at 98°C for 15 s, annealing at 60°C for 15 s, and extension at 72°C for 30 s; and then final extension at 72°C for 5 min. The average insert size for the final cDNA library was 300 ± 50 base pairs (bp). At last, a 2 × 150-bp paired-end sequencing was performed on an illumina Novaseq™ 6,000 (LC-BioTechnology Co., Ltd., Hangzhou, China) following the vendor’s recommended protocol.

Fastp software was used to remove low-quality bases, undetermined bases, and reads containing adaptor contamination. Hisat2 software was applied to align reads to the reference genome of chicken. Mapped reads of IP and Input libraries were provided for R package exomePeak for m⁶A peak identified. MEME and HOMER were used for de novo and known motif finding, followed by localization of the motif with respect to peak summit. Called peaks were annotated by intersection with gene architecture by using R package ChIPseeker. Then, StringTie was applied to determine the expression level for all mRNAs from Input libraries by calculating FPKM. Differentially expressed mRNAs were selected with log2 (fold change) > 1 or log2 (fold change) < −1 and p-value < 0.05 by R package edge R.

Data were presented as the mean ± SE. Statistical analysis was done using SPSS 22.0 and GraphPad Prism 8.0 softwares, Paired Student’s t-tests were performed between LPS-treated and control samples. Difference with p < 0.05 was defined as the threshold for significance.

The plasma triglycerides’ level 2 h after LPS treatment was increased compared to the control group (Supplementary Table S2), which indicated a rapid effect of LPS-induced inflammation on the lipid metabolism. Thus, the following tests were carried out on samples obtained 2 h after LPS treatment. Antibody-mediated capture of m⁶A and subsequent colorimetric analysis revealed that global m⁶A modification level in chicken liver was increased 2 h after the LPS treatment (Figure 1A). The expressions of m⁶A associated factors were then detected. The expression of the m⁶A writer METTL14 (p < 0.05) and reader YTHDF3 (p < 0.01) increased at the mRNA and protein levels after LPS treatment. Meanwhile, expression of FTO was increased at the protein level (p < 0.05) but not at the mRNA level.

MeRIP-Seq identified 7,815 m⁶A peaks representing 5,066 genes of chicken liver in the control group. In the LPS-treated group, 5,940 peaks were identified and represented the transcripts of 4,014 genes (Supplementary Table S3). According to the differences and overlaps in m⁶A modified genes, 321 genes were uniquely modified in the LPS-treated group, and 3,693 of them were commonly distributed in both groups. The maximum m⁶A number in all transcripts was 11, and the majority of genes (4,499/5,066) had one or two m⁶A modified sites (Figure 2D). The top 15 genes that contain the most m⁶A peaks were shown in Supplementary Table S4.

The m⁶A methylomes were further mapped by HOMER software. The top consensus motifs in the 7,815 peaks in the control group were AAACC and GGACG (Figure 2E). The distribution of the m⁶A peaks in the transcriptome of control and LPS group was also analyzed. Transcripts were divided into 5′UTR, start code, CDS, stop codon, and 3′UTR. m6A peak distribution showed similar patterns in control and LPS groups. In general, m6A peaks were highly distributed in the CDS and 3′UTR and less abundant in 5′UTR (Figure 2F). The m⁶A peak was most correlated with two regions: in the vicinity of stop codon and the start code region. Approximately 20% of CDS m⁶A peaks were distributed in the first exon, and the percentage of m⁶A peaks was stable along the transcript length but increases dramatically at the end of the CDS at threefold higher than at the beginning. In the 3′UTR, the peaks were enriched near the stop codon, then dramatically decreased, and maintained a low level along the length of 3′UTR (Figure 2G).

The comparison of the abundance of the m⁶A peaks between LPS group and control revealed 1,200 differentially methylated m⁶A peaks distributed in 1,039 genes (fold change > 2, p < 0.05) (Supplementary Table S3). A total of 1,101 hypo- and 99 hyper-methylated m⁶A peaks were found in the LPS group (Figure 3A). Analysis of abnormal m⁶A peaks distribution showed that the hyper-methylated sites were mostly in the CDS region, and the hypo-methylated sites preferentially in the region nearing the stop codon and the 3′UTR (Figure 3B).

GO and KEGG analyses of hypo- and hyper-methylated peaks were conducted using the DAVID 6.8 database (http://david.ncifcrf.gov/), and the results were considered statistically significant when p < 0.05. For hypo- and hyper-methylated genes, 32 and 9 significantly enriched GO items were identified, respectively. The top 10 hypo-methylated gene oncology items and all the hyper-methylated gene oncology items are shown in Figures 3C,D, respectively. The hyper-methylated genes were mostly enriched in the immune reaction-associated items, such as TLR2 signaling pathway, T cell proliferation, innate immunity, and TLR9 signaling pathway. The hypo-methylated genes were enriched in the items related to cell development and gene transcription. KEGG analysis of the hypo-methylated genes revealed eight pathways that were significantly enriched, namely, steroid biosynthesis, TGFβ signaling pathway, Notch signaling pathway, lysine degradation, mTOR signaling pathway, PPAR signaling pathway, insulin signaling pathway, and Wnt signaling pathway (Figure 3E). No significantly enriched pathway was identified when analyzing the hyper-methylated genes possibly due to the small number of genes.

In the RNA sequencing dataset (m⁶A-Seq Input library), 2,854 differently expressed genes were detected during the acute liver inflammation induced by LPS. Among them, 917 genes were upregulated and 1,948 genes were downregulated (fold change > 2, p < 0.05). Data scatter plot o was shown in Figure 4A. Gene ontology and pathway analysis revealed that the upregulated genes in the LPS-treated group were significantly enriched in biological process and pathways closely related to immune response (Figures 4B,E). This finding was consistent with the acute phase response of chicken under LPS stimulation. The downregulated genes were mostly enriched in metabolic process and DNA replication and transcription (Figures 4C,F). Moreover, the lipid metabolism-associated process was specially enriched, such as cholesterol biosynthetic process and fatty acid biosynthetic process. Real-time PCR validated the expression of several genes, and the results agreed with the sequencing data.

Correlation between the differentially methylated genes and the corresponding mRNA expression levels was analyzed by combining m⁶A-seq and RNA data. The m⁶A methylation level was significantly positively correlated with the gene expression level (Figure 5A). As most of the differently regulated peaks distributed in the CDS, STOP, and 3′UTR region, we analyzed the possible relations between different peaks in these regions and the corresponding mRNA expressions. The results showed that the methylation level in these regions was also positively related to the mRNA expression (Supplementary Figure S1). A total of 406 genes were differentially methylated and differentially expressed (p < 0.05, fold change >2). Most of these genes were hypo-methylated and downregulated (351/406), and many were hyper-methylated and up-expressed (47/406). Not one gene was hyper-methylated and down-expressed. Only eight genes showed hypo-methylation and up-expression (Figure 5B).

The 398 hypo-down or hyper-up genes were subjected to GO analysis to find their related biological process. As shown in Figure 5C, the hyper-up transcripts were involved in several immunity-related processes. This finding is consistent with the result in Figure 3D, showing that the hyper-methylated genes were correlated with the immune response under LPS treatment. Meanwhile, significant differences in enriched biological processes were observed between hypo-down and hypo-methylated gene (Figure 3C). Here, the process related osteoblast and osteoclast were remarkably enriched. Fatty acid biosynthetic process was also enriched (Figure 5C).

The prediction of protein–protein interaction (PPI) of the differentially methylated and expressed genes was carried out by using the STRING online database, and the results were downloaded and analyzed using Cytoscape software. The following top six hub genes (degree ≥ 10) and their corresponding first-stage node were screened: FEN1, PLK2, MIB2, AGPAT2, SREBF1, and ITGB2. The network was clustered into four subparts around the hub genes. Genes interacting with the hub gene ITGB2 (also called CD18), such as NLRC5, LY9, P2RX7, and PTRRJ, were mostly related to immune response. Genes interacting with SREBF1 were closely involved in the lipid metabolism. The other two subparts were largely associated with protein ubiquitination and DNA replication, respectively (Figure 5D). By using MCODE, several functional modules were screened, and modules enriched in lipid metabolism, DNA replication, and protein ubiquitination were re-emerged (Figure 5E). These results indicated that the m⁶A modification may participate in these biological processes in the acute stage after LPS stimulation.