Abstract Background Few studies have been reported the potential role of N6-methyladenosine (m^6A) modification in osteoarthritis (OA). We investigated the patterns of m^6A modification in the immune microenvironment of OA. Methods We evaluated the m^6A modification patterns based on 22 m^6A regulators in 139 OA samples and systematically associated these modification patterns with immune cell infiltration characteristics. The function of m^6A phenotype-related differentially expressed genes (DEGs) was investigated using gene enrichment analysis. An m^6A score model was constructed using principal component analysis (PCA), and an OA prediction model was established based on the key m^6A regulators. We used real-time PCR analysis to detect the changes of gene expression in the cell model of OA. Results Healthy and OA samples showed significant differences in the expression of m^6A regulators. Nine key m^6A regulators, two m^6A modification patterns, m^6A-related genes and two gene clusters were identified. Some m^6A regulators had a strong correlation with each other. Gene clusters and m^6A clusters have high similarity, and cluster A corresponds to a high m^6A score. Immunocytes infiltration differed significantly between the two clusters, with the m^6A cluster B and gene cluster B having more types of infiltrating immunocytes than cluster A. The predictive model can also predict the progression of OA through m^6A regulators expression. The results of real-time PCR analysis showed that the gene expression in the cell model of OA is similar to that of the m^6A cluster B. Conclusions Our study reveals for the first time the potential regulatory mechanism of m^6A modification in the immune microenvironment of OA. This study also sheds new light on the pathogenesis of OA. Keywords: m^6A regulator, immune microenvironment, osteoarthritis, RNA methylation, immunocytes Introduction Osteoarthritis (OA) is the most common joint disease worldwide, impacting 240 million individuals ([41]1, [42]2). OA is a chronic joint disease, that is mainly characterized by pain, stiffness, joint deformity and limited joint activity ([43]3). With the trend of the aging population and the obesity epidemic, this widespread disease and the resulting disability have a great impact on individuals and society ([44]4). The progression of OA is driven by a series of factors, such as gene regulation, biochemical cascades, inflammation and cellular immunity ([45]5, [46]6). However, the etiology and disease progression mechanisms of OA are still unclear, which limits the development of effective treatment ([47]7). Many studies have revealed the significance of epigenetics in disease progression ([48]8, [49]9). Epigenetic modifications such as DNA methylation, RNA modification, histone modification and noncoding RNA modification have also been widely reported ([50]10). Over 100 different types of RNA modifications have been discovered, including N1-methyladenosine (m^1A), N6-methyladenosine (m^6A), 5-methylcytosine (m^5C), and N7-methylguanosine (m^7G) ([51]11). m^6A RNA methylation occurs on approximately 20%-40% of all transcripts encoded by mammalian cells, and it is the most common type of dynamic and reversible mRNA modification ([52]12, [53]13). Abnormal m^6A methylation levels are strongly linked to the progression of cancer, musculoskeletal disorders and other diseases ([54]10, [55]14). The level of m^6A methylation is primarily determined by the role of the m^6A regulator ([56]15). Methyltransferases, demethylases, and binding proteins all play a role in m^6A modification ([57]16). The m^6A methyltransferases (writers) include WTAP, METTL3, CBLL1, and RBM15B, while demethylases (erasers) consist of ALKBH5 and FTO. YTHDC1, YTHDF1, IGF2BP1, and other binding proteins (readers) can bind to m^6A and mediate its regulatory function ([58]12, [59]17). Previous studies have mostly concentrated on the role of m^6A methylation regulators in tumor development and treatment ([60]18–[61]20). However, the number of studies of m^6A regulators in nonneoplastic diseases is also increasing ([62]21–[63]23). There are only a few studies on the mechanism of m^6A regulators and OA at present. Most of these studies clarify the mechanism of FTO and METTL3 in OA, but they are controversial. The mechanism of m^6A regulators in OA is still unclear. Recent research has revealed that m^6A regulation can mediate some potential immune regulation mechanisms and has a significant impact on adaptive immunity ([64]24, [65]25). More and more studies have focused on the effects of the immune microenvironment on diseases ([66]26, [67]27). All types of immune cells are involved in cartilage injury and repair ([68]28). However, the interaction of m^6A methylation regulators with immune cells in OA is poorly understood ([69] Figure 1 ). Figure 1. [70]Figure 1 [71]Open in a new tab The dynamic regulation of m^6A RNA methylation modification, which regulated by ‘writers’, ‘erasers’ and ‘readers’ in OA and their potential biological functions for RNA. This study aimed to systematically evaluate the mechanism of m^6A regulators in OA. We analyzed the gene expression profile of OA by bioinformatics analysis. Subsequently, to further investigate the implication of m^6A regulators on the immune microenvironment, we investigated the correlations among clustering subgroups, risk mode, and immune cell infiltration. These findings can provide a theoretical basis for the progress and treatment of OA. Materials and methods Datasets preprocess The [72]GSE48556 datasets were obtained from the Gene Expression Omnibus(GEO) database. Total mRNA was extracted from peripheral blood mononuclear cells and detected using the Illumina HumanHT-12 V3.0 expression beadchip ([73]29). To preprocess the expression value, the “Normalize Between Arrays” function from the “limma” package was utilized. Analysis of m^6A regulators between OA patients and healthy controls We collated 22 recognized m^6A methylation regulators from published literature ([74]11, [75]12, [76]15, [77]30, [78]31). The following genes were screened: m^6A readers (YTHDC1/2, YTHDF1/2/3, HNRNPC, FMR1, LRPPRC, IGFBP1/2/3, RBMX, ELAVL1, IGF2BP1), m^6A writers (METTL3, WTAP, RBM15/15B, CBLL1, KIAA1429), and m^6A erasers (FTO, ALKBH5). The Wilcoxon test was used to compare the expression of 22 m^6A regulators in OA patients and healthy controls, and the differentially expressed m^6A regulators were screened with P value < 0.05. The R package “Random Forest” and gene importance plots were used to show the score of differentially expressed m^6A regulators. A nomogram was used to predict the possibility of OA in patients based on screened m^6A regulators. Correlation of m^6A RNA methylation regulators in OA patients The correlation between m^6A regulators in OA patients was investigated by the R package “corrplot” and Spearman’s correlation analysis. The R packages “ggMarginal” and “ggplot” are used to draw the correlation plot of significantly differentially expressed m^6A regulators. m^6A modification pattern identification Unsupervised clustering analysis was used to detect diverse m^6A clusters based on the 22 distinct m^6A regulators’ expression. The number and feasibility of modification patterns were identified by the consensus clustering algorithm. To categorize OA patients into distinct subtypes, we utilized the “ConsensusClusterPlus” software (1,000 iterations and an 80% resampling rate). The m^6A-expression pattern was assessed by principal component analysis (PCA). The m^6A regulators’ differential expression in different m^6A clusters is shown in the box plot and heatmap. We analyzed the difference in the levels of cytokines among m^6A clusters. Immune cell infiltration analysis We use the ssGSEA algorithm to calculate the enrichment score of immunocytes infiltration in each sample. The gene set of infiltrating immune cells was obtained from previous studies, which included activated B cells (ABCs), type 1/2/17 T helper (Th1/2/17) cells, dendritic cells (DCs), natural killer (NK) cells and other 23 human immune cell subtypes ([79]32). The combined ssGSEA score was used to compare the levels of immunocytes infiltration in distinct m^6A clusters. The determining criterion of a significant difference was P < 0.05. The correlation between major m^6A regulators and immunocytes infiltration was determined by Spearman correlation analysis. According to the expression of m ^6A regulators, which are strongly related to immune cell infiltration, OA patients were divided into two groups. Then, the difference in immunocytes infiltration levels in different subgroups was calculated (P value<0.05). Biological enrichment analysis for distinct m^6A modification patterns We identified the DEGs among distinct m^6A subgroups and the overlapping genes were extracted (|log FC|>1, adjusted P value<0.05). GO functional analysis and KEGG pathway enrichment analysis were performed to analyze the DEGs via R packages (enrichplot, circlize, RColorBrewer, dplyr, ComplexHeatmap, and so on). A P value of 0.05 was used as the cutoff. Identification and analysis of gene clusters We used the extracted overlapping genes and the clustering algorithm to determine cluster numbers and stability. We determined the OA gene cluster based on the extracted overlapping genes and unsupervised clustering analysis. We used a box plot diagram to show the m^6A regulators as differentially expressed in different gene clusters. The degree of immunocytes infiltration and expression of interleukin-associated factors were compared between gene clusters. The method was consistent with the immune infiltration analysis of m^6A clusters. Construction of m^6A gene signature The m^6A clusters and gene clusters were identified by previous methods. Then, we constructed the m^6A gene signature by PCA. As feature scores, principal components 1 and 2 were retrieved. The score was based only on the most significantly correlated genes, while untracked genes’ contributions to other set members were weighed. The calculation method for determining the m^6A gene signature score was based on previous studies (m^6Ascore =∑ (PC1[i] + PC2[i]), i = the expression of m^6A phenotypic related genes) ([80]33, [81]34). The difference in the m^6A score in m^6A clusters and gene clusters was analyzed. Cell and cell culture The human chondrosarcoma cell line SW135 (Pricella, China) was maintained in Dulbecco’s modified Eagle’s medium(DMEM)-high glucose (Gibco, United States) containing 10% Fetal Bovine Serum (FBS, Gibco, United States) and 1%penicillin/streptomycin. The chondrogenic ATDC5 cell line (Riken Cell Bank, Japan) was cultured in DMEM/F12 (Keygen, China) containing 10% FBS and 1%penicillin/streptomycin. Before the following experiments, all the cells were maintained under standard adherent conditions at 37°C under 5% CO[2] and humidified atmosphere. Real-time PCR analysis SW1353 cells and ATDC5 cells were treated with 10 ng/ml recombinant human IL-1β (Proteintech, Rosemont, IL) for 48 h ([82]35). Total RNA was extracted from cells using RNAiso plus reagent (Takara, Japan) according to the manufacturer’s instructions. RNA (1 µg) was reverse-transcribed to complementary DNA using cDNA Synthesis Super Mix (Trans Gen Biotech, China) according to the protocol of the manufacturer. The primers used for amplification are listed in the table below ([83] Table 1 ). Green qPCR Super Mix (Trans Gen Biotech, China) was used for real-time PCR using the CFX Connect Real-Time PCR Detection System (Biorad, United States). The thermal cycling conditions were 95°C for 30 s and 42 cycles at 95°C for 5 s and 60°C for 30 s. Table 1. The sequences of primer pairs used in the study. Gene Forward Reverse H:GAPDH 5’-GGAAGCTTGTCATCAATGGAAATC-3’ 5’-TGATGACCCTTTTGGCTCCC-3’ H:Mettl3 5’-CAGCACAGCTTCAGCAGTTCC-3’ 5’-CGTGGAGATGGCAAGACAGA-3’ H:Wtap 5’-CCAAGAAGGTTCGATTGAGTGA-3’ 5’-CAGACTCCTGCTGTTGTTGCTTT-3’ H:RBM15B 5’-CACAGCGTATCTGAGGTGGAG-3’ 5’-GTTCTGGAACTTGAGGAAGGCATA-3’ H:Rbm15 5’-GTGAGCGGAGCAAGAAGTTAGG-3’ 5’-CTATAACTATGCAAGCGGCTACTG-3’ H:YTHDC1 5’-GAGGGAATTTCATAACATGGGAC-3’ 5’-ATGGTGCTGATAGTAAGGATGGTGT-3’ H:Fto 5’-GCCAGGTGCCAGTCACGAAT-3’ 5’-TGTGAGGTCAAACGGCAGAG-3’ H:HNRNPC 5’-CGCTCCATGAACTCCCGTGT-3’ 5’-GTTCTGTTACTGACCCGTACATCTC-3’ H:IGFBP1 5’-GCACGGAGATAACTGAGGAGGA-3’ 5’-TCTTGTTGCAGTTTGGCAGGTA-3’ H:IGFBP3 5’-CTCAGAGCACAGATACCCAGAAC-3’ 5’-AGGCTGCCCATACTTATCCAC-3’ M:GAPDH 5’-CCTCGTCCCGTAGACAAAATG-3’ 5’-TGAGGTCAATGAAGGGGTCGT-3’ M:Mettl3 5’-AGGACTCTGGGCACTTGGATT-3’ 5’-ATGGCAAGACGGATGGAAAC-3’ M:Wtap 5’-GCAAGATGACCAACGAAGAACC-3’ 5’-TTAACTCATCCCGTGCCATAAC-3’ M:RBM15B 5’-CACAGTGTTTCTGAGGTGGAGC-3’ 5’-GCCTTGCCGTAGCCTATCTTT-3’ M:Rbm15 5’-GTTCAAACGCTTCGGTGATGTA-3’ 5’-CACAAAGGCTACCCGCTCAT-3’ M:YTHDC1 5’-CCGATACCAGGAAGTGGACAGAC-3’ 5’-TGGTGCTGGTAGTAAGGATGGTG-3’ M:Fto 5’-TGAGGTGGAGTTTGAGTGGCTG-3’ 5’-AGAATCTCACTCCTTTGTTCCACC-3’ M:HNRNPC 5’-ACGGTTCCTCATTTGACTTGG-3’ 5’-AACACGCTGACGTTTGGAAGG-3’ M:IGFBP1 5’-GCCACGAGCACCTTGTTCAG-3’ 5’-GCAGGGCTCCTTCCATTTCTT-3’ M:IGFBP3 5’-CTCAAAGCACAGACACCCAGAA-3’ 5’-GGCGGCACTGCTTCTTCTTAT-3’ [84]Open in a new tab (H, Human; Used in the SW1353. M, Mouse; Used in the ATDC5). Statistical analysis R 4.1.3 ([85]https://www.rproject.org/), Perl 5.32.1 ([86]https://www.perl.org) and R Bioconductor packages were used to analyze the data. Statistical analyses of real-time PCR were performed using GraphPad Prism software (version 9.0). Statistical significance was assessed using the student t-test. All experiments were performed independently at least three times. All the statistical P values were bilateral, with P value less than 0.05 considered statistically significant. Results Expression of m^6A regulators in OA The profile expression data consisted of 139 samples, including 106 genetics osteoarthritis and progression (GARP) samples and 33 normal samples. The positions of 22 selected m^6A regulators on the chromosomes are marked in [87]Figure 2A . The expression levels of 22 m^6A regulators in healthy patients and OA patients are shown in [88]Figure 2B . There were significant differences in 9 m^6A regulators (METTL3, WTAP, RBM15/15B, YTHDC1, HNRNPC, IGFBP1/3, and FTO) between normal and OA samples. In OA, four m^6A regulators (YTHDC1, HNRNPC, METTL3, and WTAP) were downregulated and five m^6A regulators (RBM15, RBM15B, IGFBP1, IGFBP3, and FTO) were upregulated ([89] Figure 2C ). IGFBP3 expression was dramatically higher in the peripheral blood of OA patients (P<0.001). METTL3 and HNRNPC levels in OA patients were considerably lower than those in healthy controls (P<0.01). Thus, these m^6A regulators are critical in OA. Figure 2. [90]Figure 2 [91]Open in a new tab The expression levels of screened m^6A regulators in healthy controls and OA patients. (A) The location of 22 RNA methylation m^6A regulators on human chromosomes. (B) The expression levels of 22 m^6A RNA methylation regulators between healthy controls and OA patients (*P < 0.05, **P < 0.01, ***P < 0.001). (C) Heatmap of expression levels of 9 key m^6A RNA methylation regulators in healthy controls and OA patients (*P < 0.05, **P < 0.01, ***P < 0.001). We further analyzed the correlations among the 22 m^6A regulators in the OA group ([92] Figure 3A ). Some m^6A regulators had a strong correlation with each other ([93] Figures 3B–F : writers and readers, [94]Figure 3G : erasers and writers, [95]Figures 3H, I : erasers and readers). Figure 3. [96]Figure 3 [97]Open in a new tab Correlation of m^6A regulators in OA patients. (A) The relationship between m^6A RNA methylation regulators in OA (A fork indicates P>0.05). (B–F) The correlation among writers and readers with a strong correlation coefficient in OA. (G) The correlation among erasers and writers in OA. (H, I) The correlation among erasers and readers in OA (B–I) Comparison between regulators with strong correlation coefficient). m^6A methylation modification patterns mediated by regulators According to m^6A regulator expression, the clustering stability was analyzed with cluster numbers from 2 to 9 ([98] Figure 4A and [99]Figure 1S ). The optimal k value was determined, and k = 2 was eventually selected as the optimal cutoff ([100] Figure 4B ). PCA revealed that OA patients could be classified into two m^6A clusters that did not intersect. Thus, two is an appropriate cluster number ([101] Figure 4C ). There were 45 OA patients in cluster A and 61 in cluster B. IGFBP1, IGFBP3, RBM15B, RBM15 and WTAP were highly expressed in cluster A. FTO, HNRNPC, METTL3 and YTHDC1 were expressed at lower levels in cluster A than in cluster B ([102] Figures 4D, E ). Figure 4. [103]Figure 4 [104]Open in a new tab Identification of 2 distinct m^6A modification patterns in OA. (A) Consensus clustering matrix (CCM) with K=2. (B) Change in the area under the cumulative distribution function (CDF) curve for k from 2 to 9. (C) PCA showed that two m^6A clusters did not intersect (2 modification patterns in OA are an appropriate choice). (D) Heatmap of the expression of 9 key m^6A RNA methylation regulators in 2 m^6A clusters. (E) Box plot of the expression of 9 key m^6A regulators in the 2 m^6A modification patterns (*P < 0.05, **P < 0.01, ***P < 0.001). Immune microenvironment features related to two m^6A modification patterns The expression of immune cells was examined to demonstrate the changes in immune microenvironment features in distinct m^6A modification patterns. Between the two m^6A clusters, there were substantial variations in virtually all immune cells observed ([105] Figure 5A ). Compared with m^6A cluster A, the cluster B had a higher level of ABCs, activated CD8^+ T cells, activated DCs, eosinophils, immature B cells, immature DCs, macrophages, plasmacytoid DCs, regulatory T cells, Th1 cells and Th2 cells. Activated CD4^+ T cells, MDSCs, mast cells, monocytes, NK T cells, NK cells, neutrophils, T follicular helper cells and Th17 cells were enriched in cluster A. More immune cells are enriched into the cluster B.The correlation between the screened m^6A regulators and immunocytes was investigated. The infiltration of immunocytes was substantially linked with the levels of IGFBP1 and RBM15B ([106] Figure 5B ). OA patients were separated into two groups based on the difference in the expression of IGFBP1 or RBM15B. Although 9 types of immunocytes were enriched in the 2 IGFBP subgroup, the high IGFBP expression group may have a higher level of immunocytes infiltration ([107] Figure 5C ). In addition, the high RBM15B expression group had fewer kinds of immune cells ([108] Figure 5D ). Figure 5. [109]Figure 5 [110]Open in a new tab Immune microenvironment features in two m^6A modification patterns. (A) Levels of infiltrating immunocytes in two m^6A clusters in OA (*P < 0.05, **P < 0.01, ***P < 0.001). (B) Correlation between 9 key m^6A regulators and immune cells. (C) The effect of high and low expression of IGFBP1 on immune infiltration (*P < 0.05, **P < 0.01, ***P < 0.001). (D) The effect of high and low RBM15B expression on immune infiltration ns, P>0.05, No statistical difference; (*P < 0.05, **P < 0.01, ***P < 0.001). Biological characteristics of m^6A modification patterns Between m^6A clusters A and B, 303 genes (m^6A phenotype-related genes) were differentially expressed. GO enrichment analysis was used for analysis of enriched biological processes, molecular functions, cellular components. These genes are mostly associated with ameboidal-type cell migration (biological process), basal plasma membrane (cell component), and lyase activity (molecular function) ([111] Figures 6A, B ). To investigate the relevant activities and pathways of m^6A phenotype-related DEGs, we employed KEGG pathway enrichment analysis ([112] Figure 6C ). Figure 6. [113]Figure 6 [114]Open in a new tab Biological characteristics related to m^6A modification patterns. (A) GO analysis of m6A phenotype-related DEGs from three perspectives: biological process, cellular composition and molecular function. (B) The circle diagram of gene enrichment numbers for each GO item. (C) KEGG pathway enrichment analysis of m^6A phenotype-related DEGs. m^6A phenotype-related DEGs in OA The 303 m^6A phenotype-associated DEGs were analyzed by unsupervised clustering analysis. It was most appropriate to divide the sample into two gene clusters ([115] Figure 2S ). There were considerable discrepancies between the two gene clusters, and only one m^6A methylation-related gene was highly expressed in cluster B ([116] Figure 7A ). The evident difference in m^6A phenotype-related gene expression between the two gene clusters was significant in controlling the establishment of immune cell infiltration ([117] Figure 7B ). The m^6A regulators also had significantly different expression in different gene clusters ([118] Figure 7C ). However, these assays were unable to predict the m^6A methylation pattern in particular individuals. We developed an m^6A scoring system to quantify the pattern of m^6A modifications in each OA patient. m^6A cluster A had a considerably higher m^6A score than m^6A cluster B. Furthermore, the m^6A score was considerably higher in gene cluster A. The m^6A score in gene cluster A was significantly higher than that in gene cluster B ([119] Figure 7D ). An alluvial plot was used to display the characteristic variations in each patient with OA. The results of m^6A regulators typing were similar to those of genotyping ([120] Figure 7E ). We also analyzed the differential expression of immune regulatory genes and immune checkpoints in m^6A clusters and gene clusters. There are significant differences in the expression of immune regulatory genes and immune checkpoints between A and B ([121] Figures 7F, G ). Figure 7. [122]Figure 7 [123]Open in a new tab Construction of m^6A signatures. (A) Heatmap of m^6A phenotype-related in different gene clusters. (B) The immunocytes infiltrating levels in distinct gene clusters (*P < 0.05, **P < 0.01, ***P < 0.001). (C) The expression levels of 9 key m^6A regulators in different gene clusters (*P < 0.05, **P < 0.01, ***P < 0.001). (D) The m^6A score in different gene clusters and m^6A clusters. (E) Alluvial diagram showing the changes in m^6A clusters, gene clusters, and m^6A score. (F) Differences in the expression of immune regulatory genes in different gene clusters and m^6A clusters (*P < 0.05, **P < 0.01, ***P < 0.001). (G) Differences in the expression of immune checkpoints in different gene clusters and m^6A clusters (*P < 0.05, **P < 0.01, ***P < 0.001). Predictive model in OA We compared the advantages of the random forest model (RFM) and support vector machine model (SVMM). The RFM was better than the SVMM ([124] Figures 8A–C ). The optimal cutoff point was selected for analysis according to the random forest tree model. The importance score of m^6A regulators was obtained ([125] Figures 8D, E ). An OA prediction model was established and evaluated. The calibration curve, clinical influence curve and decision curve all showed that the model was accurate ([126] Figures 8F–H ). The predictive model can predict the incidence of OA disease by score ([127] Figure 8I ). Figure 8. [128]Figure 8 [129]Open in a new tab Predictive model in OA. (A) Boxplot of RFM and SVMM residuals (the smaller the number, the better the model). (B) Reverse cumulative distribution map of RFM and SVMM residuals (The smaller the number is, the better the model is). (C) Evaluation of the accuracy of RFM and SVMM by receiver operating characteristic (ROC) diagram (the larger the area is, the higher the accuracy is). (D) The number of optimal trees in the verification error analysis of the random forest tree map (selecting the tree value with the lowest error). (E) The importance score of m^6A regulators. (F) The correction curve of the prediction model (the closer to the ideal dotted line, the more reliable the result). (G) Clinical impact curve to judge the accuracy of the prediction model (red line for the model to predicted high-risk patients, blue dotted line for actual high-risk patients). (H) The accuracy of the decision curve detection model (the farther the end point of the red line is from the gray line, the higher the accuracy is). (I) OA prediction model based on m^6A regulators. Expression of m^6A regulator in two chondrocyte lines after the treatment of IL-1β After we found that m^6A regulators play an important role in OA, we established an in vitro cell model of OA by treating SW1353 cells and ATDC5 cells with IL-1β. The differences in m^6A regulator gene expression before and after cell treatment were assessed. After treatment of IL-1β, three m^6A regulators (IGFBP3, RBM15, and WTAP) were down-regulated and two m^6A regulators (METTL3 and YTHDC1) were up-regulated in ATDC5. In SW1353 cells after treatment with IL-1β, the expression of FTO, HNRNPC, and IGFBP1 was down-regulated, and the expression of WTAP was up-regulated ([130] Figures 9A, B ). It appeared that the gene expression in the OA model established with cell lines is similar to that of the m^6A cluster B. These may be related to the accumulation of macrophages in the m^6A cluster B, and IL-1 is mainly secreted by macrophages. Figure 9. [131]Figure 9 [132]Open in a new tab Expression of m^6A regulator in cell model of OA. (A) Column diagram of m^6A regulator in ATDC5 before and after treatment by IL-1β (*P<0.05, **P<0.01, ***P<0.001). (B) Column diagram of m^6A regulator in SW1353 before and after treatment by IL-1β ns, P<0.05, No statistical difference,*P>0.05, **P>0.01, ***P>0.001. Discussion OA is a prevalent degenerative joint disease that can be painful and uncomfortable ([133]36). Increasing evidence shows that m^6A modification is critical for inflammation and innate immunity by interacting with a variety of m^6A regulators. Abnormal modification of the m^6A gene may lead to disorders of important genes and dynamic balance, resulting in disease ([134]37). In recent years, there have been many discussions about the relationship between inflammatory factors and the pathogenesis and progression of OA. OA is now regarded as an inflammatory illness defined by inflammatory factors rather than a degenerative disease ([135]38). The inflammatory environment is critical to the progression of OA ([136]39). At present, researchers are interested in the association between m^6A modifications and OA. Thus, future studies will focus on the link between epigenetic regulation and inflammatory factors in OA. We systematically investigated the mechanism of the modification mode of m^6A in the immune microenvironment of OA. First, we discovered that the expression of the majority of m^6A regulators changed between healthy controls and OA patients, indicating that m^6A regulators were implicated in the occurrence and progression of OA. Previous studies have indicated that METTL3 plays a significant role in the pathogenesis of osteoarthritis. Liu’s study demonstrated that METTL3 regulates OA pathogenesis by promoting inflammation and extracellular matrix (ECM) synthesis ([137]7). Jiangdong Ren showed that METTL3-mediated LINC00680 accelerates OA ([138]40). However, Sang. W showed that overexpression of METTL3 leads to reduced inflammatory cytokine levels and regulates the TIMP/MMP balance ([139]14). There are a few studies on the relationship between FTO and OA. However, there is also controversy between these studies ([140]41, [141]42). In our study, both FTO and METTL3 were found to be key m^6A regulators in OA. Our study shows that METTL3 is downregulated and FTO is upregulated in OA. We further classified OA patients according to the screening results for 9 key m^6A regulators. Furthermore, several m^6A regulators exhibited expression correlations with each other, which revealed the m^6A modification regulatory network. Third, the impact of m^6A modification patterns on immunocytes infiltration was determined to strengthen the knowledge of the interaction between m^6A RNA and the immunological response. We found that IGFBP1 and RBM15B were strongly correlated with infiltrating immunocytes in OA. Lange-Brokaar’s study showed that the immunocytes participating in cartilage injury and repair mostly comprise T cells, B cells, NK cells, DCs and macrophages ([142]43). Our study also showed that these immunocytes play a significant role in distinct m^6A modification patterns. The m^6A cluster B had more kinds of infiltrating immunocytes than cluster A, and the degree of immune cell infiltration could not be judged. However, these validated the accuracy of our immunophenotypic categorization through key m^6A regulators. It also showed that there is a certain correlation between the m^6A regulator and immune infiltration. Different immune microenvironments has different effects on osteoarthritis, which also provides a basis for immune research in osteoarthritis. Furthermore, the m^6A phenotype-related DEGs and their biological functions were identified. Our study confirmed that m^6A phenotype-related DEGs were significantly related to immune cell infiltration. We further found that the different expression levels of immune regulatory and immune checkpoints in different m^6A clusters and gene clusters. This finding may provide new ideas for the study on treatment of OA. We set up an m^6A score model to assess m^6A modifications in individual OA patients. We also found that the m^6A clusters and gene clusters contained almost the same sets of patients. Then, we established a prediction model that can predict the occurrence and progression of OA. We can score patients according to the expression of the m^6A regulator, and then get the probability of OA by score. Finally, we tried to use the OA cell model to verify whether the changes in the m^6A gene were consistent with our bioinformatics analysis. The results of real-time PCR analysis were not consistent with the differences in m^6A regulator expression between OA patients and healthy patients. This may be due to the fact that the cell model was different from the human body and that the samples came from the peripheral blood of the patient. But the change of m^6A regulators’ expression proved that these genes play a role in the progression of OA. More importantly, we found the change in gene expression in the OA cell model was almost consistent with that of the m^6A cluster B. We believed that this phenomenon is due to the fact that we used IL-1β for cell modeling, and the macrophages that mainly secrete IL-1 were mainly enriched in the B cluster. The study of m^6A regulators is widely used in the field of oncology, and better treatment can be achieved through molecular subtyping ([143]44, [144]45). In the field of OA, however, few studies have focused on the mechanism of m^6A regulators. We are the first to systematically analyze the mechanism of m^6A in OA. We proved that m^6A modification has a role in regulating the OA immune microenvironment. Limitations of the study These results will provide new inspiration for the study of the pathogenesis and treatment of OA. However, there are still some deficiencies in our research. First, our study is based on bioinformatics analysis, and although our results are accurate and reliable in theory, they need to be verified by more experiments. In the future, we hope that single-cell RNA seq can be performed to better understand how m^6A modification affects immune cell infiltration in OA patients. In addition, the expression profile dataset ([145]GSE48556) contains only 139 samples, and a larger sample size would be more beneficial for bioinformatics analysis. Therefore, we hope that more OA expression profile data and studies with larger sample sizes will become available. Third, because OA is a nonneoplastic disease, its survival curve is rarely studied, so it is impossible to establish a model considering the m^6A score and survival. However, the severity of knee joint damage can be studied, and more clinical data will establish the relationship between the m^6A score and disease progression. In this way, m^6A score will be critical in predicting the progression of the disease. But the accuracy of m^6A scoring model need validated by a large number of patient samples. Finally, there were few genes in the GO and KEGG analyses of m^6A phenotype-related DEGs, which could be due to the limited sample size. However, some of the enriched biological functions mentioned in our study have been proven to participate in OA by other studies, which indicates that our results are worth using for reference. We hope that the expression profile datasets with a larger sample size become available. Conclusion In summary, our study reveals for the first time the potential regulatory mechanism of m^6A methylation modification in the immune microenvironment of OA. This study also sheds new light on the pathogenesis of OA. The difference in m^6A modification mode is a significant contributor to the complexity of the OA microenvironment. Individual OA patients can utilize the m^6A score to assess their m^6A modification pattern and immune cell infiltration features. The prediction model we have established is helpful in predicting disease development in OA patients. Additional studies are needed to fully clarify the molecular mechanisms of m^6A regulation and biological function during OA. Data availability statement Publicly available datasets were analyzed in this study. This data can be found here: The [146]GSE48556 datasets were obtained from the GEOdatabase([147]https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE 48556). Author contributions GS and YO designed this study. TW, HF and SC conducted literature searches and management. YO, YT and HM were responsible for data management and statistical analysis. GZ, ZW, WY interpreted the findings. YO performed the manuscript writing. The manuscript was revised by all authors, and the final version was approved by all. All authors contributed to the article and approved the submitted version. Funding This work was supported by the National Natural Science Foundation of China (NO.81960881), which was hosted by GS. The funders were not involved in the study design, data collection, statistical analysis or writing of the manuscript. The [148]GSE48556 datasets were obtained from the GEO database ([149]https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE48556). Conflict of interest The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. Publisher’s note All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher. Supplementary material The Supplementary Material for this article can be found online at: [150]https://www.frontiersin.org/articles/10.3389/fimmu.2022.1018701/fu ll#supplementary-material Figure 1S The best m6A cluster number was evaluated by the CCM. (A–G) K=3-9 CCM. [151]Click here for additional data file.^ (527.6KB, tif) Figure 2S The best gene cluster number was evaluated by the CCM. (A–I) K=2-9 CCM. (J) Relative change in the area under the CDF curve for k from 2 to 9. [152]Click here for additional data file.^ (418.5KB, tif) References