This article is mentioned in:

Introduction

Rheumatoid arthritis (RA) is a chronic inflammatory disease, characterized by synovial cell proliferation and excessive production of pro-inflammatory cytokines, leading to the destruction of diarthrodial joints cartilage and bones (1,2). RA is considered an autoimmune disease (3) that can cause severe disability and premature mortality (4), with ~1% of individuals afflicted worldwide (5). Although the etiology of RA has not been fully elucidated, numerous studies have demonstrated that it is a multifactorial disease that results from interactions between genetic and environmental factors (6).

Over the past 15 years, rheumatologists have developed a number of therapeutic strategies aiming to treat RA (7). Moreover, researchers have shown that the characterization of the cytokine signaling pathways involved in RA provides important opportunities for identifying pro-inflammatory cytokines that can be targeted in the context of novel therapeutic interventions (8–10). However, these agents are not effective in all patients (11). Furthermore, it has long been recognized that patients with RA have an increased risk for certain types of cancer (12,13). For example, it was demonstrated that RA-associated interstitial lung disease accounts for mortality of up to 20% of RA patients (14). Therefore, there is an urgent need to investigate the molecular mechanisms underlying the processes involved in RA.

High-throughput mRNA sequencing (RNA-seq), which allows simultaneous identification of transcripts and estimation of their abundance (15,16), is a recently developed transcriptome profiling approach that uses deep-sequencing technologies (17). It has fostered numerous advances in the characterization and quantification of transcripts, since it allows a nearly complete characterization of transcriptomes (18). A recent study making use of RNA-seq showed that the activation and proliferation of RA synovial fibroblasts relate to the pathogenesis of RA (19). Here, we performed a comprehensive meta-analysis of the transcriptome data from based on the data of Heruth et al (19), which were derived from RNA of healthy and RA patients. Our findings may enhance the understanding of the mechanisms underlying the process of joint destruction, and allow a more selective and specific application of therapeutic agents that target pro-inflammatory cytokines and thus, a more effective treatment of patients with RA and other inflammatory disorders.

Materials and methods

RNA-seq data

The RNA-seq data were downloaded from the Sequence Read Archive (SRA). Two samples of synovial fibroblasts were available for meta-analysis: One was from individuals with RA, and the control sample (CT) was from healthy individuals (Table I). The downloaded data come from Illumina RNA-seq (19) of paired-end (2 × 100 bp) cDNA libraries prepared from each RNA sample, with sequencing performed as in (20,21).

Table I Information on raw RNA-seq data from rheumatoid arthritis (RA) and control (CT) samples, downloaded from the Sequence Read Archive (SRA).

Table I

Information on raw RNA-seq data from rheumatoid arthritis (RA) and control (CT) samples, downloaded from the Sequence Read Archive (SRA).

Sample name	Type	Library	Data size (bp)	SRA no.
RA	Rheumatoid arthritis	Paired-end	9.1 G	SRR364313
CT	Healthy	Paired-end	8.2 G	SRR364315

Quality analysis of raw RNA-seq reads

SolexaQA is a software that allows investigation and trimming of sequences based on their base quality scores. To exclude analytical error from sequencing errors and low data quality, the SolexaQA software (22) was used to process the raw data, based on the read-level quality of each sequence. Raw sequences with quality scores <20 and a length <25 bp were removed from the final dataset that was used for downstream analyses.

Calculation of expression values and identification of differentially expressed genes (DEGs)

First, the TopHat software version 2.0.8 (23) was used to map the sequencing reads to the reference human genome (version hg19, UCSC Genome Browser). TopHat is an efficient read-mapping algorithm designed to align reads from an RNA-seq experiment to a reference genome; we used the following settings: maximum 20 multiple hits per read and maximum 2 mismatches allowed. Second, Cufflinks software package (version 2.0.2) was used to assemble transcripts and calculate their relative expression levels, expressed as fragments per kilobase of exon per million fragments mapped (FPKM) using the default parameters and an average insert size of 200 bp (24). Finally, the Cuffdiff module in Cufflinks was used to estimate differential expression, expressed as a ratio of RA to control expression for each transcript, along with the statistical significance of the observed differences (25). Genes with log2 fold change (FC) ≥1 and p-value ≤0.05 were selected as DEGs. In addition, the expression values of the differentially expressed genes were normalized by z-score transformation prior to visualization with the heatmap function available in the R statistical package (26).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis

The GOstat software (27) was used to conduct GO functional annotation and enrichment analysis of DEGs (significantly at P<0.05). The KAAS annotation server (28) was used to identify the KEGG metabolic pathways in which the identified DEGs are involved. Finally, the KOBAS server 2.0 (29), integrating data from a number of human disease databases, was used to identify DEGs related to RA.

Results

Identification of DEGs

Quality control statistics on the RNA-seq are shown in Table II. A total of 293 genes were identified as differentially expressed between the RA and CT samples based on our criteria (log2 FC ≥1; p-value ≤0.05). Fig. 1 shows the expression patterns of DEGs in the RA and CT samples. The relative expression (FKPM), log2 FC, and associated p-values of the top 10 DEGs (in terms of FKPM at the RA sample) are shown in Table III.

Figure 1 Expression values of genes differentially expressed between rheumatoid arthritis (RA) and healthy control (CT) samples. Values are expressed as fragments per kilobase of exon per million fragments mapped (FPKM). Red represents high expression, black represents median expression, and green represents low expression.

Table II Quality control statistics of mRNA sequence reads.

Table II

Quality control statistics of mRNA sequence reads.

Reads	Healthy control	Percentage	Rheumatoid arthritis	Percentage
Original	80,782,262	-	89,757,726	-
Low-quality	9,605,518	11.89	6,707,374	7.47
Remaining	71,176,744	88.11	83,050,352	92.53

Table III The top 10 significantly differentially expressed genes.

Table III

The top 10 significantly differentially expressed genes.

Gene	Chr	RA FPKM	CT FPKM	Log2 FC	p-value
CHI3L1	Chr1	5,847.720	0.211	−14.753	3.02E-14
MMP-1	Chr11	78.271	0.017	−12.166	9.28E-11
SMOC2	Chr6	36.802	0.010	−11.719	1.95E-07
Ror2	Chr9	8.913	0.009	−9.995	1.23E-05
VIT	Chr2	9.313	0.013	−9.522	1.34E-04
HEYL	Chr1	0.043	28.459	9.355	1.29E-14
MCAM	Chr11	0.155	153.674	9.955	0
EFHD1	Chr2	0.054	103.317	10.893	5.70E-12
FOXE1	Chr9	0.009	18.429	11.069	2.92E-12
WFDC1	Chr16	0.020	76.095	11.882	5.09E-13

[i] Chr, human chromosome; RA, rheumatoid arthritis; FPKM, fragments per kilobase of exon per million fragments mapped; FC, fold change; CHI3L1, chitinase 3-like 1; MMP-1, matrix metalloproteinase-1; SMOC2, SPARC related modular calcium binding 2; Ror2, receptor tyrosine kinase-like orphan receptor 2; VIT, vitrin; HEYL, hes-related family bHLH transcription factor with YRPW motif-like; MCAM, melanoma cell adhesion molecule; EFHD1, EF-hand domain family, member D1; FOXE1, forkhead box E1; WFDC1, WAP four-disulfide core domain 1.

GO analysis

The GOstat tool retrieves GO annotations and investigates over-representation of these in a given gene list (27). We used GOstat to functionally annotate the identified DEGs based on the cut-off point of P<0.05. The 293 DEGs were enriched for 309 GO terms. The top 10 GO terms in terms of p-value are shown in Table IV. This analysis indicated that the identified DEGs are mainly involved in anatomical structure development, cell membrane formation and stability, and biological adhesion.

Table IV The top 10 significantly enriched Gene Ontology (GO) terms among differentially expressed genes.

Table IV

The top 10 significantly enriched Gene Ontology (GO) terms among differentially expressed genes.

GO ID	Name	Description	p-value
0048856	Anatomical structure development	The biological process whose specific outcome is the progression of an anatomical structure from an initial condition to its mature state. This process begins with the formation of the structure and ends with the mature structure, whatever form that may be including its natural destruction. An anatomical structure is any biological entity that occupies space and is distinguished from its surroundings. Anatomical structures can be macroscopic such as a carpel, or microscopic such as an acrosome.	8.76e-44
0032502	Developmental process	A biological process whose specific outcome is the progression of an integrated living unit: An anatomical structure (which may be a subcellular structure, cell, tissue, or organ), or organism over time from an initial condition to a later condition.	1.19e-36
0007275	Multicellular organismal development	The biological process whose specific outcome is the progression of a multicellular organism over time from an initial condition (e.g. a zygote or a young adult) to a later condition (e.g. a multicellular animal or an aged adult)	1.36e-35
0048731	System development	The process whose specific outcome is the progression of an organismal system over time, from its formation to the mature structure. A system is a regularly interacting or interdependent group of organs or tissues that work together to carry out a given biological process.	2.12e-33
0032501	Multicellular organismal process	Any biological process, occurring at the level of a multicellular organism, pertinent to its function.	3.89e-32
0009653	Anatomical structure morphogenesis	The process in which anatomical structures are generated and organized. Morphogenesis pertains to the creation of form.	1.62e-28
0005886	Plasma membrane	The membrane surrounding a cell that separates the cell from its external environment. It consists of a phospholipid bilayer and associated proteins.	2.66e-28
0031226	Intrinsic component of plasma membrane	The component of the plasma membrane consisting of gene products and protein complexes that have some covalently attached part (e.g. peptide sequence or GPI anchor), which spans or is embedded in one or both leaflets the membrane.	8.71e-25
0005887	Integral component of plasma membrane	The component of the plasma membrane consisting of gene products and protein complexes that have some part that penetrates at least one leaflet of the membrane bilayer. This component includes gene products that are buried in the bilayer with no exposure outside the bilayer.	8.01e-24
0022610	Biological adhesion	The attachment of a cell or organism to a substrate or other organism.	1.01e-21

KEGG pathway analysis

The 293 DEGs were found to be involved in 131 pathways, including Wnt signaling, cell adhesion molecules (CAMs), neuroactive ligand-receptor interaction, PI3K-Akt signaling, cytokine-cytokine receptors interactions, calcium signaling, regulation of actin cytoskeleton and focal adhesion (Table V). By combining the functional annotation of DEGs with data from disease databases, we found that 16 genes among the DEGs have been previously associated with the occurrence of RA (Table VI).

Table V KEGG pathways that DEGs are mainly involved in.

Table V

KEGG pathways that DEGs are mainly involved in.

Pathway id	Description	Count
4080	Neuroactive ligand-receptor interaction	10
4151	PI3K-Akt signaling pathway	9
4060	Cytokine-cytokine receptor interaction	9
4020	Calcium signaling pathway	7
4512	ECM-receptor interaction	7
4514	Cell adhesion molecules	7
4810	Regulation of actin cytoskeleton	7
4510	Focal adhesion	7
4974	Protein digestion and absorption	7
5410	Hypertrophic cardiomyopathy	7

[i] KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes; PI3K, phosphoinositide-3-kinase; ECM, extracellular matrix.

Table VI Genbank information on differentially expressed genes (n=16) that have been associated with rheumatoid arthritis.

Table VI

Genbank information on differentially expressed genes (n=16) that have been associated with rheumatoid arthritis.

Acc. no.	Gene symbol	Name	Expression
NM_003881	WISP2	WNT1 inducible signaling pathway protein 2	Down
NM_001066	TNFRSF1B	Tumor necrosis factor receptor superfamily, member 1B	Down
NM_000877	IL1R1	Interleukin 1 receptor, type 1	Down
NM_002996	CX3CL1	Chemokine (C-X3-C motif) ligand 1	Up
NM_007365	PADI2	Peptidyl arginine deiminase, type II	Up
NM_001276	CHI3L1	Chitinase 3-like 1	Down
NM_004864	GDF15	Growth differentiation factor 15	Down
NM_000214	JAG1	Jagged 1	Up
NM_000612	IGF2	Insulin-like growth factor 2	Up
NM_004878	PTGES	Prostaglandin E synthase	Down
NM_001145938	MMP-1	Matrix metalloproteinase-1	Down
NM_019111	HLA-DRA	Major histocompatibility complex, class II, DR α	Up
NM_000396	CTSK	Cathepsin K	Down
NM_005118	TNFSF15	Tumor necrosis factor superfamily, member 15 transcript variant 1	Up
NM_006379	SEMA3C	Sema domain, Ig, short basic domain, secreted, semaphorin 3C	Down
NM_000692	ALDH1B1	Aldehyde dehydrogenase 1 family, member B1	Up

Discussion

RA is a systemic inflammatory disorder that commonly affects the diarthrodial joints (30). The pathogenesis of RA is characterized by the influx of immune system cells (31), which induce the production of pro-inflammatory cytokines, decreased synthesis of anti-inflammatory cytokines and the subsequent activation and proliferation of synovial fibroblasts (32). Remission in RA is an increasingly attainable goal, but there is no widely used definition of remission that is stringent yet achievable, and that could be uniformly applied as a criterion of clinical outcome (33). Therefore, the identification and characterization of genes related to RA is important for the understanding of the pathogenesis of this disease and the identification of novel anti-inflammatory therapeutic targets.

In the present study, we identified genes differentially expressed between RA and healthy individuals. A total of 293 genes were identified as DEGs, including these encoding the receptor tyrosine kinase-like orphan receptor 2 (Ror2), chitinase-3-like 1 (CHI3L1), matrix metalloproteinases (MMPs), interleukin (IL)-26, and v-maf musculoaponeurotic fibrosarcoma oncogene homolog B (MafB). Ror2 has been associated with RA. A previous study indicated that the Wnt5a-Ror2 pathway is crucial for osteoclastogenesis in physiological and pathological environments, and may represent a therapeutic target for bone diseases, including RA (34). Sonomoto et al (35) demonstrated that IL-1β can effectively and rapidly induce the differentiation of human mesenchymal stem cells into osteoblasts, as well as mineralization, mainly through the non-canonical Wnt5a-Ror2 pathway.

The CHI3L1 or YKL40 gene encodes for the human cartilage glycoprotein 39 (HC-gp39), which is secreted by synovial fibroblasts, macrophages, neutrophil granulocytes, and chondrocytes. Its expression is regulated by NF-κB (36). It was recently suggested that the YKL-40 protein might be implicated in the pathogenesis of RA and that its level may indicate the degree of joint inflammation (37). The carbohydrate-binding motif in YKL-40 specifically activates the Akt signaling pathway in colonic epithelial cells (38). The level of YKL-40 in the serum varies depending on the RA status of patients, and is increased in 54% of patients with clinically active disease (39). YKL-40 may thus be suitable for assessing the disease activity and pathophysiology of RA.

MMPs are members of an enzyme family that contain a zinc ion on their active site, which is required for their catalytic activity. MMPs are critical for maintaining tissue allostasis. These enzymes are active at neutral pH, and can therefore catalyze the physiological turnover of extracellular matrix (ECM) macromolecules (40). A previous study showed that the levels of MMP-8 and -9 in the systemic circulation are representative of the levels of these enzymes in the inflamed joint, and suggested that MMP-9 may be involved in degradation of the joint collagen. The study further confirmed the hypothesis of an MMP/TIMP imbalance in RA (41). In addition, inhibition of MMP-13 was shown to reduce cartilage erosions in two out of three tested animal models of RA, strongly supporting the development of drugs targeting MMPs in order to reduce or halt joint destruction in patients with RA (42).

IL-26, a member of the IL-10 cytokine family, induces the production of pro-inflammatory cytokines by epithelial cells. It was recently demonstrated that IL-26 is constitutively produced by fibroblast-like cells known as synoviocytes in RA, and that it induces the secretion of pro-inflammatory cytokines by myeloid cells and favors the formation of T helper cells producing interleukin 17 (Th17). The authors suggested that IL-26 is a pro-inflammatory cytokine located upstream of the pro-inflammatory cascade, and may constitute a promising target to treat RA and chronic inflammatory disorders (43). Although synoviocytes are present in all joints, only synoviocytes from RA patients can produce IL-26 (44).

MafB was another gene found as differentially expressed in our study. The MafB protein is a putative tumor suppressor in the myeloid lineage, with a key role in monocytopoiesis (45), as well as in monocyte-dendritic cell differentiation (46).

The GO functional annotation analysis showed that DEGs are enriched for a total of 309 GO terms. The top 10 GO terms mainly referred to anatomical structure development, membrane formation and stability, and biological adhesion. It is notable that alterations in biological adhesion are involved in RA. For example, COL5A1 (collagen, type V, α1) was expressed at significantly higher levels in chondrocytes from the damaged region of osteoarthritic cartilage than in those from the intact region (47).

KEGG signaling pathway analysis revealed that the DEGs are predicted to be involved in a total of 131 pathways, including Wnt and calcium signaling, as well as CAM-related pathways. The Wnt signaling pathway plays a key role in cell renewal. A few studies showed that it is also involved in RA pathogenesis (48,49). In the synovial membrane of patients with RA, the Wnt and Fz genes are expressed at higher levels compared to those observed in patients without RA (50). The Wnt proteins are glycoproteins that bind to the Fz receptors on the cell surface, thereby affecting a number of important biological processes, such as cell differentiation, embryonic development, limb development, and joint formation (51). Enhanced knowledge of the role(s) of the Wnt signaling pathway in RA is expected to improve our understanding of the different RA clinical features and improve prognosis. Both calcium signaling and CAM-related pathways, upregulated in RA patients showing a traditional Chinese medicine heat pattern, have been suggested to be important for T-lymphocyte interactions, and thus, constitute candidate targets for RA therapy (52).

Cross-referencing with disease databases revealed that 16 of the identified DEGs have been previously associated with RA; these include genes encoding MMP-1, interleukin-1 receptor type 1 (IL1R1) and PADI2. The aggressive phenotype of synovial fibroblasts in RA is characterised by the increased expression of MMP-1 (53). Collagenases MMP-1 and -13 play an important role in collagen degradation in RA and osteoarthritis, while gelatinases MMP-2 and -9 may be involved in arthritis by degrading non-collagen matrix components in the joints (54). The IL1R1 receptor, once activated upon IL-1 binding, activates NF-κB, which is a modulator of expression of inflammatory and immune genes (55). Nevertheless, further experiments are needed to study the roles that these genes may play in the occurrence of RA.

Overall, our study provided a list of candidate molecular targets for RA treatment and further research. In recent years, gene-targeted therapy of RA has become a popular approach, and related ongoing research appears promising. Despite significant therapeutic advances, RA treatment remains an unsolved medical issue. Therefore, the molecular mechanism(s) underlying this disease need to be further investigated.

Acknowledgements

This study was supported by the National Natural Science Foundation of China of the year 2011 (81171734) and the Applied Basic Research Project of Yunnan Province (2011FZ313).

References