We modeled (Kianmeher et al., 2016; Liu et al., 2021) two groups of female 6-week-old BALB/C mice: asthma group and blank control group (n = 3 per group). They had the same feeding conditions and growth environment. Immunization solution: dissolve 20 mg ovalbumin (OVA) in 1 ml normal saline (NS); after OVA is completely dissolved, dilute 0.4–10 ml and mix well; and then mix it with the same volume of liquid aluminum adjuvant and place on a shaking table at 4°C for 30 min. Challenge solution: add 0.5 g OVA in 10 ml NS, fully dissolve it, and shake it on a shaking table at 4°C for 30 min. Immunization: mice were injected intraperitoneally on days 0 and 12, each with 0.2 ml, and the control group was treated with an equal volume of normal saline. Challenge: on days 18–23, the mice were atomized by ultrasound in a closed container at a dose of 10 ml once a day for 20 min. Lung tissue was taken 24 h after last atomization and immediately stored in liquid nitrogen. All experimental procedures used in this study were approved and conducted according to the Guidelines for the Care and Use of Laboratory Animal Management and Ethics Committee of Zhejiang Chinese Medical University.

RNA extraction and MeRIP-seq were conducted by the Beijing Genomics Institute. In brief, total RNA was extracted using the TRIzol reagent (Invitrogen), and the concentration and integrity of total RNA were measured by using the Qubit RNA HS Assay kit and an Agilent 2100 Bioanalyzer (Agilent Technology), respectively. For MeRIP experiment, about 10 ug of total RNA from each sample was fragmented using 10X RNA fragmentation buffer (Invitrogen, AM8740) by incubating in a preheated thermal cycler for 10 min at 70°C. The fragmented RNA was purified by using a RNA Clean & Concentrator™ kit (Zymo, R1018). Protein A and protein G magnetic beads (Invitrogen, 10002D, 10004D) were washed twice with IP buffer (150 mM NaCl, 10 mM Tris–HCl pH 7.5, and 0.1% IGEPAL CA-630 in nuclease-free water) before incubating with 5 ug m6A antibody (Synaptic Systems) at RT for 10 min. After two washes with IP buffer, antibody–bead complexes were resuspended in 500 μl of the IP reaction mixture including fragmented total RNA, and incubated for 4 h at 4°C. The immunoprecipitated m6A RNA with protein A/G magnetic beads was then washed three times with IP buffer for 10 min each at 4°C. Then, the beads complexes were resuspended in 500 μl reagent and collected by using the RNA Clean & Concentrator™ kit. For library preparation, the MeRIP libraries comprising eluted RNA were constructed using the SMARTer Stranded Total RNA-Seq Kit version 2 (Takara/Clontech), according to the manufacturer’s protocol. In brief, the eluted m6A RNA and input RNA were directly used for first-strand cDNA synthesis without fragmentation. After that, all the following steps were based on the manufacturer of the SMARTer Stranded Total RNA-Seq Kit version 2. Libraries for IP RNA were PCR amplified for less than 16 cycles and input libraries for less than 12 cycles. All libraries were analyzed by using an Agilent 2100 Bioanalyzer (Agilent Technologies) and quantified by using real-time PCR. Finally, the different libraries were pooled according to effective concentration and target downstream data volume and then sequenced on the HiSeq platform with the PE150 sequencing strategy.

After obtaining the original offline data, we performed bioinformatics analysis. First, we used Trim Galore software (version: 0.5.0) to carry out quality control steps such as removing connector sequences and low-quality bases from the original data. The parameters used were: --stringency 4 --quality 22—clip R1 13 –clip R2 13 --length 30 --paired. Hisat2 software (version: 2.1.0) (Kim et al., 2015) was used for reading alignments, and the reference genome was from the UCSC RefSeq database. Second, the m6A-modified area was detected by exomePeak software (version: 2.1.2) (Meng et al., 2014), and the parameters used were window width = 200, sliding step = 30, fragment length = 150, and fold enrichment = 1.5. Third, the differential m6A modification was identified by the R package exomePeak; it combined the peaks of the samples to be compared, calculated the cumulative number of reads in the combined peak of each sample, standardized these reads, made two groups of samples at a comparable level, and then tested whether there was a significant difference in the number of reads between the two groups of samples within the peak. The parameters used were window width = 200, sliding step = 30, fragment length = 150, diff peak abs fold change = 2, fold enrichment = 1.5, and fragment length = 200; the differential m6A modification was analyzed by annotation, distribution statistics, and motif identification. Fourth, we used StringTie software (Pertea et al., 2016) to calculate the expression of mRNA and displayed it with TPM, where TPM = Ri/Li*(1e6/sum (Ri/Li)). For lncRNAs, we used previously published methods (Yang et al., 2017); in brief, transcripts were first assembled using StringTie (Pertea et al., 2016) to obtain all transcripts in each sample and then labeled. The coding ability of candidate lncRNAs was then predicted using CNCI (Sun et al., 2013) and CPC (Kong et al., 2007) software, and the expression of lncRNAs was standardized using TPM. For circRNA, we used find_circ (Memczak et al., 2013) and CIRCexplorer2 (Zhang et al., 2014) for identification and calculated its expression using SRPBM (Zheng et al., 2016), where SRPBM = Ri/(Rtotal*Li). After obtaining the gene expression of each sample, differential gene analysis was carried out by edgeR software, and then the real differential genes were screened by threshold. Finally, the overlap of differential m6A-associated genes and differentially expressed genes was analyzed.

Data for validation were obtained from U-BIOPRED (Unbiased Biomarkers for the Prediction of Respiratory Disease Outcomes) (Bigler et al., 2017), a multicenter prospective cohort study involving 16 clinical centers in 11 European countries, and downloaded from Gene Expression Omnibus datasets (https://www.ncbi.nlm.nih.gov/geo/). The serial number is GSE69683, and the sample type is blood. A total of 87 healthy controls and 411 asthma cases were selected. The platform used was the GPL13158 [HT_HG-U133_Plus_PM] Affymetrix HT HG-U133 + PM Array Plate. R software and annotation packages were used to obtain gene symbols of the dataset. The differential expression of genes between asthma cases and healthy controls was analyzed by using the Wilcoxon test, and the up/down/unchanged genes were visualized using the R package “ggplot2”. The potential m6A-modified genes in patients with asthma were identified by univariate logistic regression and were cut off by p < 0.05. The least absolute shrinkage and selection operator (LASSO) Cox regression was used for feature selection and dimension reduction (McEligot et al., 2020), and the risk scores of potential asthma-related genes were calculated for verification [we used the rms package (Zhang J. A. et al., 2020) for the nomogram plotting and the nomogramFormula package (Bi et al., 2020) to calculate total points and probabilities of nomogram]. Receiver operating characteristic (ROC) curve analysis was used to evaluate distinguishing performance. The m6A methylation peak was displayed by IGV software according to the TDF file of sequencing samples, and the minimum value of data range was set to 0 to remove those non-specific peaks and low enrichment peaks (Thorvaldsdottir et al., 2013).

The original offline data includes the splice sequences introduced during library preparation and bases with low quality. These factors will lead to fewer reads to the genome, resulting in less information, so it needs to be filtered. We used Trim Galore to control the quality of the original offline data. Q20 and Q30 in the quality control results were calculated according to the correct rate of base recognition during sequencing and were the key indicators of base quality [Calculation formula: Qphred = −10log10P (error)]. It was found that the base ratio of Q20 in the asthma group and the control group was higher than 95%, and the proportion of Q30 bases was higher than 90%, which proves that the sequencing quality of this data is good and reliable. The results are shown in Table 1. Next, we used the software hisat2 to align the clean data to the reference genome. The results are shown in Table 2.

MeRIP-Seq enriched and sequenced the m6A-modified region on RNA; therefore, in the m6A-modified region, the number of reads covered by IP will be significantly higher than that of the input, thus forming a “peak.” The location of m6A modification on RNA can be obtained by detecting the location of these peaks. We used exomePeak software for peak detection. The number, total length, and average length of peaks in each group were counted. The results are shown in Table 3. The overlapping of peaks between samples was analyzed (Figures 1A,B). According to the results, the number of overlapped peaks in the asthma group was 13,481 and the number of overlapped peaks in the control group was 12,444, accounting for the majority of total peaks, which proved that the consistency of m6A modification was high. At the same time, the overlap of peak modification genes in the sample was visualized (Figures 1C,D). According to the results, the number of overlapped peak modification genes in the asthma group was 4,198, and that in the control group was 3,990.

The differential m6A modification (i.e., differential peak) was identified by using the R package exomePeak, and then the differential peak was counted by the R package Chipseeker (version: 2.16.0). A total of 568 peaks were detected, including 416 m6A methylation peaks that were significantly upregulated and 152 significantly downregulated (log2-fold change was used in this experiment to represent the ratio of standardized reads between the asthma group and the control group; the value >0 indicates that m6A modification in the asthma group is higher than that in the control group and vice versa). The number, total length, and average length of differential peaks between the asthma group and the control group are shown in Table 3 and Supplementary Table S1. The results of the top 20 differently expressed peaks are shown in Table 4. Next, we counted the distribution of the m6A peak. First, we counted the distribution of differential m6A modifications on chromosomes. The statistical method is to calculate the number of differential peak coverage of each base on the chromosome and draw the figure with the statistical file. (Figure 2A). Second, we analyzed the distribution patterns of differential peaks on mRNA; we found that the differential peaks were mainly distributed on CDS and 3′ UTR, and the highest peak of distribution was at the junction of CDS and 3′ UTR (Figure 2B). Finally, to understand its specific distribution on mRNA, the number of differential peaks distributed on each gene element was counted based on the location of differential m6A. This statistic helps to understand whether the distribution of differential m6A modifications has a preference for gene elements. According to the results, we found that the differential peaks had the highest percentage of 3′ UTR distribution with 42.61% (Figure 2C).

A motif refers to a short DNA sequence with a specific pattern. These sequences are likely to be important regulatory regions and related to biological functions. The m6A-modified sequence may have some characteristics. Therefore, the motif of the peak region was identified. We used motif analysis software HOMER (Hansen et al., 2016) to search motifs with high reliability in the peak area, and the width, p-value, and general position information of each motif in each peak sequence were obtained. The results are shown in Figure 3A. The canonical m6A motif—GGACU—was identified in the motif of the m6A-modified region in all samples of the asthma group [m6A mainly occurs on the motif of RRACH, where R = G/A; H = A/C/U (Wei et al., 2021)].

To investigate the biological functions of differentially m6A-modified genes, we analyzed them with Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). GO is divided into three ontologies: molecular functions (MFs), cellular components (CCs), and biological processes (BPs). Figure 3B shows the results of the GO analysis, and the enrichment table is shown in Supplementary Table S2. The most significantly enriched MF was small GTPase binding, and the most significantly enriched BP was cost chromatin modification. KEGG is a simulation of the biological system, including the molecular wiring diagram of interaction, reaction, and relationship network composed of molecular structural units of genes, proteins, and compounds, as well as the information of diseases and drugs. Figure 3C shows the results of KEGG enrichment, and the enrichment table is shown in Supplementary Table S3. We found that the signal pathway with the most significant enrichment of differentially m6A-modified genes was the Rap1 (Ras-proximate-1) signaling pathway. Rap1 is a small GTPase, a small cytosolic protein that acts as a cellular switch and is essential for efficient signal transduction (Burbach et al., 2007). The latest study found that it is associated with seasonal allergies and asthma symptoms in children (Tiwari et al., 2021). In a separate study (Hu et al., 2015), asthmatic HASM (human airway smooth muscle) cells were also found to exhibit increased constitutive direct binding of the small Rap1 GTPase-activating protein Rap1GAP to the Gi protein alpha subunit, which contributes to the synergistic promotion of Ras activation. This suggests that small GTPase binding and Rap1 signaling pathway have important roles in asthma.

m6A-seq input library is equivalent to RNA-seq library, which can be used to analyze gene expression and identify differentially expressed genes. Therefore, we used it to analyze the gene expression levels of mRNA, lncRNA, and circRNA. We used the software StringTie to calculate gene expression of mRNA and then standardized it with TPM and displayed it with a density map (Figure 4A). For lncRNA, we identified 13,208 lncRNA transcripts by combining the results of two prediction software (CNCI and CPC) (Figure 4B) and then standardized the expression of lncRNA with TPM and displayed it with a density map (Figure 4C). For circRNA, due to the high false-positive rate of circRNA identification, we used two software applications (find_circ and Circexplorer2) for circRNA identification, took the intersection of the results of the two software as the final prediction result, and identified a total of six circRNA transcripts (Figure 4D). Then, SRPBM was used to calculate the expression of circRNA and display it with a density map (Figure 4E). Finally, the TPM of mRNA was analyzed by a principal component analysis (PCA) (Figure 4F), and the correlation between two samples was calculated (Figure 4G). According to the results, it can be seen that the asthma group and the control group were significantly separated.

After obtaining the gene expression of each sample, the differential gene analysis was carried out between the asthma group and the control group. Because mRNA, lncRNA, and circRNA come from the same library, we combined these three RNAs for analysis when calculating differential expression, and there were 68,555 genes after combination. Next, the differentially expressed genes (DEGs) were identified by edgeR software [If the p-value < 0.05 and the absolute value of log2 (fold change) > 1, it is considered to be a DEG]. We found that compared with the control group, there were 2,505 significantly upregulated genes and 4,715 significantly downregulated genes in the asthma group (Supplementary Table S4). DEGs were displayed using MA map, volcano map, and heat map (Figures 5A–C). Table 5 shows the top 20 genes with the most differences in mRNA. Third, the DEGs were analyzed by GO and KEGG enrichment (Supplementary Tables S5, S6). The top 30 enrichment results are shown in Figures 5D,E. The most significantly enriched CC was an immunoglobulin complex, the most significantly enriched MF was immunoglobulin receptor binding, and the most significantly enriched BP was the antigen receptor-mediated signaling pathway. As for KEGG, the most significantly enriched item was the cytokine–cytokine receptor interaction. It is suggested that the DEGs are closely related to immune function.

To understand the overlapping relationship between genes associated with differential m6A peak and differentially expressed genes, we carried out a series of statistics. First, the differential peaks associated with each gene (including significant and insignificant) were obtained. One gene may match with multiple peaks, and in this case, this gene will appear multiple times; there may also be no matching differential peak, in which case the gene will not appear. Second, the overlapping relationship between the genes obtained in the previous step and the differentially expressed genes was identified. We found that compared with the control group, in the asthma group, one gene with a 2-fold increase in the m6A modification level and a 2-fold increase in the expression level was shown in red, four genes with a 2-fold decrease in the m6A modification level and a 2-fold decrease in the expression level were shown in blue, and 14 genes with a 2-fold increase in the m6A modification level and a 2-fold decrease in the expression level were shown in orange. The peaks matched by these genes were all significantly different. The results are shown in Figure 5F and Table 6.

To evaluate the clinical significance of gene expression regulated by m6A modification, the GEO database was explored. We selected 14 genes with differentially methylated m6A peaks and simultaneous differential expression according to the correlation analysis of differential m6A peaks and DEGs (Wang et al., 2019; Zhang et al., 2021). Because there was very large interference inducible genes among these 14 genes, taking the intersection of MeRIP-seq and RNA-seq and GEO, it was found that there were eight co-expressed genes, of which three were differential, namely, BCL11A, MATK, and CD300A (Figures 6A–C). Next, we used univariate logistic regression to analyze these three genes and found that they were all potential m6A-regulated genes of asthma (Figure 7A). Third, LASSO Cox regression was used for feature selection and dimension reduction, and it was found that these three genes were all important for asthma (Figures 7B,C). At the same time, the risk scores of the three genes were compared between the asthma group and the healthy control group. It was found that the risk scores of the three genes in the asthma group were significantly higher than that in the healthy control group (p = 4e-13) (Figures 7D,E); the ROC curve also illustrated that the three genes possess a good performance in classifying patients with asthma and healthy controls (Figure 7F). Finally, the methylation peaks of BCL11A, MATK, and CD300A were visualized according to the TDF file. It was found that the methylation peaks of three key m6A-regulated genes in asthma were different from those in the control group, and they were mainly distributed in exons and enriched in 3’ UTR (Figures 8A–C).