Identification of differentially expressed miRNAs and mRNAs in synovial of osteoarthritis via RNA-sequencing

Background Osteoarthritis (OA) is the most common form of arthritis and a leading cause of disability. This study attempted to investigate the key mRNAs and miRNAs related to OA. Patients and methods From April 17th, 2018 to May 17th, 2018, five patients with OA and three normal controls were enrolled in this present study. To identify the differentially expressed mRNAs (DEmRNAs) and miRNAs (DEmiRNAs) between patients with OA and normal controls, RNA-sequencing was performed. Then, DEmiRNA-target DEmRNAs analysis and functional annotation of DEmiRNA-target DEmRNAs were performed. To validate the RNA-sequencing results, quantitative real time-PCR (RT-PCR) and western blot analysis were performed as well. Results A total of 1068 DEmRNAs, 21 DEmiRNAs and 395 DEmiRNA-DEmRNA pairs were identified in synovial tissues of patients with OA. The functional annotation of DEmiRNA-target DEmRNAs revealed that Pathways in cancer and PI3K-Akt signaling pathway were significantly enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. QRT-PCR and western blot results revealed that except for TLR7, the expression level of the others was consistent with the RNA-sequencing results, generally. Conclusion The findings of this present study may provide new clues for the roles of DEmRNAs and DEmiRNAs in the pathogenesis of OA.


Introduction
Osteoarthritis (OA) is the most frequent musculoskeletal disease and leads to functional decline and loss in quality of life [1]. Currently, there is no effective treatment to prevent the initiation and progression of the disease while the severity of OA often worsens with age [2]. Reportedly, it estimates that among population over 60 years old, there was 9.6% of men and 18% of women suffer from symptomatic OA in the worldwide [3]. Although the disease is with a trait of destruction of articular cartilage, pathological changes in subchondral bone and associated synovitis, the pathogenesis is poor understanding which is thought to be multifactorial [4]. Inability to diagnose early and poorly understanding of the pathophysiology makes early diagnosis to be a key factor in the prevention and management of disease [5]. Recently, molecular biology has been reported to play a vital role in explaining its disease pathophysiology, and gene regulation has been demonstrated to be involved in driving an imbalance between the expression of catabolic and anabolic factors, leading eventually to osteoarthritic cartilage degeneration [6]. There is increasing attention on the influence of dysregulation at a molecular level on the pathogenic process.
MicroRNAs (miRNAs) have been detected widely in eukaryotes to play roles in regulating growth, development, differentiation, and metabolism in model organisms which are short (approximately 22 nucleotides), non-coding, RNA regulators of gene expression [7]. Little was known of the function of miRNAs at the stage when it is first identified in the early 1990's by Lee et al. (in Caenorhabditis elegans) [8]. Due to the link between alterations in miRNA expression levels and multifarious disease processes have been linked, aberrant expression in pathological states has drawn public attention [9]. Studies based on highly specific patterns of miRNA expression correlate with development and several diseases have revealed the potential for therapeutic manipulation of miRNAs [10,11]. Recently, a large number of studies have been done to explore miRNAs and genes associated with OA [5,6,[12][13][14]. Even so, studies on biomarkers for OA are still urgent to be performed.
In this study, differentially expressed miRNAs (DEmiRNA) and mRNAs (DEmRNAs) in synovial tissues of patients with OA were identified by RNAsequencing. DEmiRNA-target DEmRNAs analysis and functional annotation of DEmiRNA-target DEmRNAs were performed. Quantitative real time-PCR (RT-PCR) and western blot analysis were performed to validate the RNA-sequencing results. In doing so, we hope this study could represent a new avenue to understand the pathogenesis and develop potential biomarkers for OA.

Patients and samples
Five patients with OA and three normal controls were recruited in this study from April 17th, 2018 to May 17th, 2018. According to the criteria of the American College of Rheumatology, OA were diagnosed [15]. Healthy employees with no symptoms or signs of OA, or any other type of arthritis, or any painful condition of the joints, were included as controls. Beside, participants with no personal or family history of OA were selected as control subjects. The participants with history of joint diseases, including inflammatory arthritis (rheumatoid arthritis or any other autoimmune disease), post-traumatic or post-septic arthritis, poliomyelitis, skeletal dysplasia, were excluded. Radiographic evaluation of all participants were performed. Table 1 displayed the detailed information of all these participants. The written informed consent for use of their samples from every participant were provided in the present study. Ethical approval for this study was granted by the ethics committee of People's Hospital of Deyang City (2017-043). Synovial tissues of every participant were obtained.

RNA isolation and sequencing
Following the manufacturer's protocol, we used TRIzol reagent (Invitrogen, Carlsbad, CA, USA) to isolate total RNA from samples. The concentration and purity of RNA was determined with Nanodrop ND-2000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA), and the integrity of RNA was confirmed via a 2% agarose gel. With an Agilent 2100 Bioanalyzer (Agilent, Palo Alto, CA, USA), the RNA integrity number (RIN) value was obtained. With QiaQuick PCR Purification Kit, the mRNA library was constructed. The 18-30 nt RNA was obtained from the total RNA. By using Tru-seqTM Small RNA Sample Prep Kit, adapter ligation and reverse transcription polymerase chain reaction (PCR) were performed to obtain the cDNA. Sequencing was performed based on HiSeq x-ten platform (Illumina) and SE50, BGIseq, respectively.

Functional annotation of DEmRNAs between patients with OA and normal controls
To further research the biological function of DEmR-NAs, Gene Ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of the DEmRNAs between patients with OA and normal controls were performed by DA-VID 6.8 (https://david.ncifcrf.gov/). A value of p < 0.05 was defined as the criteria of statistical significance.

Identification of DEmiRNAs in patients with OA compared with normal controls
With Bowtie (bowtie-bio.sourceforge.net), the alignment between the cleaned miRNA sequencing reads was aligned to the human genome (GRCh38.p7 assembly) based on Genome human UCSC reference annotation. With miRDeep2 (https://www.mdc-berlin.de/8551903/ en/), the transcription abundance of miRNAs was determined. Based on the read count of each sample, the DEmiRNAs in OA compared to normal controls were calculated with an R-bioconductor package, DESeq2 (http://www.bioconductor.org/packages/release/bioc/ html/DESeq2.html). For DEmiRNAs in OA compared to normal controls, the threshold was defined as base Mean > 100, p < 0.05 and |log 2 FC| > 2. With R package "pheatmap", hierarchical clustering analysis of DEmR-NAs was conducted.

DEmiRNA-target DEmRNAs analysis
Given miRNAs tend to decrease the expression of their target mRNAs, we selected target genes from DEmRNAs that expressed inversely with that of miRNA for further research. Firstly, the putative targeted DEmRNAs of DEmiRNAs were predicted by six bioinformatic algorithms (RNA22, miRanda, miRDB, miRWalk, PICTAR2 and Targetscan). Then, with miRWalk, the confirmed targeted DEmRNAs of DEmiRNAs were obtained. Thirdly, the confirmed DEmiRNA-DEmRNA pairs were derived from miRWalk and the DEmiRNA-DEmRNA pairs recorded by ≥4 algorithms. Based on the obtained DEmiRNA-DEmRNA pairs, DEmiRNA-DEmRNA interaction networks between OA and normal controls were constructed by using Cytoscape software (http://www. cytoscape.org/).

Functional annotation of DEmiRNA targets
To further research the biological function of the target DEmRNAs of DEmiRNAs, GO analysis was performed t by using Gorilla (http://cbl-gorilla.cs.technion.ac.il/) with a p-value < 0.01, and KEGG pathway analysis was performed by using webgestalt (http://www.webgestalt. org/) with a p-value < 0.05.

Quantitative real time-PCR (RT-PCR) and western blot analysis
Eight synovial tissues were obtained from four patients with OA and four normal controls. The written informed consent for use of their samples from every participant were provided in the present study. Ethical approval for this study was granted by the ethics committee of People's Hospital of Deyang City (2017-043).
Total RNA was isolated with the Trizol reagent (Invitrogen, USA). The qRT-PCR reactions were performed based on SuperReal PreMix Plus (Invitrogen, USA) in ABI 7300 Real-time PCR Detection System. With 2 -ΔΔCt method, relative gene expression was determined. The human GAPDH and ACTB were used as endogenous controls for mRNA expression, and the human hsa-U6 was used as endogenous controls for miRNA expression in analysis.
Synovial tissues were lysed on ice with RIPA lysis buffer, and then the supernatants were collected by centrifugation at 12,000 rpm at 4°C for 30 min. The protein concentrations were detected with BCA protein assay. The protein extracts were separated by 10% SDS-PAGE, transferred onto PVDF membranes, and probed using primary and then secondary antibodies. The primary antibodies were as follows: rabbit anti-GAPDH, TIMP3 Antibody, CTSS Antibody and TLR7 Antibody. The blots were visualized with ECL reagent.

DEmRNAs and DEmiRNAs between patients with OA and normal controls
Expression profiling of mRNA and miRNA extracted from five patients with OA and three normal controls identified 1068 DEmRNAs (516 up-regulated and 552 down-regulated DEmRNAs) and 21 DEmiRNAs (6 upregulated and 15 down-regulated DEmiRNAs). The top 10 up-and down-regulated DEmRNAs and all DEmiR-NAs were showed in Table 2 and Table 3, respectively. Hierarchical clustering analysis of top 50 up-and downregulated DEmRNAs and DEmiRNAs was displayed in Fig. 1 and Fig. 2, respectively. The raw-data have been uploaded to Gene Expression Omnibus (GEO) (GSE143514, https://www.ncbi.nlm.nih.gov/geo/query/ acc.cgi?acc=GSE143514).

Functional annotation of DEmiRNA targets
Base on GO enrichment analysis, RNA processing (p = 5.80E-04), response to activity (p = 9.04E-04), nucleus (p = 3.56E-04) and alpha-tubulin binding (p = 5.80E-04) were significantly enriched GO terms in OA (Fig. 5a-c). According to the KEGG pathway enrichment analysis, the DEmiRNA-target DEmRNAs were significantly enriched in Pathways in cancer (p = 3.36E-02) and PI3K-Akt signaling pathway (p = 4.67E-02) (Fig. 5d-e) [16].  The results of qRT-PCR and western blot Three DEmRNAs (including TLR7, CTSS and TIMP3) and two DEmiRNAs (including hsa-miR-17-5p and hsa-miR-20b-5p) were selected to use for qRT-PCR validation. Based on the results of RNA sequencing, TLR7, CTSS and TIMP3 were up-regulated while hsa-miR-17-5p and hsa-miR-20b-5p were downregulated in OA. Except for TLR7, expression of the others in the qRT-PCR results was consistent with that in this present study, generally (Fig. 6a). Western blot results revealed that the protein level of CTSS and TIMP3 were up-regulated while the protein level of TLR7 was down-regulated in OA (Fig. 6b).

Discussion
OA is a progressive disease with a long silent period, which shows signs of cartilage degradation, mild-tomoderate synovial inflammation and altered bone structure, resulting in severe destruction and impaired function of the affected joints [17]. This study performed RNA-sequencing to attempt to obtain the key miRNAs and mRNAs associated with OA. Here, we discussed two DEmiRNAs (including hsa-miR-17-5p and hsa-miR-20b-5p), which were the top two DEmiRNAs that covered most DEmRNAs, and three DEmRNAs (including TLR7, CTSS and TIMP3), which were target genes of the two DEmiRNAs mentioned above. Toll like receptor 7 (TLR7), a member of the Tolllike receptor (TLR) family, encoded by TLR7, plays a fundamental role in pathogen recognition and activation of innate immunity [18]. From Drosophila to humans, TLRs are highly conserved and share structural and functional similarities [18]. In humans, the TLR family, a set of 10 type I transmembrane receptors, have specificity for different types of pathogenassociated molecular patterns (PAMPs) [19]. PAMPs are expressed on infectious agents, and mediate the production of cytokines necessary for the development of effective immunity [19]. In humans, ten TLRs have been identified, and TLR-1, TLR-2, TLR-4, TLR-5, TLR-6, and TLR-10 are expressed on the cell surface while TLR-3, TLR-7, TLR-8, and TLR-9 are expressed in endosomes [20]. TLR-7 recognizes single-stranded RNA (ssRNA) in endosomes and is activated by the synthetic antiviral compound imiquimod [21]. TLR activation has been demonstrated to be involved in the    [24]. Due to TLR7 was found to be up-regulated in OA in our study, we speculated that TLR7 may involve in OA. However, down-regulated TLR7 was found in the results of qRT-PCR and western blot, its exact role in OA need to be confirmed in the further study.
CTSS (cathepsin S), encoded by CTSS, is a member of the peptidase C1 family [25]. In addition, CTSS is a lysosomal cysteine proteinase that may participate in the degradation of antigenic proteins to peptides for presentation on MHC class II molecules [25]. Appleton et al. identified increased expression of CTSS in the OA model [26]. The study of Lambert et al. showed that CTSS was significantly up-regulated in inflamed synovial biopsy samples which was confirmed at the protein level by using immunohistochemistry, and it was also significantly up-regulated in cartilage catabolism pathway based on their results [27]. In agreement with previous studies, we observed that CTSS was up-regulated in patients with OA which may imply the importance of CTSS in OA.
TIMP3 (full name, TIMP metallopeptidase inhibitor 3), is a member of TIMP family which is inhibitors of the matrix metalloproteinases, a group of peptidases implicated in degradation of the extracellular matrix (ECM) [28]. In bone, TIMP-3 plays the role as a local cytokine, regulating bone metabolism through suppressing osteoblast differentiation and inducing osteoblast apoptosis [29,30]. TIMP3 was reported to be associated with RA and arthritis based on a search of the U. S. National Library of Medicine database (MED-LINE) (http://www.ncbi.nlm.nih.gov/IEB/Research/Ace mbly/av.cgi) [31]. A case-control study in a Chinese Han population linked TIMP3 polymorphism with severe knee OA [32]. In this study, TIMP3 was detected to be dysregulated in patients with OA. Hence, we hypothesized that TIMP3 was involved with OA.
Hsa-miR-17-92 polycistron encodes a cluster of seven miRNAs derived from the c-myc regulated c13orf25 locus at chromosome 13q31.3, including miR-17-5p [33]. MiR-17-5p was reported to have a role as a tumor suppressor in breast cancer cells and as a key regulator of the G1/S phase cell cycle transition [7,34]. MiR-17-5p may act as an oncogene or a tumor suppressor in different cellular contexts [7]. Although no report linked hsa-miR-17-5p with OA, it was a down-regulated miRNA that covered most DEmRNAs in this study. In addition, TLR7, CTSS and TIMP3 were detected to be target genes of hsa-miR-17-5p. Given the results of our analysis, we hypothesized that hsa-miR-17-5p may participate in OA by regulating TLR7, CTSS and TIMP3. MiR-20b-5p is transcribed from the miR-106a~363 clusters which is reported to be involved in several process [35]. Ma et al. suggested that miR-20b-5p may play a vital role in multiple sclerosis (MS) pathogenesis [36]. In the study of Luo et al., miR-20b-5p was found to play role in promoting myoblast differentiation and repressing myoblast proliferation [35]. Interestingly, hsa-miR-20b-5p was down-regulated in patients with OA. In addition, TLR7 and CTSS were targeted by hsa-miR-20b-5p which indicated that miR-20b-5p may implicate in OA. To our best knowledge, we are the first to report the role of hsa-miR-20b-5p in OA.

Conclusions
In conclusion, we identified 1068 DEmRNAs (516 upregulated and 552 down-regulated DEmRNAs), 21 DEmiRNAs (6 up-regulated and 15 down-regulated DEmiRNAs) and 395 DEmiRNA-DEmRNA pairs in synovial tissues of patients with OA, and emphasized MF, molecular function. d-e KEGG pathways. d. PI3K-Akt signaling pathway. e. Pathways in cancer. The red and green rectangles represented the particles which regulated by the up-and down-regulated target DEmRNAs of DEmiRNAs, respectively. Appropriate copyright permission to use the signalling pathways was obtained the importance of several mRNAs and miRNAs which may implicate in OA. These findings may provide new avenue to understand the mechanism of OA. A limitation of present study is the small sample size for RNA sequencing. The exact function of these mRNAs and miRNAs in OA need to be determined with further large sample research.
Additional file 1 : Figure S1 Flow chart of the analyses. Ethics approval and consent to participate This study was approved by the ethical committee of People's Hospital of Deyang City. Written informed consent for use of their samples from every participant were provided in the present study.

Consent for publication
Not applicable.