Skip to main content

Identification and analysis of RNA-5-methylcytosine-related key genes in osteoarthritis

Abstract

Background

5-methylcytosine (m5C) modification is widely associated with many biological and pathological processes. However, knowledge of m5C modification in osteoarthritis (OA) remains lacking. Thus, our study aimed to identify common m5C features in OA.

Results

In the present study, we identified 1395 differentially methylated genes (DMGs) and 1673 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 133 genes was significantly affected by m5C methylation. A protein–protein interaction network of the 133 genes was constructed using the STRING database, and the cytoHubba plug-in of Cytoscape was used to hub genes were screen out 11 hub genes, including MMP14, VTN, COL15A1, COL6A2, SPARC, COL5A1, COL6A3, COL6A1, COL8A2, ADAMTS2 and COL7A1. The Pathway enrichment analysis by the ClueGO and CluePedia plugins in Cytoscape showed that the hub genes were significantly enriched in collagen degradation and extracellular matrix degradation.

Conclusions

Our study indicated that m5C modification might play an important role in OA pathogenesis, and the present study provides worthwhile insight into identifying m5C-related therapeutic targets in OA.

Peer Review reports

Introduction

Osteoarthritis (OA) is the most common type of arthritis and one of the leading causes of pain and disability worldwide. It is estimated that approximately 27–31 million people in the US have symptomatic OA [1]. The economic costs related to OA are predicted to be between 1% and 2.5% of the gross domestic product (GDP) of developed countries [2]. Despite marked progress in OA treatment, there is currently no effective disease-modifying therapy available to relieve the pain and inhibit the progression of OA other than joint arthroplasty. This is largely due to the incomplete understanding of the pathogenesis of OA and the deficiency of effective therapeutic targets. Thus, the further elucidation of the underlying molecular mechanism of OA is urgently needed for the development of novel therapies.

Epigenetic modification plays critical regulatory roles in various biological processes of eukaryotic cells. Hundreds of modifications have been discovered in RNA, and among these abundant modifications, 5-methylcytosine (m5C) methylation has received great attention. The m5C modification participates in multiple biological processes. It was found that ALYREF promoted mRNA export in an m5C-dependent manner, and the knockdown of NSUN2, an m5C methyltransferase, significantly reduced the cytoplasmic-to-nuclear ratios of mRNAs [3]. Yang et al. demonstrated that RNA m5C modification regulated mRNA stabilization during the maternal-to-zygotic transition in zebrafish [4]. Chen et al. found that DNA damage could induce m5C modification in mRNAs at sites of DNA damage, and this is an important mechanism of DNA repair [5]. Aberrant mRNA m5C modifications are also widely involved in the pathogenesis of a variety of diseases. The m5C modification level significantly affected the proliferation and migration of HEK293 cells [6]. Chen et al. demonstrated that m5C was significantly hypermethylated in bladder cancer and primarily enriched in oncogenic pathways. The m5C modification participated in bladder cancer pathogenesis by regulating mRNA stability [7]. A study by Luo et al. showed that decreased m5C methylation (caused by the inhibition of NSUN2) can reverse atherosclerosis and endothelial inflammation after aortic transplantation [8]. Similar to N6-methyladenosine (m6A), which is the most prevalent internal mRNA modification, m5C methylation is dynamically regulated by methyltransferase or “writers” (which catalyze the deposition of m5C), demethylase or “erasers” (which catalyze the removal of m5C), and binding protein or “readers” (which recognize the modified nucleotides). While the list of m5C readers and erasers is still under debate, the m5C writers have been acknowledged and mainly include the NOL1/NOP2/sun (NSUN) family and DNA methyltransferase 2 (DNMT2) [9]. Studies have indicated that the m5C regulators modulate the expression of a variety of oncogenes in an m5C-dependent manner [10]. An illustration of the regulatory role of m5C modification in the pathogenesis of multiple diseases might facilitate progressions in treatment. A study by Lyko et al. showed that the reduced m5C modification of tRNA by DNMT2 inhibition suppressed metabolic activity in the human cancer cell line [11]. NSUN2 promoted gastric cancer development was promoted by CDKN1C repression, and the knockdown of NSUN2 suppressed tumor growth [12].

The study of the functions of m5C modification in diverse pathological processes is gaining increasing attention and has significant implications. However, the regulatory role of m5C modification in OA pathogenesis remains unknown. With advances in high-throughput technologies, researchers have successfully used m5C-RNA immunoprecipitation (IP) to characterize m5C methylation sites in several types of cancers, such as hepatocellular carcinoma [13] and ovarian cancer [14], thereby laying a foundation for the further elucidation of the mechanism. In the present study, we collected three OA knee cartilage tissues and three normal knee cartilage tissues to obtain the first transcriptome-mRNA m5C modification profile by methylated RNA immunoprecipitation next-generation sequencing (MeRIP-seq). Significant differences in the numbers and distributions of m5C peaks were noted in OA and normal cartilage tissues. Differentially methylated genes (DMGs) were identified, and gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed. We identified 11 m5C-modified hub genes via bioinformatics analyses. The pathway enrichment analysis of hub genes was performed using the ClueGO and CluePedia plugins in Cytoscape.

Results

Flow chart

A flow chart of our study is shown in Fig. 1.

Fig. 1
figure 1

Flow chart of identification and analysis of m5C-related hub genes in OA

General features of m5C methylation map in OA and normal knee cartilage

In general, we found 33,272 m5C peaks and mapped up to 11,714 annotated genes in OA samples (Fig. 2A), and we found 16,145 m5C peaks and mapped up to 7502 annotated genes in normal samples (Fig. 2B); Peak information for each sample has been listed in Additional file 1: Table S1.

Fig. 2
figure 2

Overview of m5C peaks in mRNAs and differential m5C methylation analysis. A Number of m5C peaks in OA and normal groups. B Number of m5C methylated genes in OA and normal groups. C Preferential location of m5C in mRNA transcripts. Each transcript was divided into three parts: 5’untranslated region, coding DNA sequence, and 3’untranslated region. D The number of m5C peaks on each mRNA. Most mRNAs had only one m5C peak in both OA and normal samples. E Visualization of m5C at the chromosome level in OA and normal samples. F Number of genes with hyper- and hypomethylated m5C sites in OA cartilage

To identify the preferential location of m5C in the mRNA, we analyzed the metagene profiles of the m5C peaks in the mRNA transcriptome. The m5C peaks were distributed in all regions of the mRNA, but mostly in the CDS region. The distribution in OA samples was similar to that in the normal samples at the 5′UTR and 3′UTR. The m5C peaks in OA samples were more distributed near the stop codon, and the m5C peaks in the normal samples were more distributed near the start codon (Fig. 2C). We also conducted statistical analyses of the number of peaks on each methylated mRNA. When the number of peaks detected for the same gene in different samples was inconsistent, we selected the maximum count for analysis. In the OA group, the number of genes with one, two, three, four, five, and more than five peaks are 7263, 2577, 703, 351, 176, and 234, respectively. In the normal group, the number of genes with one, two, three, four, five, and more than five peaks are 4014, 1913, 600, 375, 225, and 338, respectively. We found that most mRNAs in the OA and normal samples had only one m5C peak (Fig. 2D). We used Circos software to determine the distribution of m5C at the chromosome level. It was found that the numbers and distributions of m5C peaks on each chromosome were different between the OA and normal knee cartilage tissues (Fig. 2E).

We used DiffReps software to identify 1607 upregulated m5C methylation peaks in 1215 mRNAs and 326 downregulated m5C methylation peaks in 286 mRNAs. Among these genes, 1109 genes and 180 genes had only upregulated m5C peaks or downregulated m5C peaks, respectively, and a total of 106 genes showed both up- and downregulated m5C peaks (Fig. 2F, Additional file 2: Table S2).

GO and KEGG pathway analysis

GO analysis indicated that the mRNAs with upregulated m5C peaks in OA cartilage were mainly enriched in the cellular component organization or biogenesis (GO term: BP) (Fig. 3A), intracellular and intracellular organelle (GO term: CC) (Fig. 3B), and protein binding and enzyme binding (GO term: MF) (Fig. 3C) categories. The KEGG analysis results showed that these mRNAs were primarily involved in focal adhesion and proteoglycans in cancer (Fig. 3D).

Fig. 3
figure 3

GO functional enrichment analysis and KEGG pathway analysis of mRNAs with upregulated m5C modification. A Biological process annotation diagram. B Cellular component annotation diagram. C Molecular function annotation diagram. D KEGG annotation

GO analysis indicated that the mRNAs with downregulated m5C peaks were mainly enriched in the mRNA metabolic process and mRNA processing (GO term: BP) (Fig. 4A), membrane-bounded organelle and intracellular (GO term: CC) (Fig. 4B), and heterocyclic compound binding and organic cyclic compound binding (GO term: MF) (Fig. 4C) categories. The KEGG analysis results showed that these mRNAs were significantly enriched in the PI3K-Akt signaling pathway and spliceosome (Fig. 4D).

Fig. 4
figure 4

GO functional enrichment analysis and KEGG pathway analysis of mRNAs with downregulated m5C modification. A Biological process annotation diagram. B Cellular component annotation diagram. C Molecular function annotation diagram. D KEGG annotation

Association analysis of MeRIP-seq and RNA-seq data

In total, 1673 differentially expressed genes (DEGs) were identified in the RNA-Seq data (Additional file 3: Table S3). To explore the effect of m5C methylation on transcriptional expression in OA, we conducted a conjoint analysis of MeRIP-Seq data and RNA-Seq data. The expressions of 133 genes were affected by m5C methylation, i.e., 89 hypermethylated and upregulated genes, 24 hypermethylated and downregulated genes, 10 hypomethylated and upregulated genes, and 10 hypomethylated and downregulated genes (Fig. 5, Additional file 4: Table S4).

Fig. 5
figure 5

The association analysis between DEGs and DMGs. 133 genes were identified, including 89 hypermethylated and upregulated genes, 24 hypermethylated and downregulated genes, 10 hypomethylated and upregulated genes, and 10 hypomethylated and downregulated genes

Identification and the pathway enrichment analysis of hub genes

The STRING database was used to construct a PPI network of the 133 selected genes (Fig. 6). We used the cytoHubba plug-in of Cytoscape to further analyze the PPI network. Node scores were calculated using the three algorithms of cytoHubba, and the top 15 genes were selected as key genes (Fig. 7A-C). Taking the intersections of the key genes obtained by these algorithms, 11 genes, i.e., MMP14, VTN, COL15A1, COL6A2, SPARC, COL5A1, COL6A3, COL6A1, COL8A2, ADAMTS2 and COL7A1, were identified as hub genes (Fig. 7D). All 11 genes were upregulated and hypermethylated in OA cartilage. The m5C-modified sites of these genes were shown in Table 1.

Fig. 6
figure 6

A protein–protein interaction network of the 133 genes constructed using the STRING database

Fig. 7
figure 7

Identification of the m5C-modified hub genes in OA. A Top 15 genes ranked by Closeness algorithm. B Top 15 genes ranked by Degree algorithm. C Top 15 genes ranked by MNC algorithm. D Venn diagram of the identified genes by the three algorithms

The Reactome function enrichment analysis using the ClueGO and CluePedia plug-in of Cytoscape showed that these hub genes were mainly enriched in collagen degradation and extracellular matrix (ECM) proteoglycans (Fig. 8).

Fig. 8
figure 8

The function enrichment of the 11 hub genes by the ClueGO and CluePedia plugins in Cytoscape

Table 1 The m5C-methylated peaks in 11 hub genes

Discussion

In the present study, we have given an overview of the distinct m5C modification profiles of mRNAs in OA and normal cartilage for the first time. We found abundant m5C methylation sites on the mRNAs of knee cartilage tissues, especially in degenerative cartilage, indicating an overall upregulated m5C methylation level in OA. This is consistent with previous studies in multiple tumor tissues, which also found rich depositions of m5C in transcriptomes [7, 14, 15]. However, a study by Legrand et al. only found a few m5C methylated mRNAs in the mouse transcriptome, and this suggested that mRNAs were sparsely methylated [16]. Huang et al. collected seven types of human tissue and 11 types of mouse tissue to examine the landscape of m5C in mammals and found a wide range of methylation levels in different tissues [17]. Thus, it can be assumed that the m5C modification levels in mRNAs vary across species and tissues. Our study has indicated that tissue specificity endows OA cartilage with rich m5C site distribution in mRNAs, which is an important basis for m5C to play a critical role in OA pathogenesis. On the other hand, it is still challenging to detect the transcriptome-wide m5C profile accurately and systematically. Several detection methods have been reported for identifying the mRNA m5C map, including MeRIP-Seq, Aza-IP [18], miCLIP [19], and RNA bisulfite sequencing (BS-seq) [20]. These methods have their own benefits and limitations, and no consensus has been reached on the best option. The variable detection results may be partly due to the method selection. These controversial findings prompt the development of more robust techniques for accurate transcriptome-wide m5C detection in mRNAs.

In addition to mRNAs, m5C modification is also widely present in other regulatory non-coding RNAs. Several studies have revealed substantially elevated m5C levels of lncRNAs and circRNAs in carcinoma tissues compared with adjacent non-tumor tissues [15, 21]. In fact, most of the m5C-modified RNAs are ribosomal RNAs (rRNAs) and transfer RNA (tRNAs), and the deposition of m5C in these RNAs imparts important biological functions [22]. For instance, it was found that m5C modification was essential to stabilizing the tRNA secondary structure and maintaining metabolic stability [23]. Changed m5C levels in rRNA could modulate the respective lifespans in yeast, worms and flies [24]. The unmethylated status of 28 S rRNA caused by the loss of an RNA methyltransferase for m5C increased survival under cellular stress [25]. Nevertheless, m5C-related studies on noncoding RNAs in OA pathogenesis are still lacking, and further exploration should be conducted in this area.

KEGG analysis revealed that genes with downregulated m5C peaks were enriched in the PI3K-Akt signaling pathway and the hypoxia-inducible factor-1 (HIF-1) signaling pathway. These two signaling pathways are crucial to OA pathogenesis. Articular cartilage is maintained in a low-oxygen environment in the lifetime of an individual, and HIF-1α is essential for the adaptation of chondrocytes to hypoxia [26]. Thus, HIF-1α is an important protective factor in chondrocyte homeostasis. The PI3K-Akt pathway is also critical in both the normal metabolism of cartilage and OA pathogenesis. Previous studies showed that the PI3K-Akt pathway extensively participates in the inflammatory response, chondrocyte proliferation, apoptosis, and autophagy [27, 28]. The multiple functions of this axis make it difficult to verify the general effects of this signaling pathway in OA, thereby setting up obstacles to the exploration of related therapeutic targets. The dominant effects or the cross-talk among the diverse effects of PI3K-Akt in OA are the prominent issues that need to be addressed. Our study has indicated that m5C modification might participate in the regulation of these signaling pathways in OA and might provide a new perspective for solving the above issues.

We identified 11 m5C-modified hub genes in OA cartilage via bioinformatic analysis, and these hub genes were primarily enriched in collagen degradation and ECM degradation. Progressive ECM degradation is one of the hallmarks of OA. Our results indicate that m5C modification might play an important role in this catabolic process. COL15A1, COL6A2, COL5A1, COL6A3, COL6A1, COL8A2, and COL7A1 are collagen-encoding genes. Multiple studies have shown that these genes are markedly upregulated in OA pathogenesis [29,30,31]. Type II collagen is a major component of normal articular cartilage, and the elevated expression of ECM turnover-related genes is indicative of ECM component remodeling and aberrant collagen accumulation. It was shown that aberrant collagen-encoding gene expression might contribute to fibrosis [32]. It was also indicated that the ECM turnover-related genes were downstream of the TGF-β signaling pathway in OA-related synovial fibrosis, and OA cartilage damage could stimulate elevated TGF-β expression, thereby facilitating fibrosis [33]. Our results have shown that m5C modification might be involved in the activation of aberrant collagen-encoding genes and suggest another regulatory mechanism of ECM component remodeling in OA cartilage.

MMP14 and ADAMTS2 are well-known metalloproteinases. Their elevated expression considerably reduces some of the dominating components of the ECM in cartilage, such as Type II collagen and aggrecan, thereby promoting OA progression [34]. Due to the important role of metalloproteinases in OA pathogenesis, many metalloproteinase inhibitors have been developed as disease-modifying therapeutics, and some of them have entered clinical trials [35, 36]. Among them, only doxycycline has been approved by the FDA. Nevertheless, the effect of doxycycline on reducing OA symptoms is still controversial. Hsien-Tsung Lu et al. found that injectable hyaluronic acid doxycycline significantly ameliorated the progression of OA [37]. However, a triple-blinded randomized controlled trial showed that doxycycline was ineffective in OA symptom relief in the short and medium term [38]. Further studies on the regulatory mechanism of metalloproteinases in OA pathogenesis are needed for the development of more effective and specific inhibitors. It has been shown that alterations in epigenetic modification can affect the vitality and expression of metalloproteinases. The ubiquitination of MMP14 by FBXO6, a ubiquitin ligase subunit, decreased its expression, leading to inhibited OA development [39]. ADAMTS2 was also recognized as a crucial ECM-related regulator in pancreatic cancer. The knockdown of the m6A demethylase FTO significantly increased the m6A level and decreased the mRNA expression of ADAMTS2 [40]. As the m5C modification effects of these ECM-related protease have never been clarified, we can only make some suggestions based on our findings that the elevated m5C modification of these proteases might contribute to their increased expression and aberrant activity in OA pathogenesis.

It should be acknowledged that our study had several limitations. First, the small sample size may affect the accuracy of the results. Variables other than OA may have an impact on the results, and larger sample sizes can help minimize the bias introduced by these variables. Second, our results were mainly based on high-throughput sequencing and bioinformatics analysis, and further molecular experiments are needed to verify our hypotheses. Despite these defects, we believe our study offers a valuable reference for further investigation of m5C-related therapeutic targets in OA.

Conclusion

In this study, we provided the first overview of the patterns and characteristics of m5C modifications in OA and normal cartilage tissue. We obtained a rich set of m5C modification sites in mRNA and observed significant differences between the two groups. We believe that our study provides valuable reference for further research on therapeutic targets related to m5C in OA.

Materials and methods

Sample collection and ethical approval

Knee cartilage samples of medial condyles were collected from three patients who underwent knee arthroplasty due to advanced OA and three patients who underwent thigh amputation due to trauma. Clinical characteristics of included patients were shown in Table 2. Fresh samples were immediately frozen in liquid nitrogen and stored at -80 ℃ for detection. This study was approved by the institutional ethics board of the First affiliated Hospital of Zhengzhou University.

Table 2 Clinical characteristics of included patients

RNA extraction and preparation

The total RNA from each sample was extracted using TRIzol reagent (Invitrogen Corporation, Carlsbad, CA). The content of rRNA in the total RNA was reduced using the Ribo-Zero rRNA Removal Kit (Illumina, Inc., CA, USA). The concentration of total RNA was measured using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). RNA with an OD260/OD280 ratio range of 1.8–2.1 was marked as acceptable.

MeRIP-Seq construction and sequencing

The m5C-IP-Seq service was provided by CloudSeq Inc. (Shanghai, China). Total RNA was subjected to immunoprecipitation with the GenSeq® m5C-IP Kit (GenSeq Inc.) by following the manufacturer’s instructions. RNA libraries for IP and input samples were then constructed with the GenSeq® Low Input Whole RNA Library Prep Kit (GenSeq, Inc.) by following the manufacturer’s instructions. Libraries were qualified using the Agilent 2100 bioanalyzer and then sequenced on a NovaSeq platform.

Sequencing data analysis

Briefly, paired-end reads were harvested from the Illumina Novaseq 6000 sequencer, and Q30 was used for quality control. After 3’ adaptor-trimming and low-quality read removal via cutadapt software (v1.9.3), clean reads of all libraries were aligned to the reference genome (UCSC MM10) using the Hisat2 software (v2.0.4). Methylated sites on RNAs (peaks) were identified using MACS software. Differentially methylated sites were identified by diffReps with a fold change threshold > 2 and P-value < 0.0001. The differentially expressed genes were identified using the DESeq2 R package with a P-value < 0.05 and fold change > 2 as the cutoff criteria. The Gene Ontology (GO) project contains three parts: biological processes (BP), molecular functions (MF), and cellular components (CC). Differentially m5C-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) [41]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.

Identification and function enrichment analysis of hub genes

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database was used to analyze the PPI network of the selected genes, and Cytoscape (v3.9.1) was used to identify the network modules in the resulting file. The degree, closeness and MNC algorithms of the cytoHubba plugin in Cytoscape (v3.9.1) were used to calculate scores of nodes in the PPI network. The hub genes were determined by overlapping the top 15 genes identified by these algorithms.

The function enrichment analysis of m5C-modified hub genes was performed with the ClueGO [42] and CluePedia [43] plugins of Cytoscape (v.3.9.1). The analysis type is set to REACTOME pathways. Pathways with P < 0.05 were shown.

Data Availability

The data used to support the findings of this study are available in the NCBI repository. The data is accessible via NCBI GEO submission ID: GSE235610. It can be viewed at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE235610.

References

  1. Cisternas MG, Murphy L, Sacks JJ, Solomon DH, Pasta DJ, Helmick CG. Alternative methods for defining Osteoarthritis and the impact on estimating prevalence in a US Population-Based survey. Arthritis Care Res. 2016;68(5):574–80.

    Google Scholar 

  2. Guermazi A, Niu J, Hayashi D, Roemer FW, Englund M, Neogi T, Aliabadi P, McLennan CE, Felson DT. Prevalence of abnormalities in knees detected by MRI in adults without knee osteoarthritis: population based observational study (Framingham Osteoarthritis Study). BMJ (Clinical Research ed). 2012;345:e5339.

    PubMed  Google Scholar 

  3. Yang X, Yang Y, Sun BF, Chen YS, Xu JW, Lai WY, Li A, Wang X, Bhattarai DP, Xiao W, et al. 5-methylcytosine promotes mRNA export - NSUN2 as the methyltransferase and ALYREF as an m(5)C reader. Cell Res. 2017;27(5):606–25.

    CAS  PubMed  PubMed Central  Google Scholar 

  4. Yang Y, Wang L, Han X, Yang WL, Zhang M, Ma HL, Sun BF, Li A, Xia J, Chen J, et al. RNA 5-Methylcytosine facilitates the maternal-to-zygotic transition by preventing maternal mRNA decay. Mol Cell. 2019;75(6):1188–1202e1111.

    CAS  PubMed  Google Scholar 

  5. Chen H, Yang H, Zhu X, Yadav T, Ouyang J, Truesdell SS, Tan J, Wang Y, Duan M, Wei L, et al. M(5)C modification of mRNA serves a DNA damage code to promote homologous recombination. Nat Commun. 2020;11(1):2834.

    CAS  PubMed  PubMed Central  Google Scholar 

  6. Xue S, Xu H, Sun Z, Shen H, Chen S, Ouyang J, Zhou Q, Hu X, Cui H. Depletion of TRDMT1 affects 5-methylcytosine modification of mRNA and inhibits HEK293 cell proliferation and migration. Biochem Biophys Res Commun. 2019;520(1):60–6.

    CAS  PubMed  Google Scholar 

  7. Chen X, Li A, Sun BF, Yang Y, Han YN, Yuan X, Chen RX, Wei WS, Liu Y, Gao CC, et al. 5-methylcytosine promotes pathogenesis of bladder cancer through stabilizing mRNAs. Nat Cell Biol. 2019;21(8):978–90.

    CAS  PubMed  Google Scholar 

  8. Luo Y, Feng J, Xu Q, Wang W, Wang X. NSun2 Deficiency protects endothelium from inflammation via mRNA methylation of ICAM-1. Circ Res. 2016;118(6):944–56.

    CAS  PubMed  Google Scholar 

  9. Bohnsack KE, Höbartner C, Bohnsack MT. Eukaryotic 5-methylcytosine (mµC) RNA methyltransferases: mechanisms, Cellular Functions, and links to Disease. Genes 2019, 10(2).

  10. Yang Z, Zhang S, Xia T, Fan Y, Shan Y, Zhang K, Xiong J, Gu M, You B. RNA modifications meet tumors. Cancer Manage Res. 2022;14:3223–43.

    CAS  Google Scholar 

  11. Schaefer M, Hagemann S, Hanna K, Lyko F. Azacytidine inhibits RNA methylation at DNMT2 target sites in human cancer cell lines. Cancer Res. 2009;69(20):8127–32.

    CAS  PubMed  Google Scholar 

  12. Mei L, Shen C, Miao R, Wang JZ, Cao MD, Zhang YS, Shi LH, Zhao GH, Wang MH, Wu LS, et al. RNA methyltransferase NSUN2 promotes gastric cancer cell proliferation by repressing p57(Kip2) by an m(5)C-dependent manner. Cell Death Dis. 2020;11(4):270.

    CAS  PubMed  PubMed Central  Google Scholar 

  13. Zhang Q, Zheng Q, Yu X, He Y, Guo W. Overview of distinct 5-methylcytosine profiles of messenger RNA in human hepatocellular carcinoma and paired adjacent non-tumor tissues. J Translational Med. 2020;18(1):245.

    CAS  Google Scholar 

  14. Meng L, Zhang Q, Huang X. Comprehensive Analysis of 5-Methylcytosine profiles of Messenger RNA in Human High-Grade Serous Ovarian Cancer by MeRIP sequencing. Cancer Manage Res. 2021;13:6005–18.

    Google Scholar 

  15. He Y, Zhang Q, Zheng Q, Yu X, Guo W. Distinct 5-methylcytosine profiles of circular RNA in human hepatocellular carcinoma. Am J Translational Res. 2020;12(9):5719–29.

    CAS  Google Scholar 

  16. Legrand C, Tuorto F, Hartmann M, Liebers R, Jacob D, Helm M, Lyko F. Statistically robust methylation calling for whole-transcriptome bisulfite sequencing reveals distinct methylation patterns for mouse RNAs. Genome Res. 2017;27(9):1589–96.

    CAS  PubMed  PubMed Central  Google Scholar 

  17. Huang T, Chen W, Liu J, Gu N, Zhang R. Genome-wide identification of mRNA 5-methylcytosine in mammals. Nat Struct Mol Biol. 2019;26(5):380–8.

    CAS  PubMed  Google Scholar 

  18. Khoddami V, Cairns BR. Identification of direct targets and modified bases of RNA cytosine methyltransferases. Nat Biotechnol. 2013;31(5):458–64.

    CAS  PubMed  PubMed Central  Google Scholar 

  19. Hussain S, Sajini AA, Blanco S, Dietmann S, Lombard P, Sugimoto Y, Paramor M, Gleeson JG, Odom DT, Ule J, et al. NSun2-mediated cytosine-5 methylation of vault noncoding RNA determines its processing into regulatory small RNAs. Cell Rep. 2013;4(2):255–61.

    CAS  PubMed  PubMed Central  Google Scholar 

  20. Hussain S, Aleksic J, Blanco S, Dietmann S, Frye M. Characterizing 5-methylcytosine in the mammalian epitranscriptome. Genome Biol. 2013;14(11):215.

    PubMed  PubMed Central  Google Scholar 

  21. He Y, Shi Q, Zhang Y, Yuan X, Yu Z. Transcriptome-wide 5-Methylcytosine functional profiling of long non-coding RNA in Hepatocellular Carcinoma. Cancer Manage Res. 2020;12:6877–85.

    CAS  Google Scholar 

  22. Song H, Zhang J, Liu B, Xu J, Cai B, Yang H, Straube J, Yu X, Ma T. Biological roles of RNA m(5)C modification and its implications in Cancer immunotherapy. Biomark Res. 2022;10(1):15.

    PubMed  PubMed Central  Google Scholar 

  23. Squires JE, Patel HR, Nousch M, Sibbritt T, Humphreys DT, Parker BJ, Suter CM, Preiss T. Widespread occurrence of 5-methylcytosine in human coding and non-coding RNA. Nucleic Acids Res. 2012;40(11):5023–33.

    CAS  PubMed  PubMed Central  Google Scholar 

  24. Schosserer M, Minois N, Angerer TB, Amring M, Dellago H, Harreither E, Calle-Perez A, Pircher A, Gerstl MP, Pfeifenberger S, et al. Methylation of ribosomal RNA by NSUN5 is a conserved mechanism modulating organismal lifespan. Nat Commun. 2015;6:6158.

    CAS  PubMed  Google Scholar 

  25. Janin M, Ortiz-Barahona V, de Moura MC, Martínez-Cardús A, Llinàs-Arias P, Soler M, Nachmani D, Pelletier J, Schumann U, Calleja-Cervantes ME, et al. Epigenetic loss of RNA-methyltransferase NSUN5 in glioma targets ribosomes to drive a stress adaptive translational program. Acta Neuropathol. 2019;138(6):1053–74.

    CAS  PubMed  PubMed Central  Google Scholar 

  26. Zeng CY, Wang XF, Hua FZ. HIF-1α in Osteoarthritis: from pathogenesis to therapeutic implications. Front Pharmacol. 2022;13:927126.

    CAS  PubMed  PubMed Central  Google Scholar 

  27. Zhang Y, Vasheghani F, Li YH, Blati M, Simeone K, Fahmi H, Lussier B, Roughley P, Lagares D, Pelletier JP, et al. Cartilage-specific deletion of mTOR upregulates autophagy and protects mice from osteoarthritis. Ann Rheum Dis. 2015;74(7):1432–40.

    CAS  PubMed  Google Scholar 

  28. 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(4):400–9.

    CAS  PubMed  Google Scholar 

  29. Wang W, Liu Y, Hao J, Zheng S, Wen Y, Xiao X, He A, Fan Q, Zhang F, Liu R. Comparative analysis of gene expression profiles of hip articular cartilage between non-traumatic necrosis and osteoarthritis. Gene. 2016;591(1):43–7.

    CAS  PubMed  Google Scholar 

  30. Guo SM, Wang JX, Li J, Xu FY, Wei Q, Wang HM, Huang HQ, Zheng SL, Xie YJ, Zhang C. Identification of gene expression profiles and key genes in subchondral bone of osteoarthritis using weighted gene coexpression network analysis. J Cell Biochem. 2018;119(9):7687–95.

    CAS  PubMed  Google Scholar 

  31. Zheng L, Chen W, Xian G, Pan B, Ye Y, Gu M, Ma Y, Zhang Z, Sheng P. Identification of abnormally methylated-differentially expressed genes and pathways in osteoarthritis: a comprehensive bioinformatic study. Clin Rheumatol. 2021;40(8):3247–56.

    PubMed  Google Scholar 

  32. van der Slot-Verhoeven AJ, van Dura EA, Attema J, Blauw B, Degroot J, Huizinga TW, Zuurmond AM, Bank RA. The type of collagen cross-link determines the reversibility of experimental skin fibrosis. Biochim Biophys Acta. 2005;1740(1):60–7.

    PubMed  Google Scholar 

  33. Remst DF, Blom AB, Vitters EL, Bank RA, van den Berg WB, Blaney Davidson EN, van der Kraan PM. Gene expression analysis of murine and human osteoarthritis synovium reveals elevation of transforming growth factor β-responsive genes in osteoarthritis-related fibrosis. Arthritis Rheumatol. 2014;66(3):647–56.

    CAS  PubMed  Google Scholar 

  34. Troeberg L, Nagase H. Proteases involved in cartilage matrix degradation in osteoarthritis. Biochim Biophys Acta. 2012;1824(1):133–45.

    CAS  PubMed  Google Scholar 

  35. Shaw T, Nixon JS, Bottomley KM. Metalloproteinase inhibitors: new opportunities for the treatment of rheumatoid arthritis and osteoarthritis. Expert Opin Investig Drugs. 2000;9(7):1469–78.

    CAS  PubMed  Google Scholar 

  36. Rangasamy L, Geronimo BD, Ortín I, Coderch C, Zapico JM, Ramos A, de Pascual-Teresa B. Molecular Imaging Probes based on Matrix Metalloproteinase inhibitors (MMPIs). Molecules 2019, 24(16).

  37. Lu HT, Sheu MT, Lin YF, Lan J, Chin YP, Hsieh MS, Cheng CW, Chen CH. Injectable hyaluronic-acid-doxycycline hydrogel therapy in experimental rabbit osteoarthritis. BMC Vet Res. 2013;9:68.

    PubMed  PubMed Central  Google Scholar 

  38. Snijders GF, van den Ende CH, van Riel PL, van den Hoogen FH, den Broeder AA. The effects of doxycycline on reducing symptoms in knee osteoarthritis: results from a triple-blinded randomised controlled trial. Ann Rheum Dis. 2011;70(7):1191–6.

    CAS  PubMed  Google Scholar 

  39. Wang G, Chen S, Xie Z, Shen S, Xu W, Chen W, Li X, Wu Y, Li L, Liu B, et al. TGFβ attenuates cartilage extracellular matrix degradation via enhancing FBXO6-mediated MMP14 ubiquitination. Ann Rheum Dis. 2020;79(8):1111–20.

    CAS  PubMed  Google Scholar 

  40. Wang W, He Y, Wu L, Zhai LL, Chen LJ, Yao LC, Yu KH, Tang ZG. N(6) -methyladenosine RNA demethylase FTO regulates extracellular matrix-related genes and promotes pancreatic cancer cell migration and invasion. Cancer Med. 2023;12(3):3731–43.

    CAS  PubMed  Google Scholar 

  41. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

    CAS  PubMed  PubMed Central  Google Scholar 

  42. Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, Fridman WH, Pagès F, Trajanoski Z, Galon J. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinf (Oxford England). 2009;25(8):1091–3.

    CAS  Google Scholar 

  43. Bindea G, Galon J, Mlecnik B. CluePedia Cytoscape plugin: pathway insights using integrated experimental and in silico data. Bioinf (Oxford England). 2013;29(5):661–3.

    CAS  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This study was supported by funds from the Youth Foundation of the First Affiliated Hospital of Zhengzhou University (YNQN-70131).

Author information

Authors and Affiliations

Authors

Contributions

YY and XJZ designed the study; LST collected samples; LY and LXM contributed to data collection and statistical analysis.; YY and LST wrote this manuscript. XJZ contributed to the critical revision of the final manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Jianzhong Xu.

Ethics declarations

Ethical approval consent to participate

The study has been approved by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University (2022-KY-0854-002), and the patient’s written informed consent was obtained before the study began. We certify that the study was carried out in accordance with the ethical standards as laid down in the Declaration of Helsinki.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Yu, Y., Lu, S., Liu, X. et al. Identification and analysis of RNA-5-methylcytosine-related key genes in osteoarthritis. BMC Genomics 24, 539 (2023). https://doi.org/10.1186/s12864-023-09651-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-023-09651-4

Keywords