Expression profiling in spondyloarthropathy synovial biopsies highlights changes in expression of inflammatory genes in conjunction with tissue remodelling genes

Background In the spondyloarthropathies, the underlying molecular and cellular pathways driving disease are poorly understood. By undertaking a study in knee synovial biopsies from spondyloarthropathy (SpA) and ankylosing spondylitis (AS) patients we aimed to elucidate dysregulated genes and pathways. Methods RNA was extracted from six SpA, two AS, three osteoarthritis (OA) and four normal control knee synovial biopsies. Whole genome expression profiling was undertaken using the Illumina DASL system, which assays 24000 cDNA probes. Differentially expressed candidate genes were then validated using quantitative PCR and immunohistochemistry. Results Four hundred and sixteen differentially expressed genes were identified that clearly delineated between AS/SpA and control groups. Pathway analysis showed altered gene-expression in oxidoreductase activity, B-cell associated, matrix catabolic, and metabolic pathways. Altered "myogene" profiling was also identified. The inflammatory mediator, MMP3, was strongly upregulated (5-fold) in AS/SpA samples and the Wnt pathway inhibitors DKK3 (2.7-fold) and Kremen1 (1.5-fold) were downregulated. Conclusions Altered expression profiling in SpA and AS samples demonstrates that disease pathogenesis is associated with both systemic inflammation as well as local tissue alterations that may underlie tissue damaging modelling and remodelling outcomes. This supports the hypothesis that initial systemic inflammation in spondyloarthropathies transfers to and persists in the local joint environment, and might subsequently mediate changes in genes directly involved in the destructive tissue remodelling.


Background
The underlying processes driving disease progression in the spondyloarthropathies (SpA) are very poorly understood. The disease transitions from an initial inflammatory insult through an inflammation-driven tissue destruction phase to an osteoproliferative phase which in the worst cases results in joint fusion. SpAs mainly present in the axial skeleton and the inaccessibility of these joints and subsequent lack of sample availability together with the slow disease progression hinders research such that the dysregulated molecular and cellular mechanisms driving disease remain largely unknown.
Expression profiling studies of affected tissues in SpA offer a hypothesis free approach to elucidating underlying pathogenic mechanisms. Previously ours and other groups have focussed largely on peripheral blood samples, either from whole blood [1,2] or from total [3,4] or partial [5] PBMC fractions. These studies provide valuable infor-mation regarding the systemic immunological processes involved in SpA, they are less informative regarding local inflammatory and tissue damage processes, in particularly the mechanisms underlying joint damage and the progression from inflammation to osteoproliferation in SpA.
Until very recently, only two small-scale tissue expression profiling studies have been undertaken in SpA, in synovial biopsies [6] and sacroiliac joint fluid [7], and no comprehensive genomic profiling study had been reported in joint tissue in SpA.
Peripheral arthritis is present in significant numbers of SpA patients with estimates between 14-20% of AS patients and 18-26% of Undifferentiated SpA patients [8]. In ankylosing spondyltitis (AS) patients with both axial and peripheral inflammation, anti-TNF treatments, such as adalimumab, have shown efficacy in reducing both peripheral and axial disease [9]. This site inclusive treatment efficacy suggests similar disease processes are occurring in these different joint environments. Subsequently this provides some justification for assessment of molecular changes within affected knee joints, that are a more accessible tissue site, as a viable approach for elucidating joint specific disease processes in SpA.
In early 2013, Yeremenko et al. published a study in which they undertook a large-scale gene expression profiling study comparing knee synovial biopsies from SpA, rheumatoid arthritis (RA) and gout patients. They demonstrated that many inflammatory genes and pathways were shared across RA and SpA. However, a "myogenic" profile was evident in the SpA samples which delineated them from the RA samples [10].
We have undertaken a similar approach, comparing archived formaldehyde-fixed paraffin-embedded (FFPE) synovial biopsies from AS, SpA, normal control and osteoarthritis (OA) patients. We similarly identified an enhanced myogene signature in our AS/SpA samples. Additionally we have also identified a number of other pathways that may contribute to tissue remodelling as well as inflammatory pathways.

Patients
Fifteen knee synovial biopsy tissue samples consisting of six seronegative spondyloarthropathy (SpA), two ankylosing spondylitis (AS), three osteoarthritis (OA) and four normal control biopsies were obtained from the Synovial Tissue Bank at the Repatriation General Hospital in Adelaide, South Australia (Additional file 1: Table S1). Biopsies were taken arthroscopically under direct vision biopsying with sampling of macroscopically abnormal appearing synovium. All patients provided informed written consent. Ethical approval for this study was obtained from the Southern Adelaide Health Service/Flinders University Human Research Ethics Committee.

RNA preparation and Microarray analysis
RNA was extracted from the biopsies embedded in formaldehyde-fixed paraffin-embedded (FFPE) tissue blocks using the Arcturus Paradise © Plus Reagent System (Molecular Devices, Sunnyvale, CA) as per the manufacturer's protocol. All the biopsy was used for the RNA extraction.
200 ng of RNA was used in the Illumina Whole-Genome DASL® (cDNA-mediated Annealing, Selection, Extension, and Ligation) Gene Expression Assay according to the Illumina protocol. This technique has been specifically developed for whole-genome expression profiling of degraded RNA samples from archived tissue biopsies. RNA is first converted to cDNA through a reverse transcription reaction with biotinylated primers. The biotinylated cDNA is then annealed to assay oligonucleotide probes specific for each of the 24000 cDNAs targeted by the array. The hybridized oligonucleotides are then extended and ligated in a second-strand cDNA synthesis forming a synthetic template that is transferred to a PCR reaction containing a fluorescently labelled primer. The labelled PCR product strand is then isolated and the fluorescent products were hybridised to Illumina Ref-8 Expression BeadChips and scanned. Gene expression is then quantified by the level of fluorescent hybridization to the BeadChip. Data was processed in GenomeStudio (Illumina) and analysed using Lumi [11] and BRB ArrayTools [12] as described previously [3]. Data was transformed by variance stabilization transformation (VST) [13] then normalized by robust spline normalization (RSN) [14]. This data has been uploaded to the NCBI GEO database and assigned accession number GSE41038.
Of the 24,500 cDNAs on the DASL arrays, 20,700 were found to be expressed in at least one sample and were included in the analysis. For analysis, AS and SpA samples were grouped together and compared with a control group consisting of OA and normal samples. Differentially expressed genes were identified by unpaired t-test with multivariate permutation correction. The evaluation of which Gene Ontology (GO) classes are differentially expressed between control and affected bones was performed using a functional class scoring analysis as described previously [2]. Efron-Tibshirani's Gene Set Analysis (GSA) was used which uses "maxmean" statistics for assessing significance of pre-defined genesets. Gene Ontology analysis was performed using BRB-ArrayTools.

Quantitative PCR
Quantitative PCR validation (qPCR) was carried out as described previously [3] in nine normal and OA samples as well as in seven SpA and AS samples. Due to low RNA yields obtained from the biopsies four of the array samples lacked sufficient RNA for confirmation qPCR follow-up but an additional five control samples were obtained for the qPCR analysis generating a partially independent confirmation cohort (Additional file 1: Table S1).
Relative expression of genes of interest were determined using the ΔCT method or standard curve method. Comparisons between different patient groups were undertaken using Mann-Whitney tests.

Immunohistochemistry
For the MMP3 immunohistochemistry, three AS, five SpA, 9 normal and 24 RA biopsies were stained. Tissue sections were blocked for endogenous peroxidase before digestion with proteinase K. This was followed by incubation first with a mouse anti-human MMP3 primary antibody (Santa Cruz Biotechnology, Santa Cruz, CA) for 2 hrs at room temperature (RT) then with a donkey anti-mouse IgG secondary antibody (Jackson ImmunoResearch Labs, West Grove, PA) for 40 mins at RT. Antibody staining was visualised with an ABC kit (Vector Laboratories, Burlingame, CA) using an AEC chromagen substrate (Dako, Cambellfield, Australia) before counterstaining with haematoxylin and mounting with Aquatek (Merck, Kilsyth, Australia). Staining was quantified using NIS Elements Br 3.0 software (Nikon, Lidcombe, Australia).

Results
To maximise the power of the study we grouped the eight AS and SpA samples together (AS-SpA) and compared them with a control group consisting of seven normal and OA (a non-inflammatory arthritis) (OA-Norm) samples for the analysis. The validity of this grouping was confirmed by unsupervised clustering that showed no differences between AS and SpA nor between OA and normal samples (data not shown). However, unsupervised clustering clearly delineated between the AS-SpA and OA-Normal groups, with only one sample from each group misclassifying ( Figure 1A, sensitivity 88%, specificity 86%).
To identify differentially expressed genes we undertook a class comparison of the two groups which showed this clustering was driven by 416 differentially expressed genes (p < 0.01) ranging from a 4.7-fold up-regulation to a 4.6-fold down-regulation ( Figure 1B, Additional file 2: Table S2).
To ascertain if there was a correlation in the tissues with systemic inflammatory genes dysregulated in our previous PBMC expression profiling studies [3] we compared the genelists. Using Gene-set Enrichment Analysis (GSEA) to calculate the degree of enrichment of the synovial biopsy genelist in the transcriptome of the AS PBMCs, Efron-Tibshirani's GSA maxmean test showed the synovial geneset was enriched in the PBMC transcriptome with a p-value of 0.005. A number of immune/inflammation-associated genes were altered in the two datasets (highlighted in Additional file 2: Table S2). The upregulated genes were CD40 (a member of the TNF receptor superfamily); CLEC12A (a member of the C-type lectin/C-type lectinlike domain superfamily); and FCGR1A (a high-affinity Fc-gamma receptor). Conversely, TSC22D3, which plays a key role in the anti-inflammatory and immunosuppressive effects of glucocorticoids, was downregulated in both PBMCs and synovial biopsies.
To identify changes in pathways that might mediate disease we undertook Gene Ontology (GO) analysis. In the synovial biopsies, a number of inflammatory pathways showed altered expression including those involving oxidoreductase activity (including the cyclooxygenases which mediate prostaglandin production), B-cell activity, interferon-γ response and myeloid cell activation (Table 1).
We also specifically focused on gene expression changes that might contribute directly to the tissue remodelling seen in affected joints in SpA. The tissue remodelling inflammatory genes, matrix metalloproteinase 1 (MMP1, 3.5-fold, p = 0.001) and matrix metalloproteinase 3 (MMP3, 4.7-fold, p = 0.005) showed marked up-regulation in AS/ SpA biopsies (Table 2). Quantitative PCR confirmed these changes showing an 11-fold upregulation in MMP-3 expression ( Figure 1C). Robust MMP-3 protein expression was detected by immunohistochemistry in AS biopsies (6700 OD/mm2) with low expression in SpA (32 OD/ mm2) and RA (652 OD/mm2) samples. MMP-3 protein expression was not detected in normal control samples ( Figure 2). MMP-3 RNA levels were also higher in the two AS samples than in the SpA samples, though not significantly. The prostaglandin E receptor 4 (PTGER4) was also upregulated (1.24-fold by microarray, 1.9-fold by qPCR, p < 0.05). Gene ontology analysis identified matrix catabolic and metabolic pathway dysregulation (Figure 1).
Two Wnt pathway inhibitory genes were down-regulated in our microarray dataset ( Figure 1C, Table 2); DKK3 (2.7-fold, p = 0.003) and Kremen1 (1.5-fold, p = 0.007). Quantitative PCR data supported the array findings with DKK3 down-regulated 2.7-fold (p = 0.09, Figure 1C, Table 2); DKK3 was in fact undetectable in the AS samples with low levels of expression in the SpA samples (data not shown).
A recent study demonstrated a strong enhancement of a "myogene signature" in AS and SpA synovial biopsies [10]. We also saw alterations in a number of myocyteassociated pathways (Table 1). However when we looked specifically at the genes differentially expressed in the myogene signature in the Yeremenko study we did see not strong expression changes suggesting our myogene signature was due to a different subset of genes (data not shown).

Discussion and conclusions
Using whole genome expression profiling in archived synovial biopsies we have established changes in key pathways and genes that might mediate both the inflammatory changes and the tissue remodelling downstream of the inflammation in SpA and AS.
Estimates of the incidence of peripheral arthritis are between 20-50% in AS and SpA patients [8,15,16]. It has been proposed that the aetiopathogenesis of peripheral and axial SpA are similar [8,17]. In both cases inflammation arises close to the enthesis with the inflammatory infiltrate sharing many common features at the two sites [17]. Whether enthesitis is the underlying initiating pathology driving disease in SpA is still a subject of considerable debate [18,19].
As might be expected in inflammatory arthritidies such as SpA and AS, immune pathways are affected. Comparison of this synovial tissue dataset with our previously published PBMC dataset [3] identified a subset of inflammatory genes and pathways that were altered in both studies. Similar dysregulation in the interferon response and myeloid cell pathways was seen possibly reflecting systemic changes. Localised tissue inflammatory pathways such as the oxidoreductase pathways however are altered in synovial tissue but not PBMCs. Differentially regulated pathways potentially mediate the progression from systemic inflammation to localised inflammatory-driven tissue damage.
In synovium, a number of closely-associated inflammatory pathways involved in oxidoreductase activity were identified, which includes the monooxygenase pathways involved in nitric oxide production, and the cyclooxygenase pathways producing COX-1 and COX-2 produce prostanoids such as prostaglandins. COX-2 expression has previously been demonstrated in SpA-affected joints [20]. Cyclo-oxygenase inhibition using non-steroidal anti-inflammatory drugs is a mainstay of therapy in AS, and there is even suggestive evidence that such treatment may retard the progression of ankylosis in the disease [21][22][23]. Prostaglandin E receptor 4 (PTGER4) was also upregulated, which has been shown to be associated with AS in genomewide association studies Gene ontology (GO) analysis was undertaken using BRB Array Tools. Significantly altered genesets were identified using Efron-Tibshirani's Gene Set Analysis (GSA) test which tests which genesets contain more differentially expressed genes than would be expected by chance. The threshold of determining significant gene sets was 0.005.  [24]. This is of particular interest, as PTGER4 through its ligand PGE 2 is a good molecular candidate to link physical stress at entheses with bone formation [25], and in driving inflammation through stimulation of IL-23 production by dendritic cells [26]. Further alterations at the tissue level were seen in pathways affecting collagen metabolism and catabolism, cellmotility and extracellular matrix interactions reflecting the inflammatory joint destruction and tissue remodelling seen in SpA. These were not altered in our studies on whole blood and PBMCs [2,3].
MMP-3 was one of the most strongly upregulated genes. Members of the MMP family of stromelysins have been well documented to play roles in inflammationmediated tissue destruction. Elevated serum levels of MMP-3 have been indicated in AS as a systemic biomarker of disease progression and activity [27], and correlate well with BASDAI [28] and response to TNFblockade treatment [29,30]. In a study on SpA patients with peripheral joint involvement, high serum MMP-3 correlated closely with increased synovial fibroblast MMP-3 production supporting a local joint source for the serum levels. MMP3 levels have been suggested to be the best predictor of peripheral arthritis treatment response [31]. In fact high MMP3 production was proposed as a diagnostic biomarker for peripheral involvement rather than global inflammation in SpA. High serum MMP3 levels (presumably originating from the synovitis) differentiated those patients suffering from axial and peripheral SpA from those with only axial SpA [31]. Even though synovial inflammation in RA is usually more destructive than that in SpA, MMP3 levels are still higher in SpA suggesting a different tissue remodelling role for MMP3 in SpA.
The Wnt pathway has been identified as playing an important role in mediating bone formation and release of inhibition of this pathway has been suggested to contribute to osteoproliferation both in AS [32] and in mouse models of SpA [33]. Downregulation of Wnt inhibitors, such as DKK3 and Kremen1, as suggested by the current data, could therefore generate permissive signals for the excess bone formation seen in AS. Osteoproliferation/bone formation in the synovial joints of SpA patients has not been described however, though bone formation in the affected entheses of SpA patients has been demonstrated [18,19].
In a similar study to this one, Yerenmenko et al. undertook a large scale whole genome expression profiling study comparing SpA with RA and gout synovial biopsies rather than OA and normal samples [10]. The key finding from this study was the identification of a 296-gene "myogene" expression profile that was highly enriched for genes associated with muscle/myocyte/ myofibroblast biology. Interestingly, they did not report strong upregulation of inflammatory genes possibly due to the comparison being between two inflammatory arthritidies, although MMP1 was upregulated in the SpA samples. They also reported altered expression of genes in the Wnt pathway.
Similarly we also saw changes in "myogene" associated pathways, further supporting their proposal for fibrotic changes in the synovium of SpA patients. The specific gene changes underlying these pathways were not the same in the two studies but this may reflect the different patient cohorts and tissue processing (FFPE vs. fresh frozen). Analysis of our previous expression profiling studies in PBMCs and whole blood showed the absence of a myogene signature in these datasets suggesting it is a disease-site specific phenomenon [2,3]. Interestingly, gene ontology analysis comparing expression profiling of spines and knees in proteoglycan induced spondylitis (PGISp) mice showed a greater number of muscle-associated pathways upregulated in the knee joints suggesting this may be a unique feature of peripheral disease [33]. The significance of the myogene profile though remains to be elucidated however.
Two samples (1 OS-Normal and 1 AS-SpA) misclassified during the clustering analysis. There were no technical issues identified that might underline this so we can assume the reasons were biological. The misclassification of the sample probably reflects the compounded biological variation in SpA patients due to a combination of genetic factors and disease heterogeneity reflecting onset, severity and symptoms.
Although we identified some key pathways and genes of interest in this study it must e regarded as an exploratory study at this time. Despite some of the findings agreeing with previous studies [10], further independent validation studies are required to confirm the significance of our initial findings.
By adopting a whole genome profiling approach this study has identified gene signatures differentiating SpA from non-SpA samples and highlighting pathways that might play key pathophysiological roles in AS. Further, the candidate gene changes we have highlighted possible disease pathways that might control the progression through the inflammation and tissue destructive/osteoproliferative phases of spondyloarthropathy and provide guidance for focusing research efforts to elucidate disease mechanisms.

Additional files
Additional file 1: Table S1. Clinical data for patients and controls.
Additional file 2: Table S2. Class comparison of the Norm-OA and AS-SpA groups identified 416 differentially expressed genes (p < 0.01) ranging from a 4.7-fold up-regulation to a 4.6-fold down-regulation.