Our study is based on the classic OVA model of allergic asthma. The schematic design of the animal modeling is shown in Figure 1A . Airway hyperresponsiveness is the most typical pathological feature of allergic asthma. To explore the ability of OVA to effectively induce airway hyperresponsiveness in mice, we tested the pulmonary function of mice. In the OVA-induced asthma mouse model (hereafter abbreviated as OVA group), lung resistance was increased ( Figure 1B ) and the dynamic compliance (C _dyn) ( Figure 1C ) was significantly decreased at 3.125, 12.5, and 50 mg/ml methacholine compared with those of phosphate-buffered saline (PBS)-treated control mice (hereafter abbreviated as PBS group).

Th2-type immune response is crucial for triggering the development of airway hyperresponsiveness. Therefore, we next examined airway inflammation in the OVA group using ELISA, hematoxylin and eosin (HE) staining, and leukocyte classification and counting. In ELISA ( Figure 2A ), the levels of IgE in serum and IL-4, IL-5, and IL-13 in bronchoalveolar lavage fluid (BALF) were measured in mice in the OVA group at significantly higher concentrations than that in PBS control mice. Meanwhile, we classified and counted the leukocytes in BALF, and the results indicated distinctly increased numbers of total leukocytes (Total), neutrophils (Neu), lymphocytes (Lym), monocytes (Mon), eosinophils (Eos), and basophils (Bas) ( Figure 2B ) and rising percentages of Lym, Eos, and Bas ( Figure 2C ) in the OVA group compared with those in the PBS group. Furthermore, HE staining of lung tissue sections showed an increased inflammatory cell infiltration around the airway and blood vessels in OVA-induced mice ( Figure 2D ), and the corresponding histological inflammation score ( Figure 2E ) was higher in the OVA group compared with that of the PBS group.

Additionally, we performed periodic acid–Schiff (PAS) and Masson’s trichrome (MT) staining on the lung sections. PAS staining and positive scores revealed increased airway mucus secretion in the OVA group, which was implicated in goblet cell hyperplasia ( Figures 2F, G ), and the results of MT staining showed that the increase in collagen fiber deposition in the airways of mice in the OVA group confirmed the change of the airway structure and the occurrence of airway remodeling ( Figures 2H, I ).

After six MeRIP-seq libraries and six RNA-seq libraries were sequenced, an average of 92,112,289 raw reads (13.82 GB data) from MeRIP-seq and 79,562,999 raw reads (11.93 GB data) from RNA-seq were obtained ( Supplementary Table S1 ). Then, the reads that contained adaptor contamination, low-quality bases, and undetermined bases were filtered out and clean reads were obtained. Specifically, we collected an average of 89,544,192 valid MeRIP-seq reads (12.31 GB data) and 77,517,266 valid RNA-seq reads (10.58 GB data) ( Supplementary Table S1 ). Further analysis found that an average of 77,085,259 (93.05%) clean reads were mapped to the mouse reference genome. Of these, an average of 11,444,146 (13.84%) were multi-mapped reads and 65,641,113 (79.21%) were unique mapped reads ( Supplementary Table S2 ). Having analyzed the mapped reads in different genomic regions, we found approximately 58.41% of reads mapped to exons, 39.8% mapped to introns, and 1.79% mapped to intergenic regions ( Supplementary Table S3 ).

By comparing the distributions of the MeRIP-seq and RNA-seq reads in the lung tissues of mice in the OVA and PBS groups, we identified 27,292 m6A peaks within 14,118 transcripts in the OVA group ( Supplementary Table S4 ) and 26,102 m6A peaks within 13,564 transcripts in the PBS group ( Supplementary Table S5 ). We next analyzed the enrichment of the m6A peaks in the mRNA transcripts and found that 11,419 m6A peaks were enriched in 6,418 mRNA transcripts of the OVA group ( Supplementary Table S6 ) and 11,088 m6A peaks were enriched in 6,224 mRNA transcripts of the PBS group ( Supplementary Table S7 ). Further analysis revealed an overlap of 1,166 peaks and 5,557 m6A-modified mRNAs between the two groups. The OVA group had 10,253 unique m6A peaks and 861 unique m6A-modified mRNAs, while the PBS group had 9,922 unique m6A peaks and 667 unique m6A-modified mRNAs ( Figure 3A ). Meanwhile, the distribution characteristics of the peak frequency and width within the mRNA transcripts were analyzed in the two groups. The results showed that the majority of m6A-modified mRNAs (3,595 in the OVA group and 3,493 in the PBS group) contained only one peak site ( Figure 3B ), and the widths of the majority of peak sites (7,088 in the OVA group and 7,188 in the PBS group) were within 600 bp ( Figure 3C ).

We also investigated the genic region distributions of all m6A peak locations in asthmatic and normal lung tissues. Figure 3D shows that the peaks in both the OVA and PBS groups were dominantly distributed in the coding sequence (CDS) near the stop codon. Furthermore, the distribution patterns of the m6A peaks in genic regions of the two groups also showed similar trends. The peaks were abundant in the 3' untranslated (UTR) region (52.68% in the OVA group and 52.75% in the PBS group) and the exon region (35.71% and 35.87%, respectively), followed by the 5' UTR region (11.61% and 11.39%, respectively) ( Figure 3E ) We next conducted the motif prediction in the regions of m6A peak sites. Figure 4A shows the top five high confidential sequence motifs in the two groups. We found RRACH sequence patterns in both groups, which are typical sequence motifs proven to be related to m6A (35). Visualization of the distribution of the transcriptome-wide m6A peaks across the 20 chromosomes was performed in lung tissues. The results showed that the overlapping peaks between the two groups and the unique peaks of asthma were distributed on each chromosome, and the distribution patterns coincided with the gene content density ( Figure 4B ). Most of the hypomethylated peaks were enriched in chromosomes 2 (27peaks), 5 (23 peaks), and 7 (20 peaks), while most of the hypermethylated peaks were enriched in chromosomes 7 (76 peaks), 2 (74 peaks), and 1 (73 peaks) ( Figure 4B ). In addition, the hypomethylated peaks with the largest widths were distributed on chromosomes 15, 8, and 1, while hypermethylated peaks with the largest widths were distributed on chromosomes 16, 18, and 6 ( Figure 4B ).

The m6A methylation levels of the lung tissues in the OVA and PBS groups were compared in order to explore the changes and functions of m6A methylation in asthma. We found 855 hypermethylated peaks and 300 hypomethylated peaks within whole sets of transcripts ( Supplementary Table S8 ), and all differentially methylated m6A peak sites in the two groups most often appeared in the 3' UTR region (57.2%), followed by the exon (31.28%) and the 5' UTR region (11.52%) ( Figure 5A ). Furthermore, we screened out the methylated mRNA transcripts and found that there were 353 significantly hypermethylated peaks within 329 mRNA transcripts and 150 significantly hypomethylated peaks within 143 mRNA transcripts ( Figure 5B and Supplementary Table S9 ). These mRNA transcripts containing methylated m6A peaks are hereinafter referred to as differentially methylated mRNAs. The top 20 differentially methylated mRNAs are listed in Table 1 based on log₂(Fold change). We calculated the density of these differentially methylated peaks in each chromosome and found that they were not distributed homogeneously ( Figure 5C ). The top five chromosomes with the largest number of peak distributions were chromosomes 7, 1, 11, 2, and 5. In detail, chromosome 7 harbored the most number of differentially methylated m6A peaks (33 hyper- and 9 hypomethylated peaks), followed by chromosome 1 (32 hyper- and 9 hypomethylated peaks), chromosome 11 (28 hyper- and 10 hypomethylated peaks), chromosome 2 (26 hyper- and 11 hypomethylated peaks), and chromosome 5 (28 hyper- and 9 hypomethylated peaks). The number of hypermethylated peaks in each chromosome was distinctly larger than that of hypomethylated peaks (apart from chromosomes 12 and 19). We also found three differentially methylated peaks in two protein-coding mitochondrial transcripts (MT-ND4 and MT-ND3), and all three peak sites were hypomethylated. MT-ND4 contained two hypomethylated peaks, while MT-ND3 contained one ( Supplementary Table S9 ). Besides, we mapped the differentially methylated peaks (within mRNA transcripts) to the mouse chromosomes and visualized the peak site locations across the chromosomes using the R package RIdeogram (36) ( Figure 5D ).

To explore the biological sense of the differentially methylated mRNAs and their potential role in the pathological processes of asthma, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. Figure 6A provides an overview of the distributions of the differentially methylated mRNAs enriched in various GO categories. In the biological process (BP) category, differentially methylated mRNAs were mainly enriched in “regulation of DNA-templated transcription” (87%), “signal transduction” (81%), and “multicellular organism development” (64%). In the cellular component (CC) category, the top three enriched functions were “membrane” (91%), “cytoplasm” (76%), and “nucleus” (75%). In the molecular function (MF) category, “protein binding” (82%), “metal ion binding” (47%), and “DNA binding” (26%) were the most enriched terms. We next analyzed the significant enrichment degree of the differentially methylated mRNAs in various GO terms. The differentially methylated mRNAs were found to be significantly enriched in 120 GO terms (p < 0.05), suggesting that they are related to diverse biological functions. For instance, the differentially methylated mRNAs were significantly enriched in “G protein-coupled receptor signaling pathway”, “sensory perception of smell”, “positive regulation of cell proliferation in bone marrow”, “regulation of Wnt signaling pathway” (BP; Figure 6C ), “chromocenter”, “membrane attack complex”, “CCR4-NOT core complex” (CC; Figure 6D ), “G protein-coupled receptor activity”, “calcium ion transmembrane transporter activity”, and “bioactive lipid receptor activity” (MF; Figure 6E ). Intriguingly, we also found that these differentially methylated mRNAs were significantly enriched in GO terms related to T-cell differentiation and inflammatory response, such as “T-helper 1 cell differentiation”, “cytokine activity”, and “inflammatory response to wounding”. Detailed information on the enriched GO terms and the corresponding gene sets are shown in Supplementary Table S10 .

Likewise, we analyzed the distributions of the differentially methylated mRNAs in the six KEGG categories. The results are shown in Figure 6B . In the cellular process category, the differentially methylated mRNAs were mostly distributed in “cellular community-eukaryotes” (11%) and “transport and catabolism” (7%). The most enriched pathways in the environmental information processing category were “signaling molecules and interaction” (19%) and “signal transduction” (10%). In the genetic information processing category, “transcription” (8.33%) and “folding, sorting and degradation” (8.26%) were mainly enriched. “Substance dependence” (20.93%) and “immune disease” (15.79%) were the most enriched pathways in the human disease category. In the metabolism category, the most enriched pathways were “energy metabolism” (17.39%) and “metabolism of cofactors and vitamins” (14.29%). The most enriched pathways in the organismal systems category were “sensory system” (15%) and “immune system” (14.40%). The KEGG pathways significantly enriched by the differentially methylated mRNAs were “complement and coagulation cascades”, “calcium signaling pathway”, “herpes simplex virus 1 infection”, and “hematopoietic cell lineage” (p < 0.05) ( Figure 6F ). Among these pathways, hematopoietic cell lineage is closely related to asthma, and the mRNAs enriched in this pathway included FCER2A, GP1BA, IL-7, and CD5. Detailed information on the enriched KEGG pathways and the corresponding gene sets are shown in Supplementary Table S11 .

We predicted the RBPs that might bind to the differentially methylated peaks within the whole sets of transcripts via exploring the RMBase (v2.0) database (37). The 25 potential RBPs in the 299 hypermethylated peaks and 24 potential RBPs in the 96 hypomethylated peaks were obtained. We classified these bound peaks into nine groups according to log2(FC). Figure 7A shows the binding percentage of RBPs in the total differentially methylated m6A peaks. The results displayed that the binding frequency of Cstf2, Mbnl3, Upf1, Ago2, Taf15, Fus, Fmr1, Mbnl2, and Nova in the peaks of each group shared similar distribution patterns. The RBPs which bound the most of peaks were Cstf2 (bound 226 peaks), Mbnl3 (bound 167 peaks), Upf1 (bound 121 peaks), Ago2 (bound 116 peaks), Taf15 (bound 112 peaks), and Fus (bound 109 peaks). The binding abundance of RBPs in the hypermethylated peaks is higher than that in the hypomethylated peaks ( Figure 7B ), which indicated that these RBPs prefer to target hypermethylated peaks. Additionally, the RBPs and their binding peaks were mostly distributed in the group of differentially methylated peaks with log₂(Fold change) = 1 ( Figure 7C ).

We then performed GO and KEGG enrichment analyses on these 25 RBPs, and the results revealed that RBPs were significantly enriched in GO terms and KEGG pathways related to the biogenesis and metabolism of RNA. For instance, RBPs were significantly enriched in the GO terms “RNA splicing”, “regulation of alternative mRNA splicing via spliceosome”, “mRNA processing”, “regulation of RNA splicing”, “mRNA 3' UTR binding”, “spliceosomal complex”, “mRNA splice site selection”, and “mRNA stabilization” ( Figure 7D ). The results of KEGG enrichment showed that these RBPs were significantly enriched in “mRNA surveillance pathway”, “spliceosome”, and “RNA transport” ( Figure 7E ).

To identify the differential gene expression profiles in asthmatic lung tissues, we compared the differential gene expressions of the OVA and PBS groups and found a total of 3,766 differentially expressed genes ( Supplementary Table S12 ), among which 2,170 were upregulated and 1,596 were downregulated in asthma ( Figure 8A ). A similar trend of expression was displayed in the same group of lung tissues regarding the differentially expressed genes ( Figure 8B ). We further conducted GO and KEGG enrichment analyses for these differentially expressed genes. The results of GO enrichment showed that the function of these genes was involved in immune response and regulation, including “immune system process”, “immune response”, “antigen binding”, “positive regulation of B-cell activation”, and “chemokine-mediated signaling pathway” ( Figure 8C ). KEGG analysis identified a bunch of significantly enriched pathways related to immune regulation and diseases, such as “cytokine–cytokine receptor interaction”, “chemokine signaling pathway”, “HIF-1 signaling pathway”, “systemic lupus erythematosus”, “rheumatoid arthritis”, and “asthma” ( Figure 8D ).

We next conducted a combined analysis of the MeRIP- and RNA-seq data. Significant differential expressions were found in 170 differentially methylated mRNA transcripts ( Figure 9A and Supplementary Table S13 ). Specifically, from the 127 hypermethylated mRNAs, we obtained 107 highly and 20 low expressed mRNAs. From the 43 hypomethylated mRNAs, we obtained nine highly expressed and 34 low expressed mRNAs ( Figure 9A ). The top 20 differently methylated mRNAs with differential expressions are shown in Table 2 . To clarify the relationship between the methylation levels and the expression levels of these 170 mRNAs, we further performed Pearson’s correlation analysis and found that the methylation levels of these mRNAs were positively correlated with their expression levels (R ² = 0.2007, P = 1.01e−9) ( Figure 9B ). This suggested that m6A methylation plays a crucial role in the regulation of gene expression in asthma.

To identify the roles of these differentially methylated mRNAs with expression changes in asthma, we performed GO and KEGG enrichment analyses. For the GO enrichment results, the GO terms most significantly enriched were “integral component of plasma membrane”, “G protein-coupled receptor activity”, and “G protein-coupled receptor signaling pathway” ( Figure 9C ). We also found that many GO terms related to immune function and inflammatory response were significantly enriched, such as “inflammatory response”, “innate immune response”, “positive regulation of T-cell activation”, “positive regulation of cytokine production”, “response to interleukin-1”, “chemokine receptor activity”, and “T-cell chemotaxis” ( Supplementary Table S14 ), which indicated that m6A methylation may have an immunoregulatory effect in asthma. The KEGG enrichment results showed that these differentially methylated mRNAs with expression changes were significantly enriched in “one carbon pool by folate”, “alcoholism”, “cytosolic DNA-sensing pathway”, “hematopoietic cell lineage”, among others ( Figure 9D and Supplementary Table S15 ). Given that the pathological mechanism of asthma involves immune dysregulation and excessive inflammation, we thus explored whether there are differential m6A-modified mRNAs associated with pathological processes. Interestingly, a bunch of differentially methylated mRNAs related to immune regulation were found. Here, we enumerated five representative mRNAs (CCDC88B, PLD4, ZBP1, GPR183, and CCR2) and visualized the distribution of the peak regions in these mRNAs using IGV software ( Figure 10 ).

The m6A methylation of RNA is a reversible epigenetic modification (38), and its methylation level is tightly regulated by m6A methylase and demethylase. m6A modification is mediated and facilitated by methylases such as METTL3 (methyltransferase-like 3), METTL14, and WTAP (Wilms’ tumor 1-associated protein), whereas modification can be reversed by demethylases such as FTO (fat mass and obesity-associated protein) and ALKBH5 (alkylation repair homolog 5) (35). Thus far, only two studies have directly investigated the relationship between asthma and m6A modification. One of them showed that m6A enzymes may have a role in the pathogenesis of asthma through bioinformatics analysis of the public dataset from the GEO database (39). This study analyzed the gene expressions of m6A enzymes in the peripheral blood mononuclear cells of children (average age, 8.4 years) with asthma and found two different expression patterns of the m6A enzymes that were associated with the immune responses of T helper cells (Th1 and Th2). Besides, the degree of immune infiltration was higher in patients with an upregulated methylase WTAP gene expression compared to that in patients with a downregulated expression (39). Another intensive research reported on the potential mechanism of m6A modification in the Th2 immune response of asthma (40). The research group found that OVA-induced FTO knockout mice exhibited evident defects of motile ciliogenesis and asthma-like inflammatory phenotypes, such as the upregulated expression of the Th2 cytokine (IL-13) and the downregulated expression of the Th1 cytokine (IL-12b) mRNA in lung tissues, and significantly increased levels of IL-13 and IL-4 in BALF (40). They speculated that the decrease of ciliary cells in the airway epithelium may be responsible for the aggravation of airway Th2 immune response.

Dysfunction of airway epithelial cells underlies the pathogenesis of airway hyperresponsiveness and inflammation in asthma. A recent study investigated the role of m6A modification in human airway epithelial cells (41). The results showed that the expression level of methylase METTL3 mRNA and the m6A modification level of total RNA were significantly increased in the human bronchial epithelial (HBE) and A549 cells induced by fine particulate matter (PM_2.5). Subsequent findings suggested that the METTL3 protein could maintain the stability of the oxidative stress-induced growth inhibitor 1 (OSGIN1) mRNA by increasing the m6A methylation level and thus participate in PM_2.5-induced airway epithelial cell damage (41). Likewise, another recent study examined the role of the demethylase FTO in the differentiation of HBE cells to ciliated cells (40). The results showed that the mRNA levels of FTO and master ciliary transcription factor (FOXJ1) were significantly increased in epithelial cells cultured under the condition of air–liquid interface after 14 days. The protein levels of FTO and FOXJ1 showed the same trend after 16 days of culture. In addition, FTO knockdown HBE cells exhibited inhibition of the expression of FOXJ1 and differentiation of ciliated cells; the goblet cell numbers were instead increased (40). The researchers further observed that the multi-ciliated cells and the FOXJ1 mRNA levels in the lung tissues of FTO knockout mice were lower than those in wild-type mice. These findings revealed that m6A demethylation may affect the formation of motor cilia in the airway epithelium by regulating the expression of FOXJ1.

However, no previous study has comprehensively examined the m6A modification characteristics in strictly identified animal models of asthma or patients with asthma. The above studies neither explored the mechanism of m6A modification that contributed to the asthmatic immune response nor the role of m6A modification in the most important phenotype of asthma, airway hyperreactivity. To complement these studies, we therefore employed MeRIP-seq and RNA-seq to reveal, for the first time, the m6A methylation characteristics of asthmatic mice at the whole transcriptome level and to analyze in detail the asthma-related methylated genes and their potential functions.

Airway inflammation and hyperreactivity are two typical pathological hallmarks of asthma (42). The inflammatory process of asthma is characterized by excessive Th2 immune response caused by the imbalance of the Th1/Th2 cell ratio. The increased Th2 cells subsequently secrete large quantities of cytokines such as IL-4, IL-5, and IL-13. IL-4 promotes the differentiation of memory B cells, which produce a plethora of IgE when stimulated by antigens. IL-5 induces the infiltration of eosinophilic cells into lung tissues, and IL-13 activates a vast panel of inflammatory cells and promotes mucus secretion and bronchial constriction (43). IgE and Th2-type cytokines further aggravate airway hyperreactivity and elicit airway remodeling by activating fibroblasts (28). In the OVA-induced mice, we observed typical airway hyperresponsiveness and reduced pulmonary compliance, which represent the most important indicators to evaluate the success of the establishment of a mouse model of asthma. We further detected that the OVA-induced mice presented excessive Th2-type immune response; for instance, the levels of IL-4, IL-5, and IL-13 in serum and IgE in BALF were significantly increased. Moreover, the number and percentage of eosinophils in the BALF of the model group were significantly higher than those in the control group; so were the number of total and different subtypes of white blood cells. HE, PAS, and MT staining respectively displayed massive infiltration of inflammatory cells, mucus hypersecretion, and airway structure change in lung tissues. Overall, these results validated that OVA successfully induced airway hyperresponsiveness and inflammation, goblet cell hyperplasia, and airway remodeling, which conform to the pathological changes of asthma (44, 45), thus demonstrating the effectiveness of our model.

Our results showed that there were prominent m6A modifications in both normal and OVA-induced lung tissues. The two groups have approximately the same distribution characteristics of the m6A peaks, such as the peak frequency and width within mRNA transcripts and the distribution of m6A peaks in genic regions. Nevertheless, there were a large number of unique m6A peaks and methylated mRNAs in OVA-induced lung tissues. We speculated that OVA induced a specific pattern of m6A methylated modification. m6A modification was widely distributed in the RNA transcripts of eukaryotes, and the m6A peaks were highly abundant in different organs and tissues. In this study, we identified 27,292 m6A peaks in OVA-induced lung tissues. Most of these m6A peaks were distributed in the 3' UTR and exon regions, which is inconsistent with the previously reported peak distribution trend in mouse liver (46), sheep liver (47), rat heart (48), and human tumor tissue (49, 50). We considered that such inconsistencies may be related to different species, tissues, and conditions. In the OVA-induced lung tissues, we also found various subsets of the typical sequence motif RRACH which have been confirmed to be closely related to the enrichment of m6A peaks (35). Our results demonstrated that m6A methylation might have a role in asthma disease.

We analyzed the distribution of the differentially methylated peaks in the mRNA transcripts and found that there were 353 significantly hypermethylated and 150 significantly hypomethylated m6A peaks within 329 and 143 mRNA transcripts, respectively. To determine the biological functions of the differentially methylated mRNAs in asthma, we performed GO and KEGG enrichment analyses. The GO enrichment results showed that these differentially methylated mRNAs were significantly enriched in 120 GO terms. Of these, “regulation of Wnt signaling pathway”, “T-helper 1 cell differentiation”, “G protein-coupled receptor signaling pathway”, “G protein-coupled receptor activity”, and “cytokine activity” were closely related to immune function. Several previous studies have found that the Wnt signaling pathway have immunomodulatory effects and are closely related to lung function in asthma. For instance, the Wnt signaling pathway has been shown to be positively correlated with the airway Th2 immune response in asthmatic patients (51, 52). The Wnt/β-catenin pathway combined with transforming growth factor beta (TGF-β) can regulate airway remodeling mediated by epithelial–mesenchymal transition in asthma patients (53). Wnt1 was able to reduce airway allergic inflammation and the hyperresponsiveness of asthma by inhibiting the activation of pulmonary dendritic cells (54). Activation of Wnt-3A promoted the maturation of mast cells, thereby further aggravating the infiltration of immune cells in asthmatic airways (52). Th1 cells are considered to inhibit the function of Th2 cells and thus play an anti-asthmatic role (55). The overreaction of Th2 cells led to the imbalance of the Th1/Th2 ratio, which may increase exacerbations in asthma (56). In contrast, several researches have verified that Th1 cells promoted airway neutrophil infiltration and Th2 cell-induced eosinophil recruitment. Both Th1 and Th2 cells were able to induce airway inflammation. Moreover, the airway hyperresponsiveness induced by Th1 cells was even more severe than that induced by Th2 cells (57). A study has shown in vitro that OVA, IL-2, and IL-18 could stimulate Th1 cells to produce a vast panel of cytokines such as IFN-γ, IL-9, IL-13, GM-CSF, RANTES, and MIP-1. In vivo experiments have suggested that intranasal interventions with OVA and IL-18 in mice receiving adoptive transfer of Th1 cells could induce severe airway inflammation and hyperresponsiveness (58).

Differentially methylated mRNAs were enriched on the G protein-coupled receptor (GPR)-related GO terms. Multiple studies have confirmed that many GPR family molecules were involved in the functional regulation of immune cells. For example, lysophosphatidylserine suppresses the IL-2 secretion of CD4⁺ T cells through GPR174 (59). GPR174 inhibited the proliferation of regulator T cells in the thymus (60), and GPR84 signaling could promote macrophages to produce inflammatory cytokines and chemokines, such as TNF-α, IL-6, IL-12b, CCL2, CCL5, and CXCL1 (61). Lipopolysaccharide (LPS) induced GPR84 expression in monocytes/macrophages, and GPR84, in turn, promoted IL-12b secretion by monocytes/macrophages (62), which further facilitated the differentiation of Th1 cells and the generation of related cytokines, thus affecting the balance of Th1/Th2 (63). Medium-chain free fatty acids regulated the balance of Th1/Th2 by directly acting on GPR84 (62), and GPR21 can promote the migration of macrophages to inflammatory tissues (64). In addition, the differentially methylated gene CXCL16 enriched in the GO term “cytokine activity” could enhance Th1-type immune response (65).

The KEGG enrichment results showed that the differentially methylated mRNAs significantly enriched in the pathway “hematopoietic cell lineage” were FCER2A, GP1BA, IL-7, and CD5. Hematopoietic stem cells are well-known progenitors that can differentiate into multiple essential immune cells, including T cells, B cells, natural killer cells, dendritic cells, and macrophages (66). The differentially methylated mRNAs enriched in this pathway were involved in the regulation of immune function. For example, FCER2A (CD23) is associated with the activation of B cells (67). GP1BA, a glycoprotein component on the surface of platelets is often related to clotting function (68) and has also been considered to participate in lung inflammation (69, 70). IL-7 is considered to activate eosinophils and has a pivotal function in allergic inflammation in asthma (71), but is also essential for T- and B-cell development (72). CD5, a transmembrane protein, is expressed in both T and B cells that affects thymocyte development and negatively regulates T-cell activation (73). In our case, we speculated that these differentially methylated mRNAs may be playing a role in the pathogenesis of asthma by regulating the functions of immune cells. Additionally, these differentially methylated mRNAs were also significantly enriched in pathways including “complement and coagulation cascades”, “calcium signaling pathway”, and “herpes simplex virus 1 infection”. Although these pathways have not been reported in asthma, they deserve attention.

RBPs have substantial implications in gene expression and the occurrence and progression of diseases. RBPs affect RNA fate through binding to RNAs for regulatory alternative splicing, nuclear export, stability, localization, and translation (74). The m6A sites in RNA have been found to be specifically recognized by a number of RBPs, such as the YTH domain family proteins and HNRNP family proteins (35), which are the executors of m6A-associated biological functions. Due to certain diverse RBPs that can bind to m6A peak sites, the regulation of RNA processes by m6A methylation is highly heterogeneous (75). Interestingly, recent research has reported that the RBPs were also directly involved in the pathogenesis of asthma. The researchers adopted the method of OVA combined with complete Freund’s adjuvant induction to establish refractory neutrophilic inflammation of asthma. They observed that the RBP Mex-3B knockout mice have less neutrophilic airway inflammation and airway responsiveness than did wild-type mice (76). To identify how differentially methylated mRNAs function, we explored the RMBase (v2.0) database (37) to predict the potential RBPs that may bind to the methylated peak sites. Subsequent results showed that 395 differentially methylated sites were bound by 25 RBPs in asthmatic lung tissues. Accordingly, a study also predicted 25 RBPs that bound to 163 differentially methylated sites in liver tissues of high-fat diet-induced mice (46). This hinted at the distribution of RBPs among different tissues being conserved. In order to identify the functions of these RBPs in asthmatic lung tissues, we conducted GO and KEGG enrichment analyses. The results showed that the functions of these RBPs prevailed in mRNA biogenesis and metabolism, such as “mRNA processing”, “regulation of RNA splicing”, “mRNA binding”, “mRNA splice site selection”, “mRNA stabilization”, “mRNA surveillance pathway”, and “RNA transport”. We speculated that RBPs may affect the expressions of these differentially methylated transcripts by regulating RNA processes, thus participating in the pathological mechanism of asthma, which is worth further exploration in a future study.

m6A methylation has been reported to affect mRNA stability and degradation and relies on the binding of m6A-related RBPs, thereby regulating mRNA expression (2, 38). To investigate whether m6A methylation affects gene expression in asthma, we assessed the expressions of genes in the differentially methylated sites. Firstly, we analyzed the gene expression profiles of asthmatic lung tissues and obtained 2,170 upregulated and 1,596 downregulated genes. These differentially expressed genes were mostly related to the immune regulatory mechanism of asthma, such as “immune system process”, “innate immune response”, “B-cell receptor signaling pathway”, “immunoglobulin production”, “immunoglobulin receptor binding”, “phagocytosis”, “recognition”, “antigen binding”, “chemokine signaling pathway”, and “cytokine–cytokine receptor interaction”.

Next, we proceeded to jointly analyze the differentially expressed genes and differentially methylated peaks in asthmatic lung tissues and obtained 107 hypermethylated mRNAs with high expression, 20 hypermethylated mRNAs with low expression, 9 hypomethylated mRNAs with high expression, and 34 hypomethylated mRNAs with low expression. In addition, Pearson’s correlation analysis showed that the methylation levels of these genes were positively correlated with their expression levels. In our results, the hypermethylated mRNAs with high expression accounted for the largest proportion in asthmatic lung tissues. We concluded that, in asthma, the effect of m6A methylation in maintaining the stability of mRNAs predominates. However, other differential mRNAs, such as hypermethylated mRNAs with low expression and hypomethylated mRNAs with differential expression, may also play important roles in asthma, although they accounted for only low proportion. In this case, which m6A-related RBPs regulated these differentially methylated mRNAs with differential expression and how these mRNAs play a role in asthma will be further investigated in the future.

Furthermore, the results of the enrichment analyses showed that a portion of these mRNAs was associated with immune function such as “inflammatory response”, “innate immune response”, “positive regulation of T-cell activation”, “positive regulation of cytokine production”, “response to interleukin-1”, “chemokine receptor activity”, “T-cell chemotaxis”, and “hematopoietic cell lineage”. We selected five representative differentially methylated mRNAs with differential expressions (CCDC88B, PLD4, ZBP1, GPR183, and CCR2) from the above pathways and reviewed their roles in immunity. CCDC88B has a critical impact on the regulation of T-cell function. It has been reported that CCDC88B was highly expressed in CD4⁺ and CD8⁺ T cells and could regulate the maturation, activation, and proliferation of T cells and the production of related cytokines such as IFN-γ and TNF (77). In addition, CCDC88B was also involved in the occurrence of inflammatory diseases such as colitis (78). For example, CCDC88B expression was significantly elevated in lymphatic and myeloid cells of the colon lamina propria in patients with ulcerative colitis and Crohn’s disease. Similarly, CCDC88B-positive lymphocytes and myeloid cells were distinctly increased in the inflammatory tissues from a dextran sodium sulfate-induced mouse colitis model. Further research has shown that the deficiency of CCDC88B in immune cells can markedly diminish the inflammatory injury of the intestinal epithelial barrier and play an anti-colitis role (78). PLD4 is closely associated with auto-inflammatory and autoimmune conditions (79). A previous study reported on the inflammatory phenotypes in PLD4-deficient mice, such as splenomegaly, and increased concentrations of IFN-γ and CXCL10 in plasma. Dendritic cells in PLD4-deficient mice exhibited enhanced TLR9-mediated inflammatory responses (80). PLD4 is considered to be a 5' exonuclease that can degrade the ligand of TLR9, thereby reduce TLR9-mediated inflammatory response. When both PLD3 and PLD4 genes are lacking, neonatal mice developed severe liver inflammation with extensive infiltration of CD68⁺ myeloid cells and exhibited significantly increased levels of cytokines (IL-6, IL-10, TNF, and IFN-γ) and chemokines (CXCL10, CCL7, and CCL2) in serum (80). Besides, PDL4 is also believed to play a regulatory role in B-cell activation (81). ZBP1 is a key inflammatory mediator (82) and an innate immune sensor, which can regulate programmed cell death and inflammatory response (83). For instance, ZBP1 can mediate the necrosis or apoptosis of macrophages and induce the generation of inflammasome, thereby promoting the production of IL-1β. ZBP1 also elicited the activation of the classical inflammatory signaling pathway NF-κB (84) and may be involved in IL-17-mediated immune responses (85). The G-protein-coupled receptor GPR183 was highly expressed in immune cells and played an important role in adaptive immune response (86). GPR183 can not only promote the migration of CD4⁺ T cells to the outer T region and the differentiation of T follicular helper cells (87) but also control the migration of B cells and dendritic cells and the production of T-cell-dependent antibodies (88, 89). CCR2 belongs to a type of chemokine receptors widely expressed on monocytes, macrophages, Th1, natural killer cells, basophils, and immature dendritic cells and participates in monocyte migration and Th1-type adaptive immunity (90, 91). Evidence suggests that CCR2 can promote the release of monocytes from bone marrow into peripheral circulation and mediate the further migration of monocytes to local inflammatory tissues (92). CCR2 also plays a significant part in Th2-type immune response. For example, CCR2 knockout mice induced by diesel exhaust particles displayed complete abolishment of the recruitment of monocytes and dendritic cells in lung tissues and the disappearance of Th2 response in mediastinal lymph nodes (93). Otherwise, it has been reported that both Th1 and Th2 immune responses can promote the mRNA expression of CCR2 chemokine ligands CCL2, CCL7, and CCL12 in skin tissues (94).

Based on all of the above findings, we hypothesize that the methylation alterations, by affecting the RNA life cycle via RBPs, can regulate mRNA expression and function in immune cells, thereby participating in the allergic immune response of asthma.

Taken together, we identified a list of differentially methylated mRNAs and their potential binding proteins, which may serve as vital regulators in the pathogenesis of asthma. To our knowledge, we mapped the first mouse m6A methylomic landscape in the lung tissues following OVA-induced asthma, and our results provided a novel direction for understanding the mechanism of m6A methylation affecting allergic inflammation and screening potential therapeutic targets in asthma.

However, in this study, the expressions of the m6A-related enzymes and RBPs in the lung tissues of asthmatic mice were not examined. The interaction between the differentially methylated genes and RBPs and the mechanism by which RBPs regulate the expression levels of differentially methylated genes in asthma will be the focus of future research. Moreover, the roles of the differentially methylated mRNAs and the m6A-related enzymes and RBPs in asthmatic phenotypes such as airway inflammation, hyperresponsiveness, and remodeling deserve further investigation.

Female BALB/c mice (6–8 weeks old) were raised and maintained at the Animal Center of Fudan University School of Pharmacy under specific pathogen-free conditions with food and water ad libitum. The animal study protocol was reviewed and approved by the Experimental Animal Ethics Committee of Fudan University School of Pharmacy (approval no. 2020-12-HSYY-WY-01). The mouse model of acute allergic asthma was prepared based on the procedures reported in a previous research, with slight modifications (95). The schematic design of the animal study is shown in Figure 1A . In detail, mice were randomly assigned to the normal control group (PBS, n = 12) and the asthma group (OVA, n = 12). For allergy sensitization, on days 0 and 7, each mouse was intraperitoneally injected with 20 μg/ml OVA (grade V, Sigma-Aldrich Inc., St. Louis, MO, USA) and 2 mg/ml aluminum hydroxide adjuvants (Sigma-Aldrich Shanghai Trading Co., Ltd., Shanghai, China) in 0.2 ml PBS solution or with 0.2 ml PBS alone as a control. For allergen challenge, from days 14 to 20, mice were exposed to 3% OVA (grade II, Sigma-Aldrich Inc., St. Louis, MO, USA) in PBS solution aerosolized via ultrasonic nebulizer (Yuyue Medical Equipment Co., Ltd., Shanghai, China) once daily for 30 min or exposed to PBS as a control.

A bronchial provocation test was performed to determine airway hyperresponsiveness in mice. Here, we adopted methacholine (Sigma-Aldrich Shanghai Trading Co., Ltd., Shanghai, China) as the stimulus drug. One day after the last OVA challenge (day 21), airway resistance and dynamic compliance (C _dyn) in mice were evaluated using the FinePointe RC System (DSI Buxco Electronics, Troy, NY, USA). Specific details of the measurement were described in our previous study (44, 45). Briefly, mice were intubated by tracheotomy under anesthesia and ventilated with 2 μl nebulized methacholine of different concentrations (0, 3.125, 12.5, and 50 mg/ml). The mouse bronchus was sequentially provoked with methacholine of increasing concentrations and the extrema of airway resistance and C _dyn recorded for each provocation. After the bronchial provocation test, whole blood was collected from anesthetized mice via quickly removing an eyeball. Mice were subsequently sacrificed by cervical dislocation under anesthetic status. Then, the BALF and lung tissues were collected and subjected to subsequent assays.

Whole blood was centrifuged (3,000 rpm, 30 min) and serum was separated for detecting the level of IgE. BALF was centrifuged (500 × g, 10 min) and the resultant supernatant was used for measuring the levels of IL-4, IL-5, and IL-13. These immune and inflammation-related protein levels were measured using ELISA following instructions in the kit (MultiSciences Biotech Co., Ltd., Hangzhou, China). The precipitate from BALF was collected to calculate the numbers and percentages of the different subsets of leukocytes using Auto Hematology Analyzer (Shenzhen Mindray Bio-Medical Electronics Co., Ltd., Shenzhen, China).

Lung tissues were fixed in 4% formaldehyde for more than 24 h and made into a paraffin section. HE, PAS, and MT staining was performed on lung sections to assess the degree of inflammation, mucus secretion, and collagen deposition in the lungs. Histopathologic results were interpreted using scoring methods described in previous studies via ImageJ software (96).

MeRIP-seq and RNA-seq were performed by LC-Bio Technology Co., Ltd. (Hangzhou, China). Detailed experimental procedures were described in published literature (49). In brief, three biological replicates were used in each group. Total RNA from lung tissues was isolated using Trizol (Invitrogen, Carlsbad, CA, USA) in accordance with the protocol provided by manufacturer. We then used Nanodrop ND-1000 (NanoDrop, Wilmington, DE, USA) to determine the purity and concentration of total RNA and used Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA) for agarose electrophoresis to detect and verify the integrity of total RNA (RIN values >7.0). Ribosomal RNA was removed from total RNA using the Epicentre Ribo-Zero Gold Kit (Illumina, San Diego, CA, USA) and the remaining RNA was fragmented (86°C, 7 min) using RNA Fragmentation reagents (NEB, Ipswich, MA, USA). The fragmented RNA was incubated with m6A antibody (Synaptic Systems, Göttingen, Germany) as an immunoprecipitated RNA and was used for MeRIP-seq. The fragmented RNA without immunoprecipitation (IP) as an input RNA was used for RNA-seq. The IP and input RNAs were catalyzed by reverse transcriptase (Invitrogen, Carlsbad, CA, USA) to synthesize the first strand of complementary DNA (cDNA). Then, RNase H (NEB, Ipswich, MA, USA) was used to hydrolyze the RNA strand in the RNA–DNA hybridization molecule, and Escherichia coli DNA polymerase I (NEB, Ipswich, MA, USA) was used to synthesize the second strand of cDNA. Meanwhile, dUTP (Thermo Fisher, Waltham, MA, USA) was incorporated into the second strand. Subsequently, an A-base was added to the ends of each strand for ligating the cDNAs to the indexed adapters. After the second strand of the cDNAs containing dUTP was degraded using UDG enzyme (NEB, Ipswich, MA, USA), cDNA libraries with a fragment size of approximately 300 bp were constructed by PCR amplification. The libraries were paired-end sequenced (PE150 mode) using the Illumina Novaseq 6000 platform at LC-Bio Technology Co., Ltd. (Hangzhou, China).

Clean data in fastq format were obtained by removing adapter, repeated, and low-quality sequences from the raw reads of the IP and input samples using FASTP software (97). Clean data were aligned to the Mus musculus genome from Ensembl version 101 using HISAT2 (98) for obtaining the BAM files. For peak calling and identification of the differentially methylated peaks, the BAM files were analyzed using R package exomePeak (99). IGV software was used for the visualization of peaks (100) and R package ChIPseeker was used for the annotation of peaks (101). Homer software was used to search for motifs with a high confidence level in the peak region of each group of samples (http://homer.ucsd.edu/homer/motif). Mapped reads were assembled and the gene expression levels (in fragments per kilobase of transcript per million mapped reads) were assessed by StringTie (102). The differential expressions of the genes between the two groups were analyzed using the R package edgeR (103). In our study, the screening criteria for differentially methylated peaks and differentially expressed genes were |log₂(Fold change)| > 1 and p < 0.05. We performed GO function and KEGG pathway enrichment analyses using the OmicStudio tools (https://www.omicstudio.cn/tool). The R package RIdeogram (36) was used to visualize the distribution of the differentially methylated m6A peaks within mRNAs. Correlation was analyzed using Pearson’s correlation. Differences between the OVA and PBS groups were analyzed using unpaired and two-tailed Student’s t-test, and lung function data were analyzed using two-way mixed-design ANOVA (GraphPad Prism 8.0). The data are presented as the means ± SD and differences were considered statistically significant when p < 0.05.