Systematic Analysis of Transcriptomic Profile of Chondrocytes in Osteoarthritic Knee Using Next-Generation Sequencing and Bioinformatics

The phenotypic change of chondrocytes and the interplay between cartilage and subchondral bone in osteoarthritis (OA) has received much attention. Structural changes with nerve ingrowth and vascular penetration within OA cartilage may contribute to arthritic joint pain. The aim of this study was to identify differentially expressed genes and potential miRNA regulations in OA knee chondrocytes through next-generation sequencing and bioinformatics analysis. Results suggested the involvement of SMAD family member 3 (SMAD3) and Wnt family member 5A (WNT5A) in the growth of blood vessels and cell aggregation, representing features of cartilage damage in OA. Additionally, 26 dysregulated genes with potential miRNA–mRNA interactions were identified in OA knee chondrocytes. Myristoylated alanine rich protein kinase C substrate (MARCKS), epiregulin (EREG), leucine rich repeat containing 15 (LRRC15), and phosphodiesterase 3A (PDE3A) expression patterns were similar among related OA cartilage, subchondral bone and synovial tissue arrays in Gene Expression Omnibus database. The Ingenuity Pathway Analysis identified MARCKS to be associated with the outgrowth of neurite, and novel miRNA regulations were proposed to play critical roles in the pathogenesis of the altered OA knee joint microenvironment. The current findings suggest new perspectives in studying novel genes potentially contributing to arthritic joint pain in knee OA, which may assist in finding new targets for OA treatment.


Introduction
Osteoarthritis (OA) is one of the common musculoskeletal disorders with increased incidence in the elderly [1]. The diagnosis of OA is based on clinical symptoms and radiographic findings,

Ingenuity Pathway Analysis (IPA)
The Ingenuity Pathway Analysis (IPA) software (Ingenuity systems, Redwood City, CA, USA) provides bioinformatics analysis of gene arrays; RNA and small RNA sequencing, proteomics and other biological experiments; categorizes molecules into different biological networks; and identifies related canonical pathways, upstream regulators, diseases and functions involved in a subset of molecules. The IPA core analysis provides predicted canonical pathways based on the literature and de novo network discovery based on expression changes of uploaded experimental results [17]. Additionally, construction of causal networks based on gene interactions curated from the literature is also available, which provides mechanistic hypotheses based on expression changes of the uploaded dataset. The statistical significance of the constructed regulatory networks was determined by Fisher's exact test p-value (enrichment score) that measured the overlap of observed and predicted gene sets, and Z score that measured the match of observed and predicted regulatory patterns with directional effect [25]. In the current study, differentially expressed genes with fold changes in expression between normal and OA knee chondrocytes were uploaded into IPA for core analysis to obtain pathway and network results. For "create expression analysis" setting, we took into consideration all direct and indirect relationships identified in all data sources with experimentally observed or moderate to highly predicted confidence, without specifying tissue types or species.

Database for Annotation, Visualization and Integrated Discovery (DAVID) Bioinformatics Resources
The DAVID bioinformatics resources is an open-source bioinformatics database that provides an integrated data-mining environment and categorizes a list of genes into different biological functions using different tools such as functional annotation and gene functional classification. Users are able to have an overall concept of the biological themes that are enriched in the list of genes [18]. Differentially expressed genes of interest were input into DAVID database for functional annotation analysis, and the Expression Analysis Systematic Explorer (EASE) score was set at a default cutoff value of 0.1, representing the modified Fisher's exact p-value in the DAVID database.

Gene Expression Omnibus (GEO) Database
The GEO database provides high-throughput genomic data sets, and researchers can have access to the analysis of genes of interest and their expression values through web-based tools such as GEO2R, and perform further between-group comparisons [19]. The relevant datasets related to joint tissues of OA patients (GSE114007 for cartilages, GSE51588 for subchondral bones, GSE55235 and GSE55457 for synovial tissues) were used in this study to examine whether the expression patterns of genes of interest were in line with that observed in our expression profile of OA knee chondrocytes. Detailed information of the selected datasets was provided in .

miRmap Database
MiRmap is an open-source software library developed by Vejnar et al. that is available online. The library provides miRNA target prediction by ranking the repression strength of potential targets. The repression strength is predicted by a comprehensive approach, combining thermodynamic, evolutionary, probabilistic and sequence-based features, and is given a miRmap score for each potential target. The higher score indicates higher repression strength. In our current study, the 46 differentially expressed miRNAs were sequentially inputted into the library to obtain potential target genes. Those potential targets with miRmap scores higher than 99.0 were selected for further analysis [20,21].

Statistical Analysis
To determine differentially expressed mRNAs between normal and OA knee chondrocytes, statistical analysis was performed using Cuffdiff (Cufflinks version 2.2.1), with p-value calculation for non-grouped samples using blind mode, in which all samples were treated as replicates of a single global condition and used to blind one model for a statistic test [26]. Based on the selection of competitive null hypothesis and the use of over-representation analysis as the pathway analysis method for our pathway annotation analysis [16,27], we firstly selected from our high-throughput differential expression data those mRNAs with expression values of FPKM >0.3, a fold-change >2.0 and a p-value <0.05 between normal and OA knee chondrocytes. Further, multiple testing correction for p-value adjustments using Bonferroni adjustment method, Hochberg's method and false discovery rate (FDR) controlling method [28,29] was performed using Statistical Analysis System (SAS) PROC MULTTEST (SAS version 9.4, SAS Institute, Cary, NC, USA) for the pre-selected differential expression gene list. An adjusted p-value of <0.05 indicated statistical significance. All pre-selected differential expression genes remained statistically significant after Hochberg's method and FDR controlling method corrections. In addition, expression values of target genes in OA and non-OA patients identified from selected GEO arrays were analyzed for between-group differences, using IBM SPSS Statistics for Windows, version 19 (IBM Corp., Armonk, NY, USA). For the relatively small numbers of cases in each group, the non-parametric method with the Mann-Whitney U test was used for between-group comparison. A p-value < 0.05 was considered as a statistically significant difference.

The Sequencing Quality and Mapped Reads for RNA and Small RNA Sequencing of Chondrocytes
The results of sequencing quality of RNA sequencing for normal and OA knee chondrocytes, using FastQC, are shown in Figure S1, reporting high quality scores in per-base sequence quality and per-sequence quality scores. Similarly, the quality scores of small RNA sequencing for normal and OA knee chondrocytes indicated good sequencing quality ( Figure S2). The mapped reads for both RNA and small RNA sequencing are listed in Table S2. The raw sequencing data were uploaded to the GEO repository with the GEO accession number of GSE110606. The data will be held private before publication.

The Differentially Expressed Genes in Osteoarthritic Knee Chondrocytes were Associated with Osteoarthritis Pathway, Cell Adhesion and Extracellular Matrix Organization
The expression values of normal and OA knee chondrocytes were obtained from the NGS results of RNA-sequencing and were further analyzed. Screening for gene expressions with FPKM >0.3, the difference in FPKM performance between normal and OA knee chondrocytes was displayed as a density plot, as indicated in Figure 1A. A total of 495 differentially expressed protein-coding genes with >2.0-fold-change and >0.3 FPKM were identified, where 263 genes were up-regulated and 232 genes were down-regulated in OA knee chondrocytes, compared to normal chondrocytes. The detailed information of the 495 differentially expressed genes is provided in Table S3. To determine the biological functions and pathways possibly involved in these differentially expressed genes in chondrocytes, we use the IPA software and the functional annotation tool in the DAVID database for further analysis. The results of IPA analysis identified the osteoarthritis pathway as the most highly enriched canonical pathway among these dysregulated genes, with 23 differentially expressed genes involved in this pathway (p = 4.72 × 10 −10 , Z score = 1.342), and tumor necrosis factor (TNF) as the most highly predicted upstream regulator (p = 2.20 × 10 −22 , Z score = −0.159). The top five canonical pathways and potential upstream regulators are listed in Table 1. the GEO repository with the GEO accession number of GSE110606. The data will be held private before publication.

The Differentially Expressed Genes in Osteoarthritic Knee Chondrocytes were Associated with Osteoarthritis Pathway, Cell Adhesion and Extracellular Matrix Organization
The expression values of normal and OA knee chondrocytes were obtained from the NGS results of RNA-sequencing and were further analyzed. Screening for gene expressions with FPKM >0.3, the difference in FPKM performance between normal and OA knee chondrocytes was displayed as a density plot, as indicated in Figure 1A. A total of 495 differentially expressed protein-coding genes with > 2.0-fold-change and > 0.3 FPKM were identified, where 263 genes were up-regulated and 232 genes were down-regulated in OA knee chondrocytes, compared to normal chondrocytes. The detailed information of the 495 differentially expressed genes is provided in Table S3. To determine the biological functions and pathways possibly involved in these differentially expressed genes in chondrocytes, we use the IPA software and the functional annotation tool in the DAVID database for further analysis. The results of IPA analysis identified the osteoarthritis pathway as the most highly enriched canonical pathway among these dysregulated genes, with 23 differentially expressed genes involved in this pathway (p = 4.72 × 10 −10 , Z score = 1.342), and tumor necrosis factor (TNF) as the most highly predicted upstream regulator (p = 2.20 × 10 −22 , Z score = −0.159). The top five canonical pathways and potential upstream regulators are listed in Table 1.  3, the difference in FPKM performance between normal (HC) and OA (HC-OA) knee chondrocytes is shown in a density plot. The x-axis represents logarithm to the base 10 of FPKM, and the y-axis indicates read density. (B) The 495 differentially expressed protein-coding genes in OA knee chondrocytes were input into the DAVID database to determine related biological functions, and these dysregulated genes were related to cell adhesion (42 genes), extracellular matrix organization (24 genes), positive regulation of cell migration (19 genes), positive regulation of cell proliferation (31 genes), inflammatory response (27 genes), skeletal system development (15 genes), negative regulation of endopeptidase activity (14 genes), regulation of blood pressure (10 genes), nervous system development (21 genes), and kidney development (11 genes). The selected criteria for functional annotation analysis were Expression Analysis Systematic Explorer (EASE) score = 0.1 and fold enrichment >1.3. The −log (p value) of each biological term is indicated within the pie chart.   Through the functional annotation analysis in the DAVID database, the top 10 biological processes involved in these differentially expressed genes were identified to be related to cell adhesion (42 genes, p = 6.87 × 10 −12 , false discovery rate (FDR) = 1.  Figure 1B. The functional annotation clustering analysis of the 495 differentially expressed genes indicated the enrichment score of 2.47 for a group of genes related to biological functions of cartilage and bone. In addition, other significant biological functions (p < 0.05, fold enrichment >1.3) related to cellular, tissue and structural changes in OA joint, including regulation of cartilage development, chondrocyte differentiation and bone mineralization, endochondral ossification, and the corresponding differentially expressed genes are listed in Table 2. A total of nine up-regulated genes and six down-regulated genes in the OA knee chondrocytes were identified.

Identification of Dysregulated Genes Related to Joint Structural Damage in OA
To investigate the expression patterns of genes related to biological processes of cellular, tissue and structural changes in OA joint in clinical specimens, we searched in the GEO database for OA-related datasets. One RNA sequencing dataset consisting of 18 normal knee cartilages and 20 OA knee cartilages was found (GSE114007). One array comparing the expression profiles of subchondral bones at medial and lateral tibial plateaus of both normal and OA knee joints was also found (GSE51588). Considering the close contact of cartilage and subchondral bones and possible crosstalk between these joint tissues in the OA joint microenvironment, we simultaneously investigated the expression patterns of these genes in the subchondral bone dataset. The results showed that the expression patterns of ACVRL1, ALPL, GDF6, TMEM119, WNT5A, CYTL1, and SOX9 in OA knee cartilages were consistent with our expression profile of OA chondrocytes (Figure 2A  were significantly up-regulated, whereas CYTL1 and SOX9 were significantly down-regulated in OA knee cartilages. The expression patterns of GPM6B and SMAD3 (down-regulated) in OA knee cartilage were different from that identified in our expression profile of OA knee chondrocytes (up-regulated). * indicated p < 0.05, ** indicated p < 0.01, *** indicated p < 0.001, and n.s. indicated no statistical significance between normal and OA groups. Expression patterns of genes related to the biological processes of OA in the GEO dataset of OA knee cartilages. The expression patterns of (A) 9 up-regulated genes and (B) 6 down-regulated genes related to biological processes of cellular, tissue and structural changes in OA were analyzed in a GEO dataset of normal and OA knee cartilages (GSE114007). ACVRL1, ALPL, GDF6, TMEM119, and WNT5A were significantly up-regulated, whereas CYTL1 and SOX9 were significantly down-regulated in OA knee cartilages. The expression patterns of GPM6B and SMAD3 (down-regulated) in OA knee cartilage were different from that identified in our expression profile of OA knee chondrocytes (up-regulated). * indicated p < 0.05, ** indicated p < 0.01, *** indicated p < 0.001, and n.s. indicated no statistical significance between normal and OA groups.
dataset of normal and OA knee cartilages (GSE114007). ACVRL1, ALPL, GDF6, TMEM119, and WNT5A were significantly up-regulated, whereas CYTL1 and SOX9 were significantly down-regulated in OA knee cartilages. The expression patterns of GPM6B and SMAD3 (down-regulated) in OA knee cartilage were different from that identified in our expression profile of OA knee chondrocytes (up-regulated). * indicated p < 0.05, ** indicated p < 0.01, *** indicated p < 0.001, and n.s. indicated no statistical significance between normal and OA groups.

SMAD3 and WNT5A were Involved in Growth of Blood Vessel and Cell Aggregation
Diseases and functions related to cellular, tissue and structural changes of cartilage in OA analyzed in the IPA software were selected, including osteoarthritis (p = 5.53 × 10 −11 , Z score = 0.762), abnormality of cartilage tissue (p = 5.06 × 10 −7 , Z score = 1.000), differentiation of chondrocytes (p = 6.60 × 10 −8 , Z score = −0.265), growth of blood vessel (p = 4.29 × 10 −7 , Z score = 0), and aggregation of cells (p = 1.15 × 10 −9 , Z score = 2.262); the merged network of these diseases and functions is shown in Figure 4. Among the seven significantly dysregulated genes identified in the previous section, we found that SMAD3 was centered to the merged network, and interconnected with molecules of other networks, including those in the growth of blood vessel, a feature of cartilage damage in OA [9]. In addition, WNT5A was involved in the aggregation of cells related to the formation of chondrocyte clusters in OA [7,9,30]. The 495 differentially expressed genes identified in OA knee chondrocytes were analyzed in the Ingenuity Pathway Analysis (IPA) software, and diseases and functions related to cellular, tissue and structural changes of cartilage in OA, including osteoarthritis, abnormality of cartilage tissue, differentiation of chondrocytes, growth of blood vessel, and aggregation of cells, were selected to identify target genes involved. SMAD3 was centered to the merged network and interconnected with the molecules of other networks, as indicated in light blue lines, including those in the growth of the blood vessel. WNT5A was involved in aggregation of cells, which was possibly associated with the formation of chondrocyte clusters in OA.

Identification of Differentially Expressed miRNAs and Potential miRNA-mRNA Interactions between Normal and OA Knee Chondrocytes
To determine the potential miRNA-mRNA interactions between normal and OA knee chondrocytes, we also performed small RNA sequencing by NGS, and identified 46 differentially expressed miRNAs (>2.0-fold change and >10 RPM of either origin of chondrocytes), including nine up-regulated and 37 down-regulated miRNAs in OA knee chondrocytes compared to normal chondrocytes. The detailed information of the 46 differentially expressed miRNAs is provided in Table S4. The heat map of the 46 differentially expressed miRNAs with z-score values is shown in Figure 5A. The putative targets of the 46 miRNAs were predicted using miRmap database and selected targets with miRmap score >99.0. The analytical results yielded 274 putative targets of nine up-regulated miRNAs and 1041 putative targets of 37 down-regulated miRNAs. The 274 targets of up-regulated miRNAs and 1041 targets of down-regulated miRNAs were then matched to our 232 down-regulated genes and 263 up-regulated genes obtained from high-throughput RNA sequencing, respectively. Using the Venn diagram analysis available on the website [31], four down-regulated genes and 22 up-regulated genes with potential miRNA-mRNA interactions in OA chondrocytes were identified, as depicted in Figure 5B. The 495 differentially expressed genes identified in OA knee chondrocytes were analyzed in the Ingenuity Pathway Analysis (IPA) software, and diseases and functions related to cellular, tissue and structural changes of cartilage in OA, including osteoarthritis, abnormality of cartilage tissue, differentiation of chondrocytes, growth of blood vessel, and aggregation of cells, were selected to identify target genes involved. SMAD3 was centered to the merged network and interconnected with the molecules of other networks, as indicated in light blue lines, including those in the growth of the blood vessel. WNT5A was involved in aggregation of cells, which was possibly associated with the formation of chondrocyte clusters in OA.

Identification of Differentially Expressed miRNAs and Potential miRNA-mRNA Interactions between Normal and OA Knee Chondrocytes
To determine the potential miRNA-mRNA interactions between normal and OA knee chondrocytes, we also performed small RNA sequencing by NGS, and identified 46 differentially expressed miRNAs (>2.0-fold change and >10 RPM of either origin of chondrocytes), including nine up-regulated and 37 down-regulated miRNAs in OA knee chondrocytes compared to normal chondrocytes. The detailed information of the 46 differentially expressed miRNAs is provided in Table  S4. The heat map of the 46 differentially expressed miRNAs with z-score values is shown in Figure 5A. The putative targets of the 46 miRNAs were predicted using miRmap database and selected targets with miRmap score >99.0. The analytical results yielded 274 putative targets of nine up-regulated miRNAs and 1041 putative targets of 37 down-regulated miRNAs. The 274 targets of up-regulated miRNAs and 1041 targets of down-regulated miRNAs were then matched to our 232 down-regulated genes and 263 up-regulated genes obtained from high-throughput RNA sequencing, respectively. Using the Venn diagram analysis available on the website [31], four down-regulated genes and 22 up-regulated genes with potential miRNA-mRNA interactions in OA chondrocytes were identified, as depicted in Figure 5B.

Analysis of Candidate Genes with Potential miRNA-mRNA Interactions in Gene Expression Omnibus (GEO) Database and Identification of Potential Molecular Signatures in OA Knee Joint Microenvironment
The 26 candidate genes with potential miRNA-mRNA interactions are listed in Table 3. We searched in the GEO database for OA-related datasets to further determine the expression patterns of these candidate genes in clinical specimen. Since OA is considered a disease involving all joint tissues, we took datasets of cartilage, subchondral bone and synovium into consideration while searching for related datasets. There were two representative arrays analyzing synovial tissues of non-arthritis and OA patients (GSE55235 and GSE55457). Along with the datasets of cartilages and subchondral bones from knee OA patients used in the previous section (GSE114007 and GSE51588), the expression patterns of the 26 candidate genes were sequentially analyzed in these four representative datasets to determine genes potentially dysregulated in the OA microenvironment, and are summarized in Table 4. The up-regulated MARCKS was found up-regulated in OA knee cartilages, subchondral bones of OA tibial plateaus and one of the OA synovial tissue arrays, while the down-regulated EREG was found down-regulated in subchondral bones of OA tibial plateaus and one of the OA synovial arrays. Additionally, LRRC15 was up-regulated in the OA knee cartilages,  Table 3. We searched in the GEO database for OA-related datasets to further determine the expression patterns of these candidate genes in clinical specimen. Since OA is considered a disease involving all joint tissues, we took datasets of cartilage, subchondral bone and synovium into consideration while searching for related datasets. There were two representative arrays analyzing synovial tissues of non-arthritis and OA patients (GSE55235 and GSE55457). Along with the datasets of cartilages and subchondral bones from knee OA patients used in the previous section (GSE114007 and GSE51588), the expression patterns of the 26 candidate genes were sequentially analyzed in these four representative datasets to determine genes potentially dysregulated in the OA microenvironment, and are summarized in Table 4. The up-regulated MARCKS was found up-regulated in OA knee cartilages, subchondral bones of OA tibial plateaus and one of the OA synovial tissue arrays, while the down-regulated EREG was found down-regulated in subchondral bones of OA tibial plateaus and one of the OA synovial arrays. Additionally, LRRC15 was up-regulated in the OA knee cartilages, synovial tissues and the subchondral bones of medial, but not lateral tibial plateau. The expression patterns of the 26 candidate genes in the OA knee cartilage dataset (GSE114007) is shown in Figure 6.  The genes and their expression patterns marked in bold indicated those with more consistent expression patterns in OA related datasets. * For those genes with more than one probe within the array, the expression patterns were considered consistent only if there were at least two probes with significant dysregulation in the same direction. UP, significantly up-regulated in OA (p < 0.05); DOWN, significantly down-regulated in OA (p < 0.05); n.s., non-significant between normal and OA. -indicated no identical probes within the array. † indicated there were at least two probes for the gene and the expression of the gene was significant only in one of the probes.
expression patterns in OA related datasets. * For those genes with more than one probe within the array, the expression patterns were considered consistent only if there were at least two probes with significant dysregulation in the same direction. UP, significantly up-regulated in OA (p < 0.05); DOWN, significantly down-regulated in OA (p < 0.05); n.s., non-significant between normal and OA.
--indicated no identical probes within the array. † indicated there were at least two probes for the gene and the expression of the gene was significant only in one of the probes.

Identification of Potential miRNA-mRNA Interactions of LRRC15, MARCKS, and EREG in OA Knee Chondrocytes
The consistent expression patterns of LRRC15, MARCKS, PDE3A, and EREG found in our NGS data and datasets of OA knee cartilages, subchondral bones and synovial tissues were analyzed in the miRmap database for potential miRNA regulations of these genes. Those potential miRNA regulations of LRRC15, MARCKS, PDE3A, and EREG with miRmap score >99.0 were selected, and matched to the 46 differentially expressed miRNAs identified from our miRNA data. The putative 3'UTR binding sites and sequences of the potential miRNA regulations identified from miRmap database were further validated in two other miRNA prediction databases, TargetScan and miRDB, and the results are listed in Table 5. There was one target binding site for miR-140-3p in the 3'UTR of MARCKS ( Figure S3) and two target binding sites for miR-301a-3p in the 3'UTR of EREG ( Figure S4) validated in all three miRNA prediction databases. However, no validated target binding site for miR-140-5p in the 3'UTR of LRRC15 was found. In addition, the target binding site for miR-495-3p in the 3'UTR of PDE3A was validated in miRDB, but not in the TargetScan database.

MARCKS and EREG were Potentially Involved in the Pathogenesis of Arthritic Knee Joint Pain
The 26 candidate genes were uploaded to the DAVID database and IPA software to disclose biological functions related to these differentially expressed genes of OA knee chondrocytes. The functional annotation analysis performed in the DAVID database identified 11 biological functions and their related genes; among them, three were directly related to nerve function, including axon extension, axon guidance, and synapse assembly ( Table 6). The IPA network analysis of these genes yielded three categorized networks, as listed in Table 7. Network 1 related to cellular movement, cardiovascular system development and function, and organismal development possessed the highest score among the three networks, with 13 of the 26 genes involved in this network. The interconnection of molecules in network 1 is shown in Figure 7, with MARCKS and EREG in the same network. The overlay diseases and functions analysis identified BDNF, FGF7, MARCKS, and SEMA3A to be associated with "outgrowth of neurite", as indicated by the purple frame in Figure 7. Since MARCKS was associated with the outgrowth of neurite, we further explored the potential roles of MARCKS and EREG on the pain condition of OA. Diseases and functions involved in the differentially expressed genes of OA knee chondrocytes were categorized using the IPA software, and networks of genes related to abnormality of cartilage tissue (p = 5.06 × 10 −7 , Z score = 1.000), proliferation of connective tissue cells (p = 1.96 × 10 −9 , Z score = 0.512), growth of neurites (p = 6.40 × 10 −11 , Z score = 0.918), and allodynia (p = 1.68 × 10 −7 , Z score = 2.547) were merged to identify the connections between MARCKS, EREG and other differentially expressed genes ( Figure 8). Using the "Connect" tool in the IPA, MARCKS was connected to APOE and CXCL12 in this merged network, while CXCL12 was one of the molecules involved in allodynia. In addition, EREG was linked to WNT5A, TNC, PTGS2, and TNFAIP6. Since MARCKS was associated with the outgrowth of neurite, we further explored the potential roles of MARCKS and EREG on the pain condition of OA. Diseases and functions involved in the differentially expressed genes of OA knee chondrocytes were categorized using the IPA software, and networks of genes related to abnormality of cartilage tissue (p = 5.06 × 10 −7 , Z score = 1.000), proliferation of connective tissue cells (p = 1.96 × 10 −9 , Z score = 0.512), growth of neurites (p = 6.40 × 10 −11 , Z score = 0.918), and allodynia (p = 1.68 × 10 −7 , Z score = 2.547) were merged to identify the connections between MARCKS, EREG and other differentially expressed genes (Figure 8). Using the "Connect" tool in the IPA, MARCKS was connected to APOE and CXCL12 in this merged network, while CXCL12 was one of the molecules involved in allodynia. In addition, EREG was linked to WNT5A, TNC, PTGS2, and TNFAIP6.  The merged network of genes related to the abnormality of cartilage tissue, proliferation of connective tissue cells, growth of neurites, and allodynia categorized from differentially expressed genes of OA knee chondrocytes by IPA software is shown. MARCKS was connected to APOE and CXCL12, as indicated in light blue lines, where CXCL12 was demonstrated to be simultaneously involved in all four networks. EREG was mainly linked to WNT5A, TNC, PTGS2, and TNFAIP6, as indicated with blue lines.

Discussion
In the current study, we identified dysregulated genes related to OA joint damage from primary OA knee chondrocytes to have similar expression patterns observed in the datasets of cartilages and subchondral bones of OA knees (GSE114007 and GSE51588). Through NGS and bioinformatics analysis, we also identified 46 differentially expressed miRNAs and 26 candidate genes potentially involved in OA knee chondrocytes. Of the 26 candidate genes, MARCKS, EREG, LRRC15, and PDE3A exerted similar expression patterns as observed in the datasets of cartilages, subchondral bones and synovial tissues of OA patients. Two potential miRNA-mRNA interactions were identified, including miR-140-3p-MARCKS and miR-301a-3p-EREG, and were validated systematically. These novel miRNA regulations are proposed to play critical roles in the pathogenesis of the altered OA knee joint microenvironment and to contribute partly to arthritic joint pain.
The close anatomical contact of articular cartilage to subchondral bone, altered bone structure correlating to aggravated cartilage degeneration, increased vascularization of the bone-cartilage interface observed during OA progression, and similar presentation of epigenetic phenotypes of eroded subchondral bone and cartilage in OA suggest potential crosstalk and molecular transport between cartilage, subchondral bone and the joint microenvironment [32][33][34]. The emerging role of extracellular vesicles in cartilage homeostasis and arthritis potentiated the possible mechanism of Figure 8. Analysis of the potential connections between MARCKS, EREG and OA joint pain condition. The merged network of genes related to the abnormality of cartilage tissue, proliferation of connective tissue cells, growth of neurites, and allodynia categorized from differentially expressed genes of OA knee chondrocytes by IPA software is shown. MARCKS was connected to APOE and CXCL12, as indicated in light blue lines, where CXCL12 was demonstrated to be simultaneously involved in all four networks. EREG was mainly linked to WNT5A, TNC, PTGS2, and TNFAIP6, as indicated with blue lines.

Discussion
In the current study, we identified dysregulated genes related to OA joint damage from primary OA knee chondrocytes to have similar expression patterns observed in the datasets of cartilages and subchondral bones of OA knees (GSE114007 and GSE51588). Through NGS and bioinformatics analysis, we also identified 46 differentially expressed miRNAs and 26 candidate genes potentially involved in OA knee chondrocytes. Of the 26 candidate genes, MARCKS, EREG, LRRC15, and PDE3A exerted similar expression patterns as observed in the datasets of cartilages, subchondral bones and synovial tissues of OA patients. Two potential miRNA-mRNA interactions were identified, including miR-140-3p-MARCKS and miR-301a-3p-EREG, and were validated systematically. These novel miRNA regulations are proposed to play critical roles in the pathogenesis of the altered OA knee joint microenvironment and to contribute partly to arthritic joint pain.
The close anatomical contact of articular cartilage to subchondral bone, altered bone structure correlating to aggravated cartilage degeneration, increased vascularization of the bone-cartilage interface observed during OA progression, and similar presentation of epigenetic phenotypes of eroded subchondral bone and cartilage in OA suggest potential crosstalk and molecular transport between cartilage, subchondral bone and the joint microenvironment [32][33][34]. The emerging role of extracellular vesicles in cartilage homeostasis and arthritis potentiated the possible mechanism of cell-cell communication between joint tissues through carriage of bioactive signaling molecules, such as miRNAs by extracellular vesicles [35][36][37]. Increasing evidence suggested extracellular vesicles released from chondrocytes may participate in dysregulated cartilage mineralization and contribute to OA [35,38]. Our study results identified SMAD3 to potentially play a central role in the merged network of diseases and functions related to cellular, tissue and structural changes in OA, and serves as one of the molecules in the network of "growth of vessel" (Figure 4). In the OA joint, articular cartilage loses resistance to vascularization, and the growth of blood vessels increased particularly at the bone-cartilage interface [39]. TGF-β/SMAD2/3 signaling was shown to exert a chondroprotective effect in excessive mechanical compression by blocking chondrocyte differentiation [40]. Conversely, other studies indicated TGF-β/SMAD3 negatively regulated miR-140 in OA chondrocytes [41], and SMAD3 was significantly up-regulated in human OA cartilage tissue [42]. Additionally, TGF-β/Smad3 treatment enhanced vascular endothelial growth factor expression and inhibited vascular smooth muscle apoptosis through autocrine signaling [43]. These results suggested that TGF-β/Smad3 may have dual roles in the vascularization of OA cartilage and the development of OA.
In a review article of WNT5A signaling, the induced WNT5A activates non-canonical Wnt signaling to regulate the control of tissue polarity and cell aggregation and canonical WNT signaling for the proliferation of human synovial fibroblasts in rheumatoid arthritis (RA) [44]. In addition, the Wnt/β-catenin signaling is suggested to participate in the degradation of ECM in the degenerative joint of animal models by stimulating matrix catabolic genes and activity in the articular chondrocytes [45]. Recent studies also presented the positive correlation of the expression of Wnt5a in human articular cartilage to the progression of knee OA [46], and showed that Wnt5a promotes MMP production in human chondrocyte through non-canonical Wnt signaling [47]. The up-regulated WNT signaling in the inflammatory response of human chondrocytes was also proposed [48]. In the OA joint, chondrocytes aggregate to form clusters around the bone-cartilage interface [30]. Our current study identified the up-regulated WNT5A in OA knee chondrocytes, and bioinformatics analysis suggested the involvement of WNT5A in the positive regulation of cartilage development and the aggregation of cells, suggesting the important role of WNT5A in cartilage degradation in knee OA.
The differentially expressed miRNAs in the cartilage and bone of human knee OA with potential target genes related to inflammatory processes were identified [49]. Using the NGS and bioinformatics analysis, we identified the potential interaction of miR-140-3p-MARCKS in OA knee chondrocytes. Literature has identified miR-140 as a major contributor to cartilage homeostasis, targeting genes involved in the differentiation of chondrocytes [10,50]. Miyaki et al. reported the increased expression of miR-140 during chondrogenesis in human mesenchymal stem cells, the reduced expression of miR-140 in OA cartilage and in chondrocytes in response to IL-1β, and confirmed the dual roles of miR-140 in cartilage development and homeostasis in transgenic mouse model [51,52]. In a recent study analyzing the expressions of miR-140-3p and miR-140-5p in synovial fluid of OA and non-OA patients, both miRNAs were down-regulated in OA synovial fluid and were associated with the severity of OA [53]. The down-regulated miR-140-3p and miR-140-5p were also observed in our small RNA sequencing result of OA knee chondrocytes, indicating the potential target of miR-140 as biomarkers in knee OA.
Pain is a major symptom in patients with OA. The mechanisms involved in arthritic joint pain are complicated, while structural pathologies, neuronal mechanisms of pain, and general factors such as obesity and genetic factors shall all take part in the consequence of joint pain [54]. Central and peripheral sensitizations of the nociceptive system are extensively proposed mechanisms of neuronal causes of OA joint pain [55]. Normal articular cartilage lacks vascular supply and sensory innervation [56]; structural changes of innervation in OA may be potential contributors to joint pain, with an evident increase in nerve growth and vascular penetration to the meniscus [57,58].
Suri et al. also reported the presence of nerve fibers within vascular channels of OA cartilage and subchondral bone marrow, which may contribute to pain in knee OA [33]. In our current study, we identified the increased expression of MARCKS in OA knee chondrocytes, and observed the similar expression patterns of MARCKS in OA subchondral bones and synovial tissues. Myristoylated alanine-rich protein kinase C substrate (MARCKS) is an actin-binding protein distributed in the nervous system and is localized in axons and dendrites; its phosphorylation is involved in the release of neurotransmitters and synaptic trafficking, mechanisms related to inflammatory pain [59]. Bradykinin is a peptide involved in inflammation, with its B2 receptors expressed on various OA-related cell types, implicating its involvement in OA pathogenesis [60]. Tanabe et al. also reported a novel bradykinin signaling cascade of neurite outgrowth through MARCKS phosphorylation in a neuroblastoma cell line [61]. Studies also suggested local release of neurotrophins by non-neuronal tissues in inflammatory joint lining that facilitates nerve ingrowth into tissues [62], and released neurotrophins may exert autocrine/paracrine effects that regulate articular chondrocyte metabolism [63]. This provides clues to the connection between neurite outgrowth and the pathogenesis of arthritic joint pain. There is emerging evidence regarding the possible regulatory function of miRNAs in nociceptive pathways and it is a potential therapeutic target in OA management [12]. Our current results identified the potential role of miR-140-3p-MARCKS regulation in the OA knee joint microenvironment. Whether targeting this novel miRNA regulation assists in the management of joint pain or not, and whether it takes part in the altered structural innervation or peripheral sensitization in knee OA merits further investigation.
Another novel miRNA regulation of miR-301a-3p-EREG was identified in our current study. MiR-301a-3p is a pro-inflammatory miRNA that promotes inflammation during tumorigenesis [64]. There is not much literature reporting the role of miR-301a-3p in the regulation of arthritis. Two recent studies reported the regulation of miR-301a-3p in inflammatory response of RA and chondrocytes. In patients with RA, over-expressed miR-301a-3p was found in the peripheral blood mononuclear cells and it was positively correlated to an increased proportion of a subset of helper T cell [65]. Chen et al. studied the role of miR-301a in a mouse chondrogenic cell line in response to inflammatory injury using lipopolysaccharide (LPS), and reported the increased expression of miR-301a with LPS induction and the suppression of miR-301a alleviated LPS-induced chondrogenic cell injury by targeting Sirt1, implicating the potential therapeutic target of miR-301a for OA [66]. The down-regulated EREG in our NGS result was the putative target of miR-301a-3p. Epiregulin (EREG), an EGFR ligand, is expressed in various tissues; it plays important roles in wound healing and vascular remodeling during inflammation and regulates cell proliferation and differentiation [67]. A recent study by Martin et al. demonstrated the genetic association of EREG and EGFR with chronic pain in patients with temporo-mandibular joint disorder, indicating the role of EREG in pain processing [68]. However, the association of EREG expression in the joint microenvironment remains unclear. Only one report demonstrated the increased expression of EREG and pro-inflammatory cytokines in fibroblast-like synoviocytes, suggesting the role of EREG-mediated EGFR signaling to the pathogenesis of RA [69]. In our current study, we observed the similarly down-regulated EREG in OA knee chondrocytes, subchondral bone and synovial tissue of OA patients. Further investigation into the role of miR-301a-3p-EREG regulation in OA may provide further insight into better understanding novel miRNA regulation in the pathogenesis of knee OA.
The characteristic change of OA is cartilage breakdown, but a growing consensus has proposed OA as a disease of the whole joint, involving all joint tissues [9]. In the current study, we identified several important molecules involved in the structural damage of OA knee joint and explored novel miRNA-mRNA regulations potentially involved in the pathogenesis of OA joint pain. The schematic potential molecular signatures and miRNA regulations related to OA structural damage and joint pain are summarized in Figure 9. These identified novel miRNA-mRNA regulations may provide new insights into the potential targets for arthritic joint pain. Further validation with experimental studies on clinical OA knee cartilage specimens are necessary to consolidate the clinical significance of the current in silico results.

Conclusions
Our current study identified the critical roles of SMAD3 and WNT5A in the OA knee joint microenvironment, which is potentially linked to features of ECM degradation in knee OA. Additionally, miR-140-3p-MARCKS and miR-301a-3p-EREG regulations were potentially involved in the pathogenesis of the altered OA knee joint microenvironment. The current findings suggest new perspectives in studying novel genes potentially contributing to arthritic joint pain in knee OA, which may assist in finding new targets for OA treatment.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Table S1: Available information of selected osteoarthritis related datasets from GEO database, Table S2: RNA and small RNA sequencing summary, Table S3: The list of 495 differentially expressed genes in OA knee chondrocytes, Table S4: The list of 46 differentially expressed miRNAs in OA knee chondrocytes, Figure S1: RNA sequencing quality control report for normal adult chondrocytes (HC) and osteoarthritic knee chondrocytes (HC-OA), Figure S2: Small RNA sequencing quality control report for HC and HC-OA, Figure

Conclusions
Our current study identified the critical roles of SMAD3 and WNT5A in the OA knee joint microenvironment, which is potentially linked to features of ECM degradation in knee OA. Additionally, miR-140-3p-MARCKS and miR-301a-3p-EREG regulations were potentially involved in the pathogenesis of the altered OA knee joint microenvironment. The current findings suggest new perspectives in studying novel genes potentially contributing to arthritic joint pain in knee OA, which may assist in finding new targets for OA treatment.
Supplementary Materials: The following are available online at http://www.mdpi.com/2077-0383/7/12/535/s1, Table S1: Available information of selected osteoarthritis related datasets from GEO database, Table S2: RNA and small RNA sequencing summary, Table S3: The list of 495 differentially expressed genes in OA knee chondrocytes, Table S4: The list of 46 differentially expressed miRNAs in OA knee chondrocytes, Figure S1: RNA sequencing quality control report for normal adult chondrocytes (HC) and osteoarthritic knee chondrocytes (HC-OA), Figure S2: Small RNA sequencing quality control report for HC and HC-OA, Figure S3