Systematic Analysis of Differential Expression Profile in Rheumatoid Arthritis Chondrocytes Using Next-Generation Sequencing and Bioinformatics Approaches

Abstract

Cartilage destruction in rheumatoid arthritis (RA) occurs primarily in the pannus-cartilage interface. The close contact of the synovium-cartilage interface implicates crosstalk between synovial fibroblasts and chondrocytes. The aim of this study is to explore the differentially expressed genes and novel microRNA regulations potentially implicated in the dysregulated cartilage homeostasis in joint destruction of RA. Total RNAs were extracted from human primary cultured normal and RA chondrocytes for RNA and small RNA expression profiling using next-generation sequencing. Using systematic bioinformatics analyses, we identified 463 differentially expressed genes in RA chondrocytes were enriched in biological functions related to altered cell cycle process, inflammatory response and hypoxic stimulation. Moreover, fibroblast growth factor 9 (FGF9), kynureninase (KYNU), and regulator of cell cycle (RGCC) were among the top dysregulated genes identified to be potentially affected in the RA joint microenvironment, having similar expression patterns observed in arrays of clinical RA synovial tissues from the Gene Expression Omnibus database. Additionally, among the 31 differentially expressed microRNAs and 10 candidate genes with potential microRNA-mRNA interactions in RA chondrocytes, the novel miR-140-3p-FGF9 interaction was validated in different microRNA prediction databases, and proposed to participate in the pathogenesis of joint destruction through dysregulated cell growth in RA. The findings provide new perspectives for target genes in the management of cartilage destruction in RA.

Keywords: bioinformatics.; cell cycle; chondrocytes; next-generation sequencing; rheumatoid arthritis.

Figure 1

Flowchart of study design. The primary human chondrocytes of normal (HC) and rheumatoid arthritis (HC-RA) knee cartilages were cultured and harvasted for RNA sequencing and expression profiling. Differentially expressed genes with > 2.0 fold change and > 0.3 fragments per kilobase of transcript million (FPKM) were selected for further enrichment analyses using different bioinforatmics resources. In addition, differentially expressed microRNAs with > 2.0 fold change and > 10 reads per million (RPM) were selected for further putative targets using miRmap target prediction database. The identified potential miRNA-mRNA interactions were systematically validated in different miRNA target prediction databases. Finally, rheumatoid arthritis (RA) related arrays from clinical RA joint tissue specimen were searched in the Gene Expression Omnibus (GEO) database, and the expression patterns of candidate genes of interest in these arrays were analyzed.

Figure 2

Display of differential expression patterns of normal and rheumatoid arthritis (RA) chondrocytes from deep sequencing. (A) The RNA sequencing result of differential gene expression in normal (HC) and RA chondrocytes (HC-RA) were plotted by volcano plot. The x-axis indicated the logarithm to the base 2 of expression fold-change (HC-RA/HC) and the y-axis indicated the negative logarithm to the base 10 of the p-values. Red circular marks represented up-regulated genes in HC-RA, and green triangular marks represented down-regulated genes in HC-RA. Vertical lines reflected the filtering thresholds of 2.0-fold-change, and horizontal line reflected filtering threshold of p-value = 0.05. A total of 249 significantly up-regulated and 214 significantly down-regulated genes in HC-RA were identified. (B) The comparison of the difference in fragments per kilobase of transcript per million (FPKM) performance between HC and HC-RA after screening (threshold setting: > 0.3 FPKM and > 2.0-fold-change) were displayed as density plot. The x-axis indicated the logarithm to the base 10 of FPKM, and the y-axis indicated read density.

Figure 3

The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of differentially expressed genes in DAVID database. The 463 differentially expressed genes in RA chondrocytes were uploaded into DAVID database for enrichment analysis. The top 10 GO and KEGG pathway analysis results of these dysregulated genes in RA chondrocytes were displayed in bar chart. The bars indicated the -Log₁₀(p value) of each GO and KEGG term, and the numbers shown at the right side of each bar indicated the number of genes involved in each term.

Figure 4

The Gene Set Enrichment Analysis (GSEA) result of differentially expressed genes. The 463 differentially expressed genes in RA chondrocytes were uploaded into GSEA for enrichment analysis. The h.all.v5.1.symbols.gmt [Hallmarks] gene sets database was used as the gene set collection for analysis. GSEA performed 1000 permutations. The maximum and minimum sizes for gene sets were 500 and 15, respectively. Cutoff for significant gene sets was false discovery rate < 25%.

Figure 5

The protein-protein interaction (PPI) network analysis of differentially expressed genes using STRING database. The 463 differentially expressed genes were input into STRING database for PPI network analysis, and achieved a PPI network of 460 nodes and 1180 edges, with PPI enrichment p-value < 1.0 x 10^-16. The three primary clusters of subnetworks were analyzed by plug-in MCODE in Cytoscape, and the nodes in each cluster were input into STRING database to obtain the PPI subnetworks.

Figure 6

Expression patterns of top 10 up-regulated and 10 down-regulated genes identified from RA chondrocytes in a representative RA synovial tissue array in the Gene Expression Omnibus (GEO) database. The expression values of the top 20 dysregulated genes in RA chondrocytes were analyzed in a representative array of clinical specimen of normal and RA synovial tissues from the GEO database (GSE55235). Significant up-regulation of FGF9 and KYNU, and down-regulation of RGCC were observed to have similar expression pattern in the RA synovial tissues, compared to the normal synovium. * indicated p < 0.05, ** indicated p < 0.01, *** indicated p < 0.001, and n.s. indicated no statistical significance. (Probe ID reference: FGF7, 205782_at; FGF9, 206404_at; SERPINF1, 202283_at; KYNU_1, 217388_s_at; KYNU_2, 210663_s_at; KYNU_3, 204385_at; KYNU_4, 210662_at; VGLL3, 220327_at; IGF1_1, 209541_at; IGF1_2, 209540_at; IGF1_3, 209542_x_at; IGF1_4, 211577_s_at; PCBP3, 205663_at; COL2A1_1, 217404_s_at; COL2A1_2, 213492_at; COL9A2, 213622_at; RGCC, 218723_s_at; COL11A2_1, 213870_at; COL11A2_2, 216993_s_at; F3, 204363_at and ARHGDIB, 201288_at). Some genes did not have identical probes in the array data.

Figure 7

Merged network analysis from Ingenuity Pathway Analysis (IPA) software for associations among molecules related to rheumatoid arthritis, inflammation of joint, cell cycle and apoptosis. The merged network of “Rheumatoid arthritis”, “Inflammation of joint”, “Apoptosis of chondrocytes” and “Cell cycle progression” categorized in the IPA software was obtained from IPA software. The overlay canonical pathway of “Role of osteoblasts, osteoclasts and chondrocytes in rheumatoid arthritis” identified 7 molecules interconnected to networks of RA and cell cycle progress. Among the merged network, the down-regulated RGCC in RA chondrocytes participated simultaneously in the network of RA, inflammation of joint and cell cycle progression, with additional connection to Cyclin B1 (CCNB1), Toll-like receptor 4 (TLR4), cyclin-dependent kinase 1 (CDK1) and polo-like kinase 1 (PLK1), as shown in light blue lines.

Figure 8

Differentially expressed microRNAs with potential microRNA-mRNA interactions identified in primary RA chondrocytes. (A) A total of 31 differentially expressed microRNAs (selection criteria of > 2.0-fold change and reads per million (RPM) > 10) from next generation sequencing method were identified, and the heat map according to z-score value is illustrated. (B) Using the miRmap database for microRNA target prediction, there were 80 putative targets of 2 up-regulated microRNAs and 404 putative targets of 29 down-regulated microRNAs obtained (selection criteria of miRmap score ≥ 99.0). Matching to the 214 down-regulated genes and 249 up-regulated genes identified in the RA chondrocytes, the Venn diagram analysis identified 10 up-regulated genes with potential microRNA-mRNA interactions.

Figure 9

The putative binding site of miR-140-3p on FGF9. The sequence and putative binding sites of miR-140-3p on the 3′UTR of FGF9 at the positions of 746-752, 1908-1914 and 2783-2789 were validated in miRmap (A); TargetScan (B); and miRDB (C).

Figure 10

The schematic summary of proposed molecular signatures in RA joint microenvironment.