Identifying Key Genes and Functionally Enriched Pathways in B cells from Peripheral Blood of Rheumatoid Arthritis Patients by Bioinformatics Analysis and Validation

DOI: https://doi.org/10.21203/rs.3.rs-1382093/v1

Abstract

Rheumatoid arthritis is a common chronic inflammatory disease. B cells are critical for the development of RA, but so far, we still know little about the role of B cells in the pathogenesis of RA. In the present study, we downloaded mRNA expression profiles of B cells from the peripheral blood (BCPB) of 7 RA patients and 9 healthy people from the Gene Expression Omnibus (GEO) database. Using weighted gene co-expression network analysis (WGCNA) in conjunction with differentially expressed gene (DEGs) analysis, we identified 23 RA-trait different expression genes. In the subsequent protein protein interaction (PPI) network constructed by RA-trait DEGs, we eventually identified 5 hub genes (TLR8, CCR2, CX3CR1, GIMAP6 and C1QA). Compared with healthy people, these hub genes are significantly highly expressed in BCPB of RA patients in the dataset GSE4588. qRT-PCR results showed TLR8, CCR2, CX3CR1 were compatible with our bioinformatics analyses. Combined with other biological interpretations, our current study might provide a basis for identifying key genes and potential pathways related to the pathogenesis of BCPB of RA patients.

Introduction

Rheumatoid arthritis (RA) is a common rheumatic disease clinically characterized by synovial tissue hyperplasia, pannus formation, and bone destruction[1]. The primary feature of RA is the malfunction of the body's immune system, in which pathogenic T and B cells, autoantibodies, and inflammatory cytokines are aberrantly elevated in RA patients' bodies. B cells depletion therapy has fully demonstrated B cells are crucial in the pathogenesis of RA. Relevant data have suggested that B cells play a variety of roles in RA including as an antigen presenting cell to mediate the body's inflammatory response[2]; increasing the body's bone destruction by secreting of RANKL[3]; producing IL-6, IFN-γ and other inflammatory factors, and directly secreting RF, ACPA and other autoantibodies to trigger and aggravate the body's inflammatory response[4, 5]. New research shows that amphiregulin (AREG) produced by B cells in RA patients can promote the migration and proliferation of synovial fibroblasts (FLS), and the combination of AREG and ACPA will further increase the differentiation of osteoclasts[6]. However, so far, we have not fully understood the pathogenic role of B cells in RA.

WGCNA is a classical biology method to reveal gene association patterns between different biological samples[7]. WGCNA can establish a link between gene expression and clinical manifestations to screen out key genes or therapeutic targets[8]. According to reports, WGCNA has been widely applied for identifying hub genes and key signal pathways of autoimmune diseases such as Sjogren’s syndrome[9], systemic lupus erythematosus[10], and type I diabetes[11]. Aiming to identify hub genes linked to BCPB of RA patients and understand the pathogenesis of B cells in RA, we deeply analyze the characteristic genes of BCPB of RA patients by WGCNA and DEGs analysis. Finally, we explored potential functions of RA-trait DEGs, identified hub genes and find the gene sets associated with a single hub gene.

Methods

GEO Microarray

Gene expression profiles of BCPB of 7 RA patients and 9 controls were downloaded from the GEO database. GSE4588 was used to establish a co-expression network and perform follow-up analysis. GSE4588 was log2-transformed and normalized independently using "normalizeBetweenArrays" function in "limma" R package.

Construction of Co-Expression Network

The R package "WGCNA" was used to establish co-expression network by using the first 25% of the genes with high variation in GSE4588, which included about 5000 genes. When the correlation coefficient threshold was set at 0.9 and the minimum genes' number in each module was set at 150, we chose 7 as the soft-thresholding power. Besides, we set the cut height threshold at 0.5 to allow for the merging of potentially similar modules. Finally, significant RA-correlated genes were chosen as RA-trait module genes based on the threshold values |gene significance| > 0.5 and |module membership| > 0.7.

Identification of DEGs

Using the R package "limma", we identified the DEGs of BCPB of RA patients and healthy people in the expressing data. DEGs were characterized as fold change (FC) ≥ 2 (|log2FC|≥1) and adjusted P-values<0.05.

Functional Enrichment Analysis

The overlapping genes between RA-trait module genes and DEGs were defined as RA-trait DEGs. In order to obtain deeper understanding into the function of the DEGs and RA-trait DEGs, we performed gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis by the R package "clusterProfiler". The cut-off was set to P-value<0.05.

PPI network and Hub Genes Identification

In the Search Tool for the Retrieval of Interacting Genes (STRING) database, we analyzed the RA-trait DEGs and constructed a protein protein interaction (PPI) network[12]. In PPI network, we screened out hub genes by the cytoHubba (a Cytoscape plugin), which offers 12 distinct algorithms for evaluating core degree or importance of genes.

Validation of hub genes

BCPB of 10 RA patients and 10 healthy people were collected for qRT-PCR validation to verify the findings of the bioinformatics analysis. The Ethics Committees of Shanxi Bethune Hospital approved the study protocol (approval number: 2017LL043), and all people signed an informed consent form. After Ficoll density centrifugation, peripheral blood mononuclear cells (PBMCs) were obtained from venous blood samples and purify B cells from PBMCs using CD19 Microbeads (Miltenyi Biotec, Cat#: 130050301). Total RNA of B cells was isolated by using an RNA purification kit (EZBioscience, Cat#: B004DP). Total RNA was subsequently reverse transcribed to cDNA by using PrimeScript™ RT Master Mix (TaKaRa, Cat#: RR036A), and qRT-PCR was carried out by using TB Green® Premix Ex Taq™ II (TaKaRa, Cat#: RR820A). We selected GAPDH as the internal reference gene. The primers’ sequences were showed in Table 1. The relative expression levels of mRNA were obtained by 2−ΔΔ Ct method.

Table 1

Primer name

Primer forward sequences (5’-3’)

Primer reverse sequences (5’-3’)

GAPDH

CAGAACATCATCCCTGCCTCTAC

TTGAAGTCAGAGGAGACCACCTG

TLR8

GTTGGAACTACACGGAAACC

GGACTGGCACAAATGACATC

CCR2

GTGTGTGGAGGTCCAGGAGT

TTTCCTTTTCCACGACCATC

CX3CR1

CAACAGCAAGAAGCCCAAGAG

TGAAGAAGAAGGCGGTAGTGAA

GIMAP6

AGACGCTATCTGCCAAGCC

GGCCCAGTTGTGTCACCAG

C1QA

TCTGCACTGTACCCGGCTA

CCCTGGTAAATGTGACCCTTTT

GSEA

To research potential biological function of validated hub genes in BCPB of RA patients, we utilized GSEA for a single hub gene. In the GSE4588, according to the levels of hub genes' median expression, 7 RA samples were divided into two groups including low and high expression groups. The R software "clusterprofiler" is applied to perform GSEA and adjusted P-value < 0.05 was considered significant[13].

Statistical Analysis

According to data distribution characteristics, the differences between two groups' statistical significance were tested by t test or non-parametric test. All analysis were performed by software R4.1.1 and graphpad9.0.

Results

Co-Expression Networks

Samples of B cells in GSE4588 were normalized using "normalizeBetweenArrays" function in "limma" R package (Fig. 1A, Supplementary file 2). When the correlation coefficient threshold was set to 0.9, 7 was used as the soft-thresholding power (Fig. 1B). By WGCNA analysis, we constructed 9 different co-expression modules (Fig. 1C). The module containing the greatest number of genes was the turquoise module, followed by green and pink (Fig. 1C). In addition, these modules have little correlation with other modules (Fig. 1D).

Module-Trait Correlations of BCPB of RA patients

According to module-trait correlations, multiple modules are linked to BCPB of RA patients (Fig. 2A). The turquoise module containing 1171 genes was most significantly associated with BCPB of RA patients (Fig. 2B). Significance of all genes in the turquoise module was showed in Fig. 2C (Fig. 2C). Finally, 391 RA-trait module genes were screened out for further study (Supplementary file 3).

Identification and functional annotation of DEGs

We screened out a total of 93 DEGs (Supplementary file 4), including 48 down regulated genes and 45 upregulated genes from BCPB of RA patients compared with the control samples (Fig. 3A), with the volcano map revealing the significant difference (Fig. 3B). Figures 3C–D and Supplementary file 5 show the results of the DEG enrichment analysis. GO enrichment analysis showed DEGs were significantly enriched in regulation of interleukin-1 and interleukin-1 beta production, cytoplasmic side of plasma membrane and C-C chemokine receptor activity. Besides, significant enrichment of DEGs in several KEGG pathways including Salmonella infection, Cytokine-cytokine receptor interaction was revealed (Fig. 3D).

Functional annotation of the RA-trait DEGs

We identified 23 RA-trait DEGs (Fig. 4A, Fig. 4C). GO analysis showed that the genes in RA-trait DEGs were mainly enriched in the lymphocyte mediated immunity, external side of plasma membrane and C-C chemokine receptor activity (Fig. 4B). KEGG analysis showed RA-trait DEGs were significantly enriched in the viral protein interaction with cytokine and cytokine receptor (Fig. 4D).

PPI network construction

We constructed a PPI network in the STRING database by using RA-trait DEGs (Fig. 5A). To obtain hub genes that are crucial in related pathways, we utilized 8 different algorithms in the cytoHubba (including MCC, MNC, DMNC, EPC, Degree, Closeness, EcCentricity, and Radiality) to analyze the RA-trait DEGs in PPI network. Supplementary file 7 show the top 10 genes of each algorithm. We observed 5 genes were enriched (including TLR8, CCR2, CX3CR1, GIMAP6 and C1QA) in all 8 algorithms (Fig. 5B). In addition, these hub genes were closely related to each other (Fig. 5C) and significantly highly expressed (Fig. 5D) in the dataset GSE4588.

Validation of hub genes

qRT-PCR was carried out to further confirm expression levels of candidate genes in BCPB of 10 RA patients and 10 healthy people. Five hub genes (including TLR8, CCR2, CX3CR1, GIMAP6 and C1QA) were used to apply for qRT-PCR verification. As Fig. 5 shown, TLR8 (P < 0.05), CCR2 (P < 0.01) and CX3CR1 (P < 0.01) were obviously up-regulated in BCPB of RA patients compared with healthy people. The expression level of GIMAP6 and C1QA (P > 0.05) in BCPB was not significant difference between RA patients and healthy people, but its expression still displayed the up-regulated tendency in RA samples.

Gene Set Enrichment Analysis

By GSEA, we discovered a complete list of gene sets enriched by B cells that highly express the validated hub genes, including TLR8, CCR2 and C1QA, in the peripheral blood of RA patients (Fig. 6A, Supplementary file 8). Then, we selected immune-related gene sets from the full list for further analysis. B cells with high expression of these validated hub genes in RA patients were mainly enriched in three signal pathways, including "JAK_STAT3 response", "interferon alpha response", "interferon gamma response" (Fig. 6B).

Discussion

To our knowledge, there are no other studies using WGCNA to analyze BCPB in RA patients. In this study, we obtained 9 co-expression modules through WGCNA. The module most associated with the pathogenesis of BCPB of RA patients is the turquoise module, which contains 1171 genes. We also applied GO and KEGG analysis to clarify possible biological functions of DEGs and RA-trait DEGs. Most results of RA-trait DEGs are almost indistinguishable from those of DEGs. These findings of DEGs and RA-trait DEGs further proved the important roles of multiple cytokines and chemokines in the pathogenesis of BCPB in RA patients.

In addition, among the RA-trait DEGs, the hub gene is more important to BCPB of RA patients than other genes. Most of these hub genes are inflammatory genes. Among them, TLR8, CX3CR1, CCR2 are related to RA[14–17], and gene polymorphisms of CCR2 are related to RA susceptibility[18]. The activation of the signal pathway where TLR8 is located can cause immune cells in RA patients to produce IFN-α, IFN-γ, IL-6 and other inflammatory factors, which are involved in the maintenance and development of RA[19]. The autoantigens in RA patients can also make the body's B cells produce autoantibodies through the TLR signaling pathway, but more specific and detailed mechanisms need to be further studied [19, 20]. Recently, Aluri et al. found that gain-of-function variants in the TLR8 gene have been described in patients with immunodeficiency and bone marrow failure, some of which had positive antineutrophil antibodies[21]. This finding may reflect the growing body of evidence pointing to the role of bone marrow in the pathogenesis of RA. CC motif chemokine receptor-2 (CCR2) on B cells can participate in the body's inflammatory response by binding to monocyte chemoattractant protein-1[22]. In the RA patients' peripheral blood, ratio of CCR2+ B cells correlates positively with the level of IL-6[23]. As another important chemokine receptor, CX3CR1 has been studied extensively on T cells and monocytes in RA[15, 17]. Compared with TFH cells, CX3CR1 is more abundantly expressed in TPH cells in RA, which can promote the migration of TPH cells to the inflamed sites[15]. CX3CR1 is also involved in the differentiation of human monocytes into osteoclasts[24]. Li et al. found CX3CR1 influences B cell differentiation and the humoral immune response by regulating BCR signaling via actin remodeling[25]. Related research shows blockade of the CX3CL1-CX3CR1 axis could treat rheumatic diseases, so the role of CX3XR1 in B cells in RA is worthy of our further study[26]. In short, the role of these hub genes in B cells in RA has not yet been reported, and the relationship between these hub genes and B cells in RA remains to be studied in the future.

GSEA showed that B cells with high expression of hub gene were mainly enriched in five signal pathways, including "JAK_STAT3 response", "interferon alpha and gamma response" and "TNFα_signaling_via_NFκB". Relevant studies have shown that inhibiting the activity of Janus kinase in RA patients will weaken the proliferation and development of B cells in their bodies[27]. In addition, the intracellular part of the corresponding receptors on the B cell membrane of interferon α and interferon γ can be coupled with Janus kinase, so JAK inhibitors can block pro-inflammatory effects of interferon alpha and interferon gamma mediated by B cells in RA patients[28–30]. These fingdings indicate that JAK inhibitors may be an ideal drug for reducing inflammation caused by B cells in RA patients. A recent study also shows that "TNFα_signaling_via_NFκB" is activated in B cells and related to bone destruction in RA[31]. However, how B cells cause or aggravate the disease through the complement pathway in RA remains to be further studied.

Although our study investigates the co-expression gene networks associated with BCPB of RA patients by WGCNA analysis, this study still has some limitations. We have not further studied the effects of these hub genes on the proliferation and differentiation of B cells in RA patients. Besides, we have not studied the effects of these hub genes on the secretion of cytokines and autoantibodies by B cells in RA patients.

Conclusions

In conclusion, our research finds three key genes and functional biological pathways linked to the pathogenesis of BCPB of RA patients. These findings shed new light on the pathogenesis of B cells in RA, but the exact pathogenic mechanism of these key genes remains to be further explored.

Abbreviations

RA: Rheumatoid arthritis; BCPB: B cells from the peripheral blood; WGCNA: Weighted gene co-expression network analysis; GEO: Gene Expression Omnibus; DEGs: differentially expressed gene; PPI: protein protein interaction; FC: fold change; GO: gene ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; STRING: Retrieval of Interacting Genes; PBMCs: peripheral blood mononuclear cells; JAK: Janus kinase.

Declarations

Ethical Approval and Consent to participate

The Ethics Committees of Shanxi Bethune Hospital approved the study protocol (approval number: 2017LL043), and all people signed an informed consent form.

Consent for publication

Not applicable.

Availability of supporting data

Gene expression profiles of GSE4588 can be downloaded from the GEO database(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE4588).

Competing interests

The authors declare that they have no competing interests.

Funding

This work was supported by the National Natural Science Foundation of China [grant number 81771768] and by the applied basic research project of Shanxi Science and Technology Department [grant number 201901D111416].

Authors' contributions

Fengping Wu, Jinfang Gao and Xingxing Zhao performed the experiment and wrote the manuscript. Fengping Wu, Jinfang Gao and Xingxing Zhao have contributed equally to this work and share first authorship. Jie Kang, Xuexue Wang, Qing Niu and Jiaxi Liu analyzed the data and drew pictures. Liyun Zhang reviewed and revised the manuscript. All authors read and approved the final manuscript.

Acknowledgements

We thank all study participants, patients and all authors of GSE4588.

Authors' information

Fengping Wu1†, Jinfang Gao2†, Xingxing Zhao3†, Jie Kang4, Xuexue Wang4, Qing Niu1, Jiaxi Liu4, Liyun Zhang2*

1School of Basic Medical Sciences, Shanxi Medical University, Taiyuan, 030001, China

2Department of Rheumatology, Shanxi Bethune Hospital, Shanxi Academy of Medical Sciences, Tongji Shanxi Hospital, Third Hospital of Shanxi Medical University, Taiyuan, 030032, China

3School of Basic Medicine Sciences, Shanxi University of Chinese Medicine, Jinzhong, 030619, China

4Third Hospital of Shanxi Medical University, Shanxi Bethune Hospital, Shanxi Academy of Medical Sciences, Tongji Shanxi Hospital, Taiyuan, 030032, China

References

  1. Veale DJ, Orr C, Fearon U. Cellular and molecular perspectives in rheumatoid arthritis. Semin Immunopathol. 2017;39:343–354.
  2. Takemura S, Klimiuk PA, Braun A, Goronzy JJ, Weyand CM. T cell activation in rheumatoid synovium is B cell dependent. J Immunol. 2001;167:4710–4718.
  3. Meednu N, Zhang H, Owen T, Sun W, Wang V, Cistrone C, et al. Production of RANKL by Memory B Cells: A Link Between B Cells and Bone Erosion in Rheumatoid Arthritis. Arthritis Rheumatol. 2016;68:805–816.
  4. Barr TA, Shen P, Brown S, Lampropoulou V, Roch T, Lawrie S, et al. B cell depletion therapy ameliorates autoimmune disease through ablation of IL-6-producing B cells. J Exp Med. 2012;209:1001–1010.
  5. McInnes IB, Schett G. Cytokines in the pathogenesis of rheumatoid arthritis. Nat Rev Immunol. 2007;7:429–442.
  6. Mahendra A, Yang X, Abnouf S, Adolacion J, Park D, Soomro S, et al. Beyond Autoantibodies: Biologic Roles of Human Autoreactive B Cells in Rheumatoid Arthritis Revealed by RNA-Sequencing. Arthritis Rheumatol. 2019;71:529–541.
  7. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. Bmc Bioinformatics. 2008;9:559.
  8. Presson AP, Sobel EM, Papp JC, Suarez CJ, Whistler T, Rajeevan MS, et al. Integrated weighted gene co-expression network analysis with an application to chronic fatigue syndrome. Bmc Syst Biol. 2008;2:95.
  9. Yao Q, Song Z, Wang B, Qin Q, Zhang JA. Identifying Key Genes and Functionally Enriched Pathways in Sjogren's Syndrome by Weighted Gene Co-Expression Network Analysis. Front Genet. 2019;10:1142.
  10. Zhao X, Zhang L, Wang J, Zhang M, Song Z, Ni B, et al. Identification of key biomarkers and immune infiltration in systemic lupus erythematosus by integrated bioinformatics analysis. J Transl Med. 2021;19:35.
  11. Riquelme MI, Lubovac-Pilav Z. Gene Co-Expression Network Analysis for Identifying Modules and Functionally Enriched Pathways in Type 1 Diabetes. Plos One. 2016;11:e156006.
  12. 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-D613.
  13. Liberzon A, Birger C, Thorvaldsdottir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1:417–425.
  14. Leng RX, Di DS, Ni J, Wu XX, Zhang LL, Wang XF, et al. Identification of new susceptibility loci associated with rheumatoid arthritis. Ann Rheum Dis. 2020;79:1565–1571.
  15. Rao DA, Gurish MF, Marshall JL, Slowikowski K, Fonseka CY, Liu Y, et al. Pathologically expanded peripheral T helper cell subset drives B cells in rheumatoid arthritis. Nature. 2017;542:110–114.
  16. Lee YH, Bae SC, Song GG. Meta-analysis demonstrates association between TLR polymorphisms and rheumatoid arthritis. Genet Mol Res. 2013;12:328–334.
  17. Rana AK, Li Y, Dang Q, Yang F. Monocytes in rheumatoid arthritis: Circulating precursors of macrophages and osteoclasts and, their heterogeneity and plasticity role in RA pathogenesis. Int Immunopharmacol. 2018;65:348–359.
  18. Leng RX, Di DS, Ni J, Wu XX, Zhang LL, Wang XF, et al. Identification of new susceptibility loci associated with rheumatoid arthritis. Ann Rheum Dis. 2020;79:1565–1571.
  19. Li J, Wang X, Zhang F, Yin H. Toll-like receptors as therapeutic targets for autoimmune connective tissue diseases. Pharmacol Ther. 2013;138:441–451.
  20. Lee YH, Bae SC, Song GG. Meta-analysis demonstrates association between TLR polymorphisms and rheumatoid arthritis. Genet Mol Res. 2013;12:328–334.
  21. Aluri J, Bach A, Kaviany S, Chiquetto PL, Kitcharoensakkul M, Walkiewicz MA, et al. Immunodeficiency and bone marrow failure with mosaic and germline TLR8 gain of function. Blood. 2021;137:2450–2462.
  22. Moadab F, Khorramdelazad H, Abbasifard M. Role of CCL2/CCR2 axis in the immunopathogenesis of rheumatoid arthritis: Latest evidence and therapeutic approaches. Life Sci. 2021;269:119034.
  23. Kholodnyuk I, Kadisa A, Svirskis S, Gravelsina S, Studers P, Spaka I, et al. Proportion of the CD19-Positive and CD19-Negative Lymphocytes and Monocytes within the Peripheral Blood Mononuclear Cell Set is Characteristic for Rheumatoid Arthritis. Medicina (Kaunas). 2019;55.
  24. Muraoka S, Kaneko K, Motomura K, Nishio J, Nanki T. CX3CL1/fractalkine regulates the differentiation of human peripheral blood monocytes and monocyte-derived dendritic cells into osteoclasts. Cytokine. 2021;146:155652.
  25. Li N, Jiang P, Chen A, Luo X, Jing Y, Yang L, et al. CX3CR1 positively regulates BCR signaling coupled with cell metabolism via negatively controlling actin remodeling. Cell Mol Life Sci. 2020;77:4379–4395.
  26. Tanaka Y, Hoshino-Negishi K, Kuboi Y, Tago F, Yasuda N, Imai T. Emerging Role of Fractalkine in the Treatment of Rheumatic Diseases. Immunotargets Ther. 2020;9:241–253.
  27. Moura RA, Fonseca JE. JAK Inhibitors and Modulation of B Cell Immune Responses in Rheumatoid Arthritis. Front Med (Lausanne). 2020;7:607725.
  28. Wang SP, Iwata S, Nakayamada S, Sakata K, Yamaoka K, Tanaka Y. Tofacitinib, a JAK inhibitor, inhibits human B cell activation in vitro. Ann Rheum Dis. 2014;73:2213–2215.
  29. Sonomoto K, Yamaoka K, Kubo S, Hirata S, Fukuyo S, Maeshima K, et al. Effects of tofacitinib on lymphocytes in rheumatoid arthritis: relation to efficacy and infectious adverse events. Rheumatology (Oxford). 2014;53:914–918.
  30. Banerjee S, Biehl A, Gadina M, Hasni S, Schwartz DM. JAK-STAT Signaling as a Target for Inflammatory and Autoimmune Diseases: Current and Future Prospects. Drugs. 2017;77:521–546.
  31. Zhang LL, Xiao H, Zhang F, Wu YJ, Shu JL, Li Y, et al. BAFF, involved in B cell activation through the NF-kappaB pathway, is related to disease activity and bone destruction in rheumatoid arthritis. Acta Pharmacol Sin. 2021;42:1665–1675.