- Department of Orthopedics, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, China
Although N6-methyladenosine (m6A) modification is closely associated with the pathogenesis of osteoarthritis (OA), the mRNA profile of m6A modification in OA remains unknown. Therefore, our study aimed to identify common m6A features and novel m6A-related therapeutic targets in OA. In the present study, we identified 3962 differentially methylated genes (DMGs) and 2048 differentially expressed genes (DEGs) using methylated RNA immunoprecipitation next-generation sequencing (MeRIP-seq) and RNA-sequencing. A co-expression analysis of DMGs and DEGs showed that the expression of 805 genes was significantly affected by m6A methylation. Specifically, we obtained 28 hypermethylated and upregulated genes, 657 hypermethylated and downregulated genes, 102 hypomethylated and upregulated genes, and 18 hypomethylated and downregulated genes. The differential gene expression analysis based on GSE114007 revealed 2770 DEGs. The Weighted Gene Co-expression Network Analysis (WGCNA) based on GSE114007 identified 134 OA-related genes. By taking the intersection of these results, ten novel aberrantly expressed, m6A-modified and OA-related key genes were identified, including SKP2, SULF1, TNC, ZFP36, CEBPB, BHLHE41, SOX9, VEGFA, MKNK2 and TUBB4B. The present study may provide valuable insight into identifying m6A-related pharmacological targets in OA.
1 Introduction
Osteoarthritis (OA) is the most common type of arthritis and one of the leading causes of pain and disability worldwide. The primary features of OA include articular cartilage degeneration, synovial inflammation, and subchondral bone sclerosis (Peat and Thomas, 2021). Despite advances in OA treatment, there are still challenges in attenuating its devastating effects. One of these challenges is the lack of methods to reverse cartilage degradation and inhibit progression in advanced OA. The current effective treatment, joint replacement, causes financial, physical and psychological burden. To develop new therapies, it is urgently needed to elucidate the underlying molecular mechanism of OA pathogenesis (Quicke et al., 2022).
N6-methyladenosine (m6A) modification is the most common RNA epigenetic modification in mammalian cells. The epigenetic regulation of m6A methylation regulates gene expression and translation, therebyaffecting cell development and differentiation. Accumulating evidence indicates that m6A modification participates in the regulation of critical biological processes and pathogenesis of numerous human diseases (Zhao et al., 2020; Jiang et al., 2021). A variety of m6A-modified mRNAs have the potential to serve as novel biomarkers for diagnosis or therapeutic targets for chemotherapy (Ianniello et al., 2019; Liu K. et al., 2022; Zhang H. et al., 2022; Qiu et al., 2022). Furthermore, the importance of m6A modification is underscored by its evolutionary conservation across species and its presence in all eukaryotic organisms. The modification is also dynamically regulated by a complex interplay of enzymes, readers, and erasers, which further emphasizes its critical role in cellular processes.
Recently, the regulatory role of m6A methylation in OA has garnered increasing attention. The overall m6A modification was found to be elevated in IL-1β treated chondrocytes, and the dysregulation of m6A modification in chondrocytes led to the aberrant expression of genes involved in extracellular matrix (ECM) synthesis and degradation, which plays a critical role in the development and progression of OA (Liu et al., 2019). Moreover, several enzymes involved in m6A modification have been found to be differentially expressed in OA chondrocytes, suggesting that their dysregulation may contribute to the pathogenesis of OA. For example, m6A writer protein METTL3 was found to be upregulated in chondrocytes and associated with ECM degradation and proliferation repression (Ren et al., 2022). The eraser protein FTO was found to be downregulated in chondrocytes, resulting in increased apoptosis and ECM degradation (Yang et al., 2021). Several studies have revealed the pleiotropic roles of m6A modification in OA. M6A modification can regulate OA by mediating cellular senescence (Chen et al., 2022), apoptosis (Liu et al., 2019), inflammation response (Shi et al., 2022), and ECM degradation (Sang et al., 2021). Given the significant impact of m6A modification on chondrocyte metabolism, it is essential to search for and identify more core substrate genes of m6A modification in OA. These m6A-modified genes could be promising therapeutic targets and research focus to further elucidate the underlying mechanism of OA pathogenesis. However, until now, the transcriptome-wide m6A profile, the functions and pathways in which m6A-methylated genes are enriched, and the potential impacts of m6A modifications on mRNA expression in OA remain unknown.
In the present study, we collected three OA knee cartilage tissues and three normal knee cartilage tissues to obtain the first transcriptome-wide m6A modification profile of mRNAs in OA using methylated RNA immunoprecipitation next-generation sequencing (MeRIP-seq). Then, we performed a co-expression analysis to determine the correlation between m6A modification and mRNA expression in OA. Moreover, combining the results of weighted gene co-expression network analysis (WGCNA) based on GSE114007, we screened ten novel m6A-modified key genes in OA.
2 Materials and methods
2.1 Sample collection
Knee cartilage samples of medial condyle were collected from three patients who underwent knee arthroplasty and three patients who underwent thigh amputation due to trauma. Rheumatoid arthritis and other metabolic diseases were excluded in all patients. Clinical characteristics of included patients were shown in Table 1. Fresh samples were immediately frozen in liquid nitrogen and stored at −80°C to be detected. This study was approved by the institutional ethics board of the First affiliated Hospital of Zhengzhou University (2022-KY-0854-002).
2.2 MeRIP-sequencing
Total RNA was extracted using TRIzol™ Reagent (Invitrogen, 15596018). The concentration of total RNA was measured by Qubit RNA HS assay kit (Invitrogen, Q32852). 100μg total RNA was fragmented into 100-200 nt RNA fragments using 10X RNA Fragmentation Buffer (100 mM Tris-HCl, 100 mM ZnCl2 in nuclease-free H2O). The reaction was stopped by adding 10X EDTA (0.5M EDTA). Methylated RNA immunoprecipitation was performed using EpiTM m6A immunoprecipitation kit (Epibiotek, R1804). Briefly, the fragmented RNA was incubated with anti-m6A monoclonal antibody (Abcam, ab208577) for 3 h at 4°C and then with protein A/G magnetic beads (Invitrogen, 8880210002D/10004D) at 4°C for an additional 2 h to obtain immunoprecipitated RNA fragments. The m6A-enriched RNA was purified using TRIzol™ Reagent (Invitrogen, 15596018). The library was prepared by smart-seq method. Both the input samples without IP and the m6A IP samples were subjected to 150-bp, paired-end sequencing on an Illumina NovaSeq 6000 sequencer.
2.3 Data processing
Cutadapt (v2.5) was used to trim adapters and filter for sequences, and the remaining reads were then aligned to the human Ensemble genome GRCh38 (mouse Ensemble genome GRCm38) using Hisat2 aligner (v2.1.0). The m6A peaks on the mRNA of the three OA samples were combined to obtain the m6A peaks of the OA group, and the normal group was treated in the same way. The m6A peaks and the differential m6A peaks were identified using the exomePeak R package (v2.13.2) under the parameters “PEAK_CUTOFF_p-value = 0.05, PEAK_CUTOFF_FDR = NA, and FRAGMENT_LENGTH = 200.” The Gene Ontology (GO) project contains three parts: biological processes (BP), molecular functions (MF), and cellular components (CC). Differentially m6A-modified genes were used to perform GO functional analysis to annotate and speculate on the function of these differentially methylated genes. Pathway analysis using the Kyoto Encyclopedia of Genes and Genomes (KEGG) was conducted with differentially methylated genes for annotation and inference of the pathways they could be involved in. GO and KEGG analyses were performed using the clusterprofile R package (v3.6.0) with a p-value <0.05 were considered statistically significant. The m6A-RNA-related genomic features were visualized using the Guitar R package (v1.16.0). We identified m6A peaks with a p-value <0.05 for the de novo motif analysis using homer (v4.10.4). The differential gene expression analysis was performed using the DESeq2 R package with p-value <0.05 and fold change >2 as the cutoff criteria in the data processing of RNA-seq of input.
2.4 GEO datasets download and gene differential expression analysis
Datasets were downloaded from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo) database. The GSE114007 dataset includes 19 OA knee cartilage samples and 19 normal knee cartilage samples (Fisch et al., 2018). The differential expression analysis was performed using the wilcox. test by the limma package (v3.40.6). Differences were defined as statistically significant when the p-value <0.05 and fold change >2.
2.5 WGCNA analysis
The gene co-expression networks were performed by the WGCNA package based on the mRNA expression data of GSE114007. Briefly, a weighted correlation network was constructed by establishing a matrix of pairwise correlations between all the mRNAs across the measured samples. Then, the adjacency matrix was constructed using the power β = 13 as the soft-threshold. Next, the connectivity measure per gene was calculated based on the connection strengths with other mRNAs in the network. Subsequently, genes with strong co-expression relationships were acquired as modules by average linkage hierarchical clustering. After that, the gene expression profiles of each module were summarized by the module eigengene (ME), which was the first principal component of each module’s gene expression profile. The module with the highest correlation score was defined as the key module. Hub genes of the key module were identified by calculating the gene significance (GS) and module membership (MM). |GS| > 0.6 and |MM| > 0.8 were set as cutoff criteria. The detailed procedure of the WGCNA is shown in the WGCNA R package (Langfelder and Horvath, 2008).
3 Results
3.1 Flow chart
A flow chart of our study is shown in Figure 1.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1168365/fgene-14-1168365-HTML/image_m/fgene-14-1168365-g001.jpg)
FIGURE 1. Flow chart of identification of m6A-related key genes and construction of a diagnostic model.
3.2 Common features of the m6A modification of mRNA in OA and normal knee cartilage
In this study, six samples were detected for mRNA m6A methylation, including three OA knee cartilage samples and three normal knee cartilage samples. After removing the low-quality data, approximately 42111259 reads were acquired from every OA IP sample, and approximately 43334627 reads were acquired from every normal IP sample. In addition, the consensus sequence GGAC, the common m6A motif described in human diseases, was found to be the most significant peak using MOMER software (Figure 2A). Moreover, the m6A peaks were significantly related to different gene locations in both the OA and control samples: mostly in the 3′ UTR and the CDS region (Figure 2B).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1168365/fgene-14-1168365-HTML/image_m/fgene-14-1168365-g002.jpg)
FIGURE 2. Transcriptome-wide analysis of m6A peaks. (A) Common motif enriched in the m6A peak in the mRNA identified in OA and normal cartilage. (B) The distribution of m6A peaks in the length of mRNA between OA and normal cartilage.
3.3 Different m6A modifications between the OA sample and the control sample
We identified 8511 genes with 25295 peaks in the OA sample and 7541 genes with 17474 peaks in the control sample, indicating an overall hypermethylated level in the OA cartilage. On average, 2.97 and 2.31 m6A sites occurred per gene mutation in the OA and normal cartilage, respectively. In addition, 6172 genes were m6A methylated in both OA and control groups (Figure 3A). By comparing the differential m6A modifications in the OA and control samples, we identified 3120 differentially methylated genes (DMGs) containing hypermethylated sites and 1218 DMGs containing hypomethylated sites in the OA sample. Among them, 376 DMGs contained both hypermethylated and hypomethylated sites (Figure 3B, Supplementary Table S1). The results of the methylation heatmap and cluster analyses showed that: there were relative consistencies within the groups and marked differences between the groups (Figure 4). The top ten up- and down-methylated peaks in the OA cartilage and their related information are shown in Table 2.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1168365/fgene-14-1168365-HTML/image_m/fgene-14-1168365-g003.jpg)
FIGURE 3. Statistics for genes with m6AI modifications in OA and normal cartilage. (A) Number of genes with m6A sites in OA and normal samples. (B) Number of genes with hyper- and hypomethylated m6A sites in OA.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1168365/fgene-14-1168365-HTML/image_m/fgene-14-1168365-g004.jpg)
FIGURE 4. Cluster analysis of m6A methylation in OA and normal knee cartilage. Heatmap of differentially methylated m6A peaks in 3962 DMGs.
3.4 Gene enrichment analysis of DMGs
The GO analysis revealed that the DMGs were mainly involved in biological processes such as translational initiation, RNA catabolic processes, viral gene expression, and mRNA catabolic processes (Figure 5A). In terms of cellular component, these genes were mainly enriched in cell-substrate junctions, focal adhesions, cell-substrate adherens junctions, and ribosomes (Figure 5B). In terms of molecular function, these genes were mainly engaged in cadherin binding, cell adhesion molecule binding, ubiquitin-like protein ligase binding, and ubiquitin protein ligase binding (Figure 5C). In addition, the KEGG analysis showed that the DMGs were chiefly enriched in pathways including ribosomes, spliceosomes, RNA transport and the cell cycle (Figure 5D).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1168365/fgene-14-1168365-HTML/image_m/fgene-14-1168365-g005.jpg)
FIGURE 5. GO functional enrichment analysis and KEGG pathway analysis of DMGs. (A) Biological process annotation diagram (Top10). (B) Cellular component annotation diagram (Top10). (C) Molecular function annotation diagram (Top10). (D) KEGG annotation.
3.5 Analysis of the effects of differential m6A modification of genes on related mRNA expression
A total of 759 upregulated genes and 1289 downregulated genes were identified from the RNA-seq of input (Supplementary Table S2). The association analysis of these 2048 differentially expressed genes (DEGs) with 3962 DMGs showed that the expression of 805 genes was affected by m6A methylation, including 28 hypermethylated and upregulated genes, 657 hypermethylated and downregulated genes, 102 hypomethylated and upregulated genes, and 18 hypomethylated and downregulated genes (Figure 6, Supplementary Table S3).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1168365/fgene-14-1168365-HTML/image_m/fgene-14-1168365-g006.jpg)
FIGURE 6. The association analysis between DEGs and DMGs. 805 genes were identified, including 28 hypermethylated and upregulated genes, 657 hypermethylated and downregulated genes, 102 hypomethylated and upregulated genes, and 18 hypomethylated and downregulated genes. Genes with the top 10 m6A modification sites were annotated.
3.6 Identification of the m6A-modified key genes in OA
A total of 2770 DEGs, including 1049 downregulated and 1721 upregulated genes, were identified from GSE114007 under the thresholds of |log2FC| > 1 and p < 0.05 (Supplementary Table S4). In addition, to identify the hub genes correlated with OA, we conducted the WGCNA for the genes gained from GSE114007. When the soft thresholding was set at 13, the scale-free topology fit index reached 0.80 (Figure 7A) and 23 co-expressed gene modules were detected (Figure 7B). The module of midnight-blue showed the strongest significant correlations with OA (r = 0.88, p = 2e-13) (Figure 7C). A total of 134 genes were screened as hub genes in this module according to |MM|>0.8 and |GS|>0.6 (Supplementary Table S5). Finally, by overlapping the WGCNA results, the DEGs and the DMGs, we identified ten m6A-modified and OA-related key genes, including SKP2, SULF1, TNC, ZFP36, CEBPB, BHLHE41, SOX9, VEGFA, MKNK2 and TUBB4B (Figure 8; Table 3).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1168365/fgene-14-1168365-HTML/image_m/fgene-14-1168365-g007.jpg)
FIGURE 7. WGCNA analysis. (A) Identification of optimal soft threshold power for the co-expression network. (B) Clustering dendrogram of genes. (C) Module-trait relationships.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1168365/fgene-14-1168365-HTML/image_m/fgene-14-1168365-g008.jpg)
FIGURE 8. Venn diagram of the WGCNA, DEGs and DMGs. Ten m6A-modified key genes were identified by crossover.
4 Discussion
Despite m6A modification being a hot emerging field in OA pathogenesis with rapidly expanding interest, the full understanding of the molecular mechanisms of m6A modification in OA is still in its infancy. Most relevant studies were designed from the perspective of the m6A regulatory proteins. Among the studied regulators, the “writer” METTL3 is the most well-documented. It was shown that METTL3-mediated m6A modification attenuated the RNA stability of ATG7 and negatively regulated autophagy in the OA fibroblast-like synoviocytes (Chen et al., 2022). Another study showed that METTL3 accelerated OA by regulating the NF-kB signaling pathway and the extracellular matrix (ECM) degradation (Liu et al., 2019). METTL3 also catalyzes the m6A methylation of non-coding RNAs in OA pathogenesis. The expression of LINC00680 could be increased by METTL3 in an m6A-dependent manner in chondrocytes, facilitating OA progression (Ren et al., 2022). A study on temporomandibular joint arthritis showed that METTL3 inhibited apoptosis and autophagy of chondrocytes in inflammation through the m6A/YTHDF1/Bcl2 signaling axis, providing a promising therapeutic strategy for temporomandibular joint arthritis (He et al., 2022).
In our present study, we overviewed the transcriptome-wide m6A profile and obtained abundant m6A-modified mRNAs. These findings might offer researchers valuable options for further research. Additionally, we identified ten m6A-modified key genes through a series of bioinformatics methods, six of which are OA-related, including SULF1, CEBPB, VEGFA, SKP2, SOX9 and TNC. However, none of them has been studied from the m6A perspective in previous literature. Specifically, SULF1 is an extracellular heparan sulfate proteoglycan-specific 6-O-endosulfatasethat can increase type II collagen expression and decrease MMP-13 expression in chondrocytes (Otsuki et al., 2017). In SULF1 knockout mice, MMP-13 were elevated whereas type II collagen and aggrecan were reduced, resulting in more severe spontaneous cartilage degeneration than in wild-type mice (Otsuki et al., 2010). CEBPB is a transcription factor of the basic-leucine zipper class and plays critical roles in proliferation and differentiation of a variety of tissues (Niehrs and Calkhoven, 2020). It enhanced the promoter activity of MMP13, and knockout of CEBPB in mice caused resistance to OA with decreased MMP-13 expression (Hirata et al., 2012). Additionally, CEBPB was identified as a substrate of m6A modification in autoimmune conditions (Bechara et al., 2021) and cholangiocarcinoma (Zhu et al., 2022). VEGFA is an essential member of vascular endothelial growth factor family and a crucial regulator of vasculogenesis and angiogenesis, with multiple signaling pathways and factors involved in the regulation of its expression. M6A modification facilitated the stability of VEGFA by preventing its mRNA degradation in colorectal cancer (Zhang G. et al., 2022; Liu X. et al., 2022). Moreover, in bone mesenchymal stem cells, m6A modification participated in regulating the alternative splicing of VEGFA (Tian et al., 2019). However, how VEGFA is regulated by m6A in OA pathogenesis remains unknown, warranting further research. SKP2, also known as p45, is a kinase-associated protein that participate in cell apoptosis, proliferation, migration, and invasion in multiple types of human cancer (Wang et al., 2012; Cai et al., 2020). In OA pathogenesis, SKP2 was shown to modulate the MAPK signaling pathway and ECM degradation (Feng et al., 2020), and suppressing SKP2 by miR-337-3p could inhibit ubiquitination of DUSP1 and promote chondrocyte proliferation (Jian et al., 2021). SKP2 is an important ligase of ubiquitination, but its relationship with m6A methylation remains unclear. SOX9 is a critical transcription factor and master regulator in chondrogenesis (Song and Park, 2020), promoting the expression of ECM protein during chondrogenic differentiation (Yang et al., 2022). Elevated expression of SOX9 could alleviate the severity of OA, making it a promising therapeutic target in OA (Tong et al., 2021). Nevertheless, the molecular mechanisms regulating SOX9 expression are still largely unknown. In a testicular cancer cell line, phosphorylation of SOX9 promoted its nuclear translocation (Malki et al., 2005). In degenerative human endplate cartilage, METTL3 modulated SOX9 m6A modification and reduced its mRNA stability, thereby inhibiting the expression of the type II collagen (Xiao et al., 2022). These results suggest that the expression of SOX9 is epigenetically regulated. However, the epigenetic regulation of SOX9, including m6A modification, has not been explored in OA pathogenesis. TNC is a large molecular glycoprotein and component of the extracellular matrix with wide involvement in cartilage development and chondrogenesis (Chiquet and Fambrough, 1984). Its expression decreased along with the maturation of chondrocytes and almost disappeared in normal adult articular cartilage (Chevalier et al., 1994). TNC has been reported to induce inflammatory mediators and accelerate cartilage degradation (Patel et al., 2011; Chockalingam et al., 2013), but it also has the potential to repair cartilage defects (Matsui et al., 2018), with elevated expression in OA cartilage indicating a self-repair mechanism (Hattori et al., 2021). The dual roles of TNC may be derived from the versatile nature of this glycoprotein. Although TNC has been widely studied in OA, the regulatory role of m6A modification in TNC has never been addressed. Thus, further studies in the m6A perspective may help to elucidate the underlying mechanism of the two opposing roles of TNC in OA. In addition, four m6A-modified key genes found in our study, including MKNK2, TUBB4B, BHLHE41, and ZFP36, have never been studied in OA pathogenesis and may be valuable for further exploration as potential therapeutic targets according to our results.
Intriguingly, our study revealed that most hypermethylated mRNAs were downregulated, while most hypomethylated mRNAs were upregulated in OA cartilage. This finding may have implications for further understanding the role of m6A modification in OA. M6A modification exerts different functional mechanisms, including mRNA splicing, mRNA translation, mRNA degradation, mRNA stability, and mRNA maturation in different pathological processes (Jiang et al., 2021). Our results indicated that decreased mRNA stability and elevated mRNA decay, which generally result in reduced mRNA expression, might be essential processes involved in the regulatory metabolisms of many m6A-modified mRNAs in OA. On the other hand, m6A-binding proteins can recognize m6A structures and regulate the subsequent biological processes of mRNA. These proteins have been extensively studied in multiple disorders, including the YTH domain family (YTHDF1/2/3 and YTHDC1/2), insulin-like growth factor 2 mRNA-binding proteins (IGF2BP1-3), and heterogeneous nuclear ribonucleoproteins (HNRNPs, hnRNPC, hnRNPG, and hnRNPA2B1) (Feng et al., 2022). Among them, YTHDF2 has been the most well-recognized m6A-binding protein that destabilizes mRNAs and decreases the expression of target genes (Paris et al., 2019; Wang and Lu, 2021). However, the role of YTHDF2 in the pathogenesisof OA is still poorly understood. It could be assumed that YTHDF2 might be an important m6A-binding protein in OA, and further study of YTHDF2 in OA might be promising.
The balance between catabolic and anabolic signaling pathways is critical for maintaining cartilage homeostasis, and its disturbance contributes to OA. By GO and KEGG enrichment analyses, we found that the DMGs in OA cartilage were involved in RNA catabolic processes. This finding suggested that m6A modification affected metabolic processes that reflected the specific activities of the pathogenesis of OA. GO analysis also showed that the DMGs were correlated with the molecular functions of ubiquitin ligase. This is in accordance with previous observations in cancer. In osteosarcoma, an m6A modification disorder led to increased expression of the ubiquitin ligase RNF40, resulting in unchecked replication (Yadav et al., 2022). Ubiquitin ligases have been demonstrated to have an important impact on skeletal homeostasis by regulating the stability and function of various signaling factors (Liang et al., 2020). For instance, the overexpression of PARKIN, an E3 ubiquitin ligase, could accelerate the bone healing process (Zhang et al., 2020). PARKIN could also inhibit inflammation and aging-related metabolic changes in chondrocytes and attenuate OA progression through ubiquitination (Xu et al., 2020). Taken together, m6A RNA methylation may work in concert with other epigenetic mechanisms, such as ubiquitination, to participate in the pathogenesis of OA.
It should be acknowledged that our study has several limitations. First, the small sample size may affect the accuracy of the results. Second, our results are mainly based on high-throughput sequencing and bioinformatics analysis, and further molecular experiments are needed to verify our hypotheses.
5 Conclusion
In this study, we outlined the patterns and characteristics of m6A modification in OA and normal cartilage tissues for the first time. We obtained abundant m6A-modified sites in mRNAs and observed significant differences between the two groups. Then we identified ten m6A-modified key genes in OA pathogenesis by in-depth bioinformatics analysis. We believe our study offers a valuable reference for further investigation of m6A-related therapeutic targets in OA.
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 the First Affiliated Hospital of Zhengzhou University (2022-KY-0854-002). The patients/participants provided their written informed consent to participate in this study.
Author contributions
YY and JX designed the study; SL collected samples; YL contributed to data collection and statistical analysis; YY and SL wrote this manuscript. JX contributed to the critical revision of the final manuscript. All authors read and approved the final manuscript.
Funding
This study was supported by funds from the Youth Foundation of the First Affiliated Hospital of Zhengzhou University (YNQN-70131).
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/fgene.2023.1168365/full#supplementary-material
References
Bechara, R., Amatya, N., Bailey, R. D., Li, Y., Aggor, F. E. Y., Li, D. D., et al. (2021). The m(6)A reader IMP2 directs autoimmune inflammation through an IL-17- and TNFα-dependent C/EBP transcription factor axis. Sci. Immunol. 6 (61), eabd1287. doi:10.1126/sciimmunol.abd1287
Cai, Z., Moten, A., Peng, D., Hsu, C. C., Pan, B. S., Manne, R., et al. (2020). The Skp2 pathway: A critical target for cancer therapy. Semin. Cancer Biol. 67 (2), 16–33. doi:10.1016/j.semcancer.2020.01.013
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
Chevalier, X., Groult, N., Larget-Piet, B., Zardi, L., and Hornebeck, W. (1994). Tenascin distribution in articular cartilage from normal subjects and from patients with osteoarthritis and rheumatoid arthritis. Arthritis Rheum. 37 (7), 1013–1022. doi:10.1002/art.1780370706
Chiquet, M., and Fambrough, D. M. (1984). Chick myotendinous antigen. II. A novel extracellular glycoprotein complex consisting of large disulfide-linked subunits. J. Cell Biol. 98 (6), 1937–1946. doi:10.1083/jcb.98.6.1937
Chockalingam, P. S., Glasson, S. S., and Lohmander, L. S. (2013). Tenascin-C levels in synovial fluid are elevated after injury to the human and canine joint and correlate with markers of inflammation and matrix degradation. Osteoarthr. Cartil. 21 (2), 339–345. doi:10.1016/j.joca.2012.10.016
Feng, H., Yuan, X., Wu, S., Yuan, Y., Cui, L., Lin, D., et al. (2022). Effects of writers, erasers and readers within miRNA-related m6A modification in cancers. Cell Prolif. 56, e13340. doi:10.1111/cpr.13340
Feng, Y., Li, B., Li, S. J., Yang, X. C., Lv, T. T., Shang, H., et al. (2020). Skp2/p27 axis regulates chondrocyte proliferation under high glucose induced endoplasmic reticulum stress. Eur. Rev. Med. Pharmacol. Sci. 24 (17), 9129–9138. doi:10.26355/eurrev_202009_22860
Fisch, K. M., Gamini, R., Alvarez-Garcia, O., Akagi, R., Saito, M., Muramatsu, Y., et al. (2018). Identification of transcription factors responsible for dysregulated networks in human osteoarthritis cartilage by global gene expression analysis. Osteoarthr. Cartil. 26 (11), 1531–1538. doi:10.1016/j.joca.2018.07.012
Hattori, T., Hasegawa, M., Unno, H., Iino, T., Fukai, F., Yoshida, T., et al. (2021). TNIIIA2, the peptide of tenascin-C, as a candidate for preventing articular cartilage degeneration. Cartilage 13 (2), 1367S–1375S. doi:10.1177/1947603520912300
He, Y., Wang, W., Xu, X., Yang, B., Yu, X., Wu, Y., et al. (2022). Mettl3 inhibits the apoptosis and autophagy of chondrocytes in inflammation through mediating Bcl2 stability via Ythdf1-mediated m(6)A modification. Bone 154, 116182. doi:10.1016/j.bone.2021.116182
Hirata, M., Kugimiya, F., Fukai, A., Saito, T., Yano, F., Ikeda, T., et al. (2012). C/EBPβ and RUNX2 cooperate to degrade cartilage with MMP-13 as the target and HIF-2α as the inducer in chondrocytes. Hum. Mol. Genet. 21 (5), 1111–1123. doi:10.1093/hmg/ddr540
Ianniello, Z., Paiardini, A., and Fatica, A. (2019). N(6)-Methyladenosine (m(6)A): A promising new molecular target in acute myeloid leukemia. Front. Oncol. 9, 251. doi:10.3389/fonc.2019.00251
Jian, S., Luo, D., Wang, Y., Xu, W., Zhang, H., Zhang, L., et al. (2021). MiR-337-3p confers protective effect on facet joint osteoarthritis by targeting SKP2 to inhibit DUSP1 ubiquitination and inactivate MAPK pathway. Cell Biol. Toxicol. doi:10.1007/s10565-021-09665-2
Jiang, X., Liu, B., Nie, Z., Duan, L., Xiong, Q., Jin, Z., et al. (2021). The role of m6A modification in the biological functions and diseases. Signal Transduct. Target Ther. 6 (1), 74. doi:10.1038/s41392-020-00450-x
Langfelder, P., and Horvath, S. (2008). Wgcna: an R package for weighted correlation network analysis. BMC Bioinforma. 9, 559. doi:10.1186/1471-2105-9-559
Liang, C., Liang, G., Zheng, X., Huang, Y., Huang, S., and Yin, D. (2020). RSP5 positively regulates the osteogenic differentiation of mesenchymal stem cells by activating the K63-linked ubiquitination of akt. Stem Cells Int. 2020, 7073805. doi:10.1155/2020/7073805
Liu, K., Ouyang, Q. Y., Zhan, Y., Yin, H., Liu, B. X., Tan, L. M., et al. (2022a). Pharmacoepitranscriptomic landscape revealing m6A modification could be a drug-effect biomarker for cancer treatment. Mol. Ther. Nucleic Acids 28, 464–476. doi:10.1016/j.omtn.2022.04.001
Liu, Q., Li, M., Jiang, L., Jiang, R., and Fu, B. (2019). METTL3 promotes experimental osteoarthritis development by regulating inflammatory response and apoptosis in chondrocyte. Biochem. Biophys. Res. Commun. 516 (1), 22–27. doi:10.1016/j.bbrc.2019.05.168
Liu, X., He, H., Zhang, F., Hu, X., Bi, F., Li, K., et al. (2022b). m6A methylated EphA2 and VEGFA through IGF2BP2/3 regulation promotes vasculogenic mimicry in colorectal cancer via PI3K/AKT and ERK1/2 signaling. Cell Death Dis. 13 (5), 483. doi:10.1038/s41419-022-04950-2
Malki, S., Nef, S., Notarnicola, C., Thevenet, L., Gasca, S., Méjean, C., et al. (2005). Prostaglandin D2 induces nuclear import of the sex-determining factor SOX9 via its cAMP-PKA phosphorylation. Embo J. 24 (10), 1798–1809. doi:10.1038/sj.emboj.7600660
Matsui, Y., Hasegawa, M., Iino, T., Imanaka-Yoshida, K., Yoshida, T., and Sudo, A. (2018). Tenascin-C prevents articular cartilage degeneration in murine osteoarthritis models. Cartilage 9 (1), 80–88. doi:10.1177/1947603516681134
Niehrs, C., and Calkhoven, C. F. (2020). Emerging role of C/EBPβ and epigenetic DNA methylation in ageing. Trends Genet. 36 (2), 71–80. doi:10.1016/j.tig.2019.11.005
Otsuki, S., Hanson, S. R., Miyaki, S., Grogan, S. P., Kinoshita, M., Asahara, H., et al. (2010). Extracellular sulfatases support cartilage homeostasis by regulating BMP and FGF signaling pathways. Proc. Natl. Acad. Sci. U. S. A. 107 (22), 10202–10207. doi:10.1073/pnas.0913897107
Otsuki, S., Murakami, T., Okamoto, Y., Hoshiyama, Y., Oda, S., and Neo, M. (2017). Suppression of cartilage degeneration by intra-articular injection of heparan sulfate 6-O endosulfatase in a mouse osteoarthritis model. Histol. Histopathol. 32 (7), 725–733. doi:10.14670/hh-11-838
Paris, J., Morgan, M., Campos, J., Spencer, G. J., Shmakova, A., Ivanova, I., et al. (2019). Targeting the RNA m(6)A reader YTHDF2 selectively compromises cancer stem cells in acute myeloid leukemia. Cell Stem Cell 25 (1), 137–148. doi:10.1016/j.stem.2019.03.021
Patel, L., Sun, W., Glasson, S. S., Morris, E. A., Flannery, C. R., and Chockalingam, P. S. (2011). Tenascin-C induces inflammatory mediators and matrix degradation in osteoarthritic cartilage. BMC Musculoskelet. Disord. 12, 164. doi:10.1186/1471-2474-12-164
Peat, G., and Thomas, M. J. (2021). Osteoarthritis year in review 2020: Epidemiology and therapy. Osteoarthr. Cartil. 29 (2), 180–189. doi:10.1016/j.joca.2020.10.007
Qiu, F. S., He, J. Q., Zhong, Y. S., Guo, M. Y., and Yu, C. H. (2022). Implications of m6A methylation and microbiota interaction in non-small cell lung cancer: From basics to therapeutics. Front. Cell Infect. Microbiol. 12, 972655. doi:10.3389/fcimb.2022.972655
Quicke, J. G., Conaghan, P. G., Corp, N., and Peat, G. (2022). Osteoarthritis year in review 2021: Epidemiology and therapy. Osteoarthr. Cartil. 30 (2), 196–206. doi:10.1016/j.joca.2021.10.003
Ren, J., Li, Y., Wuermanbieke, S., Hu, S., and Huang, G. (2022). N(6)-methyladenosine (m(6)A) methyltransferase METTL3-mediated LINC00680 accelerates osteoarthritis through m(6)A/SIRT1 manner. Cell Death Discov. 8 (1), 240. doi:10.1038/s41420-022-00890-0
Sang, W., Xue, S., Jiang, Y., Lu, H., Zhu, L., Wang, C., et al. (2021). METTL3 involves the progression of osteoarthritis probably by affecting ECM degradation and regulating the inflammatory response. Life Sci. 278, 119528. doi:10.1016/j.lfs.2021.119528
Shi, L., Hu, H., Sun, P., Li, Z., Ji, L., Liu, S., et al. (2022). RPL38 knockdown inhibits the inflammation and apoptosis in chondrocytes through regulating METTL3-mediated SOCS2 m6A modification in osteoarthritis. Inflamm. Res. 71 (7-8), 977–989. doi:10.1007/s00011-022-01579-x
Song, H., and Park, K. H. (2020). Regulation and function of SOX9 during cartilage development and regeneration. Semin. Cancer Biol. 67 (1), 12–23. doi:10.1016/j.semcancer.2020.04.008
Tian, C., Huang, Y., Li, Q., Feng, Z., and Xu, Q. (2019). Mettl3 regulates osteogenic differentiation and alternative splicing of vegfa in bone marrow mesenchymal stem cells. Int. J. Mol. Sci. 20 (3), 551. doi:10.3390/ijms20030551
Tong, J., Deng, C., Sun, G., Zhou, J., Zhong, P., Wang, T., et al. (2021). Electroacupuncture upregulates HIF-1α and SOX9 expression in knee osteoarthritis. Evid. Based Complement. Altern. Med. 2021, 2047097. doi:10.1155/2021/2047097
Wang, G., Chan, C. H., Gao, Y., and Lin, H. K. (2012). Novel roles of Skp2 E3 ligase in cellular senescence, cancer progression, and metastasis. Chin. J. Cancer 31 (4), 169–177. doi:10.5732/cjc.011.10319
Wang, J. Y., and Lu, A. Q. (2021). The biological function of m6A reader YTHDF2 and its role in human disease. Cancer Cell Int. 21 (1), 109. doi:10.1186/s12935-021-01807-0
Xiao, L., Hu, B., Ding, B., Zhao, Q., Liu, C., Öner, F. C., et al. (2022). N(6)-methyladenosine RNA methyltransferase like 3 inhibits extracellular matrix synthesis of endplate chondrocytes by downregulating sex-determining region Y-Box transcription factor 9 expression under tension. Osteoarthr. Cartil. 30 (4), 613–625. doi:10.1016/j.joca.2022.01.002
Xu, L., Wu, Z., He, Y., Chen, Z., Xu, K., Yu, W., et al. (2020). MFN2 contributes to metabolic disorders and inflammation in the aging of rat chondrocytes and osteoarthritis. Osteoarthr. Cartil. 28 (8), 1079–1091. doi:10.1016/j.joca.2019.11.011
Yadav, P., Subbarayalu, P., Medina, D., Nirzhor, S., Timilsina, S., Rajamanickam, S., et al. (2022). M6A RNA methylation regulates histone ubiquitination to support cancer growth and progression. Cancer Res. 82 (10), 1872–1889. doi:10.1158/0008-5472.Can-21-2106
Yang, J., Zhang, M., Yang, D., Ma, Y., Tang, Y., Xing, M., et al. (2021). m(6)A-mediated upregulation of AC008 promotes osteoarthritis progression through the miR-328-3p‒AQP1/ANKH axis. Exp. Mol. Med. 53 (11), 1723–1734. doi:10.1038/s12276-021-00696-7
Yang, L., Ren, Z., Yan, S., Zhao, L., Liu, J., Zhao, L., et al. (2022). Nsun4 and Mettl3 mediated translational reprogramming of Sox9 promotes BMSC chondrogenic differentiation. Commun. Biol. 5 (1), 495. doi:10.1038/s42003-022-03420-x
Zhang, G., Wang, T., Huang, Z., Chen, Y., Sun, L., Xia, X., et al. (2022a). METTL3 dual regulation of the stability of LINC00662 and VEGFA RNAs promotes colorectal cancer angiogenesis. Discov. Oncol. 13 (1), 89. doi:10.1007/s12672-022-00557-3
Zhang, H., Huang, T., Ren, X., Fang, X., Chen, X., Wei, H., et al. (2022b). Integrated pan-cancer analysis of CSMD2 as a potential prognostic, diagnostic, and immune biomarker. Front. Genet. 13, 918486. doi:10.3389/fgene.2022.918486
Zhang, W., Hou, W., Chen, M., Chen, E., Xue, D., Ye, C., et al. (2020). Upregulation of parkin accelerates osteoblastic differentiation of bone marrow-derived mesenchymal stem cells and bone regeneration by enhancing autophagy and β-catenin signaling. Front. Cell Dev. Biol. 8, 576104. doi:10.3389/fcell.2020.576104
Zhao, W., Qi, X., Liu, L., Ma, S., Liu, J., and Wu, J. (2020). Epigenetic regulation of m(6)A modifications in human cancer. Mol. Ther. Nucleic Acids 19, 405–412. doi:10.1016/j.omtn.2019.11.022
Keywords: osteoarthritis, m6A methylation, methylome profile, WGCNA, merip-seq
Citation: Yu Y, Lu S, Li Y and Xu J (2023) Overview of distinct N6-Methyladenosine profiles of messenger RNA in osteoarthritis. Front. Genet. 14:1168365. doi: 10.3389/fgene.2023.1168365
Received: 17 February 2023; Accepted: 26 April 2023;
Published: 09 May 2023.
Edited by:
Wieslawa Lesniak, Polish Academy of Sciences, PolandReviewed by:
Piyush Khandelia, Birla Institute of Technology and Science, IndiaXiaoyan Shao, Nanjing University, China
Copyright © 2023 Yu, Lu, Li and Xu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jianzhong Xu, xjzzzu@yeah.net