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 information 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.

Method

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 gene-sets. 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).

Briefly, cDNA was generated from 1 μg of total RNA using the Bioline cDNA synthesis Kit (Bioline, London, UK) according to manufacturer’s instructions. Candidate genes were assayed using the predesigned TaqMan assays (DKK3 = Hs00247426_m1, PTGER4 = Hs00168761_m1, MMP3 = Hs00968308_m1). For normalisation, expression levels of the housekeeping gene, RPL32, [8] were measured by SYBR green based qRT-PCR using specific forward (5′-CCCCTTGTGAAGCCCAAGA-3′) and reverse (5′-GACTGGTGCCGGATGAACTT-3′) primers. All assays were carried out using SensiMix dT RT-PCR reagent (Quantace, Sydney, Australia) under the following conditions; 50°C for 2 min, 95°C for 10 min, and 40 cycles of 95°C for 15 s and 60°C for 60s.

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%).

Figure 1

Histological and immunohistochemical findings. (A) Diffuse lymphocytic infiltration of the dermis predominantly composed of large pleomorphic cells with irregular nuclei and prominent nucleoli (hematoxylin-eosin staining, magnification 20×). By immunoperoxidase staining, the neoplastic cells showed marked and diffuse expression of CD3 (B) and CD30 (C).

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 lectin-like 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).

Table 1

Gene ontology analysis of differentially expressed genes

Gene ontology term	GO category	Efron-Tibshirani’s GSA test p-value
Collagen/extracellular matrix association
Collagen binding	GO:0005518	< 0.005
Collagen metabolic process	GO:0032963	0.05
Collagen catabolic process	GO:0030574	0.035
Integrin complex	GO:0008305	< 0.005
Substrate-bound cell migration	GO:0006929	< 0.005
Negative regulation of cytoskeleton organization	GO:0051494	0.025
Nitric oxide/monooxygenase activity
Regulation of monooxygenase activity	GO:0032768	< 0.005
Regulation of oxidoreductase activity	GO:0051341	< 0.005
Regulation of nitric-oxide synthase activity	GO:0050999	0.005
Nitric oxide biosynthetic process	GO:0006809	0.015
Positive regulation of monooxygenase activity	GO:0032770	0.005
Nitric oxide metabolic process	GO:0046209	0.01
Regulation of nitric oxide biosynthetic process	GO:0045428	0.005
Regulation of calcidiol 1-monooxygenase activity	GO:0060558	0.005
Immune associated functions
B cell receptor signaling pathway	GO:0050853	0.015
Regulation of leukocyte mediated immunity	GO:0002703	0.02
Response to interferon-gamma	GO:0034341	0.005
Positive regulation of B cell activation	GO:0050871	0.04
Humoral immune response mediated by circulating immunoglobulin	GO:0002455	0.02
Myeloid cell activation during immune response	GO:0002275	0.01
Regulation of lymphocyte mediated immunity	GO:0002706	0.035
Muscle/myocyte/myofibroblast biology
Dystroglycan binding	GO:0002162	0.015
Myosin complex	GO:0016459	0.01
Muscle filament sliding	GO:0030049	0.02
Actin-myosin filament sliding	GO:0033275	0.02
Myosin II complex	GO:0016460	0.025
Actomyosin structure organization	GO:0031032	0.04
Regulation of myoblast differentiation	GO:0045661	0.03
Muscle myosin complex	GO:0005859	0.02
Myoblast differentiation	GO:0045445	0.02
Regulation of myotube differentiation	GO:0010830	0.015
Myotube differentiation	GO:0014902	0.005

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.

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).

Table 2

Expression levels of candidate genes on the microarray and in the qPCR confirmation study

Array data	Norm-OA (9)		AS-SpA (6)		AS-SpA/Norm-OA	p-value
MMP3	480.66		2260.32		4.7	0.005255
MMP1	199.47		701.47		3.51	0.001053
DKK3	517.49		190.98		-2.71	0.003263
KREMEN1	988.68		674.73		-1.46	6.50E-03
PTGER4	317.58		393.08		1.24	4.70E-02
qPCR	Norm-OA (9)	SD	AS-SpA (7)	SD	AS-SpA/Norm-OA	p-value
MMP3	0.00026	0.00058	0.00295	0.00412	11.35	0.07
DKK3	0.101	0.073	0.037	0.061	-2.7	0.086
PTGER4	66.6	62.9	126.7	50.1	1.9	0.042

Figure 2

Immunohistochemistry showing MMP3 protein levels in synovial biopsies. MMP3 was undetectable in normal biopsies (A) with very low levels seen in RA (C) with similarly low levels in SpA (D) samples and high expression of MMP3 in AS samples (B). Quantitation of the stain is shown in (E).

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 myocyte-associated 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-23]. Prostaglandin E receptor 4 (PTGER4) was also upregulated, which has been shown to be associated with AS in genomewide association studies [24]. This is of particular interest, as PTGER4 through its ligand PGE2 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, cell-motility 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 inflammation-mediated 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 TNF-blockade 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.

Abbreviations

AS: Ankylosing spondylitis; PGISp: Proteoglycan-induced spondylitis; IHC: Immunohistochemistry; RA: Rheumatoid arthritis; IVD: Intervertebral disc; MRI: Magnetic resonance imaging; TNF: Tumour necrosis factor; HLA: Human leukocyte antigen; LRP: Low-density lipoprotein receptor-related protein; GSK: Glycogen synthase kinase; H&E: Haematoxylin and eosin; EDTA: Ethylenediaminetetraacetic acid; PG: Proteoglycan; RNA: Ribonucleic acid; MMP: Matrix metalloproteinase.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

GPT conceived and designed study, analysed data and wrote manuscript. RD generated and analysed data and edited manuscript. ARP provided technical expertise, analysed data and wrote manuscript. HW generated and analysed data and edited manuscript. SK generated and analysed data and edited manuscript. MS provided clinical samples, expert opinion and edited the manuscript. MAB conceived study and wrote manuscript. All authors read and approved the final manuscript.

Pre-publication history

The pre-publication history for this paper can be accessed here:

http://www.biomedcentral.com/1471-2474/14/354/prepub

Additional file 1: Table S1

Clinical data for patients and controls.

Click here for file

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.

Click here for file