Significance of m6A regulatory factor in gene expression and immune function of osteoarthritis

One of the most prevalent posttranscriptional modifications of eukaryotic mRNA is the RNA N6-methyladenosine (m⁶A) regulator, which plays a significant role in various illnesses. The involvement of m⁶A regulators in osteoarthritis (OA) is not fully known. By comparing nonosteoarthritic and osteoarthritic patients, 26 important m⁶A regulators were identified from the gene expression omnibus GSE48556 dataset. Seven candidate m⁶A regulators (IGFBP3, WTAP, IGFBP1, HNRNPC, RBM15B, YTHDC1, and METTL3) were screened using a random forest model to assess the likelihood of OA. A column line graph model founded on seven m⁶A modulator candidates was created. According to decision curve analysis, patients might profit from the column line graph model. Based on chosen relevant m⁶A modifiers, a consensus clustering approach was utilized to categorize OA into two m⁶A categories (group A and group B). To measure the m⁶A pattern, a principal component analysis technique was created to generate the m⁶A score for every sample. Cluster A patients exhibited more excellent m⁶A scores than cluster B patients. Furthermore, we discovered that patients with lower and higher m⁶A scores had varied immunological responses using the m⁶A type. At last, m⁶A regulators contribute significantly to the progression of OA. Our research on m⁶A patterns might help to guide further OA immunotherapeutic techniques.

Introduction

The most prevalent musculoskeletal illness, osteoarthritis (OA), is characterized by localized joint space restriction, osteophyte, subchondral bone sclerosis, and subchondral cysts (Ariani et al., 2019). Because of the considerable advancement in OA studies, it is generally thought that OA is a complicated degenerative disease influenced by various variables, including age, obesity, inflammation, trauma, and genetics (Yucesoy et al., 2015). In recent years, there have also been studies that bone arthritis is a low-grade inflammatory state closely linked to some inflammatory mediators. Patients with OA have high amounts of Interleukin (IL)-1, IL-10, IL-17, and tumor necrosis factor (TNF) in their blood and synovial fluid (Hussein et al., 2008). Other studies also found that the leukocytes and monocytes necessary to regulate the immune response are more marvelous and proinflammatory in OA patients, and inflammation and body mass index are related to this (Loukov et al., 2018). Therefore, early screening and effective prevention of OA patients from an immunological point of view will profoundly influence the control of OA.

N6-methyladenosine (m⁶A), 5-methylcytidine (m⁵C), N6-2′-O-dimethyladenosine (m⁶Am), pseudouridine (Ψ), 1-methyladenosine (m¹A), and 5-hydroxymethylcytidine (hm⁵C) have all been studied extensively in latest days (Jonkhout et al., 2017). Among them, m⁶A is the most popular, abundant and well-characterized internal modification in messenger RNA (mRNA). The N6-methylation of altered RNA may influence its fold, durability, degradation, and cellular interaction(s), involving it in splicing, translation, export, and decay processes (Oerum et al., 2021). Recent studies have reported that m⁶A RNA methylation can regulate promoter–proximal pausing of RNA polymerase II, and m⁶A RNA methylation mediated in methyltransferase-like 3 (METTL3) can promote antitumor immunity in natural killer cells (Akhtar et al., 2021; Song et al., 2021). m⁶A has been shown to have a critical role in a variety of disorders, including cancer (Liu et al., 2018; Sun et al., 2019), asthma (Dai et al., 2021), osteoporosis (Wu et al., 2018), and Alzheimer’s disease (Han et al., 2020). And METTL3–mediated m⁶A modification of ATG7 regulates autophagy and can promote OA progression (Chen et al., 2022). Nevertheless, the specific role of the m⁶A regulatory factor in OA immunity is unclear.

This research synthesized the involvement of the m⁶A regulatory factor in the identification and subtype categorization of OA using the GSE48556 dataset from the GEO database. Based on seven candidate m⁶A regulators, we created a gene model for predicting OA susceptibility [Insulin-Like Growth Factor Binding Protein 3 (IGFBP3), Wilms tumor 1–associating protein (WTAP), heterogeneous nuclear ribonucleoproteins (HNRNPC), Insulin-Like Growth Factor Binding Protein 1 (IGFBP1), RNA-binding motif protein 15B (RBM15B), YTH domain–containing protein 1 (YTHDC1), and METTL3], and predict the prevalence of OA from the perspective of m⁶A regulation. In addition, we compared two types of immune responses, which provides a new idea for the treatment of OA (Figure 1 flowchart).

FIGURE 1

FIGURE 1. Flow chart.

Materials and methods

Data gathering

The GEO database (http://www.ncbi.nlm.nih.gov/geo/) was used to pick a dataset (GSE48556) of 106 and 33 non-OA healthy adults. Through peripheral blood gene analysis between OA patients and non-OA patients, we extracted 26 m⁶A regulators (Supplementary File S1) from the dataset. The following are the criteria we use to choose datasets: first, the sample size must be larger than 100; second, the data set must include both a treatment and a healthy control; third, in both the patient and control groups, there must be at least 10 cases in each. These regulators consisted of nine writers (METTL3, METTL14, METTL16, WTAP, VIRMA, RBM15, RBM15B, CBLL1, and ZC3H13), two erasers (FTO and ALKBH5), and 15 readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, IGFBP1, IGFBP2, IGFBP3, RBMX, ELAVL1, and IGF2BP1) (Dai et al., 2021; Jiang et al., 2022).

Creating the random forest model and the support vector machine model

As a training model, create a random forest (RF) and support vector machine (SVM) method to forecast the incidence of OA. To assess the model, plot the reverse summary statistics of residuals, boxplots of residuals, and receiver operating characteristic (ROC) curves. The RF approach is a regression tree method that permits bootstrap aggregation and predictor randomization to obtain a high level of predicted accuracy (Rigatti, 2017). In our work, an RF model was constructed using the “Random Forest” package of the R statistical program (Version 4.0.4). Potential m⁶A regulators were chosen from a pool of 26 m⁶A regulators to forecast the incidence of OA. In our investigation, ntrees and mtry were adjusted to 100 and 3, correspondingly. The relevance of 26 m⁶A regulators was then assessed, and the suitable important m⁶A regulators were chosen using 10x cross-validation. The y-axis of the 10x cross-validation curve correlates to the model’s accuracy when the number of m⁶A regulators is changed. SVM is a supervised machine learning method based on the statistical learning theory of structural risk reduction. In our investigation, each data point was represented as a point in n-dimensional space (where n is the number of m⁶A regulators). We then discovered an ideal hyperplane that effectively distinguishes between the two groups.

Establishment of nomogram model

We utilized the “RMS” tool in R to forecast the occurrence of OA using a nomogram model depending on the chosen seven candidate m⁶A regulators. The calibration graph is also used to evaluate the projected value to the actual value, followed by an analysis of the decision curve and a clinical impact curve to see whether the model is helpful for illness prediction in patients (Iasonos et al., 2008).

Recognition of m⁶A pattern

To detect distinct m⁶A patterns, the “ConsensusClusterPlus” software program in R is used (Wilkerson and Hayes, 2010), and the “limma” package in R is used to search differentially expressed genes (DEGs) between m⁶A trends. p < 0.01 was used as a screening criterion. The “ConsensusClusterPlus” software program in R software is used to do functional enrichment analysis to better understand the various mechanisms of DEGs involvement in OA, and the findings are shown in an enrichment circle diagram.

Quantification of m⁶A patterns and assessment of immune cell abundance

We utilize the principal component analysis (PCA) technique to produce the m⁶A score, which is calculated as follows: PC1i = m⁶A score, where PC1 is principal component 1, and i is the DEG expression (Zhang et al., 2020). Gene set enrichment analysis (ssGSEA) is also used to evaluate the gene expression patterns in the samples. Then these genes were searched in the input data to quantify the presence of these genes in immune cells in every sample.

Statistical analysis

Correlations among writers, erasers, and reader are discovered using linear regression analysis. Then, the residuals were tested in the linear regression via quantile–quantile plots (Q–Q plot). Kruskal–Wallis analysis is used to assess component differences. All parametric analyses are defined as two tests, with p < 0.05 considered statistically significant. R version 4.0.4 was used for all statistical studies.

Results

The landscape of the 26 m⁶A regulators in osteoarthritis

The R software package “limma” was used to compare the expression levels of m⁶A regulators in patients with and without OA. Seven important m⁶A regulators (IGFBP3, RBM15B, IGFBP1, METTL3, YTHDC1, HNRNPC, and WTAP) were screened and visualized using heatmaps and histograms. We discovered that RBM15, RBM15B, and YTHDF1 were abundantly expressed in OA patients, but the other critical m⁶A regulators showed lower or equivalent expression in OA cases in contrast to individual patients (Figures 2A, B). Use the “RCircos” software tool to see the chromosomal locations of m⁶A regulators (Figure 2C).

FIGURE 2

FIGURE 2. m⁶A modulators in osteoarthritis (OA). (A) Heat maps of the expression of seven important m⁶A regulators in patients with nonosteoarthritis and OA. (B) Histogram of the differences in the expression of 21 m⁶A regulators identified between subjects with nonosteoarthritis and OA. (C) The chromosomal position of the m⁶A regulator. *p < 0.05, **p < 0.01, and ***p < 0.001.

Correlation among the writer, eraser, and readers in osteoarthritis

To explore the relationship between writer, eraser, and reader in OA, we used linear regression analysis to examine the correlation one to one. Then, we tested the residuals in the linear regression via Q–Q plot (Supplementary Figure S1). We found only significant differences between eraser and reader and between writer and reader. HNRNPC and LRPPRC expression levels in OA patients were substantially linked with FTO. The expression levels of YTHDF2 and FMR1 are highly negatively correlated with FTO. High expression of YTHDF3, HNRNPC, LRPPRC, YTHDC1, and YTHDC2 was positively correlated with METTL3, whereas patients with high expression of FMR1 had lower levels of METTL3 expression. Writer and eraser have no significant correlation (Figure 3).

FIGURE 3

FIGURE 3. Correlation among writers, erasers, and readers in OA (A–J). Reader genes: YTHDF2, FMR1, HNRNPC, LRPPRC, YTHDF3, YTHDC1, and YTHDC2; writer and eraser genes: METTL3 and FTO.

Creating the random forest model and support vector machine model

To predict the incidence of OA, we developed RF and SVM models. The RESIDUE of the RF model is negligible, according to both the “inverse cumulative distribution of residues” (Figure 4A) and the “remnant boxline plot” (Figure 4B). The majority of the samples in the model have tiny residuals, suggesting that the model is superior. As seen by the ROC curve, RF is more reliable than SVM (Figure 4C). The correlation plot between model uncertainties and the number of selection trees (Figure 4D) revealed a consistent inaccuracy; thus, we picked 500 trees as the variable of the final model. Throughout the creation of the RF model, the variable significance of the output results (Gini coefficient technique) was examined in the context of declining accuracy and reducing mean square error. Figure 4E reveals that seven m⁶A regulators (IGFBP3, WTAP, IGFBP1, HNRNPC, RBM15B, YTHDC1, and METTL3) are more important than five points.

FIGURE 4

FIGURE 4. Random forest(RF) model and support vector machine model. (A) The residual distribution of the RF and support vector machine (SVM) models was displayed using the reverse cumulative distribution of residuals. (B) The residual distribution of the RF and SVM models was displayed using boxplots of residuals. (C) The receiver operating characteristic curves showed how accurate the RF and SVM models were. (D) The tree model demonstrates that the error has stabilized. (E) The seven m⁶A regulators’ relative relevance.

Construction of the nomogram model

To predict the occurrence of OA patients, a line graph model was built in R using the “rms” package against seven potential m⁶A regulators (Figure 5A). It can be seen from the prediction curve that the column chart model is accurate (Figure 5B). The red line stays above the gray and black lines between 0 and 1, as shown in the Decision Curve Analysis (DCA) curve, indicating that this column line chart model may be advantageous for patients with OA (Figure 5C). The clinical impact curve shows that the line plot model has excellent predictive power (Figure 5D).

FIGURE 5

FIGURE 5. The nomogram model is created. (A) The nomogram model was built using the m⁶A regulators of the seven candidates. (B) The calibration curve reveals the nomogram model's predictive capacity. The calibration curve reveals the nomogram model’s predictive capacity. (C) Patients with OA may benefit from relying on the nomogram model. (D) The clinical relevance curve is used to analyze the nomogram model’s clinical effect.

Two different m⁶A modes distinguished by important m⁶A regulators

Two m⁶A patterns (A and B) were identified using the “ConsensusClusterPlus” tool of R software for the seven m⁶A regulators. The difference in delta area is seen in Figures 6A, B because the area tends to remain stationary when k = 2; this is the best k-value to use. When k = 2, the cluster research findings are the most trustworthy, and the resultant classification is the greatest. Figure 6C displays the average values when k assumes various values. Cluster A contained 41 cases, and cluster B contained 65 cases. The differential expression values of the seven key m⁶A regulators between the two groups are then shown using heatmaps and histograms (Figure 6D). RBM15B, IGFBP1, and IGFBP3 are expressed at higher levels in group A than those in group B, whereas METTL3, WTAP, YTHDC1, and HNRNPC showed the opposite (Figure 6E). The PCA shows that seven important m⁶A regulators can fully distinguish between these two m⁶A modes (Figure 6F). A total of 639 DEGs with m⁶A correlations were chosen from the two m⁶A modes (Supplementary File S2). GO functional enrichment analysis was used to investigate the likely processes of these DEGs in osteoarthritis, and the findings were shown using an enrichment circle chart (Figure 6G). These genes were mostly enriched in the GO terms 0016829 and 0016323, which were associated with lyase activity and basolateral plasma membrane, correspondingly.

FIGURE 6

FIGURE 6. Clustering of the seven major m⁶A regulators in OA by consensus. (A) Matrix consensus of the seven major m⁶A regulators for k = 2. (B) The difference of delta area for k = 2–9. (C) The area under the CDF curve. (D) Cluster A and B expression heatmaps of the seven important m⁶A regulators. (E) Histogram showing differential expression of the seven major m⁶A regulators in clusters A and B. (F) The expression patterns of the seven key m⁶A regulators were analyzed using principal component analysis, which revealed significant differences in gene sequences between the two m⁶A patterns. (G) The influence of the 639 m⁶A-related differentially expressed genes (DEGs) on the onset and progression of OA is investigated using gene ontology analysis. *p < 0.05, **p < 0.01, and ***p < 0.001.

The quantity of immune cells in OA samples was then calculated using ssGSEA, and the association between seven key m⁶A regulators and immune cells was assessed. We discovered that the m⁶A methyltransferase RBM15B is linked to a variety of immune cells (Figure 7A). The variations in immune cell infiltration between individuals with elevated/low RBM15B expression were investigated. Immune cell infiltration has been found to be lower in subjects with elevated RBM15B expression compared with patients with low RBM15B expression. (Figure 7B). At last, we looked at the differences between the two immune cell infiltration m⁶A patterns. Cluster A was shown to be linked with B-cells, CD8 T-cells, as well as dendritic cells, whereas cluster B was discovered to be related with macrophages and neutrophils (Figure 7C).

FIGURE 7

FIGURE 7. Gene collection from a single sample and enrichment analysis. (A) Correlation between immune cells and seven key m⁶A modulators found in different species. (B) Differences in the number of infiltrative immune cells between groups with lower and higher RBM15B expression. (C) Immune cell infiltration differences between clusters A and B. *p < 0.05, **p < 0.01, and ***p < 0.001.

Differentiation of two different m⁶A–related gene patterns and characteristics

Using a consistent clustering algorithm, individuals with OA were categorized into several genetic subgroups based on m⁶A–associated DEGs. We discovered two distinct m⁶A–related gene clusters (gene group A and gene group B), which correspond to the m⁶A pattern grouping (Figure 8A). Figure 8B shows the delta area of the curve. When k = 2, the area under the curve tends to be stationary; thus, 2 is the most appropriate k-value. Figure 8C depicts the cumulative distribution function as k varies. When k = 2, the cumulative distribution function (CDF) achieves an approximate maximum, cluster analysis findings are the most trustworthy, and the resultant classification is the best. DEGs in gene clusters A and B have distinct degrees of expression (Figure 8D). Infiltration of immune cells between gene clusters A and B is comparable with that seen in the m⁶A pattern (Figures 8E, F). This confirms that the consensus clustering technique grouping is valid once again. When we compare the m⁶A scores of the two m⁶A patterns, we can observe that group A or gene group A has a higher m⁶A score than group B or group B (Figures 8G, H).

FIGURE 8

FIGURE 8. Clustering of the DEGs in two m⁶A clusters by consensus. (A) For k = 2, consensus matrices of the seven major m⁶A regulators. (B) The delta area difference for k = 2–9. (C) The CDF curve’s area under the curve. (D) Heatmap of DEG expression in gene clusters A and B. (E) Histogram showing differential expression of the seven major m⁶A regulators in gene clusters A and B. (F) Levels of immune cell infiltration in gene groups A and B. (G) M⁶A score differences between gene clusters A and B. (H) The m⁶A score differences between gene clusters A and B. *p < 0.05, **p < 0.01, and ***p < 0.001.

Role of m⁶A patterns in distinguishing osteoarthritis

The Sankey plot depicts the association between the m⁶A pattern, the m⁶A gene pattern, and the m⁶A score (Figure 9A). IL-4, IL-5, IL-13, thymic stromal lymphopoietin (TSLP), and IL-33 were used to investigate the m⁶A pattern’s relationship to OA. The findings revealed that TSLP, IL-4, IL-5, and IL-13 expression levels in cluster A were considerably greater than those in cluster B (Figures 9B, C).

FIGURE 9

FIGURE 9. Influence of m⁶A pattern in OA. (A) A Sanchi plot illustrating the link between m⁶A pattern, m⁶A gene pattern, and m⁶A score. (B) Different transcriptional activity of thymic stromal IL-4, IL-5, IL-13, lymphopoietin (TSLP), and interleukin (IL)-33 clusters A and B. (C) Different transcriptional activity of IL-4, IL-5, IL-13, TSLP, and IL-33 between gene clusters A and B. *p < 0.05, **p < 0.01, and ***p < 0.001.

Discussion

OA is a degenerative joint disease that most often affects the elderly. It is marked by gradual cartilage deterioration, as well as synovial membrane inflammation, osteophyte production, and subchondral bone sclerosis (Charlier et al., 2019). Although the m⁶A regulator has been linked to a variety of biological processes, its involvement in OA remains unknown. The purpose of our study was to see whether m⁶A regulators were involved in OA patients.

Through differential expression analysis in individuals with non-OA and OA, we initially identified seven significant m⁶A regulators out of 26 m⁶A regulators. To forecast the incidence of OA, an RF model was built, and seven prospective m⁶A regulators (IGFBP3, WTAP, IGFBP1, HNRNPC, RBM15B, YTHDC1, and METTL3) were chosen among 26 m⁶A regulators. Then, using seven candidate m⁶A regulators, a line chart model was created, and the DCA curve revealed that this column line chart model might be effective for patients with OA.

Among these seven important m⁶A regulatory factors, IGFBP1 and IGFBP3 are proteins that regulate the growth and proliferation of somatic cells which play a vital role in the secretion and role of growth hormone (Ranke, 2015). They can increase the stability of mRNA to regulate mRNA expression. It has been shown that IGFBP1 and IGFBP3 are associated with extensive extracellular matrix loss, cell death, and the formation of an inflammatory environment in terms of initiating the typical OA response pattern, providing new ideas for the assessment and treatment of OA (Schlichting et al., 2014). Another research identified chondrocyte enlargement as a detrimental event in the development of OA, which was significantly averted when chondrocytes were knocked out of IGFBP3 (Evans et al., 2015). WTAP is a nuclear protein that, like WT1, is found in speckles across the nucleoplasm and partly colocalizes with splicing factors (Little et al., 2000). WTAP may be found on nuclear plaques that are rich in pre-mRNA processing factors, which is essential for m⁶A methyltransferase catalytic activity in vivo and can also control the recruitment of the m⁶A methyltransferase complex to mRNA targets (Ping et al., 2014). HNRNPC (nuclear heterogeneous ribonucleoprotein C) is a widely expressed RNA-binding protein that aids in the shearing of introns and the proper assembly of exons during RNA splicing (Liu et al., 2015). It has been shown that the HNRNP family can cause degradation of the extracellular matrix and cartilage loss associated with it, which in turn can initiate OA response patterns (Ruedel et al., 2013). The WTAP-METTL3 m⁶A methyltransferase complex contains RBM15, which is an RNA-binding protein. WTAP interacts with RBM15B and related RBM15 to recruit complexes to target mRNAs and cleavage factors like as SF3B1 to promote selective shearing (Patil et al., 2016). YTHDC1 controls mRNA splicing by attracting and modifying pre-mRNA splicing factors, allowing them to reach specific mRNAs’ critical areas (Xiao et al., 2016). The best-known m⁶A methyltransferase, METTL3, is involved in the reversible epi-transcriptome control of m⁶A modification (Liu et al., 2020). Some studies have shown that intra-articular injection of synovially targeted METTL3 siRNA inhibits the propagation of cellular senescence in the joint and facilitates destabilisation of the medial meniscus (DMM)-induced cartilage destruction (Chen et al., 2022).

The complexity of OA creates obstacles to its effective treatment. At present, some studies have suggested that the intervention of these factors that affect the development and differentiation of Th17 cells, such as mTORC1/mTORC2, may lead to the regulation of OA (Ren et al., 2016; Dar et al., 2018). Some researchers now believe that adaptive immunity and OA are related (Geyer and Schonfeld, 2018). The adaptive immune response is the procedure through which antigen-specific T/B lymphocytes are activated, increased, and differentiated into effector cells in response to antigenic stimulation in the body, resulting in various biological consequences. T-lymphocytes (T-cells) orchestrate the elimination of infections, whereas B lymphocytes (B-cells) create antibodies (Woodell May and Sommerfeld, 2019). Type 2 helper T-cells and mastocytes both release IL-4. Chondrocytes revealed lower sensitivity to the effects of IL-4 in the course of OA, which might be responsible for the fast deterioration of the articular cartilage (Millward–Sadler et al., 1999). IL-13, similar to IL-4, is a mediator of the body’s immune response and is secreted by a variety of cells, such as CD4⁺ T-cells, CD8⁺ T-cells, mast cells, basophils, eosinophils, and natural killer cells. IL-13 can inhibit the proinflammatory effects of TNF alpha relative to fibroblasts from OA patients (Yeh et al., 1995; Hsieh et al., 2013). IL-5 modulates the transcription of genes encoding in B-cell proliferation, cell survival, maturation, and effector activities, as well as innate and acquired immunological responses and eosinophilia (Dougan et al., 2019). The cytokines TSLP and IL-33 signaling pathways are essential in various inflammatory and allergy reactions (Hong et al., 2020). We discovered two m⁶A groups based on seven critical m⁶A regulators using the consensus clustering approach. Type A shows a significant increase in natural killer cells and neutrophils. At the same time, there is also a significant increase in inflammatory factors, including IL4, Il5, and IL33. Type B showed significantly increased B-cells, CD8 T-cells, Th1, and Th2 cells. It is more suggested that in the two types, the innate immune response is more dominant in type A, and the adaptive immune response is more dominant in type B. Different immune response methods can better guide the clinical immunotherapy of OA. We subsequently used 639 m⁶A-associated DEGs to check the accuracy of these findings in the m⁶A gene pattern. At last, to measure the m⁶A design, we use the PCA approach to calculate the m⁶A score for each sample. Because OA is a progressive disease with different treatments at different stages, our study can distinguish patients with high scores of acquired immune response from gene methylation modification regulation. Can it provide a new idea for treating OA in the future?

Conclusion

In conclusion, we established the prediction model of OA using the m⁶A regulatory factor. From the perspective of m⁶A modification, OA patients are classified to identify the population with a robust immune response, which provides a new idea for the phased inflammatory immunotherapy of OA diseases.

Data availability statement

The datasets used during the present study are available from GEO (https://www.ncbi.nlm.nih.gov/geo/) database, accession number GSE48556.

Ethics statement

Written informed consent was not obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author contributions

CYW is the corresponding author of the article who contributes the most. She designed the whole study. YZ and XYX collaborated on data curation. Statistical analysis was carried out by FJ and CYW. The original draft of the manuscript was completed by XYX. And CYW revised the manuscript. All of the authors approved the final manuscript.

Funding

This work was funded by the Nanjing Municipal Science and Technology Bureau (Grant number of 2019060002). The funding body had no role in the study design, data collection, analysis, and interpretation of data.

Acknowledgments

We would like to thank the Core Facility of the First Affiliated Hospital of Nanjing Medical University for its help in the experiment.

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: https://www.frontiersin.org/articles/10.3389/fphys.2022.918270/full#supplementary-material

Supplementary Figure S1 | Q–Q plot.

Supplementary File S1 | 26 m⁶A regulators.

Supplementary File S2 | 639 m⁶A–correlated DEGs.

References

Akhtar J., Renaud Y., Albrecht S., Ghavi-Helm Y., Roignant J. Y., Silies M., et al. (2021). m(6 A RNA methylation regulates promoter- proximal pausing of RNA polymerase II. Mol. Cell. 81 (16), 3356–3367.e6. doi:10.1016/j.molcel.2021.06.023.e6

PubMed Abstract | CrossRef Full Text | Google Scholar

Ariani A., Manara M., Fioravanti A., Iannone F., Salaffi F., Ughi N., et al. (2019). The Italian Society for Rheumatology clinical practice guidelines for the diagnosis and management of knee, hip and hand osteoarthritis. Reumatismo 71 (S1), 5–21. doi:10.4081/reumatismo.2019.1188

PubMed Abstract | CrossRef Full Text | Google Scholar

Charlier E., Deroyer C., Ciregia F., Malaise O., Neuville S., Plener Z., et al. (2019). Chondrocyte dedifferentiation and osteoarthritis (OA). Biochem. Pharmacol. 165, 49–65. doi:10.1016/j.bcp.2019.02.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen X., Gong W., Shao X., Shi T., Zhang L., Dong J., et al. (2022). METTL3-mediated m(6)A modification of ATG7 regulates autophagy-GATA4 axis to promote cellular senescence and osteoarthritis progression. Ann. Rheum. Dis. 81 (1), 87–99. doi:10.1136/annrheumdis-2021-221091

PubMed Abstract | CrossRef Full Text | Google Scholar

Dai B., Sun F., Cai X., Li C., Liu H., Shang Y. (2021). Significance of RNA N6-methyladenosine regulators in the diagnosis and subtype classification of childhood asthma using the gene expression omnibus database. Front. Genet. 12, 634162. doi:10.3389/fgene.2021.634162

PubMed Abstract | CrossRef Full Text | Google Scholar

Dar H. Y., Azam Z., Anupam R., Mondal R. K., Srivastava R. K. (2018). Osteoimmunology: The Nexus between bone and immune system. Front. Biosci. 23 (3), 464–492. doi:10.2741/4600

PubMed Abstract | CrossRef Full Text | Google Scholar

Dougan M., Dranoff G., Dougan S. K. (2019). GM-CSF, IL-3, and IL-5 family of cytokines: Regulators of inflammation. Immunity 50 (4), 796–811. doi:10.1016/j.immuni.2019.03.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Evans D. S., Cailotto F., Parimi N., Valdes A. M., Castano-Betancourt M. C., Liu Y., et al. (2015). Genome-wide association and functional studies identify a role forIGFBP3 in hip osteoarthritis. Ann. Rheum. Dis. 74 (10), 1861–1867. doi:10.1136/annrheumdis-2013-205020

PubMed Abstract | CrossRef Full Text | Google Scholar

Geyer M., Schonfeld C. (2018). Novel insights into the pathogenesis of osteoarthritis. Curr. Rheumatol. Rev. 14 (2), 98–107. doi:10.2174/1573397113666170807122312

PubMed Abstract | CrossRef Full Text | Google Scholar

Han M., Liu Z., Xu Y., Liu X., Wang D., Li F., et al. (2020). Abnormality of m6A mRNA methylation is involved in Alzheimer's disease. Front. Neurosci. 14, 98. doi:10.3389/fnins.2020.00098

PubMed Abstract | CrossRef Full Text | Google Scholar

Hong H., Liao S., Chen F., Yang Q., Wang D. Y. (2020). Role of IL-25, IL-33, and TSLP in triggering united airway diseases toward type 2 inflammation. Allergy 75 (11), 2794–2804. doi:10.1111/all.14526

PubMed Abstract | CrossRef Full Text | Google Scholar

Hsieh J., Shiau A. L., Lee C. H., Yang S. J., Lee B. O., Jou I. M., et al. (2013). CD8+ T cell-induced expression of tissue inhibitor of metalloproteinses-1 exacerbated osteoarthritis. Int. J. Mol. Sci. 14 (10), 19951–19970. doi:10.3390/ijms141019951

PubMed Abstract | CrossRef Full Text | Google Scholar

Hussein M. R., Fathi N. A., El-Din A. M. E., Hassan H. I., Abdullah F., Al-Hakeem E., et al. (2008). Alterations of the CD4(+), CD8 (+) T cell subsets, interleukins-1beta, IL-10, IL-17, tumor necrosis factor-alpha and soluble intercellular adhesion molecule-1 in rheumatoid arthritis and osteoarthritis: Preliminary observations. Pathol. Oncol. Res. 14 (3), 321–328. doi:10.1007/s12253-008-9016-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Iasonos A., Schrag D., Raj G. V., Panageas K. S. (2008). How to build and interpret a nomogram for cancer prognosis. J. Clin. Oncol. 26 (8), 1364–1370. doi:10.1200/JCO.2007.12.9791

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang F., Hu Y., Liu X., Wang M., Wu C. (2022). Methylation pattern mediated by m(6)A regulator and tumor microenvironment invasion in lung adenocarcinoma. Oxid. Med. Cell. Longev. 2022, 2930310. doi:10.1155/2022/2930310

PubMed Abstract | CrossRef Full Text | Google Scholar

Jonkhout N., Tran J., Smith M. A., Schonrock N., Mattick J. S., Novoa E. M. (2017). The RNA modification landscape in human disease. RNA 23 (12), 1754–1769. doi:10.1261/rna.063503.117

PubMed Abstract | CrossRef Full Text | Google Scholar

Little N. A., Hastie N. D., Davies R. C. (2000). Identification of WTAP, a novel Wilms' tumour 1-associating protein. Hum. Mol. Genet. 9 (15), 2231–2239. doi:10.1093/oxfordjournals.hmg.a018914

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu N., Dai Q., Zheng G., He C., Parisien M., Pan T. (2015). N(6)-methyladenosine-dependent RNA structural switches regulate RNA-protein interactions. Nature 518 (7540), 560–564. doi:10.1038/nature14234

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu S., Zhuo L., Wang J., Zhang Q., Li Q., Li G., et al. (2020). METTL3 plays multiple functions in biological processes. Am. J. Cancer Res. 10 (6), 1631–1646.

PubMed Abstract | Google Scholar

Liu Z. X., Li L. M., Sun H. L., Liu S. M. (2018). Link between m6A modification and cancers. Front. Bioeng. Biotechnol. 6, 89. doi:10.3389/fbioe.2018.00089

PubMed Abstract | CrossRef Full Text | Google Scholar

Loukov D., KarampatoS S., Maly M. R., Bowdish D. M. E. (2018). Monocyte activation is elevated in women with knee-osteoarthritis and associated with inflammation, BMI and pain. Osteoarthr. Cartil. 26 (2), 255–263. doi:10.1016/j.joca.2017.10.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Millward-Sadler S. J., Wright M. O., Lee H., Nishida K., Caldwell H., Nuki G., et al. (1999). Integrin-regulated secretion of interleukin 4: A novel pathway of mechanotransduction in human articular chondrocytes. J. Cell. Biol. 145 (1), 183–189. doi:10.1083/jcb.145.1.183

PubMed Abstract | CrossRef Full Text | Google Scholar

Oerum S., Meynier V., Catala M., Tisne C. (2021). A comprehensive review of m6A/m6Am RNA methyltransferase structures. Nucleic Acids Res. 49 (13), 7239–7255. doi:10.1093/nar/gkab378

PubMed Abstract | CrossRef Full Text | Google Scholar

Patil D. P., Chen C. K., Pickering B. F., Chow A., Jackson C., Guttman M., et al. (2016). m(6 A RNA methylation promotes XIST-mediated transcriptional repression. Nature 537 (7620), 369–373. doi:10.1038/nature19342

PubMed Abstract | CrossRef Full Text | Google Scholar

Ping X., Sun B. F., Wang L., Xiao W., Yang X., Wang W. J., et al. (2014). Mammalian WTAP is a regulatory subunit of the RNA N6-methyladenosine methyltransferase. Cell. Res. 24 (2), 177–189. doi:10.1038/cr.2014.3

PubMed Abstract | CrossRef Full Text | Google Scholar

Ranke M. B. (2015). Insulin-like growth factor binding-protein-3 (IGFBP-3). Best. Pract. Res. Clin. Endocrinol. Metab. 29 (5), 701–711. doi:10.1016/j.beem.2015.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren W., Yin J., Duan J., Liu G., Tan B., Yang G., et al. (2016). mTORC1 signaling and IL-17 expression: Defining pathways and possible therapeutic targets. Eur. J. Immunol. 46 (2), 291–299. doi:10.1002/eji.201545886

PubMed Abstract | CrossRef Full Text | Google Scholar

Rigatti S. J. (2017). Random forest. J. Insur. Med. 47 (1), 31–39. doi:10.17849/insm-47-01-31-39.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruedel A., Stark K., Kaufmann S., Bauer R., Reinders J., Rovensky J., et al. (2013). N‐cadherin promoter polymorphisms and risk of osteoarthritis. FASEB J. 28 (2), 683–691. doi:10.1096/fj.13-238295

PubMed Abstract | CrossRef Full Text | Google Scholar

Schlichting N., Dehne T., Mans K., Endres M., Stuhlmuller B., Sittinger M., et al. (2014). Suitability of porcine chondrocyte micromass culture to model osteoarthritis in vitro. Mol. Pharm. 11 (7), 2092–2105. doi:10.1021/mp5000554

PubMed Abstract | CrossRef Full Text | Google Scholar

Song H., Song J., Cheng M., Zheng M., Wang T., Tian S., et al. (2021). METTL3-mediated m(6)A RNA methylation promotes the anti-tumour immunity of natural killer cells. Nat. Commun. 12 (1), 5522. doi:10.1038/s41467-021-25803-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun T., Wu R., Ming L. (2019). The role of m6A RNA methylation in cancer. Biomed. Pharmacother. 112, 108613. doi:10.1016/j.biopha.2019.108613

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilkerson M. D., Hayes D. N. (2010). ConsensusClusterPlus: A class discovery tool with confidence assessments and item tracking. Bioinformatics 26 (12), 1572–1573. doi:10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

Woodell May J. E., Sommerfeld S. D. (2019). Role of inflammation and the immune system in the progression of osteoarthritis. J. Orthop. Res. 38 (2), 253–257. doi:10.1002/jor.24457

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu Y., Xie L., Wang M., Xiong Q., Guo Y., Liang Y., et al. (2018). Mettl3-mediated m(6)A RNA methylation regulates the fate of bone marrow mesenchymal stem cells and osteoporosis. Nat. Commun. 9 (1), 4772. doi:10.1038/s41467-018-06898-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao W., Adhikari S., Dahal U., Chen Y. S., Hao Y. J., Sun B. F., et al. (2016). Nuclear m(6)A reader YTHDC1 regulates mRNA splicing. Mol. Cell. 61 (4), 507–519. doi:10.1016/j.molcel.2016.01.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Yeh L. A., Augustine A. J., Lee P., Riviere L. R., Sheldon A. (1995). Interleukin-4, an inhibitor of cartilage breakdown in bovine articular cartilage explants. J. Rheumatol. 22 (9), 1740–1746.

PubMed Abstract | Google Scholar

Yucesoy B., Charles L. E., Baker B., Burchfiel C. M. (2015). Occupational and genetic risk factors for osteoarthritis: A review. Work 50 (2), 261–273. doi:10.3233/WOR-131739

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang B., Wu Q., Li B., Wang D., Wang L., Zhou Y. L. (2020). m(6 A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol. Cancer 19 (1), 53. doi:10.1186/s12943-020-01170-0

PubMed Abstract | CrossRef Full Text | Google Scholar