MicroRNA-1202 plays a vital role in osteoarthritis via KCNQ1OT1 has-miR-1202-ETS1 regulatory pathway

Background This study aimed to explore the molecular mechanism of osteoarthritis (OA) and provide information about new genes as potential targets for OA treatment. Methods Gene expression profile of GSE105027, including 12 OA serum samples (OA group) and 12 healthy serum samples (ctrl group), was downloaded. The differentially expressed miRNAs (DEMs) as well as miRNA-mRNAs interactions were investigated, followed by function and pathway investigation. Then the protein-protein interaction (PPI) network was performed. Furthermore, the long non-coding RNA (lncRNA)-miRNA-mRNA interactions (competing endogenous RNAs, ceRNAs) were investigated. Results A total of 17 downregulated miRNAs were revealed between OA and ctrl groups. These DEMs such as has-miR-1202 were mainly enriched in GO functions like histone acetyltransferase binding and KEGG pathways like cellular senescence. The integrated PPI network analysis showed that has-miR-1202, has-miR-33b-3p, has-miR-940, has-miR-4284, and has-miR-4281 were 5 downregulated miRNAs in this network. Furthermore, the lncRNA-miRNA-mRNA interactions such as KCNQ1OT1-has-miR-1202-ETS1 were revealed in the present ceRNA network. Conclusion Key DEMs such as miR-33b-3p, miR-940, and miR-1202 may be involved in OA. miR-1202 may regulate OA development via histone acetyltransferase pathway binding function and cellular senescence pathway. Furthermore, KCNQ1OT1-has-miR-1202-ETS1 might be vital for the process of OA.


Background
Osteoarthritis (OA) is the wear-and-tear type of arthritis that causes joint pain and stiffness [1]. Worldwide, there are approximately 0.25 billion people suffering with OA [2]. The epidemiology of OA is very complex due to various factors such as genetic and biomechanical components [3], and the molecular mechanisms underlying OA are still not completely understood.
MicroRNAs (miRNAs) are commonly used to identify novel OA genes and their related inflammatory network [4]. A previous study shows that the upregulation of miR-1 controls the development of OA via targeting FZD7 of Wnt/β-catenin signaling [5]. Zhang et al. indicated that miR-320 might target matrix metalloproteinase-13 and further affected the process of OA [6]. Despite miRNAs, long non-coding RNAs (lncRNAs) are correlated with OA and are regulated by OA-associated factors [7,8]. The previous study shows that lncRNAs take part in many kinds of pathological processes in OA such as extracellular matrix (ECM) [9]. A previous bioinformatics analysis indicates that lncRNAs (such as uc.343) and predicted target homeobox gene C8 (HOXC8) are promising therapeutic targets for OA [10]. Actually, the lncRNAs-miRNAs-mRNAs (competing endogenous RNAs, ceR-NAs) regulatory network is an important tripartite axis in the regulation of the disease process [11,12]. Although the identification of molecular factors contributes to the therapeutic interventions in OA [13], the key miRNAs and the possible regulation mechanism of ceRNA in OA are still unclear.
The bioinformatics analysis of gene expression profiles provides new opportunities to uncover the potential OA-related genes such as cathepsin H (CTSH) and cathepsin S (CTSS) [14][15][16][17]. In order to explore the potential differentially expressed miRNAs (DEMs) as a novel breakthrough in the clinical treatment of OA, Ntoumou et al. identified a circulating miRNA signature for OA patients using DEMs analysis and related pathway investigation [18]. They found that miRNAs such as miR-140-3p regulated the metabolic processes of OA patients. However, the integrated regulatory mechanism involving DEMs, mRNAs, and lncRNAs on OA progression is still unclear. Based on the gene expression profile provided by Ntoumou et al., the present bioinformatics analysis revealed the DEMs between OA samples and control samples. Then, the miRNAs-mRNAs relations were explored, followed by the miRNAs-mRNAs regulatory network investigation. The protein-protein interaction (PPI) network was further constructed according to mRNAs in the miRNAs-mRNAs network. Finally, the ceRNA network was constructed and analyzed. The present study aimed to reveal the potential mechanism of OA and provide information about new genes as potential targets for OA treatment.

Data resource and preprocessing
Gene expression profile data in the dataset GSE105027 [18] were downloaded from Gene Expression Omnibus (GEO) database. This dataset was produced on GPL21575 Agilent-070156 Human_miRNA_V21.0_ Microarray 046064 platform. A total of 12 serum samples of OA patients (OA group) and 12 serum samples of healthy individuals (ctrl group) were included in the dataset.
The normalization for the downloaded original data was performed using the RAM75 [19,20] method based on Affy package (version 1.50.0) in Rstudio software (version 1.1.414) [21]. The normalization process in this study included background correction, normalization, and expression quantification. If different probes targeted to the same miRNA (miRNA symbol), the mean value of different probes was considered as the final expression value of this miRNA.

The investigation for DEMs
The P value and fold change (FC) of DEMs between the OA and ctrl group were calculated by Linear Models for Microarray Data (limma) package [22] in R (version 3.32.5) software. Then the P value < 0.05 and |log 2 FC| > 1 were selected as the thresholds for the identification of DEMs. The bidirectional hierarchical clustering for DEMs was then performed by pheatmap software (version 1.0.8).

The miRNA-target gene investigation and network construction
The predicted target genes of miRNAs were investigated using miRWalk 2.0 software (parameters, miRBase database; species, human) [23,24]. The miRWalk [25], RNA22 [26], miRanda [27], and Targetscan [28] were selected as default databases for interaction scanning. The genes that appeared in all four databases were selected as target genes. Then, we downloaded the target genes of miRNAs which had been already verified in the validated target module of miRWalk 2.0 software. Based on overlapping results of predicted miRNA-target genes and verified miRNA-target genes, the miRNA-mRNAs interactions were finally explored. Furthermore, the miRNA-mRNAs interaction relations were visualized by Cytoscape (version 3.2.1, http://apps.cytoscape.org/) [29].

Enrichment analysis of the DEGs
The clusterProfiler software [30] is an online tool that provides enrichment analyses including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). Based on the clusterProfiler software, the GO functional annotation and KEGG pathway enrichment analysis were performed on DEGs. P value (the significance threshold of the hypergeometric test) < 0.05 was chosen as the cutoff criterion for present enrichment analysis.

Construction of PPI network
The PPI research could reveal the interactions of proteins at the molecular level [31]. Search Tool for the Retrieval of Interacting Genes/Proteins (STRING, version 10.0) is a biological database of known and predicted PPI [32]. Using STRING database, PPI was constructed by the target genes that were regulated by miRNAs (required confidence (combined score) > 0.4). Then, the integrated PPI network was further constructed based on the PPI network constructed by target genes regulated by miRNA and miRNA-target interactions, followed by visualization using Cytoscape software. The network topology analyses, including degree centrality (DC) [33], betweenness centrality (BC) [34], and closeness centrality (CC) [35], were performed by the CytoNCA [36] (version 2.1.6; parameter, without weight) to reveal the hub protein [19].

Construction of ceRNA regulatory network
The miRNA-associated lncRNAs were explored based on the lnCeDB database (http://www.lncrnablog.com/) [37]. The lncRNA-miRNA regulatory relations were constructed using Cytoscape software. Then, based on the lncRNA-miRNA relations and miRNA-DEMs interactions, the lncRNA-miRNA-mRNA (ceRNA) network was further explored. Finally, the ceRNA network was visualized by Cytoscape software.

DEMs investigation between the OA group and ctrl group
With P value < 0.05 and |log 2 FC| > 1, a total of 17 downregulated DEMs were obtained between the OA and ctrl groups. However, the result of the upregulated DEMs was negative. The bidirectional hierarchical clustering for DEMs is shown in Fig. 1. The results showed that the expression value of OA samples and ctrl samples could basically separate the two groups.

ceRNA regulatory network investigation
The lnCeDB database was used to explore the lncRNAs related with 5 miRNAs mentioned above. The results showed that a total of 172 lncRNAs and 215 interactions were investigated in miRNAs-lncRNAs regulatory relations. With the binding-site number > 2, a lncRNA-miRNA regulatory network (including 3 miRNAs and 30 lncRNAs) was constructed using Cytoscape software. Then, based on the lncRNA-miRNA relations and miRNA-DEMs interactions, the lncRNA-miRNA-mRNA interactions (ceRNA) including KCNQ1 overlapping transcript 1 (KCNQ1OT1)-has-miR-1202-ETS1 was further explored. Finally, the ceRNA network was visualized by Cytoscape software (Fig. 4).

Discussion
OA is the most common form of joint disease [38]. miR-NAs are expressed in a different fashion in osteoarthritic   Fig. 3 The integrated protein-protein interaction network in present study. The red circle represents the mRNA, the yellow triangle represents the downregulated miRNA, the line without arrow represents the protein-protein interaction relation, and the line with an arrow represents miRNA-mRNA relationship. The dot size represents the correlation degree of node in PPI network; the larger the number of dots; the larger the dots, the greater the genetic correlation degree relative to nonosteoarthritic cartilage and may be useful for diagnosis or management of OA [33,34]. In the present study, a total of 17 downregulated miRNAs were revealed between the OA group and ctrl group. These DEMs such as has-miR-1202 were mainly enriched in GO functions like histone acetyltransferase binding and KEGG pathways like cellular senescence. The integrated PPI network analysis showed that has-miR-1202, has-miR-33b-3p, has-miR-940, has-miR-4284, and has-miR-4281 were 5 downregulated miRNAs in this network. Furthermore, the lncRNA-miRNA-mRNA interactions such as KCNQ1OT1-has-miR-1202-ETS1 were revealed in the present ceRNA network. miRNAs are increasingly implicated in the pathogenesis of complex diseases such as cancer and OA [35,39]. Based on bioinformatics and proteomic analysis, Dimitrios et al. indicated that the collaborative networks of miRNAs and proteins were associated with OA pathogenesis [40]. In previous studies, serum miRNA array analysis showed that miR-33b-3p may be used as a potential biomarker for OA [41]. Meanwhile, miR-940 is involved in the inflammatory response of chondrocytes in OA [42]. In this study, miR-33b-3p and miR-940 were downregulated in the OA group, suggesting the potential role of them in OA. As a member of the miRNA family, miR-1202 suppresses proliferation and induces endoplasmic reticulum stress in tumor cells, which suggest a potential therapy of miR-1202 in clinical treatment [43]. Despite that, miR-1202 is a primate-specific and brainenriched miRNA [44]. However, the function of miR-1202 in OA is still unclear. In this study, miR-1202 was an outstanding downregulated miRNA between the OA group and ctrl group, suggesting that miR-1202 might take part in the progression of OA.
In addition, complex interactions existed between miR-NAs and their multiple target genes, which may be important in gene regulation and the control of homeostatic pathways in OA [33,34]. As a target gene of miR-1202, Fig. 4 The ceRNA network in this study. The orange dot represents mRNA; the blue diamond represents lncRNA; and the yellow triangle represents the downregulated miRNA. The dot size represents the correlation degree of node in PPI network; the larger the dots, the greater the genetic correlation degree ETS1 is a member of the ETS family of transcription factors [45]. ETS1 has been proved to associate with regulation of immune cell function and with an aggressive behavior in tumors [46]. Xu et al. indicated that miR-221-3p modulated ETS1 expression in synovial fibroblasts of OA patients [47], suggesting the potential role of ETS1 in OA. Histone acetylation has a close relation with some biological effects such as gene transcription and is catalyzed by histone acetyltransferases [48]. A previous study shows that ETS1 synergistically activates the gene expression of guanylyl cyclase/ natriuretic peptide receptor A via histone acetyltransferase pathway [49]. However, the detailed mechanism of ETS1 and histone acetyltransferase pathway in OA is unclear. Moreover, ETS1 can activate the p16INK4a promoter through an ETS-binding site via cellular senescence function [50]. In this study, miR-1202 and ETS1 were enriched in histone acetyltransferase binding function and cellular senescence pathway. Thus, based on the present study, we speculated that the downregulated miR-1202 might target EST1 and further affect the OA progression via histone acetyltransferase pathway binding and cellular senescence pathway.
LncRNAs are proved to participate in a variety of biological processes as regulatory molecules [51]. The expression of lncRNAs has recently been reported in tumorigenesis and plays a pivotal role in regulating behavior cell cycle [52]. A previous study shows that lncRNA functions as a ceRNA to promote cartilage degradation in human OA [8]. KCNQ1OT1 is a nuclear transcript found in close proximity to the nucleolus in certain cell types [53]. It involves in certain imprinted gene network that may play a role in various diseases including Beckwith-Wiedemann syndrome and coronary artery lesion [54,55]. However, the effect of KCNQ1OT1 on OA is still unclear. Our ceRNA network analysis showed that KCNQ1OT1-has-miR-1202-ETS1 was one of the most outstanding ceRNA. Thus, we speculated that the lncRNA KCNQ1OT1 might play an important role in OA progression via sponging has-miR-1202-ETS1 interaction.
However, there were some limitations in this study such as small sample size and lack of verification test, which will introduce false-positive results and consequently influence the reliability of our findings. Thus, a large sample size with a wide verification analysis is needed in further investigation.
In conclusion, the present study showed that downregulated DEMs such as miR-33b-3p, miR-940, and miR-1202 may be involved in OA. miR-1202 might target EST1 to regulate the activation of histone acetyltransferase pathway binding and cellular senescence pathway and thus further affect OA development. Furthermore, KCNQ1OT1-has-miR-1202-ETS1 might play an important role in the process of OA.