Drug screening and identification of key candidate genes and pathways of rheumatoid arthritis.

Rheumatoid arthritis (RA), which normally manifests as a multi‑joint inflammatory reaction, is a common immunological disease in clinical practice. However, the pathogenesis of RA has not yet been fully elucidated. Rituximab (RTX) is an effective drug in the treatment of RA, however its therapeutic efficacy and mechanism of action require further investigation. Thus, the present study aimed to screen the candidate key regulatory genes and explain the potential mechanisms of RA. Gene chips of RA and normal joint tissues were analyzed and, gene chips of RTX before and after treatment were investigated. In the present study, strong evidence supporting the pathogenesis of RA and mechanism of action of RTX were also revealed. Differentially expressed genes (DEGs) were analyzed using the limma package of RStudio software. A total of 1,150 DEGs were detected in RA compared with normal joint tissues. The upregulated genes were enriched in 'interleukin‑12 production', 'I‑κB kinase/NF‑κB signaling', 'regulation of cytokine production involved in immune response' and 'cytokine metabolic process'. Functional enrichment analysis showed that RTX was primarily involved in the inhibition of 'adaptive immune response', 'B cell activation involved in immune response' and 'immune effector process'. Subsequently, leukocyte immunoglobulin‑like receptor subfamily B member 1 (LILRB1), a hub gene with high connectivity degree, was selected, and traditional Chinese medicine libraries were molecularly screened according to the structure of the LILRB1 protein. The results indicated that kaempferol 3‑O‑β‑D‑glucosyl‑(1→2)‑β‑D‑glucoside exhibited the highest docking score. In the present study, the DEGs and their biological functions in RA and the pharmacological mechanism of RTX action were determined. Taken together, the results suggested that LILRB1 may be used as a molecular target for RA treatment, and kaempferol 3‑O‑β‑D‑glucosyl‑(1→2)‑β‑D‑glucoside may inhibit the pathological process of RA.

Introduction

Rheumatoid arthritis (RA) is a chronic systemic disease accompanied by inflammatory synovitis that is mainly characterized by symmetrical distribution of invasive joint inflammation of the hand and foot (1,2). In addition, RA exhibits increased interstitial inflammatory cell infiltration and bone tissue destruction, resulting in joint deformity and loss of function (3). Immune function is considered to be the main aspect associated with RA; RA is characterized by the induction of innate immune disorders, including immune complex-mediated complement activation, osteoclast and chondrocyte activation and cytokine network dysregulation, which develop semi-autonomous features that contribute to disease progression (4,5). However, the exact mechanism of RA development remains elusive and further investigation is required.

General, surgical and pharmaceutical therapies are widely applied in RA treatment (6). The most commonly used pharmacological RA drugs include the administration of non-steroidal anti-inflammatory drugs, immunosuppressants, botanicals and biological agents (7). Rituximab (RTX), a chimeric monoclonal antibody against the CD20 ligand of B lymphocytes, has been reported to exhibit therapeutic activity in the clinical treatment of RA (8); however, its therapeutic mechanism needs to be further investigated. Although several drugs alleviate pain in patients with RA, their efficacy is limited (9), therefore the development of novel and effective drugs for RA is required.

The present study aimed to further elucidate the pathogenesis of RA and identify potential drugs for RA treatment. The expression profiles of normal, RA control and RTX-treated tissues were analyzed. A series of immune-related genes, including leukocyte immunoglobulin-like receptor subfamily B member 1 (LILRB1), were detected by screening the differentially expressed genes (DEGs). The results revealed that LILRB1 was associated with RA pathogenesis. LILRB1, an inhibitory receptor broadly expressed in leukocytes, has been demonstrated to regulate immune responses by binding to MHC class I molecules on antigen-presenting cells (10). Finally, Traditional Chinese Medicine (TCM) libraries were molecularly screened for this key functional gene in order to identify potential therapeutic drugs.

Materials and methods

Download of expression profile chip data and DEGs analysis

The screening of DEGs (11,12) in the synovial tissues of normal patients without RA and patients with RA (GSE55235) (13) was performed using the Gene Expression Omnibus (GEO) database (14) and differential gene analysis. In addition, DEG screening in RA and RTX-treated patients (GSE24742) (15) was assessed using the GEO database and R, version 3.6.2. Data quality was determined by calculating residual sign, residuals, weight, relative log expression, normalized unscaled standard errors and RNA degradation. Finally, the differences in RNA expression profiles between groups were analyzed using the pheatmap and limma R packages (16,17). |Log2 fold-change (FC)|≥1 and P<0.05 were set as the cutoff criteria for DEGs.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses (18,19)

The functions of DEGs were analyzed using the ClueGO plug-in application in Cytoscape 3.6.1 (https://cytoscape.org). In addition, KEGG pathway enrichment analysis (20) was carried out using ClueGO (21) and visualized using CluePedia (22). P<0.05 was set as the cutoff value.

Gene set enrichment analysis (GSEA)

GSEA analysis was performed using the GSEA software (23). In brief, the method consisted of the following steps. First, the gene data, including expression profiles and class distinctions, were listed. Given a defined set of genes, the goal of GSEA was to determine whether the members of the gene set were found at the top or bottom of the list. Subsequently, an enrichment score was calculated in order to identify the degree of the over-represented genes. Finally, a weighted enrichment statistics was carried out and the number of random permutations was set to 1,000 times.

Pathway analysis

DEGs were subjected to signal pathway enrichment analysis using the ConsensusPathDB database (http://cpdb.molgen.mpg.de) (24,25). The files containing the measured genetic data were uploaded to the database and subsequently pathway enrichment analysis was carried out using the gene set analysis function.

Protein-protein interaction (PPI) analysis

The PPI network was analyzed using STRING software (https://string-db.org) (26). The organism selected was ‘Homo sapiens’. The interacting protein complexes that functionally influence the physiological processes of a disease and the PPI enrichment analysis reflected the interaction between DEGs.

TCM database and molecular docking simulation of LILRB1

A total of 32,364 compounds were obtained from the TCM database (http://tcm.cmu.edu.tw) (27). The conformational energy of small molecules was minimized using the Maestro software (version 11.8, Schrödinger, LLC) by adding hydrogen atoms and removing counter ions and salts (28,29). The crystal structure of LILRB1 was downloaded from Protein Data Bank (structure no. 1UGN; http://www.rcsb.org) (30) and the protein structure was refined by removing crystalline water and ions. In addition, energy minimization was performed on the LILRB1 protein structure. TCM docking and the selection of candidate compounds was performed using the virtual screening workflow model of Schrodinger and Glide XP (extra precision), respectively (31–33).

Results

Identification of DEGs in RA

The gene expression profiles of the synovial tissues of patients with RA from GSE55235 were obtained from the GEO database. The microarray data from GSE55235 included synovial tissues from 10 healthy and 10 RA joints. A total of 1,150 DEGs were extracted from the expression profile data set, including 508 downregulated and 642 upregulated genes. |Log2FC|≥1 and P<0.05 were set as the cutoff criteria for DEGs. The distribution of DEGs is presented in a volcano plot (Fig. 1A). In addition, hierarchical cluster analysis was performed in order to obtain an overview of the expression profiles of normal and RA tissues (Fig. 1B). Finally, all heat maps demonstrated different adjustment directions and significant separation between normal and RA samples.

Figure 1.

Analysis of differentially expressed genes between normal and rheumatoid arthritis joint synovial tissue samples. (A) Volcanic map of genes between normal and rheumatoid arthritis joint synovial tissue samples. Red and green dots indicate the upregulated and downregulated genes, respectively. The black dots correspond to gene expression with |log2FC|<1. The x-axis represents false discovery rate and the y-axis denotes the value of log2FC. (B) Cluster analysis results based on the expression profiles of differentially expressed genes. FC, fold-change.

GO and KEGG analyses of DEGs in RA

GO and KEGG enrichment analyses of DEGs were performed using the ClueGO plug-in in Cytoscape software, and subsequently the upregulated and downregulated genes were analyzed. The enriched upregulated genes were mainly associated with ‘interleukin-12 production’, ‘I-κB kinase/NF-κB signaling’, ‘regulation of cytokine production involved in immune response’, ‘positive regulation of B cell differentiation’, ‘cytokine metabolic process’ and ‘activation of immune response’ (Fig. 2A). The results of the pathway analysis showed that RA mainly affected the ‘Fcγ R-mediated phagocytosis’, ‘natural killer cell-mediated cytotoxicity’, ‘B cell receptor signaling pathway’, ‘NF-κB signaling pathway’ and ‘leukocyte transendothelial migration’ (Fig. 2B). By contrast, the downregulated genes were mainly involved in ‘vasculogenesis’, ‘regulation of DNA binding transcription factor activity’, ‘cellular response to external stimulus’, ‘vascular process in circulatory system’, ‘blood circulation’ and ‘transcription factor binding’ (Fig. 3A). Finally, the pathway enrichment analysis indicated that the downregulated DEGs were mainly enriched in ‘tyrosine metabolism’, ‘FoxO signaling pathway’, ‘regulation of lipolysis in adipocytes’ and ‘colorectal cancer’ (Fig. 3B).

Figure 2.

Enrichment analysis of upregulated genes. (A) GO enrichment analysis of upregulated genes. (B) KEGG enrichment analysis of upregulated genes. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 3.

Enrichment analysis of downregulated genes. (A) GO enrichment analysis of downregulated genes. (B) KEGG enrichment analysis of downregulated genes. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Identification of DEGs based on RTX treatment data

A total of 54,675 genes were obtained from 12 control and 12 RTX-treated samples, and 941 DEGs (382 upregulated and 559 downregulated) were identified. The volcano plot of DEGs is presented in Fig. 4A. The red, green and black dots indicate the upregulated, downregulated and non-differentiated genes, respectively. Finally, the overview of the expression profiles of DEGs before and after RTX treatment (Fig. 4B), and the distribution of genes between groups were revealed using a hierarchical cluster analysis.

Figure 4.

Analysis of differentially expressed genes between synovial tissue samples from patients with rheumatoid arthritis before and after treatment with RTX. (A) Volcano plots of genes between synovial tissue samples from patients with rheumatoid arthritis before and after treatment with RTX. The red and green dots indicate the upregulated and downregulated genes, respectively. The black dots correspond to gene expression with |log2FC|<1. The y-axis represents false discovery rate and the x-axis denotes the value of FC. (B) Hierarchical cluster analysis of expression profiles of differentially expressed genes before and after treatment with rituximab. RTX, rituximab; Ctrl, control; FC, fold-change.

GSEA and pathway enrichment analyses of DEGs after RTX treatment

GSEA was performed in order to further confirm the selected DEGs from the GO functional enrichment analysis (Fig. 5A). GSEA revealed that RTX upregulated the expression of genes associated with the ‘positive regulation of dendritic development’ and the ‘regulation of neuron migration’. By contrast, the ‘adaptive immune response’, ‘anaphase-promoting complex-dependent catabolic process’, ‘antigen processing and presentation’, ‘B cell activation involved in immune response’ and ‘immune effector process’ functions were suppressed. Additionally, pathway analysis was performed, using the online ConsensusPathDB database, in order to analyze the functional and signaling pathway enrichment of the gene signatures (Fig. 5B).

Figure 5.

Function and pathway analyses of differentially expressed genes between synovial tissue samples from patients with rheumatoid arthritis before and after treatment with RTX. (A) Gene set enrichment analysis of rituximab in rheumatoid arthritis treatment. (B) Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis of downregulated genes after treatment with rituximab. GO, Gene Ontology; GPCR, G protein-coupled receptor; IL, interleukin; JAK, Janus kinase; PKA, protein kinase A; Th17, T helper 17 cell; TNF, tumor necrosis factor.

The pathway enrichment analysis also revealed that the downregulated DEGs were enriched in the ‘PI3K-Akt signaling pathway’, ‘type II interferon signaling’, ‘Janus kinase/STAT signaling pathway’, ‘interleukin-6 signaling pathway’, ‘T helper 17 cell differentiation’, and ‘interleukin-4 and interleukin-13 signaling’. The aforementioned results supported the conclusion that RTX may treat RA via regulating body immunity.

Identification of key candidate genes using STRING protein interaction network

A Venn diagram analysis of the upregulated and downregulated genes in RA and after RTX treatment groups, respectively, was performed (34). The analysis identified 13 key genes that were subsequently analyzed using the STRING database (http://string-db.org) for PPI network analysis (Fig. 6A). The results indicated that LILRB1 exhibited the highest interactivity confidence (Fig. 6B).

Figure 6.

Identification of key genes. (A) Wayne analysis of upregulated genes in RA and downregulated genes after RTX treatment. (B) Protein-protein interaction network analysis of key genes. RA, rheumatoid arthritis; RTX, rituximab.

Screening of candidate compounds for RA treatment

The PPI network analysis indicated that LILRB1 was a key node gene associated with the mechanisms of RA pathogenesis. Therefore, the LILRB1 gene was selected for virtual drug screening. A set of molecular recognition strategies for TCM compounds identification was designed via structure-based high-throughput virtual screening. The filtering process is presented in Fig. 7A. Briefly, a total of 872 TCM compounds demonstrating the highest docking score with LILRB1 were screened using the high-throughput method, and 43 of them were selected. Subsequently, 5 candidate TCM compounds were obtained via precise docking. Among them, kaempferol 3-O-β-D-glucosyl-(1→2)-β-D-glucoside exhibited the highest binding capacity. The interaction between LILRB1 and kaempferol 3-O-β-D-glucosyl-(1→2)-β-D-glucoside is shown in Fig. 7B. The key amino acid residues Cys144, Arg72, Pro90, Val15 and Glu136 of the kaempferol 3-O-β-D-glucosyl-(1→2)-β-D-glucoside binding site interacted with LILRB1 via hydrogen bonds (Fig. 7B). The 3D conformation of the interaction and surface binding of kaempferol 3-O-β-D-glucosyl-(1→2)-β-D-glucoside and LILRB1 with other proteins were analyzed in order to identify additional protein interactions (Fig. 8A and B). The aforementioned results indicated that kaempferol 3-O-β-D-glucosyl-(1→2)-β-D-glucoside may decrease the biological activity of LILRB1 by inhibiting its active center, thereby exhibiting therapeutic effects on RA.

Figure 7.

Binding mode of receptor-ligand interaction. (A) Protocol flowchart of LILRB1 inhibitor discovery strategy. (B) 2D combination mode of LILRB1-kaempferol 3-O-β-D-glucosyl-(1→2)-β-D-glucoside. LILRB1, leukocyte immunoglobulin-like receptor subfamily B member 1; TCM, Traditional Chinese Medicine; HTV, high throughput virtual; SP, standard precision; XP, extra precision.

Figure 8.

Binding mode of LILRB1-kaempferol 3-O-β-D-glucosyl-(1→2)-β-D-glucoside interaction. (A) Protein surface binding mode of LILRB1-kaempferol 3-O-β-D-glucosyl-(1→2)-β-D-glucoside. (B) Cartoon combination mode of LIRNB1-kaempferol 3-O-β-D-glucosyl-(1→2)-β-D-glucoside. LILRB1, leukocyte immunoglobulin-like receptor subfamily B member 1.

Discussion

Immune dysregulation has been implicated in the pathogenesis of RA through the anti-immunoglobulin G antibodies (35). Although several mechanisms regarding the role of immune regulation in RA pathology have been proposed (36), the exact mechanism should be further investigated. Therefore, in the present study, the aim was to delineate the underlying molecular mechanisms and identify the key regulatory genes involved in the developmental process of RA.

As gene expression profiles may reflect disease status (37), in the present study the expression of DEGs was compared in normal, RA control and RTX-treated tissues selected from the GO database. The DEGs functional analysis indicated that RA was associated with the induction of inflammation responses and immune activities. These observations were consistent with the clinical characteristics of the disease. However, following RTX treatment, inflammation and immune-related signaling pathways were widely suppressed. These results were consistently with the previously reported inflammation- and immune-related characteristics of RA pathogenesis (3). The upregulated and downregulated genes in RA and RTX-treated groups, respectively, were subjected to Venn analysis and subsequently 13 selected genes were screened in order to identify the key regulatory genes. Topological analysis was conducted to screen the node genes and LILRB1, a receptor expressed on the membrane of immune cells (10). The analysis demonstrated that LILRB1 exhibited a high degree of connectivity. The results of the present study indicated that LILRB1 may control inflammatory responses and cytotoxicity by inducing a targeted immune response and limiting autoreactivity. In addition, it was hypothesized that upstream regulatory signaling pathways of the downregulated LILRB1 gene could serve as a target of RTX. However its pharmacological mechanisms of action should be further investigated. Thus, the results indicated that LILRB1 was a potential target of RA treatment.

In order to clarify the target activity of LILRB1, potential drugs that target LILRB1 were screened using the TCM libraries. Kaempferol 3-O-β-D-glucosyl-(1→2)-β-D-glucoside showed the strongest targeting activity. Kaempferol, a natural flavonol, has been reported to act as an antioxidant by reducing oxidative stress during the treatment of several diseases (38). However, the precise molecular mechanisms of kaempferol need to be further studied. The results of the present study revealed that kaempferol could bind to the active pocket of LILRB1, suggesting inhibition of the immune response signal transduction via the receptor. Thus, the activity of kaempferol, a potential anti-RA drug, may be mediated via targeting the immune-related receptor. Additional studies investigating the clinical and molecular basis are required to further elucidate the effect of LILRB1 in regulating the pathogenesis of RA and to demonstrate the therapeutic value of kaempferol in RA treatment.

In conclusion, this study may enhance the understanding of RA development processes, and suggested that kaempferol could be a potential novel and effective drug in RA treatment in clinical practice.