- Open Access
RNA-seq analysis of synovial fibroblasts brings new insights into rheumatoid arthritis
Cell & Biosciencevolume 2, Article number: 43 (2012)
Rheumatoid arthritis (RA) is a chronic autoimmune-disease of unknown origin that primarily affects the joints and ultimately leads to their destruction. Growing evidence suggests that synvovial fibroblasts play important roles in the initiation and the perpetuation of RA but underlying molecular mechanisms are not understood fully. In the present study, Illumina RNA sequencing was used to profile two human normal control and two rheumatoid arthritis synvovial fibroblasts (RASFs) transcriptomes to gain insights into the roles of synvovial fibroblasts in RA.
We found that besides known inflammatory and immune responses, other novel dysregulated networks and pathways such as Cell Morphology, Cell-To-Cell Signaling and Interaction, Cellular Movement, Cellular Growth and Proliferation, and Cellular Development, may all contribute to the pathogenesis of RA. Our study identified several new genes and isoforms not previously associated with rheumatoid arthritis. 122 genes were up-regulated and 155 genes were down-regulated by at least two-fold in RASFs compared to controls. Of note, 343 known isoforms and 561 novel isoforms were up-regulated and 262 known isoforms and 520 novel isoforms were down-regulated by at least two-fold. The magnitude of difference and the number of differentially expressed known and novel gene isoforms were not detected previously by DNA microarray.
Since the activation and proliferation of RASFs has been implicated in the pathogenesis of rheumatoid arthritis, further in-depth follow-up analysis of the transcriptional regulation reported in this study may shed light on molecular pathogenic mechanisms underlying synovial fibroblasts in arthritis and provide new leads of potential therapeutic targets.
Rheumatoid arthritis [RA] is a chronic, systemic autoimmune disorder associated with both genetic and environmental factors. RA affects 1% of the world’s population, develops most commonly in adults between 40 – 70 years old, and occurs more frequently in women than in men [1–4]. xAlthough the etiology of the disease has not been elucidated fully, the pathogenesis of RA is characterized by the influx of cells from both the innate and the adaptive immune systems . These cells induce increased pro-inflammatory cytokine production, decreased synthesis of anti-inflammatory cytokines, and the subsequent activation and proliferation of synovial fibroblasts (SFs) [3, 4]. Rheumatoid arthritis synovial fibroblasts (RASFs) produce additional cytokines, chemokines and matrix-degrading enzymes which ultimately leads to the thickening and progressive destruction of joint membrane, cartilage and bone [5–7]. Characterization of the cytokine signaling pathways involved in RA has provided a significant opportunity for identifying pro-inflammatory cytokines which can be targeted for novel therapeutic intervention. The development of biological response modifiers (BRMs), particularly the TNF, IL-1, and IL-6 antagonists, have led to major advances in RA therapy [3, 7]. However, these agents are not effective in all patients, underscoring the genetic heterogeneity of the disease and the need for the development of additional BRMs . RASFs are intricately involved in the pathogenesis of RA and provide a source for the identification of new genes and pathways that can be targeted for therapeutic intervention.
With the advent of next generation DNA sequencing technologies , such as RNA sequencing (RNA-seq), a more comprehensive and accurate transcriptome analysis has become feasible and affordable. In RNA-seq, short fragments of complementary DNA (cDNA) are sequenced (reads) and then mapped onto the reference genome. RNA-seq enables not only the identification of differentially expressed genes, but also the precise quantitative determination of exon and isoform (alternative splicing) expression, along with the characterization of transcription initiation sites (TSSs) and new splicing variants . In the present study, we performed a comprehensive transcriptome analysis of RNA from RASFs from two adult female RA patients and the SF RNA from two healthy female donors, using the RNA-seq technique. We found significant differences in the expression levels of both genes and gene isoforms between normal SFs and RASFs RNA samples. These data provide broader and deeper insights, particularly with respect to isoform expression, into the effect of RA on the transcriptional regulation of synovial fibroblasts and a rich resource for further experimentation into the pathogenesis of the disease.
Human SFs RNAs from two healthy control donors and two patients with RA were purchased from Cell Applications, Inc. (San Diego, CA). Diseased samples were age and sex-matched with normal controls (Additional file 1). Paired-end cDNA libraries for each RNA sample were prepared and sequenced using the Illumina TruSeq RNA Sample Preparation Kit, as outlined previously [11, 12].
Quality analysis of RNA-seq data
Real-time analysis of the sequencing run was performed by the Illumina HiSeq Control Software. Clusters of identical sequences were generated on the Illumina cBot and the number of those clusters was reported, along with the percentage of those clusters passing an internal quality filter. Across the 4 samples, between 433,000 and 482,000 raw clusters were detected, with a median of 446,000 clusters per lane. Between 90.9% and 95.0% of those clusters passed the filter, with a median of 93.2% of the clusters passing the filter. Each lane was aligned in real-time with the phiX genome and between 0.80% and 0.84% of the clusters aligned, with a median of 0.81% aligned. Our control lane of phiX produced 290,000 clusters with 97.9% passing the filter and 99.08% aligning to the phiX genome. All these values were within the recommended limits established by Illumina.
Post-run quality analysis of RNA-seq data was carried out as described by Twine et al. . The total number of reads produced from each sample was between 80,782,262 and 89,757,726, with a mean across all samples of 84,177,268 (Table 1). The difference in the number of reads between the control samples and the RA samples was not statistically significant (Student’s t-test, p=0.27). To assess the quality of the reads, data was pulled from the TopHat log files as well as the output files. Between 0.10% and 0.15% of the reads were removed due to low quality before mapping to the reference genome began. Between 82.8% and 89.1% of the total reads mapped to the human genome. To ensure the uniform coverage across the genome, the data was visualized using a local copy of the Integrative Genomics Viewer. An example of the reads for both normal and RA patient samples mapped against chromosome 1 is shown in Figure 1. The average alignment was computed across the genome and those alignment scores were log-transformed (base 2) to better visualize the full range of the data. As expected, no reads mapped to the centromere or areas of the chromosome without genes.
Differentially expressed genes and isoforms
After mapping the sequencing reads to the reference genome with TopHat, transcripts were assembled and their relative expression levels were calculated with Cufflinks in Fragments Per Kilobase of exon per Million fragments mapped (FPKM). The sub-program, Cuffdiff was then used to calculate the differential expression on the gene and transcript level, as well as the calculation of alternative promoter usage and alternative splicing. Cufflinks calculates the differential gene expression with the ratio of the RA group to the control group for every gene and transcript along with the statistical significance of the values. Two categories of differential gene/isoform expression were identified. The first category consists of genes/isoforms expressed only in control SFs or only in RASFs. The second category consists of genes/isoforms in which expression of both samples in each group was up-regulated or down-regulated two-fold or greater between control SFs and RASFs.
Overall, there are 12,977 expressed genes in the control SFs and 13,445 expressed genes in the RASFs, which were aligned to the reference genome (Table 2). There are 214 genes, whose expressions were only detected in the normal SFs, while 682 genes whose expressions were only detected in RASFs. There are 122 up-regulated and 155 down-regulated genes in RASFs with at least two-fold change compared to the SFs (Table 2). As for known isoforms, there are 20,647 in the normal SFs and 21,102 in RASFs. Among them, there are 526 known isoforms, whose expressions were detected only in the normal SFs, while 981 known isoforms whose expressions were detected only in RASFs. There are 343 up-regulated and 262 down-regulated known isoforms in RASFs by at least two-fold change compared to the SFs (Table 2). For novel isoforms whose annotations are not known in the current reference gene or transcript database, there are 42,124 expressed in the normal SFs and 42,171 expressed in RASFs. Among them, there are 105 novel isoforms whose expressions were only detected in the normal SFs, while 152 novel isoforms were only detected in RASFs. There are 561 up-regulated and 520 down-regulated novel isoforms in RASFs by at least two-fold change compared to the SFs (Table 2).
Genes expressed only in control SFs or only in RASFs
The top 10 up- and down-regulated genes expressed only in control SFs or only in RASFs are presented in Table 3. An expanded list of the top 50 genes expressed only in either control SFs or in RASFs is presented in Additional file 2. Analysis of the genes expressed only in RASF reveals that nine of the top ten genes, including the major histocompatibility complex (MHC) genes HLA-A. –B, -C, and –E, are located on chromosome 6 (Table 3). Remarkably, 36 of the top 50 genes (Additional file 2) expressed only in RASFs are located on chromosome 6. The MHC, particularly the HLA-DRB1 alleles are strongly associated with RA [14–16]. A recent study by Plenge et al. has also identified associations of alleles lying outside the MHC on chromsome 6 with RA . Our observation that the CLIC1 gene (chloride intracellular protein) is expressed in RASFs correlates with the finding that CLIC1(-/-) mice were protected from development of serum transfer induced K/BxN arthritis . Two genes, the high mobility group box 1 (HMGA1) and the latent transforming growth factor beta binding protein 1 (LTBP1) have been reported to be elevated in RA [19, 20] however they are not expressed in the RASFs examined in this study (Table 3). Interestingly, HMGA1 is the only gene on chromosome 6 in the list of top 50 genes expressed in normal SFs but not expressed in RASFs (Additional file 2). The CD59 complement regulatory protein (CD59) is not expressed in RASFs in this study. This observation supports the finding that CD59 is protective as CD59 (-/-) knockout mice present with more severe symptoms in the murine antigen-induced arthritis model . An automated literature search using PubMatrix  reveals that eleven of the twenty genes listed in Table 3 have not yet been identified to be associated with RA (Additional file 3). These genes, which include chromosome 6 open reading frame 48 (C6orf48), the scavenger receptor class A, member 5 (SCARA5), CutA divalent cation tolerance homolog (CUTA), Leucine rich repeat containing 59 (LRRC59), and the protein phosphatase 1, regulatory (inhibitor) subunit 14A (PPP1R14A), may provide additional therapeutic targets. These potential targets include characterized genes, like the iron receptor SCARA5  and genes, such as C6orf48, that have not yet been well-studied. CutA, which is up-regulated in RASFs, interacts with BACE1 to regulate B-cleavage of the B-amyloid protein (APP) . CutA may play a role in the pathogenesis of Alzheimer’s, however, its role in rheumatoid arthritis remains to be elucidated. LRRC59 is required for the nuclear transport of the fibroblast growth factor 1 (FGF1) . The affect on FGF1 function resulting from decreased LRRC59 expression in RASFs warrants further investigation. PPP1R14A, which inhibits protein phosphatase 1 activity, is not expressed in RASFs compared to normal SFs, suggesting that PP1 activity will increase dramatically in RASFs. PP1 controls the Akt signal transduction pathway to regulate cell growth, cell survival, and cell differentiation .
Genes differentially expressed two-fold or greater between control SFs and RASFs
The top 10 up- and down-regulated genes, along with the expanded top 50 list, in which expression of both samples in each group was up-regulated or down-regulated two-fold or greater between control SFs and RASFs are presented in Table 4 and in Additional file 4, respectively. Three genes in the top 10 up-regulated list have been associated with rheumatoid arthritis (Additional file 3). Interleukin 26 (IL26) is up-regulated (80.8-fold) in RASFs compared to SFs. Corvaisier et al. has demonstrated that IL26 is over-expressed in arthritis and induces inflammatory cytokine production . The v-maf musculoaponeurotic fibrosarcoma oncogene homolog B (avian) (MAFB) gene is up-regulated (16.2-fold) in RASFs. Liu et al. identified polymorphisms in the MAFB gene associated with altered response to anti-TNF treatment in patients with RA . Expression of the adrenergic, alpha-2A-, receptor (ADRA2A) increased 14.4-fold in RASFs. The adrenergic, alpha-2A-, receptor may play a critical role in the proliferation and differentiation of synoviocytes . Although thrombospondin 4 (THSB4) has not yet been associated with arthritis (Additional file 3), thrombospondin 1 (THBS1) is over-expressed in RA tissue . Thrombospondin 1 and 4 are extracellular matrix remodeling proteins that have been associated with increased inflammation in coronary artery disease (CAD) [30, 31], and thus may provide a link between RA and CAD. Like THBS4, the remaining six genes up-regulated in RASFs (Table 4) have not yet been associated with RA but provide potential for further investigation. The solute carrier family members, SLC2A5, SLC14A1, and SLC12A8 are over-expressed in RASFs suggesting alterations in cellular metabolism. Complement Factor 1 (CF1) may represent a new target as the complement system plays a major role in the pathogenesis of rheumatoid disease . Expression of the plasminogen activator inhibitor gene, serpin peptidase inhibitor, clade B (ovalbumin), member 2 (SERPINB2) is decreased (-79.1-fold) in RASFs compared to control SFs (Table 4). The plasminogen activation pathway is dysregulated in arthritis . Aquoporin 1 (AQP1) expression has been shown to be up-regulated in the synovium from RA patients , but is down-regulated (-44.3-fold) in our samples. The coagulation factor X (F10), which may contribute to tissue injury and remodeling , is down-regulated (-27.5 fold). The hedgehog interacting protein (HHIP) inhibits the sonic hedgehog (SSH) signaling pathway. Inactivation of SSH inhibitor smoothened (Smo) blocks sonic hedgehog signaling and prevents osteophyte formation in the murine serum transfer arthritis model . Thus, the decrease (-26.6-fold) in HHIP expression observed in RASFs in this study may result in increased SSH activity resulting in advanced osteophyte formation.
Known isoforms expressed only in control SFs or only in RASFs
The top 10 isoforms expressed only in control SFs or only in RASFs are presented in Table 5. An expanded list of the top 50 up- and down-regulated known isoforms expressed only in either control SFs or in RASFs is presented in Additional file 5. The known isoforms identified in Table 5 correlate with the genes expressed only in control SFs or only in RASFs (Table 3 and Additional file 2). Single isoforms were detected for SCARA5, PLA2G2A, SPCS1, CITED2, IL13RA2, SLP1, FAM20A, NUMA1, PSAP, LRRC59, PPP1R14A, and SNHG6. Two isoforms were identified for PRG4, ACTG2, and CD59, while five and six isoforms exist for RPS24 and HMG1A, respectively (Table 5 and Additional file 5).
Known isoforms differentially expressed two-fold or greater between control SFs and RASFs
The top 10 up- and down-regulated known isoforms, along with the expanded top 50 list, in which expression of both samples in each group was up-regulated or down-regulated two-fold or greater between control SFs and RASFs are presented in Table 6 and in Additional file 6, respectively. Thirteen of the known isoforms identified in Table 6 can be found in the top 50 up-regulated and down-regulated genes presented in Table 4 and Additional file 4. A single isoform of IL26 is expressed 80.8-fold and correlates with the expression (80.8-fold) of the IL26 gene in RASFs. Seven known isoforms (ILI27, DHPS, BLCAP, LYNX1, C5orf13, APLP2, and CSRP1) are not represented in the top 50 regulated genes. One reason for this observation is differential isoform expression, as demonstrated by the two isoforms of Interferon, alpha-inducible protein 27 (ILI27). One ILI27 isoform is up-regulated 35.8-fold and one is down-regulated 216.8-fold. Two known isoforms were also identified for GCNT1, SLC2A5 and C5orf13 in the top 50 list.
Novel isoforms expressed only in control SFs or only in RASFs
The top 10 up- and down-regulated novel isoforms expressed only in control SFs or only in RASFs are presented in Table 7. An expanded list of the top 50 up- and down-regulated known isoforms expressed only in either control SFs or in RASFs is presented in Additional file 7. The list of the top 10 up-regulated novel isoforms includes transcripts for four unannotated genomic regions. The top 50 novel isoforms contains 21 transcripts from unannotated genomic regions. The list of top 10 down-regulated novel isoforms is divided into nine isoforms from annotated genes, including a novel transcript for HHIP, and one down-regulated novel isoform. There are transcripts for fourteen unannotated genomic regions in the top 50 down-regulated novel isoforms.
Novel isoforms differentially expressed two-fold or greater between control SFs and RASFs
The top 10 up- and down-regulated novel isoforms, along with the expanded top 50 list, in which expression of both samples in each group was up-regulated or down-regulated two-fold or greater between control SFs and RASFs are presented in Table 8 and in Additional file 8, respectively. A transcript of Fibrillin 1 (FBN1) is the top up-regulated novel isoform. Of note, a mutation in FBN1, which encodes an extracellular matrix glycoprotein, has been associated with the coexistence of Marfan’s Syndrome and ankylosing spondylitis . Novel isoforms from three unannotated regions of the genome were identified in the top 10 up-regulated novel isoforms. A total of 13 novel isoforms identified within unannotated regions of the genome were up-regulated in RASFs compared to SFs (Additional file 8). The list of top 10 down-regulated novel isoforms is divided into nine isoforms from annotated genes and one down-regulated novel isoform. A total of 10 novel isoforms within unannotated regions of the genome were down-regulated in RASFs compared to SFs. Interestingly, there are two novel transcripts for both HLA-DRB1 and SLC2A5 identified in this study (Additional file 8).
Network and pathway analyses of differentially expressed genes
To identify network and pathway connectivity, the differentially expressed gene lists of a two-fold or greater change in RASFs compared to SFs were submitted to Ingenuity Pathway Analysis (IPA) v9.0-3211 (Ingenuity Systems, Inc., Redwood City, CA), as described in the Material and Methods section. The networks affected by up-regulated genes and isoforms in RASFs compared to normal SFs are listed in Table 9. Consistent with the knowledge that RA is an immune disorder, the top network predicted to be affected by the up-regulated genes was Inflammatory Response, Immunological Disease, Cell Death, while the top network predicted to be affected by the up-regulated isoforms was Inflammatory Response, Cellular Movement, Cell-To-Cell Signaling and Interaction. The pathways affected by up-regulated genes and/or isoforms correlated with the pathways predicted to be affected by down-regulated gene expression and changes in isoform expression (Table 10). The top networks affected by down-regulated genes and isoforms in RASFs compared to normal SFs are Cellular Movement, Cell Death, and Tissue Development and Cellular Growth and Proliferation, Cell Death, Cellular Movement, respectively.
Canonical pathways analyses identified the pathways from the Ingenuity Pathways Analysis library of canonical pathways that were most significant to the data set. Genes with a two-fold or greater change in expression between SFs and RASFs and that were associated with a canonical pathway in Ingenuity’s Knowledge Base were considered for the analyses. The top canonical pathways affected by up-regulated genes and isoforms (Table 11) and the top canonical pathways affected by down-regulated genes and isoforms (Table 12) are in agreement with the networks (Tables 9 and 10) affected in RASFs. The top canonical pathways affected by up-regulated genes and isoforms (Table 11) are consistent with the knowledge that B cells, T cells, and macrophage cells play key roles in the inflammatory response and are involved in the activation and proliferation of RASFs [3, 4, 7]. These findings are further supported by the analysis of the pathways affected by the down-regulated genes and isoforms (Table 12). Dysregulation of the innate immune response and alterations in the number and types of cytokines and chemokines are well known features of RA [4, 7]. Altered cell cycle control of chromosomal replication and BRCA1 in DNA damage response, are in concordance with the hyperproliferation of synovial tissue and the corresponding decrease in apoptosis in RA [3, 38]. The identification of potential networks and pathways involved in arthritis may provide additional insights into the molecular and cellular mechanisms by which RASFs are involved in the pathogenesis of RA.
In the present study, we performed a comprehensive transcriptome analysis of human SF RNA isolated from healthy controls and patients with RA using the Illumina RNA-seq technique. It has revealed a complete picture of differentially expressed genes and their isoforms in RASFs and provided a global transcriptional insight into the novel roles of synovial fibroblasts in the pathogenesis of rheumatoid arthritis.
For RNA-seq, we used the Illumina HiScanSQ instrument to perform a 2 × 101 paired end run for all of our samples. The advantage of a paired end run is that both reads contain long range positional information, allowing for highly precise alignment of reads. We calculated the number of differentially expressed genes between RNA from two control SF and two RASF samples. We obtained a mean value of 84,177,268 reads per sample, which meets the criteria for sufficient sequence coverage for transcriptome profiling . Our mean rate of 86% total reads that map to the reference genome met quality standards of the RNA-seq technique . The breadth of the RNA sequencing reads covering chromosome 1 for both the RASFs and normal SFs indicates quality RNA-seq runs (Figure 1). Therefore, we are confident that our RNA-seq data provides an objective, high quality profile of the transcriptome in human RASFs and normal SFs.
The aim of this study was to provide a global glean into the transcriptional regulation in RASFs, which may provide mechanistic insights into the pathogenesis of rheumatoid arthritis. The activation and subsequent proliferation of SFs by proinflammatory cytokines produced by cells from both the innate and the adaptive immune systems plays a critical role in the pathogenesis of RA [3–5]. The production of additional cytokines, chemokines and matrix-degrading enzymes by RASFs leads ultimately to the progressive destruction of the joint that is a hallmark feature of RA [5–7]. However, the complete repertoire of active molecules, networks and pathways of differentially expressed genes and their isoforms of RASFs in this process are not characterized fully. Our study is filling this gap of knowledge. With RNA-seq, we found that 214 genes were not expressed in RASFs while 682 genes were only expressed in RASFs (Table 2). There are 122 up-regulated genes and 155 down-regulated genes by at least two-fold in RASFs compared to those in normal SFs. The majority of differentially expressed genes identified in this study (Tables 3 and 4 and Additional files 2 and 4) have not been previously reported to be altered in RASFs compared to normal SFs. One notable prowess of RNA-seq is to identify and quantify the expression of different isoforms of a gene. Gene isoforms are generated by alternative splicing or alternative promoter usage. Regulation of different gene isoform expression is a central aspect of most normal and disease processes. In this study, we detected more than 20,000 expressed known isoforms and more than 40,000 expressed novel isoforms (Table 2). Among them, there are 526 known isoforms which were not expressed in RASFs while 981 known isoforms were only expressed in RASFs. There are 343 up-regulated known isoforms and 262 down-regulated known isoforms by at least two-fold in RASFs compared to those in normal SFs. There are 105 novel isoforms which were not expressed in RASFs, while 152 novel isoforms were expressed only in RASFs. There are 561 up-regulated novel isoforms and 520 down-regulated novel isoforms by at least two-fold in RASFs compared to those in normal SFs. Network and canonical pathway analyses of differentially expressed genes and their known isoforms revealed that inflammatory response and cell death are represented strongly. Although these pathways have been predicted previously to correlate with RA, our study provided a more complete list of genes and isoforms involved in the inflammatory response and cell death pathways. We also identified other relevant novel networks and pathways, such as Antigen Presentation Pathway, Atherosclerosis Signalling, LXR/RXR Activation, and Role of BRCA1 in DNA Damage Response, whose dysregulation may each in part underlie their implication in the pathogenesis of RA.
Several microarray transcriptome analyses have been performed on RASFs [41–53]. The heterogeneous nature of RA and the different types of tissues used in these microarray studies leads to variations between the studies. The results from the present RNA-seq study both correlated and differed from previous microarray studies. The SFs used in our study were first isolated from synovial tissue either from healthy control donors or from patients with RA and cultured for two passages prior to RNA isolation. It should be noted, that this passage number is lower than what has been reported previously for gene profiling in SFs that have been cultured prior to RNA isolation. Del Rey et al.  and Masuda et al.  cultured SFs for 4 and 6 passages, respectively, before isolating RNA, while Haupl et al.  used immortalized SFs. The matrix metalloproteinases 1 (MMP1) and 3 (MMP3) are key players in the pathogenesis of RA . MMP1 and MMP3 were up-regulated 816.2- fold and 215.6-fold, respectively, in our study. Microarray analyses of RA synovial tissue in three separate studies detected increased MMP1 expression of 63.1-fold , 31.0-fold , and 36.6-fold . MMP3 expression was also increased 23.2-fold  and 18.7-fold  in these studies. Interleukin 1 beta (IL1B) and Interleukin 8 were up-regulated 3.2 and 9.3 fold, respectively, in RASFs from patients treated with prednisolone . In the present study, IL1B was decreased by 25.3-fold and IL8 was down-regulated 9.5-fold. Collagen, Type III, alpha 1 (COL3A1) was increased 1.76 fold in a microarray study  compared to a 1.3-fold decrease in the present study. Keratin 7 (KRT7) was down-regulated 0.49 by microarray analysis  and 14.6-fold by RNA-seq. The results presented in our study correlate well with what has been previously reported in the literature. Of the top 40 differentially expressed genes (Tables 3 and 4), 16 have been reported previously to be associated with RA (Additional file 3). Thus, we have identified 24 new potential gene targets among the genes listed in Tables 3 and 4 for further exploration. These findings are strengthened further by the ability of RNA-seq, as described above, to identify isoforms, both known and novel, that are expressed differentially in RA. With further improvements of next generation DNA sequencing techniques and further reductions of sequencing costs, it may be feasible to extend this study to analyze the transcriptomes of RASFs isolated from multiple patient samples at progressing stages of pathogenesis.
In summary, our first complete transcriptome analysis of synovial fibroblast RNA from patients with rheumatoid arthritis using RNA-seq has provided important insights into the transcriptional regulation of gene expression in RASFs. Further in-depth, follow-up analyses using large patient populations will be necessary to validate the alterations in transcriptional regulation reported in this study and to provide the resources necessary to elucidate the molecular mechanisms underlying the role of SFs in the pathogenesis of RA.
Human SF RNA from 2 healthy female donors and 2 adult female RA patients (Additional file 1) was purchased from Cell Applications, Inc. (San Diego, CA). SFs were first isolated from synovial tissue either from healthy control donors or from patients with RA and cultured for two passages prior to RNA isolation. Paired-end cDNA libraries were prepared for each sample and sequenced using the Illumina TruSeq RNA Sample Preparation Kit, as described previously [11, 12]. Briefly, the cDNA libraries were quantified using a Biotek EPOCH spectrophotometer and checked for quality and size using a Bio-Rad Experion DNA 1K chip. The four cDNA libraries were each diluted to 6 pM and spiked with 1% phiX control to improve base calling while sequencing. A 6 pM dilution of phiX control sample was also prepared for analysis. Following the Illumina cBot and HiSeq protocols, the four libraries and the phiX control underwent cluster generation on a HiSeq PE flow cell v3 and were then sequenced using a HiScanSQ (Illumina). A paired-end (2×101) run was performed using the SBS Kit (Illumina). Real-time analysis and base calling were performed using the HiSeq Control Software Version 1.4.5 (Illumina). The resulting basecalling (.bcl) files were converted to. FASTQ files using Illumina’s CASAVA 1.8 software. The number of reads for each sample type was analyzed using the Student’s t-test in SigmaPlot version 11.0 (Systat Software Inc., San Jose, CA). A p-value of below 0.05 was considered significant. The sequence data have been submitted to the NCBI Short Read Archive with accession number SRA048057.1.
Mapping of RNA-seq reads and transcript assembly and abundance estimation using Tuxedo Suite
Paired-end fastq sequence reads for each sample were aligned to the UCSC Homo sapiens reference genome hg19 using TopHat v1.3.0 [54, 55] integrated with Bowtie v0.12.7 , as described previously [11, 12]. The resulting aligned reads were analyzed further by Cufflinks v1.0.3 [55, 57]. The aligned reads were assembled into transcripts, either with or without a reference genome, and the expression of those transcripts were reported in Fragments Per Kilobase of exon per Million fragments mapped (FPKM). Cuffdiff analysis was performed, with use of the reference genome, to determine differential expression of known isoforms between pooled RA patient samples and pooled control samples. To detect novel isoforms, Cufflinks was run without a reference genome. The RA and control transcript files were compared to the reference genome using Cuffcompare to filter out previously discovered transcripts. To test the differential expression of these novel isoforms, Cuffdiff analyses were performed using the combined transcript files as the reference genome. Cuffdiff analyses were performed two ways: comparing the RA patient transcripts to the control transcripts, using the RA patient transcripts as the reference genome; and comparing the RA patient transcripts to the control transcripts, using the control transcripts as the reference genome.
Visualization of mapped reads
Aligned reads were visualized using a local copy of the Integrative Genomics Viewer (http://www.broadinstitute.org/igv/). The output files generated from TopHat were converted into files viewable in IGV by BEDTools  and then processed further by the “count” function in igvtools (included with the IGV software) to create an average alignment track viewable as a bar chart. The log2 of the frequency of the reads was plotted to better visualize the extensive range of the read coverage. Individual gene views were created by first merging the TopHat output files from the RA and control samples into two files using SAMTools . These merged files were processed in the same way as above with the “count” function in igvtools. The raw frequency of the reads was visualized in this case.
Automated literature search
Multiplex literature mining analysis was conducted with PubMatrix,  as described previously . We restricted our search to human symbols approved by HUGO Gene Nomenclature Committee (HGNC) for the top 10 genes and isoforms for each category. Terms “rheumatoid arthritis”, “osteoarthritis”, “arthritis” and “disease” were used for cross-referencing candidate genes.
Functional analysis of differentially expressed gene lists using ingenuity pathway analysis
The differentially expressed gene lists were submitted to Ingenuity Pathway Analysis (IPA) v9.0-3211 (Ingenuity Systems, Inc., Redwood City, CA). Genes with a two-fold or greater change in expression between the RA group and the control group were used. The settings for the core analysis were as follows: Ingenuity Knowledge Base; Endogenous Chemicals not included; Direct and Indirect relationships; molecules per pathway: 70; and networks per analysis: 25.
Lee DM, Weinblatt ME: Rheumatoid arthritis. Lancet. 2001, 358: 903-911. 10.1016/S0140-6736(01)06075-5
Huber LC, Distler O, Tarner I, Gay RE, Gay S, Pap T: Synovial fibroblasts: key players in rheumatoid arthritis. Rheumatol (Oxford). 2006, 45: 669-675. 10.1093/rheumatology/kel065. 10.1093/rheumatology/kel065
Bartok B, Firestein GS: Fibroblast-like synoviocytes: key effector cells in rheumatoid arthritis. Immunol Rev. 2010, 233: 233-255. 10.1111/j.0105-2896.2009.00859.x
Scott DL, Wolfe F, Huizinga TW: Rheumatoid arthritis. Lancet. 2010, 376: 1094-1108. 10.1016/S0140-6736(10)60826-4
Gierut A, Perlman H, Pope RM: Innate immunity and rheumatoid arthritis. Rheum Dis Clin North Am. 2010, 36: 271-296. 10.1016/j.rdc.2010.03.004
Cooles FA, Isaacs JD: Pathophysiology of rheumatoid arthritis. Curr Opin Rheumatol. 2011, 23: 233-240. 10.1097/BOR.0b013e32834518a3
Firestein GS: Evolving concepts of rheumatoid arthritis. Nature. 2003, 423: 356-361. 10.1038/nature01661
Buch MH, Emery P: New therapies in the management of rheumatoid arthritis. Curr Opin Rheumatol. 2011, 23: 245-251. 10.1097/BOR.0b013e3283454124
Davey JW, Hohenlohe PA, Etter PD, Boone JQ, Catchen JM, Blaxter ML: Genome-wide genetic marker discovery and genotyping using next-generation sequencing. Nat Rev Genet. 2011, 12: 499-510. 10.1038/nrg3012
Marioni JC, Mason CE, Mane SM, Stephens M, Gilad Y: RNA-seq: an assessment of technical reproducibility and comparison with gene expression arrays. Genome Res. 2008, 18: 1509-1517. 10.1101/gr.079558.108
Zhang LQ, Cheranova D, Gibson M, Ding S, Heruth DP, Fang D, Ye SQ: RNA-seq reveals novel transcriptome of genes and their isoforms in human pulmonary microvascular endothelial cells treated with thrombin. PLoS One. 2012, 7: e31229. 10.1371/journal.pone.0031229
Cheranova DGM, Chaudhary S, Zhang LQ, Heruth DP, Grigoryev DN, Ye SQ: RNA-seq analysis of transcriptomes in thrombin-treated and control human pulmonary microvascular endothelial cells. J Vis Exp. 2012, in press.
Twine NA, Janitz K, Wilkins MR, Janitz M: Whole transcriptome sequencing reveals gene expression and splicing differences in brain regions affected by Alzheimer's disease. PLoS One. 2011, 6: e16266. 10.1371/journal.pone.0016266
Newton JL, Harney SM, Wordsworth BP, Brown MA: A review of the MHC genetics of rheumatoid arthritis. Genes Immun. 2004, 5: 151-157. 10.1038/sj.gene.6364045
Raum D, Awdeh Z, Glass D, Kammer G, Khan MA, Coblyn JS, Weinblatt M, Holdsworth D, Strong L, Rossen RD: Extended haplotypes of chromosome 6 in adult rheumatoid arthritis. Arthritis Rheum. 1984, 27: 516-521. 10.1002/art.1780270506
Taylor KE, Criswell LA: Conditional analysis of the major histocompatibility complex in rheumatoid arthritis. BMC Proc. 2009, 7 (3): S36.
Plenge RM, Cotsapas C, Davies L, Price AL, de Bakker PI, Maller J, Pe'er I, Burtt NP, Blumenstiel B, DeFelice M: Two independent alleles at 6q23 associated with risk of rheumatoid arthritis. Nat Genet. 2007, 39: 1477-1482. 10.1038/ng.2007.27
Jiang L, Salao K, Li H, Rybicka JM, Yates RM, Luo XW, Shi XX, Kuffner T, Tsai VW, Husaini Y: Intracellular chloride channel protein CLIC1 regulates macrophage functions via modulation of phagosomal acidification. J Cell Sci. 2012, in press.
Uesugi H, Ozaki S, Sobajima J, Osakada F, Shirakawa H, Yoshida M, Nakao K: Prevalence and characterization of novel pANCA, antibodies to the high mobility group non-histone chromosomal proteins HMG1 and HMG2, in systemic rheumatic diseases. J Rheumatol. 1998, 25: 703-709.
Pohlers D, Beyer A, Koczan D, Wilhelm T, Thiesen HJ, Kinne RW: Constitutive upregulation of the transforming growth factor-beta pathway in rheumatoid arthritis synovial fibroblasts. Arthritis Res Ther. 2007, 9: R59. 10.1186/ar2217
Williams AS, Mizuno M, Richards PJ, Holt DS, Morgan BP: Deletion of the gene encoding CD59a in mice increases disease severity in a murine model of rheumatoid arthritis. Arthritis Rheum. 2004, 50: 3035-3044. 10.1002/art.20478
Becker KG, Hosack DA, Dennis G, Lempicki RA, Bright TJ, Cheadle C, Engel J: PubMatrix: a tool for multiplex literature mining. BMC Bioinforma. 2003, 4: 61-10.1186/1471-2105-4-61. 10.1186/1471-2105-4-61
Li JY, Paragas N, Ned RM, Qiu A, Viltard M, Leete T, Drexler IR, Chen X, Sanna-Cherchi S, Mohammed F: Scara5 is a ferritin receptor mediating non-transferrin iron delivery. Dev Cell. 2009, 16: 35-46. 10.1016/j.devcel.2008.12.002
Zhao Y, Wang Y, Hu J, Zhang X, Zhang YW: CutA divalent cation tolerance homolog (Escherichia coli) (CUTA) regulates beta-cleavage of beta-amyloid precursor protein (APP) through interacting with beta-site APP cleaving protein 1 (BACE1). J Biol Chem. 2012, 287: 11141-11150. 10.1074/jbc.M111.330209
Zhen Y, Sorensen V, Skjerpen CS, Haugsten EM, Jin Y, Walchli S, Olsnes S, Wiedlocha A: Nuclear import of exogenous FGF1 requires the ER-protein LRRC59 and the importins Kpnalpha1 and Kpnbeta1. Traffic (Copenhagen, Denmark). 2012, 13: 650-664.
Xiao L, Gong LL, Yuan D, Deng M, Zeng XM, Chen LL, Zhang L, Yan Q, Liu JP, Hu XH: Protein phosphatase-1 regulates Akt1 signal transduction pathway to control gene expression, cell survival and differentiation. Cell Death Differ. 2010, 17: 1448-1462. 10.1038/cdd.2010.16
Corvaisier M, Delneste Y, Jeanvoine H, Preisser L, Blanchard S, Garo E, Hoppe E, Barré B, Audran M, Bouvard B: IL-26 is overexpressed in rheumatoid arthritis and induces proinflammatory cytokine production and Th17 cell generation. PLoS Biol. 2012, 10: e1001395. 10.1371/journal.pbio.1001395
Liu C, Batliwalla F, Li W, Lee A, Roubenoff R, Beckman E, Khalili H, Damle A, Kern M, Furie R: Genome-wide association scan identifies candidate polymorphisms associated with differential response to anti-TNF treatment in rheumatoid arthritis. Mol med (Cambridge, Mass). 2008, 14: 575-581.
Mishima K, Otani H, Tanabe T, Kawasaki H, Oshiro A, Saito N, Ogawa R, Inagaki C: Molecular mechanisms for alpha2-adrenoceptor-mediated regulation of synoviocyte populations. Jpn J Pharmacol. 2001, 85: 214-226. 10.1254/jjp.85.214
Rico MC, Rough JJ, Del Carpio-Cano FE, Kunapuli SP, DeLa Cadena RA: The axis of thrombospondin-1, transforming growth factor beta and connective tissue growth factor: an emerging therapeutic target in rheumatoid arthritis. Curr Vasc Pharmacol. 2010, 8: 338-343. 10.2174/157016110791112296
Frolova EG, Sopko N, Blech L, Popovic ZB, Li J, Vasanji A, Drumm C, Krukovets I, Jain MK, Penn MS: Thrombospondin-4 regulates fibrosis and remodeling of the myocardium in response to pressure overload. FASEB journal: official publication of the Federation of American Societies for Experimental Biology. 2012, 26: 2363-2373. 10.1096/fj.11-190728. 10.1096/fj.11-190728
Happonen KE, Heinegard D, Saxne T, Blom AM: Interactions of the complement system with molecules of extracellular matrix: relevance for joint diseases. Immunobiology. 2012, 217: 1088-1096. 10.1016/j.imbio.2012.07.013
Busso N, Peclat V, So A, Sappino AP: Plasminogen activation in synovial tissues: differences between normal, osteoarthritis, and rheumatoid arthritis joints. Ann Rheum Dis. 1997, 56: 550-557. 10.1136/ard.56.9.550
Mobasheri A, Moskaluk CA, Marples D, Shakibaei M: Expression of aquaporin 1 (AQP1) in human synovitis. Ann Anat= Anatomischer Anzeiger: official organ of the Anatomische Gesellschaft. 2010, 192: 116-121. 10.1016/j.aanat.2010.01.001. 10.1016/j.aanat.2010.01.001
Krupiczojc MA, Scotton CJ, Chambers RC: Coagulation signalling following tissue injury: focus on the role of factor Xa. Int J Biochem Cell Biol. 2008, 40: 1228-1237. 10.1016/j.biocel.2008.02.026
Ruiz-Heiland G, Horn A, Zerr P, Hofstetter W, Baum W, Stock M, Distler JH, Nimmerjahn F, Schett G, Zwerina J: Blockade of the hedgehog pathway inhibits osteophyte formation in arthritis. Ann Rheum Dis. 2012, 71: 400-407. 10.1136/ard.2010.148262
Kiss C, Jonap I, Gergely P, Poor G: Coexistent Marfan's syndrome and ankylosing spondylitis. J Rheumatol. 2006, 33: 1199-1200.
Maas K, Westfall M, Pietenpol J, Olsen NJ, Aune T: Reduced p53 in peripheral blood mononuclear cells from patients with rheumatoid arthritis is associated with loss of radiation-induced apoptosis. Arthritis Rheum. 2005, 52: 1047-1057. 10.1002/art.20931
Sultan M, Schulz MH, Richard H, Magen A, Klingenhoff A, Scherf M, Seifert M, Borodina T, Soldatov A, Parkhomchuk D: A global view of gene activity and alternative splicing by deep sequencing of the human transcriptome. Science. 2008, 321: 956-960. 10.1126/science.1160342
Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B: Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008, 5: 621-628. 10.1038/nmeth.1226
Galligan CL, Baig E, Bykerk V, Keystone EC, Fish EN: Distinctive gene expression signatures in rheumatoid arthritis synovial tissue fibroblast cells: correlates with disease activity. Genes Immun. 2007, 8: 480-491. 10.1038/sj.gene.6364400
Del Rey MJ, Usategui A, Izquierdo E, Canete JD, Blanco FJ, Criado G, Pablos JL: Transcriptome analysis reveals specific changes in osteoarthritis synovial fibroblasts. Ann Rheum Dis. 2012, 71: 275-280. 10.1136/annrheumdis-2011-200281
Del Rey MJ, Izquierdo E, Usategui A, Gonzalo E, Blanco FJ, Acquadro F, Pablos JL: The transcriptional response of normal and rheumatoid arthritis synovial fibroblasts to hypoxia. Arthritis Rheum. 2010, 62: 3584-3594. 10.1002/art.27750
Watanabe N, Ando K, Yoshida S, Inuzuka S, Kobayashi M, Matsui N, Okamoto T: Gene expression profile analysis of rheumatoid synovial fibroblast cultures revealing the overexpression of genes responsible for tumor-like growth of rheumatoid synovium. Biochem Biophys Res Commun. 2002, 294: 1121-1129. 10.1016/S0006-291X(02)00608-3
Devauchelle V, Marion S, Cagnard N, Mistou S, Falgarone G, Breban M, Letourneur F, Pitaval A, Alibert O, Lucchesi C: DNA microarray allows molecular profiling of rheumatoid arthritis and identification of pathophysiological targets. Genes Immun. 2004, 5: 597-608. 10.1038/sj.gene.6364132
van der Pouw Kraan TC, van Gaalen FA, Huizinga TW, Pieterman E, Breedveld FC, Verweij CL: Discovery of distinctive gene expression profiles in rheumatoid synovium using cDNA microarray technology: evidence for the existence of multiple pathways of tissue destruction and repair. Genes Immun. 2003, 4: 187-196. 10.1038/sj.gene.6363975
Masuda K, Masuda R, Neidhart M, Simmen BR, Michel BA, Muller-Ladner U, Gay RE, Gay S: Molecular profile of synovial fibroblasts in rheumatoid arthritis depends on the stage of proliferation. Arthritis Res. 2002, 4: R8. 10.1186/ar427
Haupl T, Yahyawi M, Lubke C, Ringe J, Rohrlach T, Burmester GR, Sittinger M, Kaps C: Gene expression profiling of rheumatoid arthritis synovial cells treated with antirheumatic drugs. J Biomol Screen. 2007, 12: 328-340. 10.1177/1087057107299261
Lindberg J, Wijbrandts CA, van Baarsen LG, Nader G, Klareskog L, Catrina A, Thurlings R, Vervoordeldonk M, Lundeberg J, Tak PP: The gene expression profile in the synovium as a predictor of the clinical response to infliximab treatment in rheumatoid arthritis. PLoS One. 2010, 5: e11310. 10.1371/journal.pone.0011310
Green MJ, Gough AK, Devlin J, Smith J, Astin P, Taylor D, Emery P: Serum MMP-3 and MMP-1 and progression of joint damage in early rheumatoid arthritis. Rheumatol (Oxford). 2003, 42: 83-88. 10.1093/rheumatology/keg037. 10.1093/rheumatology/keg037
Biswas S, Manikandan J, Pushparaj PN: Decoding the differential biomarkers of Rheumatoid arthritis and Osteoarthritis: a functional genomics paradigm to design disease specific therapeutics. Bioinformation. 2011, 6: 153-157. 10.6026/97320630006153
Ungethuem U, Haeupl T, Witt H, Koczan D, Krenn V, Huber H, von Helversen TM, Drungowski M, Seyfert C, Zacher J: Molecular signatures and new candidates to target the pathogenesis of rheumatoid arthritis. Physiol Genomics. 2010, 42A: 267-282. 10.1152/physiolgenomics.00004.2010
Yarilina A, Park-Min KH, Antoniv T, Hu X, Ivashkiv LB: TNF activates an IRF1-dependent autocrine loop leading to sustained expression of chemokines and STAT1-dependent type I interferon-response genes. Nat Immunol. 2008, 9: 378-387.
Trapnell C, Pachter L, Salzberg SL: TopHat: discovering splice junctions with RNA-seq. Bioinformatics. 2009, 25: 1105-1111. 10.1093/bioinformatics/btp120
Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L: Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012, 7: 562-578.
Langmead B, Trapnell C, Pop M, Salzberg SL: Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009, 10: R25. 10.1186/gb-2009-10-3-r25
Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L: Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010, 28: 511-515. 10.1038/nbt.1621
Quinlan AR, Hall IM: BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010, 26: 841-842. 10.1093/bioinformatics/btq033
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R: The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009, 25: 2078-2079. 10.1093/bioinformatics/btp352
Grigoryev DN, Ma SF, Irizarry RA, Ye SQ, Quackenbush J, Garcia JG: Orthologous gene-expression profiling in multi-species models: search for candidate genes. Genome Biol. 2004, 5: R34. 10.1186/gb-2004-5-5-r34
This work is in part supported by the start-up fund and William R. Brown Missouri Endowment of Children’s Mercy Hospitals and Clinics, University of Missouri at Kansas City (Ye, S.Q.).
The authors declare that they have no competing interests.
DPH and MG carried out the experimental studies and data analysis. DPH, DNG and LQZ participated in the design of the study and performed the statistical analysis. DPH drafted the manuscript. SQY conceived of the study, and participated in its design and coordination and helped to draft the manuscript. All authors read and approved the final manuscript.
Electronic supplementary material
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.