Pain-Associated Transcriptome Changes in Synovium of Knee Osteoarthritis Patients

Joint pain causes significant morbidity in osteoarthritis (OA). The aetiology of joint pain in OA is not well understood. The synovial membrane as an innervated joint structure represents a potential source of peripheral pain in OA. Here we analyse, using a hypothesis-free next generation RNA sequencing, the differences in protein-coding and non-coding transcriptomes in knee synovial tissues from OA patients with high knee pain (n = 5) compared with OA patients with low knee pain (n = 5), as evaluated by visual analogue scale (VAS). We conduct Gene Ontology and pathway analyses on differentially expressed mRNA genes. We identify new protein-coding, long non-coding RNA and microRNA candidates that can be associated with OA joint pain. Top enriched genes in painful OA knees encode neuronal proteins that are known to promote neuronal survival under cellular stress or participate in calcium-dependent synaptic exocytosis and modulation of GABA(γ-aminobutyric acid)ergic activity. Our study uncovers transcriptome changes associated with pain in synovial microenvironment of OA knees. This sets a firm ground for future mechanistic studies and drug discovery to alleviate joint pain in OA.


Introduction
Primary knee osteoarthritis (OA) is the most common joint disorder, affecting more than 10% of the Western European population aged 60 years or older [1]. Joint pain and functional disability decrease quality of life in OA patients [2]. OA affects hands, feet, knees, hips and spine [3] with structural damage in the articular cartilage, subchondral bone, ligaments, tendons, menisci, muscles, synovium and nerve tissues [4]. The ethiopathogenesis of joint pain in OA is not well understood [5] and the therapeutic strategies to relieve OA pain are therefore limited. The correlation between the severity of pain and the degree of radiographic changes in OA knees is poor [6,7].
A number of studies substantiate the role of the central nervous system in the ethiopathogenesis of chronic pain in OA [8,9], nevertheless the sources of peripheral nociception are increasingly explored [6,10]. Peripheral pain can originate in any of the innervated joint structures, including the synovium [5]. Synovial inflammation in OA causes joint pain [11]. Inflammation decreases the activation threshold of local afferent nerve fibers in response to mechanical stimuli resulting in peripheral sensitisation. Additionally, interactions between a damaged joint and the sensory nervous system cause pain [11]. Most likely, there are two interconnected mechanisms with reverse causality: Joint damage causes synovial inflammation, whereas neurogenic inflammation contributes to joint damage, creating a positive feedback loop [12].
Genes involved in pain pathways in the nervous system are also expressed in non-neuronal cells of the joint [13]. Epidemiology studies identified a panel of candidate genes and functional genetic variants associated with OA pain [14]. To date, no hypothesis-free study with a genome-wide approach has been conducted to interrogate the molecular basis of OA pain.
Here we analyse by RNA sequencing the transcriptomes of knee synovial tissues from OA patients with high and low intensity knee pain as assessed by visual analogue scale (VAS). We uncover the molecular complexity of OA joint pain and reveal novel candidate genes associated with knee OA pain.

Patient Information
OA patients (Table 1) were enrolled in the study at the Uniklinik Balgrist, Zurich, Switzerland. Pain was assessed a day before surgery by a clinical investigator using a 10-point VAS. Patients with low (VAS scores: 0-3) or high (VAS 7-10) pain were enrolled. The indication for surgery in patients with low pain were a relevant deformity of the axis, mainly valgus, and/or loss of bone stock with a symptomatic instability and risk to fall. The Local Ethical Committee (KEK-ZH-Nr. 2013-0210) approved the study. All patients gave written informed consent. All experiments were performed in accordance with relevant guidelines and regulations.

Sample Collection and Total RNA Extraction
The standardized surgical excision (excision of the capsule with synovia in the recessus suprapatellaris) was performed. The biopsies (two per patient with the size of approximately 2 cm 2 ) were fresh frozen (dry ice) and RNA was extracted randomly from one biopsy using TRIzol ® Reagent (Thermo Fisher Scientific, Waltham, MA, USA). RNA was analysed by Qubit ® 1.0 Fluorometer (Thermo Fisher Scientific) and Agilent 2200 TapeStation (Agilent Technologies, Santa Clara, CA, USA); RNA with RIN (RNA Integrity Number) values 8.3-9.1 was used for sequencing.

Next Generation Sequencing
Libraries for mRNA/long noncoding RNA (lncRNA) and small RNA sequencing were prepared using the Illumina (San Diego, CA, USA) TruSeq ® Stranded total RNA Sample Preparation Kit and the NEBNext ® Small RNA Library Prep Set for Illumina ® , respectively. Libraries were quantified and quality checked using quantitative polymerase chain reaction (qPCR) (Roche, Basel, Switzerland) with Illumina adapter specific primers and Agilent 2200 TapeStation, respectively. Diluted indexed mRNA-seq (10 nM) and small RNA-seq (1 nM) libraries were pooled, used for cluster generation (Illumina TruSeq SR Cluster Kit v4-cBot-HS) and sequenced (Illumina HiSeq 2500, Illumina TruSeq SBS Kit v4-HS reagents, single read approach (1 × 125 bp) with 20-30 million reads per sample for mRNA/lncRNA, 10-20 million reads for miRNA).

Data Analysis
After quality control (FastQC, Babraham Bioinformatics, https://www.bioinformatics.babraham. ac.uk/) we quantified transcript expressions using the RSEM package (version 1.2.18) [16]. Transcript references were built using the Ensembl version 75 annotations of the GRCh37.p13 human genome assembly. The expression level for each gene was obtained by the aggregation of counts from all gene isoforms. Genes with less than 10 counts in ≥50% of samples for each condition were excluded from the analysis. Small RNA reads were processed using the ncPRO-seq [17] and annotated using miRbase version 20 [18], RepeatMasker and Rfam [19].
Differentially expressed genes were identified using the R/bioconductor package edgeR (v3.14.0) [20]. The Benjamini-Hochberg procedure [21] was used to adjust p values and genes with an adjusted p value ≤ 0.11 were considered differentially expressed. Gene ontology (GO) enrichment analysis was conducted using the Overrepresentation Test from PANTHER web portal (release 2015043074) [22] on biological process ontologies, using Fisher's exact test option and all genes in the database as a background. Redundant GO terms were removed using the REVIGO portal [23]. Pathway analysis was performed using the gene set enrichment analysis tool using MSigDB curated gene sets (C2, release 5) [24,25]. GO terms and pathways with adjusted p value ≤0.05 were considered enriched. Small RNAs with an absolute log2 fold change ≥1 and p value ≤0.05 were defined as differentially expressed.
Deconvolution of synovial tissue into immune cell types was evaluated using an mRNA machine learning approach of CIBERSORT [26], running the relative and absolute modes and using 500 permutations. Input matrix contained fragments per kilobase million (FPKM) for all expressed genes. CIBERSORT uses the reference dataset LM22 that consist of 22 distinct immune cell types and was constructed from gene expression profiles of those cell types measured on the Affymetrix U133A/Plus2 (Santa Clara, CA, United States and Illumina Expression BeadChip (HumanHT-12 v4) platforms.
The raw and processed RNA-seq data have been deposited under Gene Expression Omnibus (GEO) accession number GSE99662.

Protein-Coding Transcriptome and Gene Enrichment Analysis
We identified 32 genes that were differentially expressed in the synovium of OA patients with low and high pain ( Table 2). Among these 32 mRNAs, 30 mRNAs were up-regulated in OA synovium from patients with high pain compared with patients experiencing low pain. The top up-regulated genes were the Stress Responsive DNAJB4 Interacting Membrane Protein (SDIM1), Carboxypeptidase E (CPE) and Otoferlin (OTOF). Ankyrin Repeat Domain 20 Family Member A4 (ANKRD20A4) and pentraxin 3 (PTX3) were the only down-regulated genes in OA synovium from patients experiencing high pain (Table 2).
To identify biological processes enriched in the synovium of OA patients with high knee pain, we performed the GO enrichment analysis and gene set enrichment analysis (GSEA) on the 32 differentially expressed genes. Six GO sets were enriched in patients with high pain (Table 3). Only two of the enriched GO term have a relative small size (<200 genes), Peripheral nervous system development (GO:0007422) and Second-messenger-mediated signalling (GO:00019932). The characteristics of the 18 genes that were enriched in these GO terms are described in the Table 4 and Table S1. We identified two enriched Reactome pathways, specifically Nerve growth factor (NGF) signalling via Tropomyosin receptor kinase A (TRKA) from the plasma membrane and Signalling by NGF (Table 5). These two pathways were redundant with four significantly enriched genes including adenylate cyclase 3 (ADCY3), Src homology 2 domain containing (SHC) transforming protein 2 (SHC2), neurotrophic receptor tyrosine kinase 2 (NTRK2) and dual specificity phosphatase 4 (DUSP4).   Table 4. Gene Ontologies (GO) in the category of biological process and 18 enriched genes among 32 differentially expressed in synovial tissue of patients with contrast pain phenotypes. All enriched genes are up regulated in painful knee ( Table 2).

Noncoding Transcriptome
Two lncRNAs (RN7SL3 and RP11-195E2.1) and thirty-five miRNAs were identified as differentially expressed in the synovium of OA patients with contrast pain phenotypes (Tables 1 and 6). Twenty-one miRNAs were down-and fourteen were up-regulated in patients with high compared with low knee pain, but these differences were not significant when correcting for false discovery rate (FDR) ( Table 6). Table 6. List of small RNAs significantly differentially expressed in synovial tissue of patients with contrast pain phenotypes; raw p value < 0.05 and absolute log2 fold change (logFC) ≥ 1; FDR: false discovery rate).

No.
Small

Deconvolution Analysis of Harvested Synovial Tissue
There were no major differences in immune cell composition (Kruskal-Wallis test, adjusted p value ≤ 0.05) nor in the absolute content of immune cells (Kruskal-Wallis test, p = 0.17) in the synovium between patients ( Figure 1A). The expression level of CD45 (PTPRC), a leukocyte marker ( Figure 1B) did not differ between synovial tissues from patients with high and low pain (adjusted p = 0.78). In addition, fibroblast markers CLU, COL1A2, COL3A1, identifying for different subpopulations of synovial fibroblasts [27] were expressed at high levels in synovial tissues (log2 FPKM > 9) and their expression did not differ between patients with contrast pain phenotypes ( Figure 1B).

Discussion
Pain is a subjective and complex sensory experience and evaluation of pain is rather difficult. Among several methods for pain assessment [28], we used VAS since it is a clinical routine, well accepted in scientific pain literature and it influences functional scores [29]. Patients with contrast pain phenotypes were matched for gender (females), smoking status (non-smokers) and ethnic origin (Europeans) which can confound experiencing the pain [6]. Additionally, patient groups with contrast pain phenotypes did not differ significantly in age (70.6 ± 6.5 vs. 71.0 ± 6.0 in high vs. low pain group), Body Mass Index (BMI) (31.1 ± 8.5 vs. 29.0 ± 5.1 in high vs. low pain group) or the extent of joint damage (3 and 4 in Kellgren and Lawrence grade system). The use of analgesics could decrease pain in the low pain group patients 021 and 022, but appeared inefficient in controlling the pain in the high pain group patients 016, 018, 019.
We uncovered a signature of genes that are differentially expressed in the synovial tissues from knees of OA patients with contrast pain phenotypes. The synovial inflammatory mediators that are known to induce and/or respond to pain [12] were not enriched in the transcriptome profiles in synovial tissues from patients with high compared to patients with low knee pain. The deconvolution analysis of RNA-seq data showed no differences in the composition of immune cells or synovial fibroblasts between patients. Instead, the three top differentially expressed genes in patients with high pain are neuronal proteins. Two of them, specifically SDIM1 and CPE, are known for their roles in cellular responses to stress. This suggests a potential molecular connection between chronic pain and stress, as also recently demonstrated by Descalzi et al. [30]. SDIM1 is involved in cellular response to stress and inhibition of cell death [31]. In NT2 neurons SDIM1 exhibits a bi-phasic response to cell death-inducing injuries; initial down-regulation of SDIM1 is followed by an up-regulation in surviving cells [31]. Overexpression of SDIM1 improves survival of neuro-progenitor cells after injury, substantiating the pro-survival effects of SDIM1 under stress conditions [31]. We show an increased expression of SDIM1 in knee synovial tissues from OA patients with high pain. Our data suggest that acting as a stressor, the chronic pain might increase the expression of SDIM1 in the synovial microenvironment of OA knees to activate its protective functions.
Besides SDIM1, CPE can steer cellular responses to stress and we show that also CPE is significantly up-regulated in the synovium of OA patients with high pain. CPE acts as a trophic factor and promotes neuronal survival upon various stressors via increasing anti-apoptotic protein BCL-2 Bcl-2 in a MEK/ERK and/or PI3-K/AKT-dependent manner [32]. Stress-driven increase of CPE in brain has pro-survival effects in neurons [32]. Considering the sensitivity of CPE in responding to different stressors, pain as a stressor might induce CPE expression in OA knees. CPE acts through a similar molecular mechanism like brain derived neurotrophic factor (BDNF) [33] which is expressed in knee synovium of OA patients [13], its expression however does not differ between OA patients with contrast pain phenotypes. This suggests that CPE and BDNF may not function together to promote cell survival in painful OA joints; a similar scenario proposed for hippocampal neurons [32]. Instead, we show increased expression of IGF1, which is known to have neuroprotective and pro-survival effects by activating ERK and Akt [34].
Collectively, these results suggest that the cell protective functions of SDIM1 and CPE might be needed within the stress microenvironment of the OA synovium in patients with high pain. CPE and SDIM1 appear striking molecular candidates for future studies due to its potential joint protective effects.
OTOF was the second of the top three up-regulated genes in the synovium of OA patients with high pain. OTOF is a transmembrane protein required for calcium-dependent synaptic exocytosis in cochlear sensory cells [35] and mutations in OTOF cause deafness [35]. Besides, OTOF can modulate the GABAergic activity in the GAD65-dependent manner in neuronal and non-neuronal cells [36]. Future studies might reveal whether OTOF can change the secretion of GABA in the OA synovium in response to pain, thereby modulating the pain transmission.
In the gene networks enrichment analysis using GO terms and Reactome pathways, the networks that include SDIM1 and OTOF did not appear enriched, possibly due to a rather limited published data on the function of these genes. Meanwhile, CPE was linked to the anatomical structure morphogenesis. Several other up-regulated genes in the synovium of OA patients with high pain did not appear in enriched pathways and the knowledge about most of these genes is rather limited. Nevertheless, some of these genes participate in nervous system development (CLSTN2, TUBB2B) [37,38], inflammation including arthritis (AOC3) [39] and stress-induced responses (C10orf10) [40], making them promising molecular candidates in OA-driven pain.
On the other hand, the neurotrophic tyrosine kinase receptor type 2 (NTRK2), also known as TrkB was present in every enriched GO term and in both enriched Reactome pathways. This might suggest that NTRK2 is one of the key drivers, involved in the pain-related transcriptional changes in the OA synovium. The expression of the neurotrophic factors that function via TrkB including BDNF, nerve growth factor (NGF), was not altered in OA patients with high knee pain. This suggests that other molecular candidates might function via TrkB to activate downstream signaling pathways in the synovium of painful OA joint.
To date no studies have explored the global alterations in the synovial noncoding RNA expression in relation to OA joint pain. We identified two differentially expressed lncRNAs and 35 differentially expressed miRNAs in the synovium of OA patients with contrast pain phenotypes, but the differences in the miRNA expression were not significant when correcting for FDR. In all this suggests that pain-related transcriptional changes primarily affect the protein coding transcriptome. Pain-associated changes in miRNAs in OA appear less frequent or smaller, thus larger sample sizes might be required for miRNA studies.
The two identified lncRNAs (RN7SL3 and RP11-195E2.1), which were strongly down-regulated in OA patients with high knee pain, have largely unknown function. In contrast, miR-146a-3p, the top differentially expressed miRNA in the synovium of OA patient with contrast pain phenotypes, has already been linked to the pain-related pathophysiology of knee OA [41]. MiR-146a-3p, which regulates the cell repair responses to tissue damage. It was also shown to be up-regulated in the cartilage and synovium of OA patients compared with healthy controls, which directly links this miRNA to OA pathogenesis [41]. We show that miR-146a-3p is down-regulated in the synovium of OA patients with high pain. This is in line with the observed down-regulation of miR-146a-3p in peripheral and central neurons in an animal model of OA-pain [41], making miR-146a-3p a prominent molecular candidate in pain-related pathology in OA.

Conclusions
In summary, the molecular dynamics in the synovium of painful OA joints indicate that the synovium is an important tissue to investigate in the discovery of novel therapeutics to decrease OA joint pain. The increased expression of stress-inducible molecules with protective roles in cell survival emerges as a common theme in OA patients with high intensity knee pain. Future studies of stress-induced protective mechanisms should further deepen the understanding of joint pain in OA.