Runrun Zhang1,2†
Yehua Jin1,3†
Cen Chang1,3
Lingxia Xu1,3
Yanqin Bian4,5
Yu Shen4,5
Yang Sun4,5
Songtao Sun6
Steven J. Schrodi7
Shicheng Guo7*
Dongyi He1,5*- 1Department of Rheumatology, Shanghai Guanghua Hospital, Shanghai University of Traditional Chinese Medicine, Shanghai, China
- 2Department of Rheumatology, The Second Affiliated Hospital of Shandong University of Traditional Chinese Medicine, Jinan, China
- 3Guanghua Clinical Medical College, Shanghai University of Traditional Chinese Medicine, Shanghai, China
- 4Guanghua Integrative Medicine Hospital, Shanghai, China
- 5Arthritis Institute of Integrated Traditional and Western Medicine, Shanghai Chinese Medicine Research Institute, Shanghai, China
- 6Department of Orthopedics, Shanghai Guanghua Hospital, Shanghai University of Traditional Chinese Medicine, Shanghai, China
- 7Department of Medical Genetics, School of Medicine and Public Health, University of Wisconsin-Madison, Madison, WI, United States
Purpose: This study aimed to provide a comprehensive understanding of the genome-wide expression patterns in the synovial tissue samples of patients with rheumatoid arthritis (RA) to investigate the potential mechanisms regulating RA occurrence and development.
Methods: Transcription profiles of the synovial tissue samples from nine patients with RA and 15 patients with osteoarthritis (OA) (control) from the East Asian population were generated using RNA sequencing (RNA-seq). Gene set enrichment analysis (GSEA) was used to analyze all the detected genes and the differentially expressed genes (DEGs) were identified using DESeq. To further analyze the DEGs, the Gene Ontology (GO) functional enrichment and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed. The protein–protein interaction (PPI) network of the DEGs was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) and the hub genes were identified by topology clustering with the Molecular Complex Detection (MCODE)-Cytoscape. The most important hub genes were validated using quantitative real-time PCR (qRT-PCR).
Results: Of the 17,736 genes detected, 851 genes were identified as the DEGs (474 upregulated and 377 downregulated genes) using the false discovery rate (FDR) approach. GSEA revealed that the significantly enriched gene sets that positively correlated with RA were CD40 signaling overactivation, Th1 cytotoxic module, overactivation of the immune response, adaptive immune response, effective vs. memory CD8+ T cells (upregulated), and naïve vs. effective CD8+ T cells (downregulated). Biological process enrichment analysis showed that the DEGs were significantly enriched for signal transduction (P = 3.01 × 10−6), immune response (P = 1.65 × 10−24), and inflammatory response (P = 5.76 × 10−10). Molecule function enrichment analysis revealed that the DEGs were enriched in calcium ion binding (P = 1.26 × 10−5), receptor binding (P = 1.26 × 10−5), and cytokine activity (P = 2.01 × 10−3). Cellular component enrichment analysis revealed that the DEGs were significantly enriched in the plasma membrane (P = 1.91 × 10−31), an integral component of the membrane (P = 7.39 × 10−13), and extracellular region (P = 7.63 × 10−11). The KEGG pathway analysis showed that the DEGs were enriched in the cytokine–cytokine receptor interaction (P = 3.05 × 10−17), chemokine signaling (P = 3.50 × 10−7), T-cell receptor signaling (P = 5.17 × 10−4), and RA (P = 5.17 × 10−4) pathways. We confirmed that RA was correlated with the upregulation of the PPI network hub genes, such as CXCL13, CXCL6, CCR5, CXCR5, CCR2, CXCL3, and CXCL10, and the downregulation of the PPI network hub gene such as SSTR1.
Conclusion: This study identified and validated the DEGs in the synovial tissue samples of patients with RA, which highlighted the activity of a subset of chemokine genes, thereby providing novel insights into the molecular mechanisms of RA pathogenesis and identifying potential diagnostic and therapeutic targets for RA.
Introduction
Rheumatoid arthritis (RA) is an autoimmune disease characterized by synovial inflammation, hyperplasia, and cartilage and bone destruction. Clinical manifestations of RA include joint pain, swelling, stiffness, and deformation (1, 2). RA pathogenesis is thought to be related to genetic and environmental factors, obesity, diet, and gut microbiota composition (3, 4). Osteoarthritis (OA) is a joint disease characterized by the degeneration of the synovial joint and loss of articular cartilage, with primary clinical features including pain and loss of mobility (5). Genetic factors, diet, estrogen level, obesity, bone density, and joint laxity play a role in the pathogenesis of OA (6). As both the RA and OA share common physiological targets, biomarkers present in the synovial tissue that could discriminate between these diseases should be determined (7, 8).
Transcriptomics is tissue specific and it offers an avenue for the investigation of the effects of the disease at the cellular level that is likely to play an important role in the etiology of the disease (9). RNA sequencing (RNA-seq) technology has become the primary step in transcriptomic studies for the characterization of gene expression within cells and tissues. Therefore, identification of differential gene expression in the synovial tissues of patients with RA and OA using RNA-seq may provide new insights into the molecular pathophysiology of these diseases.
In this study, to better understand the functional differences at the transcriptome level between RA and OA, we analyzed the whole genes detected by gene set enrichment analysis (GSEA), identified the differentially expressed genes (DEGs) in the synovial tissue samples of patients with RA and OA using RNA-seq, analyzed the DEGs by the Gene Ontology (GO) functional enrichment and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses, constructed the protein–protein interaction (PPI) network, and screened and verified the hub genes. The validated hub genes may serve as critical molecular markers for identifying differences in the synovial tissues of patients with RA and OA due to their central role in gene expression networks.
Materials and Methods
Patient Information and Tissue Collection
In this study, we included nine patients with RA, who were diagnosed based on the 2010 American College of Rheumatology (ACR)/European League Against Rheumatism (EULAR) classification criteria for RA (10) and 15 patients with OA, who were diagnosed according to the ACR OA classification criteria (11). The synovial tissue samples of patients with RA and OA were obtained from the Guanghua Hospital, Shanghai, China. All of the involved patients underwent a knee replacement. After removing excess fat and vascular tissue, the synovial tissue samples were placed in liquid nitrogen till further use. The demographic information of patients with RA is given in Table 1 and clinical information of patients with RA and OA is given in Supplementary Table 1. This study was approved by the Ethics Committee of Guanghua Hospital of Integrated Traditional Chinese and Western Medicine (approval number: 2018-K-12) and a written consent was obtained from all the patients prior to knee replacement.
Table 1. Demographic information of patients with rheumatoid arthritis (RA).
Ribonucleic Acid Isolation and Library Preparation
Total RNA was extracted from the synovial tissue samples using TRIzol Reagent (Thermo Fisher Scientific, Waltham, Massachusetts, USA) according to the manufacturer's protocol. A NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific) was used to evaluate RNA quality and quantify each RNA sample. RNA integrity was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, California, USA). Total RNA with a standard RNA integrity number (RIN) ≥ 7.0 and 28S/18S ≥ 0.7 was subjected to RNA-seq. RNA libraries were constructed using the TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, San Diego, California, USA) according to the manufacturer's instructions.
Ribonucleic Acid Sequencing and Identification of the Differentially Expressed Genes
The libraries were sequenced on an Illumina HiSeq × 10 platform. Raw data (raw reads) in FASTQ format were first processed using Trimmomatic (12) and low-quality reads were removed to obtain clean reads. Clean reads were then mapped to the human genome (GRCh38) using HISAT2 (13). The fragments per kilobase of transcript per millions mapped reads (FPKM) of each gene were calculated using Cufflinks (14, 15) and the read counts of each gene were obtained using HTSeq-Count (16). Differential expression analysis was performed using the DESeq R package (17). An adjusted P < 0.05 and |log2FoldChange| ≥ 1.5 were set as the threshold for significant differential expression; P-value was adjusted using the false discovery rate (FDR). The expression profiling data were obtained from https://github.com/dongyihe/rheumatoidarthritis.
Gene Set Enrichment Analysis Based on RNA Sequencing Effectively Detected Genes
Gene set enrichment analysis was performed using a defined set of genes to determine statistically significant differences between the RA and OA groups, using R software (https://www.r-project.org) and the data set was obtained from the Molecular Signatures Database version 7.2 (MSigDB; GSEA-MSigDB website). MSigDB is a database of gene sets used for GSEA (18). |Normalized enrichment score (NES)| ≥ 1, P ≤ 0.05, and FDR ≤ 0.25 were selected as the cutoff criteria for statistically significant differences. ES is enrichment score, NES is the normalized ES value after correction, P indicates the confidence of enrichment results, and FDR is an estimate of the probability of false-positive results for NES, so the smaller the FDR, the more significant the enrichment (19, 20). The MSigDB gene set includes nine major collections (H:C8). C2 (curated gene sets), C5 (ontology gene sets), and C7 (immunological signature gene sets) were the target datasets for this study. The R code was obtained from https://github.com/dongyihe/rheumatoidarthritis.
Gene Ontology Functional Enrichment and the Kyoto Encyclopedia of Genes and Genomes Pathway Analysis
The DEGs were annotated using the GO functional enrichment analysis, which included biological process (BP), molecular function (MF), and cellular component (CC) domains, and the KEGG pathway analysis (21). The KEGG is a database that provides information regarding gene functions at the molecular and higher levels, including biochemical pathways (22). Annotation and visualization were performed using the clusterProfiler package (23) (an R package for comparing biological themes among gene clusters). Enrichment analysis was performed using the hypergeometric test. Adjusted P < 0.05 was selected as the cutoff criterion indicating a statistically significant difference. P was adjusted using FDR.
PPI Network Construction and Identification of Hub Genes
The protein–protein interaction network was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING), a database that provides all the exposed PPI (24). The minimum required interaction score had the highest confidence (0.900). Hub genes were screened and visualized using the Molecular Complex Detection (MCODE) and CytoHubba plugins in Cytoscape version 3.7.2 for the visualization, modeling, and analyses of the molecular and genetic interaction networks (25). The MCODE and CytoHubba plugins can identify hub genes from complex interaction networks and help to lock the hub genes in a computationally efficient manner (26).
Validation of Hub Gene Expression by Quantitative Real-Time PCR
To validate the reliability of RNA-seq analysis in identifying the DEGs and to determine the expression levels of the 10 selected hub genes, quantitative real-time PCR (qRT-PCR) was conducted. Total RNA was extracted from 9 RA and 15 OA synovial tissue samples using TRIzol Reagent (Thermo Fisher Scientific) and reverse-transcribed to complementary DNA (cDNA) using the PrimeScript™ RT Master Mix (Perfect Real Time) (Takara Bio Incorporation, Beijing, China). A qRT-PCR was then performed using the TB Green® Premix Ex Taq™ (Tli RNase H Plus) (Takara Bio Incorporation). β-actin gene was used as the internal reference. The relative messenger RNA (mRNA) expression was calculated using the 2−ΔΔCt method. The Mann–Whitney U test was used for statistical analysis and P < 0.05 indicated a significant difference.
Results
Identification of the Differentially Expressed Genes in the Rheumatoid Arthritis Synovial Tissue Samples Using RNA Sequencing
Sequencing data comprised 17,736 genes, of which 851 genes were identified as the DEGs, as they met the following threshold criteria: adjusted P < 0.05 and |log2FoldChange| ≥ 1.5. We identified 474 upregulated and 377 downregulated genes. The principal component analysis (PCA) plot of the samples is shown in Figure 1A and the volcano plot and heat map of the DEGs are shown in Figures 1B,C, respectively. Details of the top 30 DEGs are given in Table 2, details of the DEGs encoding chemokines are given in Table 3, and the information about all the DEGs is shown in Supplementary Table 2.
Figure 1. Principal component analysis (PCA) plot of samples and volcano plot and heat map of the differentially expressed genes (DEGs). (A) PCA plot of samples and (B) Volcano plot of DEGs. The red dots represent upregulated genes and the blue dots represent downregulated genes. The tagged genes are the top 30 genes with a false discovery rate (FDR) < 0.05. (C) Heat map of DEGs.
Table 2. Information of the top 30 differentially expressed genes (DEGs).
Table 3. Information of the differentially expressed chemokine genes.
Identification of RA-Related Gene Signatures Using GSEA
The first part is the enrichment score line graph: the score at the highest peak is the ES value of the gene set. In the second part, the black lines represent gene positions in the sorted gene table. The leading edge subset is the part of the genes corresponding to the origin to the peak ES of the green curve, indicating the genes that have a major contribution to the enrichment. The third part is the distribution of the rank values of all the genes after sorting. The genes corresponding to the red part of the heat map are highly expressed in RA, the genes corresponding to the blue part are highly expressed in OA, and the signal-to-noise ratio corresponding to each gene is shown in the gray area map. In C3 (curated gene sets), the significantly enriched gene sets positively correlated with the RA group were CD40 signaling up (NES = 2.38, FDR = 6.16 × 10−9) and Th1 cytotoxic module (NES = 2.50, FDR = 6.16 × 10−9) (Figure 2A). In C5 (the Ontology Gene sets), the significantly enriched gene sets positively correlated with the RA group were activation of immune response (NES = 2.15, FDR = 1.23 × 10−9) and adaptive immune response (NES = 2.62, FDR = 1.23 × 10−9) (Figure 2B). In C7, the upregulation of the gene set of effective vs. memory CD8+ T cell is related to RA-related genes (NES = 2.13, FDR = 3.17 × 10−9). The downregulation of the gene set of naïve vs. effective CD8+ T cell is related to RA-related genes (NES = 2.39, FDR = 3.17 × 10−9) (Figure 2C).
Figure 2. GSEA plot of all the detected genes. (A) C2, curated gene sets; (B) C5, the Ontology Gene sets; and (C) C7, immunologic signature gene sets. GSEA, gene set enrichment analysis; ES, enrichment score; NES, normalized enrichment score. The green line means enrichment profile. The red part of the heat map represents the genes that are highly expressed in rheumatoid arthritis (RA), the blue part of the heat map represents the genes that are highly expressed in osteoarthritis (OA), and the gray area of the heat map represents the signal-to-noise ratio of each gene.
Gene Ontology Functional Enrichment Analysis and the Kyoto Encyclopedia of Genes and Genomes Pathway Analysis
The gene ontology functional enrichment analysis revealed that the DEGs were enriched in the signal transduction (P = 3.01 × 10−6), immune response (P = 1.65 × 10−24), and inflammatory response (P = 5.76 × 10−10) functions of the BP domain. In the MF domain, the DEGs were enriched in the calcium ion binding (P = 1.26 × 10−5), receptor binding (P = 1.26 × 10−5), and cytokine activity (P = 2.01 × 10−3) functions. In the CC domain, the DEGs were enriched in the plasma membrane (P = 1.91 × 10−31), integral component of membrane (P = 7.39 × 10−13), and extracellular region (P = 7.63 × 10−11) functions (Figures 3A,B). Detailed information is given in Supplementary Table 3. The KEGG pathway analysis revealed that the DEGs were enriched in the cytokine–cytokine receptor interaction (P = 3.05 × 10−17), chemokine signaling (P = 3.50 × 10−7), T-cell receptor signaling (P = 5.17 × 10−4), and rheumatoid arthritis (P = 5.17 × 10−4) pathways (Figures 3C,D). The important pathways involved in the cytokine–cytokine receptor interaction and RA are shown in Figure 4. Detailed information is given in Supplementary Table 4.
Figure 3. The Gene Ontology (GO) functional enrichment and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of the DEGs. (A) Dot plot of the GO functional enrichment analysis (the top five terms of each domain); (B) bar plot of the GO functional enrichment analysis (the top five terms of each domain); (C) dot plot of the KEGG pathway analysis; and (D) bar plot of the KEGG pathway analysis. The size of the dots and the height of the histogram represent the number of enriched genes and the color represents the P-value. The bigger the dot, the higher the column, and the more genes are enriched. The darker the color, the smaller the P-value.
Figure 4. The two most important pathways identified through the KEGG pathway analysis. (A) The cytokine–cytokine receptor interaction pathway and (B) the rheumatoid arthritis pathway. The red boxes represent upregulated genes and the dark green boxes represent downregulated genes.
PPI Network Development and Identification of Hub Genes in Rheumatoid Arthritis
Using the STRING database, our analysis produced 833 nodes and 1,639 edges and the PPI enrichment P-value was 1.0 × 10−16. Using the MCODE plugin in Cytoscape, 30 modules were identified. The important five modules are given in Figure 5. Using both the MCODE and cytoHubba plugins in Cytoscape, we identified C-X-C motif chemokine ligand (CXCL) 13, CXCL6, CXCL3, CXCL10, C-C motif chemokine receptor (CCR) 5, CCR2, CCR7, C-X-C motif chemokine receptor 5 (CXCR5), somatostatin receptor (SSTR) 1, and SSTR3 genes as the hub genes.
Figure 5. The PPI network and hub genes (top five modules). The significant modules identified by the Molecular Complex Detection (MCODE) and the cytoHubba plugin of Cytoscape. The size and color of circles represent the degree of the hub gene. The darker the color, the bigger the circle, and the higher the degree of the core gene. DEG, differentially expressed gene; PPI, protein–protein interaction.
Validation of Rheumatoid Arthritis-Related Hub Genes
The expression levels of the 10 selected hub genes in the synovial tissue samples from patients with RA and OA were validated by qRT-PCR. The primer sequences used for this experiment are given in Table 4. Statistical analysis showed that the expression levels of CXCL13 (P < 0.0001), CXCL6 (P = 0.0252), CCR5 (P = 0.0002), CXCR5 (P = 0.0033), CCR2 (P = 0.0073), CXCL3 (P = 0.0314), and CXCL10 (P < 0.0001) in the RA synovial tissue samples were significantly higher than those in the OA synovial tissue samples, while the expression level of SSTR1 (P = 0.0486) was significantly higher in the OA synovial tissue samples than in the RA synovial tissue samples (Figure 6).
Table 4. Gene primer sequences.
Figure 6. Validation of the 10 hub DEGs identified in the synovial tissue samples of patients with RA and OA by quantitative real-time PCR (qRT-PCR). The relative expression levels of each gene were calculated using the 2−ΔΔCt method. ****, P < 0.0001, ***, P < 0.005, **, P < 0.01, *, P < 0.05.
Discussion
Rheumatoid arthritis and OA are two common types of arthritis that present with inflammation but have distinct etiologies, clinical trajectories, and treatments. The pathogenesis and manifestations of these two diseases are complex, with clinical heterogeneity in presentation and disease course (8, 27). Distinguishing between RA and OA is critically important for early diagnosis, appropriate treatment, and elucidation of the underlying pathophysiology of these disorders. Previous studies have demonstrated that synovial tissue plays an important role in the occurrence and development of RA and OA. Previous studies mostly used microarray technology to study the differential expression profile data of RA and OA, which played an important role in the study of RA. However, microarray technology can only detect known sequences, while sequencing technology can detect unknown sequences and discover unknown genes. Compared with first-generation sequencing technology, second-generation sequencing technology has high throughput and high sensitivity and may discover new disease-causing genes. In this study, second-generation sequencing was used to analyze the synovial tissue samples obtained from patients with RA and OA. We identified the DEGs between the two samples, analyzed the functions and pathways of the DEGs, and validated the hub DEGs by qRT-PCR.
In previous studies, synovial tissue datasets in the Gene Expression Omnibus (GEO) database, such as GSE55235 and GSE12021, were used for bioinformatics analysis of RA. These studies were based on different datasets and identified genes (28, 29). Li et al. (30) found that vascular endothelial growth factor A and epidermal growth factor receptor may have essential roles in the development of RA and can be used as potential biomarkers of RA. Ren et al. (29) suggested that a set of eight genes (CCR5, CCL5, CXCL9, CXCL10, CXCL13, PNOC, TLR8, and CD52) can be used to diagnose RA with excellent specificity and sensitivity. To further analyze the transcriptome of the synovial tissue samples of patients with RA, we collected the samples from patients with RA and OA present at the same rheumatology hospital, performed RNA-seq, and investigated the pathways, gene networks, and hub genes. In this study, OA was used as the control group to study RA biomarkers. We performed transcriptomic analysis of the synovial tissue samples of patients with RA and OA using RNA-seq and determined that the identified DEGs can be used as biomarkers for diagnosing the two diseases. Further study using these biomarkers should be conducted.
In this study, 17,736 genes were identified. GSEA is sensitive in detecting genes with relatively small fold changes (31). The significantly enriched curated gene sets that positively correlated with RA were CD40 signaling and Th1 cytotoxic module. CD40 signaling is associated with the production of human rheumatoid factor (32) and the CD40/nuclear factor-kappa B (NF-kB) signaling pathway plays an important role in RA pathogenesis (33). The Th1 cytotoxic module has not been reported to be related to RA, but Th1 cytotoxicity is reportedly associated with the tumor microenvironment (34). The significantly enriched Ontology Gene sets that positively correlated with the RA group were involved in the activation of the immune response and adaptive immune response. RA is an autoimmune disease that affects both innate and adaptive immunities (35). The significantly enriched immunologic signature gene sets that positively correlated with the RA group were effective vs. memory CD8+ T cells (upregulated) and naïve vs. effective CD8+ T cells (downregulated). CD8+ T cells are involved in the pathogenesis of many autoimmune diseases, mainly because of their self-reactive cytotoxic inflammatory behavior (36). Effective CD8+ T cells have proliferative and cytotoxic properties and induce the death of infected cells and effective memory CD8+ T cells have a lower ability to induce cytotoxicity than effective CD8+ T cells (36, 37).
In this study, 851 DEGs were identified, of which 474 DEGs were upregulated and 377 DEGs were downregulated. The GO functional enrichment analysis revealed that the DEGs were enriched in signal transduction, immune response, and inflammatory response (BP domain); in calcium ion binding, receptor binding, and chemokine activity (MF domain); and in the plasma membrane, an integral component of membrane, and extracellular region (CC domain). The KEGG pathway analysis showed that the DEGs were enriched in the cytokine–cytokine receptor interactions, chemokine signaling, T-cell receptor signaling, and RA pathways. The DEGs were mainly concentrated in immune and inflammation-related pathways.
Ten DEGs were identified as hub genes using the MCODE and the cytoHubba plugin of Cytoscape. According to the qRT-PCR validation, the expression levels of CXCL13, CXCL6, CCR5, CXCR5, CCR2, CXCL3, and CXCL10 in the RA synovial tissue samples were higher than that in OA synovial tissue samples, while the expression of SSTR1 showed the opposite trend. The expression levels of CCR7 and SSTR3 did not differ between the RA and OA synovial tissue samples. CXCL13, CXCL10, CXCL6, and CXCL3 are the main members of the chemokine subfamily CXC. CXCL13, a B-cell chemokine, interacts with its receptor CXCR5 to promote the migration and aggregation of B lymphocytes (38). The expression level of CXCL13 in the serum of patients with RA is positively correlated with the level of rheumatoid factor and with disease activity and treatment response in early RA (39–41). CXCL10 is a ligand for CXCR3, which may stimulate the migration of monocytes, natural killer cells, and T cells (42). The expression of CXCL10 has been detected in the serum, synovial fluid, and synovial tissue of patients with RA (43, 44). Therefore, CXCL10 could act as a disease activity marker in early RA because of its high level in the plasma of untreated early patients with RA and its association with clinical disease activity (45). This study confirmed the high expression level of CXCL10 in the RA synovial tissue samples, which was significantly higher than that in the OA synovial tissue samples. CXCL3 is associated with the invasion and metastasis of various cancers (46–48) and CXCL3 and CXCL6 are involved in the invasion and migration of various cancers (49–51). The differential expression of CXCL3 and CXCL6 in the RA and OA synovial tissue samples is not yet reported. CCR7, CCR5, and CCR2 are chemokine receptors. CCR5 is expressed in RA synovial tissue and in T-helper cell type 1 inflammatory infiltrates. The Delta32 allelic variant of CCR5 has been reported to have a protective effect on RA susceptibility (52); however, the effect of CCR5 inhibitors on RA remains controversial (53–55). CCR2 has been widely considered as a potential therapeutic target for RA and CCR2 blocking agents have been developed (56). Monocyte chemoattractant protein 1 (CCL2) and its high-affinity receptor CCR2 are central to the development of pain associated with knee OA. Thus, CCR2 plays an important role in both the RA and OA. This study found that the expression levels of CCR2 in the RA and OA synovial tissue samples were different, which are likely related to its different functions in RA and OA pathogenesis. Somatostatins can regulate diverse cellular functions such as neurotransmission and cell proliferation. SSTR1 is associated with various cancers, such as prostate cancer (57) and gastric cancer (58). The role of SSTR1 in RA and OA has not yet been studied and the present results may provide a basis for future study on arthritis.
Osteoarthritis referred to degenerative joint disease and RA referred to joint disease caused by immune disorders. OA was used as the control group for RA, which had advantages but also limitations. The main limitation of this study was that inflammation was not properly investigated. RA is characterized by persistent synovitis and systemic inflammation and in the course of the development of OA, synovial inflammation is also observed. Although study on systemic inflammation in OA remains controversial, RA and OA have different mechanisms of inflammation as elucidated by the present results. These different mechanisms will be the focus of our future study.
In this study, RNA-seq technology was used to supplement the previous microarray technology and qRT-PCR technology was used to supplement and verify the previous conclusions. At the same time, bioinformatics technology was combined with experimental technology to make the results more reliable and provide a reliable preliminary basis for future study.
Conclusion
Ribonucleic acid sequencing was used to detect differential gene expression in the RA and OA synovial tissue samples. Using bioinformatics, the DEGs were identified in the RA and OA synovial tissue samples and the GO functional and the KEGG pathway enrichment analyses of the DEGs were performed. The hub DEGs such as SSTR1, CXCR5, CXCL6, CXCL3, CXCL13, CXCL10, CCR7, and CCR2 were validated by qRT-PCR. This study enriched the expression profile data of the DEGs in the synovial tissue of patients with RA and OA and provides novel insights into the differences between RA and OA. The candidate DEG pathways might be therapeutic targets and biomarkers for RA or OA.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://github.com/dongyihe/rheumatoidarthritis.
Ethics Statement
The studies involving human participants were reviewed and approved by Shanghai Guanghua Hospital. The patients/participants provided their written informed consent to participate in this study.
Author Contributions
DH and SG contributed to the conception, design, and final approval of the manuscript. CC, LX, and YB contributed to sequencing data and statistical analyses. YSh, YSu, SS, and SJS collected the samples, helped with statistical analysis, and drafted the manuscript. The final manuscript was written by RZ and YJ. All the authors have read and approved the final version of the manuscript.
Funding
This work was funded by the National Natural Science Funds of China (82074234 and 82004166); Shanghai Chinese Medicine Development Office, National Administration of Traditional Chinese Medicine, Regional Chinese Medicine (Specialist) Diagnosis and Treatment Center Construction Project-Rheumatology; State Administration of Traditional Chinese Medicine, National TCM Evidence-Based Medicine Research and Construction Project, Basic TCM Evidence-Based Capacity Development Program; Shanghai Municipal Health Commission, East China Region based Chinese and Western Medicine Joint Disease Specialist Alliance; National Key Research and Development Project (2018YFC1705203).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher's Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmed.2022.799440/full#supplementary-material
References
1. McInnes IB. The pathogenesis of rheumatoid arthritis. N Engl J Med. (2011) 365:2205–19. doi: 10.1056/NEJMra1004965
3. Croia C, Bursi R, Sutera D, Petrelli F, Alunno A, Puxeddu I. Review one year in review 2019: pathogenesis of rheumatoid arthritis. Clin Exp Rheumatol. (2019) 37:347–57.
4. Guo S, Xu L, Chang C, Zhang R, Jin Y, He D. Epigenetic regulation mediated by methylation in the pathogenesis and precision medicine of rheumatoid arthritis. Front Genet. (2020) 11:811. doi: 10.3389/fgene.2020.00811
5. Buckwalter J, Martin J. Osteoarthritis. Adv Drug Deliv Rev. (2006) 58:150–67. doi: 10.1016/j.addr.2006.01.006
6. Felson DT. Osteoarthritis: new insights. Part 1: the disease and its risk factors. Ann Intern Med. (2000) 133:635. doi: 10.7326/0003-4819-133-8-200010170-00016
7. Scott DL, Wolfe F, Huizinga TW. Rheumatoid arthritis. Lancet. (2010) 376:1094–108. doi: 10.1016/S0140-6736(10)60826-4
8. Glyn-Jones S, Palmer AJR, Agricola R, Price AJ, Vincent TL, Weinans H, et al. Osteoarthritis. Lancet. (2015) 386:376–87. doi: 10.1016/S0140-6736(14)60802-3
9. Song X, Lin Q. Genomics, transcriptomics and proteomics to elucidate the pathogenesis of rheumatoid arthritis. Rheumatol Int. (2017) 37:1257–65. doi: 10.1007/s00296-017-3732-3
10. Aletaha D, Neogi T, Silman AJ, Funovits J, Felson DT, Bingham CO, et al. 2010 Rheumatoid arthritis classification criteria: an American College of Rheumatology/European League Against Rheumatism collaborative initiative. Ann Rheum Dis. (2010) 69:1580–8. doi: 10.1136/ard.2010.138461
11. Altman R, Asch E, Bloch D, Bole G, Borenstein D, Brandt K, et al. Development of criteria for the classification and reporting of osteoarthritis: classification of osteoarthritis of the knee. Arthritis Rheum. (1986) 29:1039–49. doi: 10.1002/art.1780290816
12. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. (2014) 30:2114–20. doi: 10.1093/bioinformatics/btu170
13. Kim D, Langmead B, Salzberg SL. HISAT a fast spliced aligner with low memory requirements. Nat Methods. (2015) 12:357–60. doi: 10.1038/nmeth.3317
14. Roberts A, Trapnell C, Donaghey J, Rinn JL, Pachter L. Improving RNA-Seq expression estimates by correcting for fragment bias. Genome Biol. (2011) 12:R22. doi: 10.1186/gb-2011-12-3-r22
15. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. (2010) 28:511–5. doi: 10.1038/nbt.1621
16. Anders S, Pyl PT, Huber W. HTSeq – a python framework to work with high-throughput sequencing data. Bioinformatics. (2015) 31:166–9. doi: 10.1093/bioinformatics/btu638
17. Simon Anders WH. Differential expression of RNA-Seq data at the gene level – the DESeq package. EMBL. (2012).
18. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Tamayo P. The molecular signatures database (MSigDB) hallmark gene set collection. Cell Syst. (2015) 1:417–25. doi: 10.1016/j.cels.2015.12.004
19. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci. (2005) 102:15545–50. doi: 10.1073/pnas.0506580102
20. Mootha VK, Lindgren CM, Eriksson K-F, Subramanian A, Sihag S, Lehar J, et al. PGC-1α-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet. (2003) 34:267–73. doi: 10.1038/ng1180
21. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene Ontology: tool for the unification of biology. Nat Genet. (2000) 25:25–9. doi: 10.1038/75556
22. Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K. KEGG new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. (2017) 45:D353–61. doi: 10.1093/nar/gkw1092
23. Yu G, Wang L-G, Han Y, He Q-Y. Clusterprofiler: an R package for comparing biological themes among gene clusters. OMICS J Integr Biol. (2012) 16:284–7. doi: 10.1089/omi.2011.0118
24. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. (2019) 47:D607–13. doi: 10.1093/nar/gky1131
25. Shannon P. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303
26. Chin C-H, Chen S-H, Wu H-H, Ho C-W, Ko M-T, Lin C-Y. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. (2014) 8:S11. doi: 10.1186/1752-0509-8-S4-S11
27. Zhao J, Guo S, Schrodi SJ, He D. Molecular and cellular heterogeneity in rheumatoid arthritis: mechanisms and clinical implications. Front Immunol. (2021) 12:790122. doi: 10.3389/fimmu.2021.790122
28. Li Z, Xu M, Li R, Zhu Z, Liu Y, Du Z, et al. Identification of biomarkers associated with synovitis in rheumatoid arthritis by bioinformatics analyses. Biosci Rep. (2020) 40:BSR20201713. doi: 10.1042/BSR20201713
29. Ren C, Li M, Du W, Lü J, Zheng Y, Xu H, et al. Comprehensive bioinformatics analysis reveals hub genes and inflammation state of rheumatoid arthritis. BioMed Res Int. (2020) 2020:1–13. doi: 10.1155/2020/6943103
30. Li H, Yang HH, Sun ZG, Tang HB, Min JK. Whole-transcriptome sequencing of knee joint cartilage from osteoarthritis patients. Bone Jt Res. (2019) 8:290–303. doi: 10.1302/2046-3758.87.BJR-2018-0297.R1
31. Zeng M, Liu J, Yang W, Zhang S, Liu F, Dong Z, et al. Multiple-microarray analysis for identification of hub genes involved in tubulointerstial injury in diabetic nephropathy. J Cell Physiol. (2019) 234:16447–62. doi: 10.1002/jcp.28313
32. Kyburz D, Corr M, Brinson DC, Von Damm A, Tighe H, Carson DA. Human rheumatoid factor production is dependent on CD40 signaling and autoantigen. J Immunol Baltim Md 1950. (1999) 163:3116–22.
33. Criswell LA. Gene discovery in rheumatoid arthritis highlights the CD40/ NF-jB signaling pathway in disease pathogenesis. Immunol Rev. (2010) 233:55–61. doi: 10.1111/j.0105-2896.2009.00862.x
34. Tosolini M, Kirilovsky A, Mlecnik B, Fredriksen T, Mauger S, Bindea G, et al. Clinical impact of different classes of infiltrating T cytotoxic and helper cells (Th1, Th2, Treg, Th17) in patients with colorectal cancer. Cancer Res. (2011) 71:1263–71. doi: 10.1158/0008-5472.CAN-10-2907
35. O'Neil LJ, Kaplan MJ. Neutrophils in rheumatoid arthritis: breaking immune tolerance and fueling disease. Trends Mol Med. (2019) 25:215–27. doi: 10.1016/j.molmed.2018.12.008
36. Carvalheiro H, da Silva JAP, Souto-Carneiro MM. Potential roles for CD8+ T cells in rheumatoid arthritis. Autoimmun Rev. (2013) 12:401–9. doi: 10.1016/j.autrev.2012.07.011
37. Tomiyama H, Matsuda T, Takiguchi M. Differentiation of human CD8 + T cells from a memory to memory/effector phenotype. J Immunol. (2002) 168:5538–50. doi: 10.4049/jimmunol.168.11.5538
38. Bao Y-Q, Wang J-P, Dai Z-W, Mao Y-M, Wu J, Guo H-S, et al. Increased circulating CXCL13 levels in systemic lupus erythematosus and rheumatoid arthritis: a meta-analysis. Clin Rheumatol. (2020) 39:281–90. doi: 10.1007/s10067-019-04775-z
39. Greisen SR, Schelde KK, Rasmussen TK, Kragstrup TW, Stengaard-Pedersen K, Hetland ML, et al. CXCL13 predicts disease activity in early rheumatoid arthritis and could be an indicator of the therapeutic ‘window of opportunity’. Arthritis Res Ther. (2014) 16:434. doi: 10.1186/s13075-014-0434-z
40. Jones JD, Hamilton B, Challener GJ, de Brum-Fernandes AJ, Cossette P, Liang P, et al. Serum C-X-C motif chemokine 13 is elevated in early and established rheumatoid arthritis and correlates with rheumatoid factor levels. Arthritis Res Ther. (2014) 16:R103. doi: 10.1186/ar4552
41. Allam SI, Sallam RA, Elghannam DM, El-Ghaweet AI. Clinical significance of serum B cell chemokine (CXCL13) in early rheumatoid arthritis patients. Egypt Rheumatol. (2019) 41:11–4. doi: 10.1016/j.ejr.2018.04.003
42. Karin N, Razon H. Chemokines beyond chemo-attraction: CXCL10 and its significant role in cancer and autoimmunity. Cytokine. (2018) 109:24–8. doi: 10.1016/j.cyto.2018.02.012
43. Patel DD, Zachariah JP, Whichard LP. CXCR3 and CCR5 ligands in rheumatoid arthritis synovium. Clin Immunol. (2001) 98:39–45. doi: 10.1006/clim.2000.4957
44. Hanaoka R, Kasama T, Muramatsu M, Yajima N, Shiozawa F, Miwa Y, et al. A novel mechanism for the regulation of IFN-γ inducible protein-10 expression in rheumatoid arthritis. Arthritis Res Ther. (2003) 5:R74. doi: 10.1186/ar616
45. Pandya JM, Lundell A-C, Andersson K, Nordström I, Theander E, Rudin A. Blood chemokine profile in untreated early rheumatoid arthritis: CXCL10 as a disease activity marker. Arthritis Res Ther. (2017) 19:20. doi: 10.1186/s13075-017-1224-1
46. Wang H, Wang T, Dai L, Cao W, Ye L, Gao L, et al. Effects of CXCL3 on migration, invasion, proliferation and tube formation of trophoblast cells. Placenta. (2018) 66:47–56. doi: 10.1016/j.placenta.2018.05.004
47. Xin H, Cao Y, Shao M-L, Zhang W, Zhang C-B, Wang J-T, et al. Chemokine CXCL3 mediates prostate cancer cells proliferation, migration and gene expression changes in an autocrine/paracrine fashion. Int Urol Nephrol. (2018) 50:861–8. doi: 10.1007/s11255-018-1818-9
48. Zhao Q-Q, Jiang C, Gao Q, Zhang Y-Y, Wang G, Chen X-P, et al. Gene expression and methylation profiles identified CXCL3 and CXCL8 as key genes for diagnosis and prognosis of colon adenocarcinoma. J Cell Physiol. (2020) 235:4902–12. doi: 10.1002/jcp.29368
49. Li J, Tang Z, Wang H, Wu W, Zhou F, Ke H, et al. CXCL6 promotes non-small cell lung cancer cell survival and metastasis via down-regulation of miR-515-5p. Biomed Pharmacother Biomedecine Pharmacother. (2018) 97:1182–8. doi: 10.1016/j.biopha.2017.11.004
50. Liu G, An L, Zhang H, Du P, Sheng Y. Activation of CXCL6/CXCR1/2 axis promotes the growth and metastasis of osteosarcoma cells in vitro and in vivo. Front Pharmacol. (2019) 10:307. doi: 10.3389/fphar.2019.00307
51. Zheng S, Shen T, Liu Q, Liu T, Tuerxun A, Zhang Q, et al. CXCL6 fuels the growth and metastases of esophageal squamous cell carcinoma cells both in vitro and in vivo through upregulation of PD-L1 via activation of STAT3 pathway. J Cell Physiol. (2021) 236:5373–86 doi: 10.1002/jcp.30236
52. Pokorny V. Evidence for negative association of the chemokine receptor CCR5 d32 polymorphism with rheumatoid arthritis. Ann Rheum Dis. (2004) 64:487–90. doi: 10.1136/ard.2004.023333
53. van Kuijk AWR, Vergunst CE, Gerlag DM, Bresnihan B, Gomez-Reino JJ, Rouzier R, et al. CCR5 blockade in rheumatoid arthritis: a randomised, double-blind, placebo-controlled clinical trial. Ann Rheum Dis. (2010) 69:2013–6. doi: 10.1136/ard.2010.131235
54. Takeuchi T, Kameda H. What is the future of CCR5 antagonists in rheumatoid arthritis? Arthritis Res Ther. (2012) 14:114. doi: 10.1186/ar3775
55. Lan Y, Wang Y, Liu Y. CCR5 silencing reduces inflammatory response, inhibits viability, and promotes apoptosis of synovial cells in rat models of rheumatoid arthritis through the MAPK signaling pathway. J Cell Physiol. (2019) 234:18748–62. doi: 10.1002/jcp.28514
56. Quinones MP, Estrada CA, Kalkonde Y, Ahuja SK, Kuziel WA, Mack M, et al. The complex role of the chemokine receptor CCR2 in collagen-induced arthritis: implications for therapeutic targeting of CCR2 in rheumatoid arthritis. J Mol Med. (2005) 83:672–81. doi: 10.1007/s00109-005-0637-5
57. Pedraza-Arévalo S, Hormaechea-Agulla D, Gómez-Gómez E, Requena MJ, Selth LA, Gahete MD, et al. Somatostatin receptor subtype 1 as a potential diagnostic marker and therapeutic target in prostate cancer. Prostate. (2017) 77:1499–511. doi: 10.1002/pros.23426
Keywords: rheumatoid arthritis, osteoarthritis, synovial tissue, RNA-seq, differential gene expression
Citation: Zhang R, Jin Y, Chang C, Xu L, Bian Y, Shen Y, Sun Y, Sun S, Schrodi SJ, Guo S and He D (2022) RNA-seq and Network Analysis Reveal Unique Chemokine Activity Signatures in the Synovial Tissue of Patients With Rheumatoid Arthritis. Front. Med. 9:799440. doi: 10.3389/fmed.2022.799440
Received: 22 October 2021; Accepted: 16 March 2022;
Published: 04 May 2022.
Edited by:
Carl Kieran Orr, Saint Vincent's University Hospital, IrelandReviewed by:
Tejpal Gill, Oregon Health and Science University, United StatesSaikat Majumder, University of Pittsburgh, United States
Copyright © 2022 Zhang, Jin, Chang, Xu, Bian, Shen, Sun, Sun, Schrodi, Guo and He. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Shicheng Guo, Shicheng.Guo@wisc.edu; Dongyi He, hedongyi1967@shutcm.edu.cn
†These authors have contributed equally to this work