m6A regulator‐mediated RNA methylation modification patterns and immune microenvironment infiltration characterization in severe asthma

Abstract N6‐methyladenosine (m6A) modification is one of the most prevalent RNA modification forms of eukaryotic mRNA and is an important post‐transcriptional mechanism for regulating genes. However, the role of m6A modification in the regulation of severe asthma has never been reported. Thus, we aimed to investigate the m6A regulator‐mediated RNA methylation modification patterns and immune microenvironment infiltration characterization in severe asthma. In this study, 87 healthy controls and 344 severe asthma cases from the U‐BIOPRED (Unbiased Biomarkers for the Prediction of Respiratory Disease Outcomes) programme were used to systematically evaluate the m6A modification patterns mediated by 27 m6A regulators and to investigate the effects of m6A modification on immune microenvironment characteristics. We found that 16 m6A regulators were abnormal and identified two key m6A regulators (YTHDF3 and YTHDC1) and three m6A modification patterns. The study of infiltration characteristics of immune microenvironment found that pattern 2 had more infiltrating immune cells and more active immune response. Besides, it was found that the eosinophils which are very important for severe asthma were affected by YTHDF3 and EIF3B. We also verified key m6A regulators with merip‐seq and found that they were mainly distributed in exons and enriched in 3′UTR. In conclusion, our findings suggested that m6A modification plays a key role in severe asthma, and may be able to guide the future strategy of immunotherapy.

which seriously affects the quality of life of patients. Besides, its frequent attacks lead to high medical costs, which has become a serious social and economic burden, and also the main cause of disability and death of asthma. 4 Therefore, to improve the understanding of the pathogenesis of severe asthma is important to ameliorate its attack and prognosis, as well as to reduce medical costs. Furthermore, it has been found that the incidence of severe asthma is familial aggregation, and there are multiple susceptibility genes. The incidence of people with susceptibility genes is greatly affected by environmental factors. 5 So, in-depth study of geneenvironment interaction will help to reveal the genetic mechanism of severe asthma.
As one of the most important discoveries and research hotspots in recent years, epigenetics plays an important role in elucidating the interaction between genes and environment and changing the course of disease. It provides instructions for when, where and how to apply genetic information. 6 The in-depth study on the epigenetic mechanism of severe asthma is conducive to investigating the relationship between genes and environmental factors, and to formulating effective treatment strategies for severe asthma. 7 Epigenetics mainly includes DNA methylation, RNA modification and histone modification. N6-methyladenosine (m6A) modification is the most prevalent form of RNA modification in eukaryotic mRNA and even viral RNA. 8,9 M6A modification has been reported since the 1970s, but the overall distribution of the modification in RNA and its effect on gene expression regulation have been poorly understood. In 2011, the first real RNA demethylase fat mass-and obesity-associated gene (FTO) was reported, and the methylation modification of m6A was proved to be reversible, which made the study of mRNA methylation come into the eyes of scientists again. 10 The regulatory proteins of m6A were composed of methyltransferases (writers), demethylases (erasers) and methylated reading proteins (readers). 11 Methyltransferases include methyltransferase-like 3 (METTL3), methyltransferase-like 3 (METTL4), RNA-binding motif protein15 (RBM15), Cbl protooncogene-like 1 (CBLL1), WT1-associated protein (WTAP) and so on. Their main function is to catalyse the m6A modification of mRNA. 12 On the contrary, the demethylases is to demethylation of bases that have undergone m6A modification, it includes FTO and Human AlkB homolog H5 (ALKBH5). 13 Methylated reading proteins are mainly proteins of the YT521-B homology (YTH) domain family, including YTH Domain Containing 1(YTHDC1), YTH Domain Containing 2(YTHDC2), YTH m6A RNA-binding protein 1 (YTHDF1), YTH m6A RNA-binding protein 2 (YTHDF2) and YTH m6A RNA-binding protein 3 (YTHDF3); its main function is to identify the bases that undergo m6A modification and thus activate downstream regulatory pathways such as RNA degradation and miRNA processing. 14 Abnormalities of these regulators can affect mRNA in many aspects, including structure, splicing, translation, stability and so on, leading to the occurrence of disease. 15 It can even explain the potential mechanism of immune regulation in some diseases. Although more and more evidence shows the regulatory role of m6A in immune response 16 and also proves the diagnostic significance of m6A regulatory pattern in childhood asthma, 17 there is still a gap in the research on severe asthma and m6A. Therefore, In this study, we systematically evaluated the m6A regulator-mediated RNA methylation modification patterns and immune microenvironment infiltration characterization in severe asthma.

| Data download and processing
Data for this study were from U-BIOPRED 18

| Difference analysis of m6A regulators
The differential expression of 27 m6A RNA methylation regulators between severe asthma cases and healthy controls was analysed by Wilcox.test and the up/down/unchanged genes were visualized with R package 'ggplot2'. The R package 'corrplot' was used to identify the correlation between the methylation regulators. The potential m6A RNA methylation regulators in patients with severe 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, and the risk scores of potential severe asthma related genes were calculated for verification. Receiver operating characteristic (ROC) curve analysis was used to evaluate the distinguishing performance.

| Non-negative matrix factorization consensus clustering
According to the expression of 27 m6A regulators, Non-negative matrix factorization (NMF) was performed to identify different m6A modification patterns. NMF was to identify potential features in gene expression profile by resolving the original matrix into k non-negative matrices. 24 The deposition was repeated, and the results were aggregated to obtain consistent clustering. According to the co-occurrence coefficient, dispersion and profile coefficient, the most suitable number of subtypes was determined. NMF was performed by an R package called 'NMF'.
Principal component analysis (PCA) was used to further verify the expression patterns of 27 m6A regulators in different modification patterns.

| Biological function analysis of m6A modification patterns
In order to study the biological functions of genes in m6A modification patterns, we compared their Kyoto Encyclopedia of Genes (KEGG) pathway enrichment and gene ontology (GO) analysis, and analysed their biological functions by Gene Set Enrichment Analysis (GSEA) with R 'clusterprofiler' package. 25 In short, n random gene sets of size k were randomly sampled in a straightforward way. The 'ont' of gsego was set to 'all', which includes molecular function (MF), biological process (BP) and cellular component (CC) categories. All terms were ranked according to the enrichment score, and the top four results were taken.

| Analysis of infiltration characteristics of immune microenvironment
Single-sample gene-set enrichment analysis (ssgsea) was used to estimate the number and reactivity of specific infiltrating immune cells, which defined an enrichment score by calculating the enrichment degree of samples in a given data set. Wilcox test was used to compare the enrichment scores among the groups. HLA is important in immune microenvironment, 23 so we also analysed the association between HLA-related genes and m6A regulators. The correlation of m6A regulators with immune reaction activity and HLA gene expression was determined by Spearman correlation analysis with R package 'corrplot'.

| Identification of genes mediated by m6A regulators
The R software package 'WGCNA' was used to construct coexpression module networks and to identify the genes mediated by m6A regulators. In briefly, there were two steps. 26 The first step was to calculate the weighted value of the correlation coefficient, that is, to take the n-th power of the gene correlation coefficient to construct a scale-free network. In the second step, hierarchical clustering tree was constructed by Person Coefficient. Based on the weighted correlation coefficient of genes, genes were classified according to expression patterns, and the adjacency matrix was transformed into topological overlap matrix (TOM). Through dynamic tree cutting, modules were formed by the differences based on TOM.
Here, we set the minimum module size to 100 and the cutting height to 0.25. Finally, the correlation between key m6A regulators and related modules was calculated, and the related module was enriched by Kyoto Encyclopedia of Genes (KEGG) pathway and analysed by Gene Ontology (GO) using The Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8; the enrichment pathway was cut off by p < 0.05.

| Mouse model and determination of m6A methylation peak
We modelled two groups of female 6-week-old BALB/C mice 27 : severe asthma group and blank control group (n = 3 per group).
They had the same feeding conditions and growth environment.

| Identification of crucial m6A regulators in severe asthma
In order to investigate the m6A regulators' contribution to the pathogenesis of severe asthma, a series of bioinformatics analysis were used. We examined the expression data of the 16 differenceexpressed genes using univariate logistic regression to identify the potential m6A regulators, and two regulators (YTHDF3 and YTHDC1) were found to be associated with severe asthma (Figure 2A). Next, LASSO Cox regression was used for feature selection and dimension reduction, and it was found that these two regulators were all important for severe asthma ( Figure 2B,C). At the same time, the risk scores of the two regulators were compared between the severe asthma group and the healthy control group. It was found that the risk scores of the two regulators in the severe asthma group were significantly higher than that in the healthy control group (p = 6.

| Three patterns of m6A RNA methylation mediated by 27 regulators
According to the value of copynetic calculated by 'NMF' package, the optimal value of k is three ( Figure 3A showed that all the regulators had significant differences (p < 0.001) ( Figure 3F). We also found that the expression of these three clusters was different in the groups of smoking, gender and race, which indicated that the methylation of m6A RNA might be related to the clinical phenotype of severe asthma ( Figure 3G,H).

| Biological function analysis of three m6A modification patterns
In order to understand the biological reactions in the three m6A modification patterns, we compared the GO analysis and KEGG pathway between them. Compared with cluster 1 and cluster 2, GO analysis results mainly involved in immune process. According to the enrichment score, the top four items were as follows: immunoglobulin complex, MHC class II protein complex, immunoglobulin complex and MHC protein complex ( Figure 4A). The results of KEGG pathway include asthma, primary immunity, systemic lupus erythematosus and hematopoietic cell lineage ( Figure 4B). The item of "asthma" ranks first, it is proved that pattern 1 and pattern 2 of m6A methylation The risk distribution between healthy and severe asthma cases, where severe asthma cases have a much higher risk score than healthy controls. (F) The discrimination ability for healthy and severe asthma cases by m6A regulators was analysed by ROC curve and evaluates by AUC value modification may affect the immune progress of asthma. Compared with cluster 2 and cluster 3, GO analysis results showed that these genes could affect acetyl CoA biosynthetic process from pyramid, DNA double-strand break processing, central complex assembly and tRNA export from nucleus ( Figure 4C). KEGG pathway results showed that these genes could affect circadian rhythms, basic transcription factors, ribosome biogenesis in eukaryotes and ribosome ( Figure 4D).
It is suggested that the m6A methylation modification in pattern 2 and pattern 3 may affect the post-transcriptional translation of RNA in severe asthma. Compared with cluster1 and cluster3, the results were still related to immunity. The first four items of GO analysis results were immunoglobulin complex, circulating, TFIID class transcription factor complex binding, mitochondrial RNA processing and immunoglobulin complex ( Figure 4E). The results of KEGG pathway showed that the first four items were primary immunity, ribosome biogenesis in eukaryotes, ribosome and aminoacyl tRNA biosynthesis ( Figure 4F).

| Infiltration characteristics of immune microenvironment
In order to investigate the differences of immune microenvironment characteristics between severe asthma and health control and between different m6A modification patterns, we evaluated the expression of infiltrating immune cell genes and HLA genes.
The results showed that the infiltration level of infiltrating immune cells in severe asthma group was higher than that in healthy con-  Figure 5F). We also found that many immune cells are co-regulated by the m6A regulator, such as eosinophils, which are closely related to severe asthma, the abundance of eosinophils is positively correlated with YTHDF3 and negatively correlated with EIF3B, which proves that the expression of YTHDF3 and EIF3B has a close influence on eosinophils in the pathogenesis of severe asthma.
Similarly, we also analysed the correlation between HLA genes and m6A DEGs, and found that HLA_DRB1-FTO was the most positively correlated pair, the most negatively correlated HLA-m6A pair was HLA_B-YTHDC2 ( Figure 5G).

| Gene expression regulatory network mediated by m6A regulators
By weighted gene co-expression network analysis (WGCNA), the weighted value of the correlation coefficient was calculated to be four, and twenty-four gene modules were identified ( Figure 6A-C).
Next, we analysed the correlation of these twenty-four gene modules with different m6A modification patterns and key m6A genes.
The results showed that cluster 2 (r = −0.58, p = 2e-38), YTHDF3 These results suggest that these genes regulated by m6A may play an important role in severe asthma by affecting these biological processes, cell components, molecular functions and pathways.

| m6A methylation peaks of YTHDF3 and YTHDC1
In order to verify the above results, we sequenced the whole methylated RNA species in the lung tissues of severe asthma mice and control mice using the previously described merip-seq method. 28 The methylation peaks of YTHDF3 and YTHDC1 were visualized according to the TDF file. It was found that the methylation peaks of key m6A regulators in severe asthma were different from those in the control group, and were mainly distributed in exons ( Figure 7A,B). In addition, it was found that both YTHDF3 and YTHDC1 transcripts had m6A enrichment in 3′UTR ( Figure 7A,B).

| DISCUSS ION
Asthma is a chronic inflammatory respiratory disease, which involves many inflammatory cells, immune cells and cell components. 1 It affects about 300 million people around the world and is expected to increase by another one-third by 2025. Severe asthma refers to the severe state of asthma, half of the cases are accompanied by First of all, we found that 16 m6A regulators were changed in severe asthma group compared with healthy control group, which indicated that m6A regulators were involved in the process of severe asthma. Through correlation analysis, it was found that they cooperated with each other significantly, which indicated that writers, regulators, which were also verified by risk score and merip-seq, the results showed that the methylation peaks of YTHDF3 and YTHDC1 in severe asthma were different from those in the control group, and were mainly distributed in exons and enriched in 3′UTR. Secondly, we identified three m6A modification patterns, and 27 m6A regulators could distinguish each pattern well, which proves the importance of m6A regulator in severe asthma again.
We also found that the expression of these three clusters was different in the groups of smoking, gender and race, which indicated that the methylation of m6A RNA might be related to the clinical phenotype of severe asthma. Thirdly, we studied the differences of immune microenvironment characteristics between severe asthma group and health control group and between different m6A modification patterns, the results showed that the immune response of severe asthma group was more active, and the pattern 2 had more infiltrating immune cells, while pattern 3 had the lowest immune response. According to the immune characteristics of different models, it can provide a theoretical basis for the classification of immune subtypes of severe asthma. This study can make us more in-depth understand the regulatory mechanism of immune microenvironment of severe asthma, which is more conducive to the clinical precise treatment. Recently, a study used this method to divide gastric cancer into three kinds of m6A modification patterns, which further understood the immune microenvironment of gastric cancer and provided the basis for precise treatment of gastric cancer. 36 Therefore, we believe that this study is meaningful for the future precise treatment of severe asthma. Eosinophilic asthma is one of the important phenotypes of severe asthma. 37 Eosinophil degranulation is involved in airway epithelial damage, and about 50%-60% of asthma airway inflammation belongs to eosinophilic type. 38,39 Therefore, we also studied the relationship between eosinophil infiltration abundance and m6A regulator, and found that it is regulated by many m6A regulators, among which there is a significant positive correlation with YTHDF3 and a significant negative correlation with EIF3B. It was proved that these two regulators played an important role in the regulation of eosinophils in severe asthma. Finally, we identified the genes under the m6A modification patterns and key m6A regulators; the expression of these genes was affected by the m6A regulators. We found that two key genes of m6A and cluster 2 were closely related to turquoise module, which indicated that the expression of genes in this module was deeply affected by the m6A regulators. Therefore, we conducted functional enrichment analysis of these genes and found that they all affected biological processes, cell components, molecular functions and pathways in varying degrees, and then played an important role in severe asthma. In conclusion, this secondary study from the U-BIOPRED programme systematically analysed the relationship between m6A regulator modification and immune microenvironment.
Combined with previous research methods 23,36,[40][41][42] and m6A detection technology-merip-seq, it has generated abundance results. Our study is the first one to analyse the m6A regulatormediated RNA methylation modification patterns and immune microenvironment infiltration characterization in severe asthma, which opens up a new direction for the research on the pathogenesis of immune-related severe asthma with m6A modification mechanism, supplements the gap in the epigenetics of severe asthma, and encourages more scholars to carry out more research in this field. Nevertheless, all our findings confirmed the effect of m6A modification on the immune characteristics of severe asthma and provide new insights into the pathogenesis of severe asthma.

ACK N OWLED G EM ENTS
We thank Dr. Jianming Zeng(University of Macau) and all the members of his bioinformatics team, biotrainee, for generously sharing their exp. This work was funded by National Key R&D Program of China (2018YFC2002500) and National Natural Science Foundation of China (82174302).

CO N FLI C T S O F I NTE R E S T
The authors confirm that there are no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of our study are openly available in the Gene Expression Omnibus (GSE69683) at https://www.ncbi. nlm.nih.gov/geo/.