Comprehensive analysis of m6A RNA methylation modification patterns and the immune microenvironment in osteoarthritis

Background: Osteoarthritis (OA) is the most common joint degenerative disease, and so far, there is no effective therapy to prevent or delay its development. Considerable attention is now being given to the impact of m6A RNA methylation modification on the disease immune regulation. However, much remains unknown about the function of m6A modification in OA.

Methods: A total of 63 OA and 59 healthy samples were applied to comprehensively examine the m6A regulators mediated RNA methylation modification pattern in OA, and evaluate the impacts of distinct patterns on the characteristics of OA immune microenvironment, including immune infiltration cells, immune responses and human leukocyte antigen (HLAs) genes expression. In addition, we screened out the m6A phenotype-related genes and further explored their potential biological functions. At last, we verified the expression of key m6A regulators and their associations with immune cells, in vitro.

Results: Most of m6A regulators was differentially expressed in OA samples compared to the normal tissues. Based on six hub-m6A regulators identified as abnormally expressed in OA samples, we developed a classifier to distinguish OA patients from healthy individuals. We noted that immune characteristics of OA were correlated with m6A regulators. For instance, YTHDF2 had a strongest significantly positive correlation with regulatory T cells (Tregs) and IGFBP2 was strongest negatively associated with dendritic cells (DCs), which were confirmed by the immunohistochemistry (IHC) staining. Two distinct m6A modification patterns were determined: pattern B had higher infiltrating immunocytes and more active immune responses than pattern A, and two patterns differed in the expression of HLA genes. We also identified 1,592 m6A phenotype-related genes that could mediate the OA synovitis and cartilage degradation by the PI3K-Akt signaling pathway. Quantitative real-time polymerase chain reaction (qRT-PCR) results indicated that IGFBP2 was significantly overexpressed, while YTHDF2 mRNA expression was decreased in OA samples, which was consistent with our findings.

Conclusion: Our research proves the essential impact of m6A RNA methylation modification on the OA immune microenvironment, and helps to explain the regulatory mechanism behind it, which may open up a new direction for more precise immunotherapy of osteoarthritis.

1 Introduction

Osteoarthritis (OA) is the most prevalent joint degenerative disease, accompanied by pain, joint deformity and disability (1). It is the main cause of disability and reduction in life expectancy among the elderly (2, 3). Although joint replacement surgery is an option, there is as yet no effective therapy to prevent the OA progression (4). There are plenty of risk factors for its development, including genetic predisposition, aging, sex, obesity, trauma, anatomical abnormality, low-grade systemic inflammation and the failure of the oxidation-antioxidant balance (5–7). Among them, the low-grade chronic inflammation plays a central role. It has been reported that immune systems are involved in the low-grade inflammation related to OA, and the inflammation of OA tissues is the results of interaction between immune system, local tissue damage, metabolic dysfunction and other factors (8). Recent researches have demonstrated that the pathogenesis of OA has a strongly close correlation with the immune microenvironment, and the immune pathways as well as the infiltrating immunocytes in OA synovial tissue are linked with the pain and progression of OA (9). Woodell-May JE and Sommerfeld SD have revealed the huge potentials of T cells and macrophages in the synovial membrane to regulate the OA inflammatory response (10), and the relevant studies on the chondrocyte apoptosis and cartilage matrix proteolysis mediated by the immune responses and immune cells have also been reported (11). Further elucidating the immune regulation mechanism of OA may explain its potential pathogenesis and generate novel therapeutic methods that may modify its course.

RNA methylation is thought to be a key regulator of disease progression, and it has important impacts on many aspects of RNA metabolism, such as RNA splicing and subcellular localization, RNA degradation, mRNA stability and translation (12–14). N6-methyladenosine (m6A) modification is the most abundant internal modifications found in eukaryotes, which recognized by m6A binding proteins (readers), installed by methyltransferases (writers), and removed by demethylases (erasers) (15, 16). m6A RNA modification actively participates in several physiological processes, including DNA damage response, tumorigenesis, cell cycle progression, differentiation and circadian rhythm maintenance (17, 18). Recent studies have revealed that m6A RNA methylation shows a strong association with inflammation responses. Dubey et al. have demonstrated that the knockdown of m6A demethylase FTO contributes to hypermethylation of pro-inflammatory cytokines, and regulation of m6A−RNA methylation could effectively suppress myocardial inflammation and dysfunction during endotoxemia in mice (19). METTL3 is highly expressed in rheumatoid arthritis and leads to fibroblast-like synoviocytes (FLSs) activation and inflammatory response via the NF-κB signaling pathway (20). Inhibition of METTL3 obviously alleviates renal injury and inflammation through IGF2BP2-dependent mechanisms, which suggests that targeting METTL3 may be a promising therapeutic approach for acute kidney injury (AKI) (21). In hematopoietic stem cell (HSC), m6A reader YTHDF2 acts as a key repressor of proinflammatory pathways, and the crucial significance of m6A mRNA modification in mediating the inflammatory pathways of HSC has also been proved (22). Several m6A regulators such as ALKBH5, METTL3, and FTO have been reported to be involved in the progression of osteoarthritis and attenuate the inflammatory response, which provides a potential strategy for OA treatment (23–25). Evidences have shown that aberrant m6A RNA methylation modification is connected with the immune regulation, and m6A modification is the main regulator of post-transcriptional immune response in cells (26). m6A RNA methylation could mediate the tumor immunity by regulating the expression of antioncogenes and oncogenes, which may provide a highly targeted mechanism for tumor immunotherapy. Served as promising predictors, m6A RNA methylation regulators are able to identify pancreatic cancer (PAAD) patients who may really benefit from treatment with immunotherapy and enhance the efficacy of immunotherapy for PAAD, which would allow to refine the therapeutic options and to better adjust the treatment strategies (27). In addition, m6A RNA methylation modification pattern take part in the regulation of immune microenvironment of rheumatoid arthritis, suggesting it may have a potential role in the immunomodulation and immunotherapy for bone and joint diseases (28). However, there is still a gap in the study on the function of m6A in osteoarthritis, and the correlation of m6A modification with the OA immune microenvironment requires further elucidation.

Nowadays, the study on the epigenetic regulation mechanism of biological processes has contributed to the understanding of the pathogenesis of various human diseases caused by abnormal epi-transcriptomic and epigenomic modification. Specifically, the increasing findings have determined the connection between the pathogenesis of diseases and aberrant modifications. In our research, we comprehensively explored the m6A regulators mediated RNA methylation modification patterns in osteoarthritis, and determined the potential regulation mechanism of m6A RNA methylation modification on the OA immune microenvironment. Our findings have showed that m6A regulators could effectively differentiate OA from healthy samples and were closely related to infiltrating immunocytes and immune pathways in OA, indicating that m6A regulators could mediate the regulation of immune microenvironment in OA. Based on 18 m6A regulators, we clustered OA samples and screened out two m6A RNA methylation modification subtypes, with distinct biological functions and immune characteristics. Additionally, 1,592 m6A phenotype-related genes and their underlying biological functions as well as molecular mechanism were observed. In vitro, we verified the expression of key m6A regulators and their associations with immune cells. Our study suggests that m6A RNA methylation modification pattern is of vital importance in the immune microenvironment of OA, and this classification of immune subtypes based on m6A regulators may provide a theoretical basis for precise OA immunotherapy.

2 Materials and methods

2.1 Clinical samples

We collected 10 OA tissues from patients who undergoing total knee arthroplasty, and 10 normal samples were gained from patients with ligament injury who underwent ligament reconstruction surgery at Qilu Hospital of Shandong University, China. Patients were randomly selected and all signed informed consent. Our research was approved by the Ethics Committee of Qilu Hospital of Shandong University.

2.2 Data pre-processing

The datasets of gene expression profile were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/): GSE82107, GSE29746, GSE55235, GSE55457, GSE117999 and GSE98918. A total of 63 OA samples and 59 healthy individuals were gained. The sva R package was applied to normalize the six gene sets (Figure S1). Detailed clinical features and platforms of patients could be found in Supplementary Files 1. 18 m6A regulators were annotated in the metadata sets: IGFBP2, YTHDF2, HNRNPA2B1, METTL3, YTHDC2, FTO, RBM15, FMR1, LRPPRC, IGFBP3, YTHDF1, HNRNPC, YTHDC1, ZC3H13, RBM15B, YTHDF3, IGFBP1 and WTAP.

2.3 Difference analysis of m6A RNA methylation regulators between OA and healthy patients

The protein-protein interaction (PPI) network among 18 m6A RNA methylation regulators was constructed using the STRING database. The correlations among m6A regulators expression in all samples or OA samples were determined by Spearman’s correlation. The limma R package was employed to screen out the differentially expressed genes (DEGs) between OA and normal tissues, and the cut-off criterion was P value < 0.05. Wilcox test was utilized to compare the differences in m6A regulators expression between OA and healthy donors. With a cut-off of P < 0.05, the univariate logistic regression was applied to find the m6A regulators linked with OA occurrence. Feature selection and dimension reduction were performed by the least absolute shrinkage and selection operator (LASSO) Cox regression (10-fold). Subsequently, the results of the univariate logistic and LASSO regression overlapped, and the corresponding coefficients of m6A RNA methylation regulators were computed using multivariate logistic regression. The performance of risk score in distinguishing OA from healthy patients was evaluated by ROC curve analysis via the pROC package in R. Values of the area under the curve (AUC) > 0.7 and P < 0.05 suggested that it had high predictive ability.

2.4 Correlation analysis of m6A RNA methylation regulators with immune characteristics in OA

The reactivity and number of immune infiltration cells as well as the activity of the immune responses were detected using the single-sample gene-set enrichment analysis (GSEA), via the gsva package in R. Wilcox tests was applied to compare the infiltrating immunocytes abundance and the enrichment fraction of immune response pathways between OA and normal patients. Using the corrplot R package, connections between hub m6A regulators, infiltrating immunocytes and immune response pathways were examined by Spearman correlation test.

2.5 Identification of m6A RNA methylation modification patterns in OA

The unsupervised consensus cluster analysis was implemented to identify the m6A modification pattern in OA samples, based on the 18 m6A RNA methylation regulators. A consensus clustering algorithm was employed to assess the number and robustness of cluster, and the consul cluster plus package in R for 1,000 iterations to ensure classification robustness. The expression pattern of 18 m6A RNA methylation regulators in different m6A modification pattern was confirmed via the principal component analysis (PCA). Next, 18 m6A regulators expression, infiltrating immunocytes abundance, immune responses and HLA genes expression of the two different m6A subtypes were compared by the Kruskal test.

2.6 Functional enrichment analysis of distinct m6A RNA methylation modification patterns

The biological function of distinct m6A RNA methylation modification pattern was analyzed via the clusterProfiler R package, including GO, KEGG and Hallmark enrichment analysis. The c5.go.v7.4.symbols, c2.cp.kegg.v7.4.symbols and h.all.v7.4.symbols gene sets were employed to reflect alterations of the biological signaling pathways. The GSVA algorithm was applied to convert the expression matrix into a scoring matrix, and the limma package in R was utlized to compare the scores of biological signaling pathways among these groups. P < 0.05 was selected as the threshold of difference analysis.

2.7 Identification of genes mediated by m6A RNA methylation regulators in OA

All samples in different m6A RNA methylation modification pattern was determined using the empirical Bayesian method to screen out DEGs between the two m6A subtypes of OA. The adjusted P < 0.0001 was applied to find the significant DEGs. m6A modification pattern-related genes and gene modules were constructed by weighted gene co-expression network analysis (WGCNA). The topology was calculated with a soft threshold of 1 to 20, and 10 was set as the optimal soft threshold. According to the topological overlap matrix (TOM), the related modules were classified. Each module contains at least 20 genes, and then similar modules were combined. The correlation between the combined modules and distinct m6A modification patterns was computed using the Pearson approach.

2.8 Immunohistochemistry staining

The correlations between the expression of YTHDF2 and Treg marker FOXP3, IGFBP2 and DCs marker CD1a were verified by IHC staining. The sections of OA and normal tissues were added into 10 mM sodium citrate solution (pH 6.0) for deparaffinization and antigen retrieval, and blocked with PBS containing 5% goat serum for 30 min RT. The primary antibody diluted in 1:200 was added into the slides and incubated overnight at 4°C, followed by incubation with the HRP-conjugated secondary antibody and stained with the 3.30-diaminobenzidine (DAB) substrate. The tissue score (IHC score) was calculated based on the percentage of stained cells and tissue staining intensity, and the intensity scoring standards were as follows: strong (+++), moderate (++), weak (+), and negative (–). YTHDF2 polyclonal antibody (24744-1-AP, Proteintech Group, Inc.) was used to detect YTHDF2. IGFBP2 polyclonal antibody (11065-3-AP, Proteintech Group, Inc.) was utilized to detect IGFBP2. FOXP3 monoclonal antibody (YM0286, ImmunoWay Biotechnology, Inc.) was applied to examine Treg marker FOXP3. CD1 monoclonal antibody (YM3068, ImmunoWay Biotechnology, Inc.) was employed to examine DCs marker CD1a. Mouse IgG (Santa Cruz, sc-2025) and rabbit IgG (Proteintech, B900610) were used as negative controls. CoraLite488 conjugated Affinipure goat anti-mouse IgG (H+L) (SA00013-1, Proteintech Group, Inc.) and CoraLite594 conjugated goat anti-rabbit IgG (H+L) (SA00013-4, Proteintech Group, Inc.) were used as secondary antibodies.

2.9 RNA extraction and qRT-PCR

Total RNAs was gained by applying 1 ml TRIzol reagent (Sigma-Aldrich) from OA and normal tissues, then reverse transcribed into cDNA using a reverse transcription kit (TaKaRa). IGFBP2 and YTHDF2 mRNA expression were detected by qRT-PCR with SYBR Green PCR Master Mixes (Thermo Fisher). The 2−ΔΔCt method was applied to calculate the relative quantification and select β-Actin as the internal control. The sequences of YTHDF2 were as follows: AGCCTCTTGGAGCAGTACAAA (sense) and GCATTATTGGGCCTTGCCTG (antisense). The sequences of IGFBP2 were as follows: GAGCAGGTTGCAGACAATGGC (sense) and TACCCGACTTGAGGGGCTT (antisense). The sequences of ß-actin were as follows: CACCATTGGCAATGAGCGGTTC (sense) and AGGTCTTTGCGGATGTCCACGT (antisense).

3 Results

3.1 Expression landscape of m6A RNA methylation regulators in OA and healthy samples

A total of 18 m6A regulators were screened out, including five writers (WTAP, ZC3H13, RBM15B, RBM15 and METTL3), 12 readers (FMR1, LRPPRC, YTHDC1, YTHDC2, HNRNPA2B1, HNRNPC, YTHDF1, YTHDF2, YTHDF3, IGFBP1, IGFBP2 and IGFBP3) and one eraser (FTO). The protein-protein interaction (PPI) network described the close connections among 18 m6A regulators, which suggested that the cross-talk among them may affect the formation of m6A modification pattern and participate in the occurrence and development of OA (Figure 1A). Next, we studied the relationship between transcriptomes by computing the pairwise relation between 18 m6A regulators expression, and noticed that positive connections were more frequent than negative ones, and the strongest positive correlation in all samples was between YTHDF3 and LRPPRC (r = 0.54, Figure 1B). In the OA samples, YTHDF3 and LRPPRC also showed the strongest positive correlation (r = 0.59, Figure 1C), indicating their important role in the m6A RNA methylation. The differential expression analysis revealed five m6A regulators (YTHDF2, HNRNFA2B1, METTL3, YTHDC2 and IGFBP2) with altered expressions in OA samples, compared with the healthy control (Figure 1D). Amongst these, the multiple change in IGFBP2 was the largest (Figure 1E). In addition, Wilcox test showed the significant difference in m6A regulators expression between OA and normal tissues, including YTHDF2, HNRNFA2B1, METTL3, YTHDC2 and IGFBP2 (Figure 1F). The above findings demonstrated that m6A RNA methylation regulators expression was highly heterogeneous in OA patients and healthy individuals, suggesting the potential functions of the abnormal m6A regulators expression patterns in the progression of OA.

FIGURE 1

Figure 1 Expression landscape of m6A regulators in OA and healthy samples. (A) The protein-protein interaction (PPI) network of 18 m6A regulators. (B) Connection analysis between 18 m6A regulators expression in all samples. (C) Connection analysis between 18 m6A regulators expression in OA samples. (D, E) Heat map and volcano map of differences in hub m6A regulators expression between health and OA patients. (F) Box plot showing 18 m6A regulators expression in OA samples compared with the control. **p < 0.01, and ***p < 0.001.

3.2 m6A RNA methylation regulators were involved in the OA progression

To determine the function of m6A RNA methylation regulators in the pathogenesis of OA, the univariate logistic regression analysis of 18 m6A regulators was carried out. This suggested that the above five expressed differentially m6A regulators were related most closely to OA occurrence (P < 0.05, Figure 2A). LASSO regression was utilized for feature selection and dimension reduction to eliminate redundant genes (Figures 2B, C), and the results indicated these five regulators were all important for OA. We then developed a classifier composed of five m6A modulators that are able to distinguish OA samples from healthy samples using the multivariate logistic regression based on the risk scores, where the risk scores of OA samples are much higher than those of healthy group (Figures 2D, E). Risk scores were calculated as expression of HNRNPA2B1*-0.5548 + expression of IGFBP2*0.8090 + expression of METTL3*-0.5232 + expression of YTHDC2*-0.8723 + expression of YTHDF2*-1.7990. Analysis of ROC curve suggested that this classifier showed a good predictive performance on the diagnosis of osteoarthritis (AUC = 0.785, Figure 2F). We concluded that m6A RNA methylation regulators could be identified as potential indicators and participated in the process of OA.

FIGURE 2

Figure 2 m6A regulators were linked with the OA process. (A) Five m6A regulators showed close correlations with the occurrence of OA through univariate logistic regression (P < 0.05). (B) LASSO coefficient profiles of the five OA-related m6A regulators. (C) In LASSO regression 10-fold cross-validation applied to turn parameter selection. (D) Multivariate logistic regression to establish the distinguishing signatures of five m6A regulators and calculate risk scores for OA. (E) Risk profiles of OA and healthy samples. (F) The ROC curve analysis of hub m6A regulators for OA diagnosis. ***p < 0.001.

3.3 Correlation analysis between m6A regulators and immune characteristics in OA

The connections between dysregulated m6A RNA methylation regulators and immune infiltration cells as well as immune response pathways were analyzed, to study its contribution to the OA immune microenvironment. Difference analysis showed that the abundance of most infiltrating immunocytes in OA samples was significantly higher than that in health, such as dendritic cells (DCs), macrophages, mast cells, TIL and Th2 cells (Figure 3A). And we found the five hub m6A regulators had close connections to several immune cells in OA samples using the correlation analysis (Figure 4A). For example, YTHDF2 showed the strongest positive relation to regulatory T cells (Tregs, r = 0.429), and the strongest negative correlation was between IGFBP2 and DCs (r = -0.325). Subsequently, we collected 10 OA samples and 10 normal tissues, and the results of IHC staining indicated that the expression of IGFBP2 in OA tissues was higher than that of the normal, while YTHDF2 expression was decreased in OA group. In OA tissues, correlation analysis suggested that IGFBP2 was negatively linked to DCs marker CD1a (Figures 5A, C), and the positive correlation between YTHDF2 and Treg marker FOXP3 was also confirmed (Figures 5B, D). The immune response of OA tissues was more active than that of healthy samples, including HLA, Cytolytic activity, Type I and Type Π IFN inflammatory responses, which suggest that these pathways might participate in the inflammatory process of OA (Figure 3B). We further analyzed the differences in other inflammatory response pathways between OA and normal samples (Figure S2). In addition, YTHDF2 had the strongest positive relation to the HLA pathway (r = 0.369), HNRNPA2B1 and Type I IFN response showed the strongest negative correlation (r = -0.355, Figure 4B). These results indicated that m6A RNA methylation regulators were of essential importance in the immune microenvironment of OA and could mediate the infiltration of immune cells and the responses to immunity in OA.

FIGURE 3

Figure 3 Differential expressions of immunocytes and immune-related response pathways in OA. (A) Differential expressions of immunocytes in OA compared to the healthy samples. (B) Differential expressions of immune response pathways in OA compared to the healthy samples. *p < 0.05, **p < 0.01, and ***p < 0.001.

FIGURE 4

Figure 4 Relationship analysis of m6A regulators with immune characteristics in OA. (A) Dot plot showing the relationship of infiltrating immune cells with hub m6A regulators in the OA immune microenvironment. The most positively related immunocyte-m6A regulator pair is YTHDF2-Treg cell (r = 0.429), and the most negative is IGFBP2-DC cell (r = -0.325), with expression and fraction status presented in the upper panel. (B) Dot plot showing the relationship of immune response pathways with hub m6A regulators in the OA immune microenvironment. The most positively associated pair is the YTHDF2-HLA pathway (r = 0.369), the most negatively correlated is HNRNPA2B1-Type I IFN response (r = -0.355), with the violin plot in the bottom panel showing expression or activity.

FIGURE 5

Figure 5 Correlations between the expression of IGFBP2 and DCs marker CD1a, YTHDF2 and Treg marker FOXP3 in OA and normal tissues. (A) Representative images for the IHC staining of IGFBP2 and CD1a in OA and normal tissues (n = 10). (B) Representative images for the IHC staining of YTHDF2 and FOXP3 in OA and normal tissues (n = 10). (C) Relation between protein levels of IGFBP2 and CD1a in OA and normal tissues via Spearman correlation. (D) Relation between protein levels of YTHDF2 and FOXP3 in OA and normal tissues via Spearman correlation.

3.4 Identification of m6A RNA methylation modification patterns

On the basis of 18 m6A regulators expression, the unsupervised consensus cluster analysis was carried out to detect m6A RNA methylation modification patterns in OA patients, and two different subtypes were found in our study (Figure 6A). PCA analysis showed that m6A regulators expression was qualitatively different, and m6A regulators effectively enabled classification of OA patients into two groups, with 35 samples in subtype A and 28 samples in subtype B (Figure 6B). Most of m6A regulators expression was significantly different between two m6A subtypes of OA, with FTO, IGFBP2, IGFBP1, IGFBP3 and YTHDF1 being highly expressed in subtype B (Figures 6C, D). We further investigated the underlying biological processes and molecular mechanisms in the different m6A RNA methylation modification pattern of OA, and KEGG pathway analysis suggested that pattern B showed more active pathways than pattern A, such as FOCAL ADHESION and ECM RECEPTOR INTERACTION (Figure 7A). And GO analysis indicated that pattern B also had more abundant biological processes, compared to the pattern A, including PURINE NUCLEOSIDE BIOSYNTHETIC PROCESS and REGULATION OF INTERGRIN MEDIATED SIGNALING PATHWAY (Figure 7B). In the Hallmark analysis, PROTEIN SECRETION and EPITHELIAL MESENCHYMAL TRANSITION were enriched in pattern B (Figure 7C). These findings have revealed the presence of diverse m6A modification patterns in OA, which performed different biological functions in the progression of OA.

FIGURE 6

Figure 6 Distinct m6A RNA methylation modification pattern of OA. (A) Two different m6A subtypes of OA by the unsupervised consensus cluster analysis. (B) Principal component analysis of the transcriptome profiles indicated the two subtypes were significantly different. (C, D) Heatmap and box plot of 18 m6A regulators expression between the two m6A modification subtypes of OA. *p < 0.05, **p < 0.01, and ***p < 0.001.

FIGURE 7

Figure 7 Functional analysis of different m6A RNA methylation modification patterns in OA. (A) Differences in KEGG pathway enrichment scores. (B) Differences in GO pathway enrichment scores. (C) Differences in Hallmark analysis scores.

3.5 Characteristic of the immune microenvironment in distinct m6A modification pattern

Next, we explored the differences in the characteristics of the immune microenvironment between two m6A modification patterns in OA, and examined the differences in the infiltrating immunocytes abundance immune response pathways and HLA genes expression. Our results demonstrated that pattern B had a relatively high infiltrating immunocytes compared with pattern A, and higher infiltration levels of macrophages, NK cells, mast cells and Tfh, whereas DCs were mainly involved in the pattern A (Figure 8A). The immune response of pattern B was also more active than that of pattern A, with regard to Cytolytic activity, HLA and Type I IFN response (Figure 8B). HLA genes expression was significantly different in two different m6A modification subtypes of OA (Figure 8C). We concluded that different m6A modification pattern mediated immune response in OA was distinguishing, pattern B mediated an active immune response while pattern A led to a mild in the OA immune microenvironment. It was verified that m6A RNA methylation modification played an essential regulatory role in the formation of different immune microenvironments of OA.

FIGURE 8

Figure 8 Characteristics of the immune microenvironment in two different m6A subtypes of OA. (A) Differences in the infiltrating immunocytes in the immune microenvironment of distinct m6A subtypes in OA. (B) Differences in the immune pathways in the immune microenvironment of distinct m6A subtypes in OA. (C) Differences in the HLA genes expression in the immune microenvironment of distinct m6A subtypes in OA. *p < 0.05, and **p < 0.01.

3.6 Identification of m6A phenotype-related genes in OA

To determine the potential mechanism of genes involvement in the regulation mediated by m6A regulators, we screened out m6A phenotype-related genes, amounting to a total of 1,592 genes. On the basis of the expressions of these genes, we performed the unsupervised consensus cluster analysis in OA samples, and three gene subtypes were identified (Figure 9A). PCA analysis showed three groups of OA patients, with 13 samples in subtype A, 16 in subtype B and 34 in subtype C (Figure 9B). We noticed the significant differences between three subtypes in 18 m6A regulators expression, immune infiltrating cells, immune pathways and HLA genes expression. Most of m6A regulators was expressed differentially among the three subtypes (Figure 9C). Compared to the patterns B and C, pattern A showed the fewest abundance of infiltrated immune cells. Pattern B had higher infiltrated levels of DCs, iDCs, pDCs, TIL and Th1 cells, while macrophages and mast cells were enriched in pattern C (Figure 9D). Patterns B and C had more active immune reactions than pattern A, especially pattern B. CCR, Check-point, Inflammation-promoting, Para-inflammation and Type I IFN response were very active in pattern B, as was HLA in pattern C (Figure 9E). Moreover, except for HLA-DOB, the expression of HLA genes was higher in pattern C than in pattern A and pattern B (Figure 9F).

FIGURE 9

Figure 9 Characteristics of the immune microenvironment in m6A phenotype-related genes modification patterns of OA. (A) Three different m6A phenotype-related genes modification subtypes of OA by the unsupervised consensus cluster analysis. (B) Principal component analysis of the transcriptome profiles indicated the three subtypes were significantly different. (C-F) Differential expressions of 18 m6A regulators, immunocytes, immune response pathways and HLA genes in the immune microenvironment of m6A phenotype-related genes modification patterns in OA. *p < 0.05, **p < 0.01, and ***p < 0.001.

We then performed functional enrichment analysis on the m6A phenotype-related genes, and the results of GO analysis showed them to be mainly involved in the extracellular matrix organization (Figure 10A). In the KEGG analysis, these genes were mostly enriched in the PI3K-Akt signaling pathway (Figure 10B). Using WGCNA to screen out gene-gene modules related to the distinct m6A modifications in OA (Figures 10C, D), and we identified two gene modules in which different m6A modification patterns were matched with their related genes (Figure 10E). These findings may shed light on the regulatory network of gene expression mediated by m6A modification patterns.

FIGURE 10

Figure 10 Functional analysis of m6A phenotype-related genes in OA. (A) GO enrichment analysis. (B) KEGG enrichment analysis. (C) Sample dendrogram and trait heatmap using WGCNA. (D) Hierarchical clustering tree shows each module. (E) Heatmap shows the relationship of m6A RNA methylation modification pattern with module eigengenes.

3.7 qRT-PCR validation of key m6A regulators

Radiological evaluations of knee joints were performed in distinct groups of patients to validate the typical imaging manifestations of OA. Our study noticed that preoperative knee X-ray and intraoperative macroscopic images of OA patients showed joint space stenosis, cartilage degeneration, subcartilage bone hyperplasia, and osteophyte formation in comparison to those of the healthy donors. (Figure 11A). The results of qRT-PCR showed that compared the normal group, the expression of IGFBP2 mRNA was significantly up-regulated, while YTHDF2 mRNA expression was decreased in OA samples, which was consistent with our findings (P < 0.01, Figures 11B, C).

FIGURE 11

Figure 11 qRT-PCR validation of key m6A regulators. (A) X ray images and intraoperative macroscopic views (arthroplasty and arthroscopy image) of knee joint from OA patients and healthy donors. (B, C) Validation of IGFBP2 and YTHDF2 mRNA expression by qRT-PCR between OA (n = 10) and normal group (n = 10). **p < 0.01, and ***p < 0.001.

4 Discussion

Osteoarthritis (OA) is considered as a degenerative inflammatory disease of the joint cartilage, with inflammatory mediators released by synovium, bone and cartilage (29). The inflammatory response of OA has been implicated in its pathogenesis, and elaborating on the inflammatory pathophysiology of osteoarthritis may represent new avenues for OA therapy (30). Extensive studies indicated that immunocytes and immune pathways effectively modulated the inflammatory response of OA and correlated with OA disease pain and progression (9). One of the m6A regulators, demethylase METTL3 could contribute to the development of OA by mediating the inflammatory response, ECM synthesis and chondrocyte apoptosis, which suggested that targeting of METTL3 or m6A RNA methylation may be a promising preventive and therapeutic approach for OA (31). In addition, m6A RNA methylation modification reportedly was involved in the immune responses of numerous diseases (32), and its functions in the immune microenvironment of various tumors had also been proved (33–35). Thus, we had reason to believe that m6A modification was connected with the regulation of OA immune microenvironment. We combined a variety of bioinformatics methods and experiments to study the m6A modification pattern in the OA immune microenvironment, and determined how m6A regulators mediated the abundance of immunocytes and immune responses in OA.

PPI networks showed that 18 m6A regulators significantly interacted with each other, suggesting that readers, erasers and writers do function in collaboration, in osteoarthritis, not in isolation. Correlation analysis has revealed that the strongest positive correlation among 18 m6A regulators was YTHDF3 and LRPPRC in both OA and all samples, which reminded us their important role in the m6A RNA methylation. Compared to the healthy samples, we identified five m6A regulators with altered expressions in OA patients, suggesting that these regulators may be related to the occurrence and development of OA. Based on the hub OA-related m6A regulators, we have created a novel classifier that could effectively distinguish OA from normal samples. These results further verified the crucial function of m6A regulators in osteoarthritis, and it could be identified as an independent diagnostic indicator of OA. We conducted qRT-PCR to validate the expression of key m6A regulators and the results proved that compared with the normal tissues, IGFBP2 mRNA was highly expressed while the mRNA expression of YTHDF2 was significantly down-regulated in OA patients, which was consistent with our findings. Subsequently, our study examined the differences in the characteristics of immune microenvironment between OA and healthy patients, including infiltrating immunocytes and immune response, and explored the potential relation of hub m6A regulators with immune characteristics in OA. We noted that the abundance of immune infiltration cells in OA tissues was higher than that in healthy individuals, and the immune response of OA group was also more active. In addition, we found YTHDF2 had the strongest positive connection with Treg cells, and IGFBP2 was strongest negatively associated with DCs. We collected samples from OA patients, and performed the IHC staining to confirm these correlations. The results showed that in OA tissues, YTHDF2 was positively linked with Treg marker FOXP3, and the negative correlation between IGFBP2 and DCs marker CD1a was also proved. Shao et al. have shown that YTHDF2 expression was positively associated with Treg markers FOXP3, CCR8, TGFB1 and STAT5B, and the activation of Treg cells by YTHDF2 led to a poor prognosis of liver hepatocellular carcinoma (36). IGFBP2 significantly facilitated the invasion and metastasis of pancreatic ductal adenocarcinoma via the NF-κB signaling pathway (37). IGFBP2 was correlated with M2 Macrophages, CD8+ T cells, which enhanced the immunosuppression of glioblastoma and it might be an effective immunotherapeutic strategy for glioblastoma therapy (38); However, there are no similar reports relating to IGFBP2 and DCs. Further studies were required to examine the potential mechanism of m6A regulators and Tregs as well as DCs in the immune regulation of OA. In terms of immune response, YTHDF2 was strongest positively connected with the HLA pathway, and HNRNPA2B1 was strongest negatively linked to the Type I IFN response, indicating the involvement of YTHDF2 and HNRNPA2B1 in the OA progression. The above findings revealed that the cross-talk among m6A regulators had the crucial function in the occurrence of OA, and might influence on regulating the immune microenvironment of OA.

On the base of 18 m6A regulators expression, the unsupervised consensus cluster analysis was carried out and two m6A modification patterns with different immune phenotypes in OA were identified. In the pattern B, the infiltration abundance of DCs, Tfh, NK cells, macrophages and mast cells were higher than that of pattern A. The immune response of pattern B was also more active, compared with the pattern A, which activated in the HLA and Type I IFN signaling pathways, reported in many chronic inflammatory diseases. And HLA genes expression significantly differed between the two m6A subtypes of OA. This classification of immune subtypes has provided insights into the regulation mechanism of the immune microenvironment in OA, which may be conducive to the formulation of more appropriate immunotherapy strategies for OA. This was similar to the research classifying gastric cancer into three m6A RNA methylation modification subtypes, which contributed to understand the regulation mechanism of tumor microenvironment and made more effective therapeutic methods for gastric cancer (39). In the periodontitis, different m6A RNA methylation modification pattern was linked with the phenotypic features and had a profound influence on the complexity and heterogeneity of the periodontitis immune microenvironment (40). We compared the biological functions of the two different m6A subtypes of OA, and enrichment analysis showed that pattern B had more biological processes and abundant pathways than pattern A, including ECM receptor interaction, purine nucleoside biosynthetic process and protein secretion. The above findings could clarify how m6A modification patterns mediated the regulation network of gene expression, which may elucidate the OA pathogenesis from m6A RNA methylation modification point of view. Subsequently, we identified 1,592 m6A phenotype-related genes that were mostly enriched in the PI3K-Akt signaling pathway. Recently studies have shown that the PI3K/Akt signaling pathway had a vital function in the subchondral bone dysfunction, cartilage degradation and synovial inflammation of OA, and affected the development of OA (41, 42). Xue et al. suggested that blocking the PI3K-AKT-mTOR pathway could effectively alleviate inflammatory reaction and induce chondrocytes to autophagy in OA (43). Having divided these genes into three m6A phenotype-related genes modification patterns, we realized the significant differences in m6A regulators expression, immune infiltrating cells, immune responses and HLA genes expression among these subtypes. By WGCNA analysis, we screened out the gene-gene modules linked with the distinct m6A RNA methylation modifications, and noticed that pattern A was closely associated with the grey module, indicating that the gene expression in this module was greatly influenced by m6A RNA methylation regulators. These results unraveled the regulatory network of gene expression mediated by m6A modification patterns, and contributed to understanding the pathologic mechanisms of OA.

We comprehensively examined the potential relation of m6A modification pattern with the OA immune microenvironment, confirming that m6A RNA methylation modification was closely connected with the OA process and had a regulatory impact on the immune microenvironment of OA. Our study filled the gap in the field of OA epigenetics, and raised a novel direction for the study about m6A modification of the immune-related pathogenesis in OA. Our studies of the genes associated with m6A modification patterns deepened our comprehension of the functional mechanism of m6A RNA methylation modification in OA, which may offer clues to appropriate therapeutic interventions and play critical functions in alleviating or treating OA. Nevertheless, our research had some limitations. The size of samples was limited, so our research results required verification through further experiments. All the data was from the mRNA expression values of genes and not directly reflected their protein expression levels, which might lead to ineffective evaluation of the immune response pathways based on protein expression. Although such limitations should not be ignored, our study did confirm the significant influence of m6A RNA methylation modification on the OA immune microenvironment and helped explain the underlying pathological mechanism.

5 Conclusions

In conclusion, our study found that m6A RNA methylation regulators showed a strong influence on the immune characteristics of OA, and uncovered the prospective regulatory mechanism of m6A modification in the OA immune microenvironment, which helped to better understand the pathological mechanism of OA, with implications for more precise immunotherapy of OA in the future.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

The studies involving human participants were reviewed and approved by the Ethics Committee of Qilu Hospital of Shandong University. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author contributions

HY and YZ designed this study and reviewed the manuscript. ZL and HL contributed to bioinformatics analysis. DL and LM collected clinical samples and did statistical analysis. TL and HS prepared the manuscript and performed the experiments. All authors contributed to the article and approved the submitted version.

Funding

Our study was supported by the Shandong Province Key R&D Program Project (Grant number: 2019GSF108205).

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/fimmu.2023.1128459/full#supplementary-material

Supplementary Materials 1 | Datasets of GSE82107, GSE29746, GSE55235, GSE55457, GSE117999 and GSE98918 downloaded from GEO database.

Supplementary Figure 1 | The datasets were combined, then normalized the original data by the sva package in R.

Supplementary Figure 2 | Differences in other inflammatory response pathways between OA and normal samples.

References

1. Bijlsma JWJ, Berenbaum F, Lafeber FPJG. Osteoarthritis: An update with relevance for clinical practice. Lancet (2011) 377:2115–26. doi: 10.1016/S0140-6736(11)60243-2

CrossRef Full Text | Google Scholar

2. Hunter DJ, Bierma-Zeinstra S. Osteoarthritis. Lancet (2019) 393:1745–59. doi: 10.1016/S0140-6736(19)30417-9

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

4. Postler A, Goronzy J, Günther KP, Lange T, Redeker I, Schmitt J, et al. Which disease-related factors influence patients’ and physicians’ willingness to consider joint replacement in hip and knee OA? results of a questionnaire survey linked to claims data. BMC Musculoskelet Disord (2020) 21:352. doi: 10.1186/s12891-020-03368-1

CrossRef Full Text | Google Scholar

5. Sandell LJ. Etiology of osteoarthritis: Genetics and synovial joint development. Nat Rev Rheumatol (2012) 8:77–89. doi: 10.1038/nrrheum.2011.199

CrossRef Full Text | Google Scholar

6. Castañeda S, Roman-Blas JA, Largo R, Herrero-Beaumont G. Osteoarthritis: A progressive disease with changing phenotypes. Rheumatology (2014) 53:1–3. doi: 10.1093/rheumatology/ket247

CrossRef Full Text | Google Scholar

7. Portal-Núñez S, Esbrit P, Alcaraz MJ, Largo R. Oxidative stress, autophagy, epigenetic changes and regulation by miRNAs as potential therapeutic targets in osteoarthritis. Biochem Pharmacol (2016) 108:1–10. doi: 10.1016/j.bcp.2015.12.012

CrossRef Full Text | Google Scholar

8. Zhuo Q, Yang W, Chen J, Wang Y. Metabolic syndrome meets osteoarthritis. Nat Rev Rheumatol (2012) 8:729–37. doi: 10.1038/nrrheum.2012.135

CrossRef Full Text | Google Scholar

9. Lopes EBP, Filiberti A, Husain SA, Humphrey MB. Immune contributions to osteoarthritis. Curr Osteoporos Rep (2017) 15:593–600. doi: 10.1007/s11914-017-0411-y

CrossRef Full Text | Google Scholar

10. Woodell-May JE, Sommerfeld SD. Role of inflammation and the immune system in the progression of osteoarthritis. J Orthop Res (2020) 38(2):253–7. doi: 10.1002/jor.24457

CrossRef Full Text | Google Scholar

11. Utomo L, Bastiaansen-Jenniskens YM, Verhaar JAN, Van Osch GJVM. Cartilage inflammation and degeneration is enhanced by pro-inflammatory (M1) macrophages in vitro, but not inhibited directly by anti-inflammatory (M2) macrophages. Osteoarthritis Cartilage (2016) 24:2162–70. doi: 10.1016/j.joca.2016.07.018

CrossRef Full Text | Google Scholar

12. Sui X, Hu Y, Ren C, Cao Q, Zhou S, Cao Y, et al. METTL3-mediated mA is required for murine oocyte maturation and maternal-to-zygotic transition. Cell Cycle (2020) 19:391–404. doi: 10.1080/15384101.2019.1711324

CrossRef Full Text | Google Scholar

13. Liu L, Wang J, Sun G, Wu Q, Ma J, Zhang X, et al. m6A mRNA methylation regulates CTNNB1 to promote the proliferation of hepatoblastoma. Mol Cancer (2019) 18:188. doi: 10.1186/s12943-019-1119-7

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

15. Hu B, Wang X, Gu X, Zou C, Gao Z, Zhang H, et al. N-methyladenosine (mA) RNA modification in gastrointestinal tract cancers: roles, mechanisms, and applications. Mol Cancer (2019) 18:178. doi: 10.1186/s12943-019-1099-7

CrossRef Full Text | Google Scholar

16. Chen X, Zhang J, Zhu J. The role of m6A RNA methylation in human cancer. Mol Cancer (2019) 18:103. doi: 10.1186/s12943-019-1033-z

CrossRef Full Text | Google Scholar

17. Lin S, Choe J, Du P, Triboulet R, Gregory RI. The m (6)A methyltransferase METTL3 promotes translation in human cancer cells. Mol Cell (2016) 62:335–45. doi: 10.1016/j.molcel.2016.03.021

CrossRef Full Text | Google Scholar

18. Wu Y, Zhou C, Yuan Q. Role of DNA and RNA N6-adenine methylation in regulating stem cell fate. Curr Stem Cell Res Ther (2018) 13:31–8. doi: 10.2174/1574888X12666170621125457

CrossRef Full Text | Google Scholar

19. Dubey PK, Patil M, Singh S, Dubey S, Ahuja P, Verma SK, et al. Increased m6A-RNA methylation and FTO suppression is associated with myocardial inflammation and dysfunction during endotoxemia in mice. Mol Cell Biochem (2022) 477(1):129–41. doi: 10.1007/s11010-021-04267-2

CrossRef Full Text | Google Scholar

20. Shi W, Zheng Y, Luo S, Li X, Zhang Y, Meng X, et al. METTL3 promotes activation and inflammation of FLSs through the NF-κB signaling pathway in rheumatoid arthritis. Front Med (Lausanne) (2021) 8:607585. doi: 10.3389/fmed.2021.607585

CrossRef Full Text | Google Scholar

21. Wang J, Wang F, Ke J, Li Z, Xu C, Yang Q, et al. Inhibition of METTL3 attenuates renal injury and inflammation by alleviating TAB3 m6A modifications via IGF2BP2-dependent mechanisms. Sci Transl Med (2022) 14(640):eabk2709. doi: 10.1126/scitranslmed.abk2709

CrossRef Full Text | Google Scholar

22. Mapperley C, Lagemaat LN, Lawson H, Tavosanis A, Paris J, Campos J, et al. The mRNA m6A reader YTHDF2 suppresses proinflammatory pathways and sustains hematopoietic stem cell function. J Exp Med (2021) 218(3):e20200829. doi: 10.1084/jem.20200829

CrossRef Full Text | Google Scholar

23. Yu B, Zeng A, Liu H, Yang Z, Fu M. MiR-654-3p, reduced by the excessive ALKBH5, alleviated the inflammation in OA by targeting TNFRSF9, the trigger of the NF-κB pathway. Biochem Biophys Res Commun (2022) 634:30–9. doi: 10.1016/j.bbrc.2022.09.103

CrossRef Full Text | Google Scholar

24. Sang W, Xue S, Jiang Y, Lu H, Zhu L, Wang C, et al. METTL3 involves the progression of osteoarthritis probably by affecting ECM degradation and regulating the inflammatory response. Life Sci (2021) 278:119528. doi: 10.1016/j.lfs.2021.119528

CrossRef Full Text | Google Scholar

25. Cai D, Zhang J, Yang J, Lv Q, Zhong C. Overexpression of FTO alleviates osteoarthritis by regulating the processing of miR-515-5p and the TLR4/MyD88/NF-κB axis. Int Immunopharmacol (2023) 114:109524. doi: 10.1016/j.intimp.2022.109524

CrossRef Full Text | Google Scholar

26. Liu C, Yang Z, Li R, Wu Y, Chi M, Gao S, et al. Potential roles of N6-methyladenosine (m6A) in immune cells. J Transl Med (2021) 19:251. doi: 10.1186/s12967-021-02918-y

CrossRef Full Text | Google Scholar

27. Guo Y, Wang R, Li J, Song Y, Min J, Zhao T, et al. Comprehensive analysis of m6A RNA methylation regulators and the immune microenvironment to aid immunotherapy in pancreatic cancer. Front Immunol (2021) 12:769425. doi: 10.3389/fimmu.2021.769425

CrossRef Full Text | Google Scholar

28. Song S, Zhao R, Qiao J, Liu J, Cheng T, Zhang S, et al. Predictive value of drug efficacy by M6A modification patterns in rheumatoid arthritis patients. Front Immunol (2022) 13:940918. doi: 10.3389/fimmu.2022.940918

CrossRef Full Text | Google Scholar

29. Colletti A, Cicero AFG. Nutraceutical approach to chronic osteoarthritis: From molecular research to clinical evidence. Int J Mol Sci (2021) 22:12920. doi: 10.3390/ijms222312920

CrossRef Full Text | Google Scholar

30. Saklatvala J. Inflammatory signaling in cartilage: MAPK and NF-kappaB pathways in chondrocytes and the use of inhibitors for research into pathogenesis and therapy of osteoarthritis. Curr Drug Targets (2007) 8(2):305–13. doi: 10.2174/138945007779940115

CrossRef Full Text | Google Scholar

31. Liu Q, Li M, Jiang L, Jiang R, Fu B. METTL3 promotes experimental osteoarthritis development by regulating inflammatory response and apoptosis in chondrocyte. Biochem Biophys Res Commun (2019) 516(1):22–7. doi: 10.1016/j.bbrc.2019.05.168

CrossRef Full Text | Google Scholar

32. Zheng Q, Hou J, Zhou Y, Li Z, Cao X. The RNA helicase DDX46 inhibits innate immunity by entrapping m (6)A- demethylated antiviral transcripts in the nucleus. Nat Immunol (2017) 18:1094–1103. doi: 10.1038/ni.3830

CrossRef Full Text | Google Scholar

33. Li J, Rao B, Yang J, Liu L, Huang M, Liu X, et al. Dysregulated m6A-related regulators are associated with tumor metastasis and poor prognosis in osteosarcoma. Front Oncol (2020) 10:769. doi: 10.3389/fonc.2020.00769

CrossRef Full Text | Google Scholar

34. Tang R, Zhang Y, Liang C, Xu J, Meng Q, Hua J, et al. The role of m6A-related genes in the prognosis and immune microenvironment of pancreatic adenocarcinoma. PeerJ (2020) 8:e9602. doi: 10.7717/peerj.9602

CrossRef Full Text | Google Scholar

35. Xu S, Tang L, Dai G, Luo C, Liu Z. Expression of m6A regulators correlated with immune microenvironment predicts therapeutic efficacy and prognosis in gliomas. Front Cell Dev Biol (2020) 8:594112. doi: 10.3389/fcell.2020.594112

CrossRef Full Text | Google Scholar

36. Shao X, Dong J, Zhang H, Wu Y, Zheng L. Systematic analyses of the role of the reader protein of n 6-methyladenosine RNA methylation, YTH domain family 2, in liver hepatocellular carcinoma. Front Mol Biosci (2020) 7:577460. doi: 10.3389/fmolb.2020.577460

CrossRef Full Text | Google Scholar

37. Gao S, Sun Y, Zhang X, Hu L, Liu Y, Chua C, et al. IGFBP2 activates the NF-κB pathway to drive epithelial mesenchymal transition and invasive character in pancreatic ductal adenocarcinoma. Cancer Res (2016) 76:6543–54. doi: 10.1158/0008-5472.CAN-16-0438

CrossRef Full Text | Google Scholar

38. Liu Y, Song C, Shen F, Zhang J, Song S. IGFBP2 promotes immunosuppression associated with its mesenchymal induction and FcγRIIB phosphorylation in glioblastoma. PloS One (2019) 14(9):e0222999. doi: 10.1371/journal.pone.0222999

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

40. Zhang X, Zhang S, Yan X, Shan Y, Liu L, Zhou J, et al. m6A regulator-mediated RNA methylation modification patterns are involved in immune microenvironment regulation of periodontitis. J Cell Mol Med (2021) 25(7):3634–45. doi: 10.1111/jcmm.16469

CrossRef Full Text | Google Scholar

41. Sun K, Luo J, Guo J, Yao X, Jing X, Guo F. The PI3K/AKT/mTOR signaling pathway in osteoarthritis: a narrative review. Osteoarthritis Cartilage (2020) 28:400–9. doi: 10.1016/j.joca.2020.02.027

CrossRef Full Text | Google Scholar

42. Lin C, Shao Y, Zeng C, Zhao C, Fang H, Wang L, et al. Blocking PI3K/AKT signaling inhibits bone sclerosis in subchondral bone and attenuates post-traumatic osteoarthritis. J Cell Physiol (2018) 233:6135–47. doi: 10.1002/jcp.26460

CrossRef Full Text | Google Scholar

43. Xue JF, Shi ZM, Zou J, Li XL. Inhibition of PI3K/AKT/mTOR signaling pathway promotes autophagy of articular chondrocytes and attenuates inflammatory response in rats with osteoarthritis. BioMed Pharmacother (2017) 89:1252–61. doi: 10.1016/j.biopha.2017.01.130

CrossRef Full Text | Google Scholar