Identification of HBEGF+ fibroblasts in the remission of rheumatoid arthritis by integrating single-cell RNA sequencing datasets and bulk RNA sequencing datasets

Abstract

Background: Fibroblasts are important structural cells in synovium and play key roles in maintaining the synovial homeostasis. By single-cell RNA sequencing (scRNA-seq), subpopulation of synovium-resident cells has been reported to protect intra-articular structures from chronic inflammation and promote tissue repair. However, a significant number of researchers have concentrated on the role of fibroblasts in the progress of rheumatoid arthritis (RA) while few reports had described the contribution of distinct fibroblast subsets in the RA remission. It is helpful to understand the role of fibroblast subpopulations in the RA process to provide predictive biomarkers and address RA remission mechanisms. Here, we found HBEGF+ fibroblasts that contributed to RA remission by integrating scRNA-seq datasets and bulk RNA sequencing (bulk RNA-seq) datasets.

Method: Three single-cell RNA datasets of cells harvested from RA patients were processed and integrated by Seurat and Harmony R packages. After identifying cell types by classic marker genes, the integrated dataset was used to run CellChat for analysis of cell-cell communication. Specially, EGF signaling pathway was found and HBEGF+ fibroblasts were identified based on HBEGF expression. Differential expressed genes of HBEGF+ were shown in heatmap and volcano plot and used to run gene ontology (GO) enrichment analysis. Next, bulk RNA-seq datasets of synovium under different conditions (health, osteoarthritis (OA), rheumatoid arthritis, before and after classical treatment) were compared to show expression change of HBEGF and gene markers that are mainly expressed by HBEGF+ fibroblasts such as CLIC5, PDGFD, BDH2, and ENPP1. Finally, two single-cell RNA sequencing datasets of synovial cells from mice were integrated to identify Hbegf+ fibroblasts and calculate the population of Hbegf+ fibroblasts under different joint conditions (health, K/BxN serum transfer arthritis (STA), and remission of STA).

Result: After integrating three single-cell RNA sequencing datasets, we identified 11 clusters of synovial cells, such as fibroblasts, mural cells, endothelial cells, CD4+ T cells, CD8+ T cells, natural killer cells, synovium macrophage, peripheral blood macrophages, plasma cells, B cells, and STMN1+ cells. We found fibroblasts had an extensive communication network with other clusters and interacted with synovial macrophages through EGF signaling pathway via analysis of cell-cell communication between synovial cells. HBEGF, ligand to EGF signaling pathway, was highly expressed by a subset of fibroblasts and macrophages, and EGFR, receptor to EGF signaling pathway, was highly expressed by fibroblasts and meniscus cells. Moreover, HBEGF was downregulated under RA state and it had an increase after classical treatment. We then defined fibroblasts with high expression of HBEGF as HBEGF+ fibroblasts. In addition, we also found that HBEGF+ fibroblasts highly expressed CRTAC1, ITGB8, SCARA5, THBS4, and ITGBL1, genes relative to encoding extracellular matrix proteins and engaged in positive regulation of cell migration and motility, cellular component movement, and cell growth by GO enrichment analysis. We eventually identified HBEGF+ fibroblasts specially expressed CLIC5, PDGFD, BDH2, and ENPP1, which positively correlated with the expression of HBEGF. Moreover, the expression of CLIC5, PDGFD, BDH2, and ENPP1 was downregulated under RA state and elevated by classical therapy. On the contrary, the HBEGF+ macrophages specially expressed SLAMF8, GK, L1RN, and JAK2, which negatively correlated with the expression of HBEGF. The expression was upregulated in SLAMF8, GK, L1RN, and JAK2 under the RA state, whereas it was decreased after classical treatment. In mice, the number of Hbegf+ fibroblasts was reduced in the RA synovium but increased in the RA remitting synovium.

Conclusions: HBEGF+ fibroblasts play a role in the remission of rheumatoid arthritis, and HBEGF has potential to become a novel biomarker for prediction of RA progress.

Fig. 1

Cross-talk analysis in RA synovium. a UMAP of single-cell RNA-seq data of 29,382 cells from three human datasets. Eleven clusters at UMAP of integrated dataset. b Dot plot showing the average expression level of canonical marker genes of each cluster. c Cross-talk analysis between each cluster in RA synovium

Fig. 2

EGF signaling pathway in RA synovium. a Heatmap showing the EGF signaling interaction scores between each cluster in RA synovium. b Ligand (HBEGF) and receptor (EGFR) of EGF signaling pathway are shown in dot plots. c Other ligands (AREG, BTC, EGF, EREG, and TGFA) of EGF signaling pathway shown in dot plots. d 6708 cells from OA meniscus and cartilage at UMAP and receptor (EGFR) of EGF signaling pathway are shown in dot plot (e). f Boxplot showing expression of ligands (HBEGF, AREG, BTC, EGF, EREG, and TGFA) using bulk RNA-seq profiles of healthy joint synovium (n = 26), RA joint synovium (n = 53), and OA joint synovium (n = 33). Significance determined by Student’s t test (P = 0.0003). e Boxplot showing HBEGF expression of bulk RNA-seq profiles from RA synovium before and after triple DMARD treatment (n=19). Significance determined by Student’s paired t test (P = 0.05433)

Fig. 3

Bioinformation of HBEGF+ fibroblasts. a HBEGF+ fibroblasts of RA synovium are displayed at UMAP. b Differential gene expression of HBEGF+ fibroblasts and HBEGF− fibroblasts are shown in a volcano plot. c Heatmap shows top 10 expressed gene markers of each cluster. d HBEGF expression in synovial fibroblasts (n = 25) and macrophages (n = 12) by bulk RNA sequencing (P = 0.7911) are plotted as log2 count + 1. e Differential expression genes between high HBEGF group and low HBEGF group are highlighted on the HBEGF+ fibroblast versus HBEGF− fibroblast plot from b. f, g GO enrichment analysis of HBEGF+ fibroblasts (f) and HBEGF− fibroblasts (g)

Fig. 4

Population change of HBEGF+ fibroblasts. a Expression of CLIC5, PDGFD, BDH2, ENPP1, SLAMF8, IL1RN, GK, and JAK2 are shown in Vin plot using single-cell RNA-seq profiles. b The Pearson correlation between HBEGF and gene markers such as CLIC5 (Cor = 0.5917), PDGFD (Cor = 0.4276), BDH2 (Cor = 0.4415), ENPP1 (Cor = 0.4097), SLAMF8 (Cor = −0.5497), L1RN (Cor = −0.2266), GK (Cor = −0.3257), and JAK2 (Cor = −0.4875) is displayed on a dot plot. c Expression of CLIC5, PDGFD, BDH2, ENPP1, SLAMF8, IL1RN, GK, and JAK2 are displayed in boxplot using bulk RNA-seq profiles of healthy joint synovium (n =26), RA joint synovium (n = 53), and OA joint synovium (n = 33). Significance determined by Student’s t test (CLIC5, P = 0.0018, PDGFD, P = 0.0238, BDH2, P = 0.0001, ENPP1, P = 0.0399, SLAMF8, P = 0.0001, L1RN, P = 0.0602, GK, P =0.1461, JAK2, P = 0.0001). d Expression of CLIC5, PDGFD, BDH2, ENPP1, SLAMF8, IL1RN, GK, and JAK2 displayed in a boxplot using bulk RNA-seq profiles from RA synovium before and after triple DMARD treatment (n=19). Significance determined by Student’s paired t test (CLIC5, P = 0.0355, PDGFD, P = 0.0022, BDH2, P = 0.0010, ENPP1, P = 0.0124, SLAMF8, P = 0.0038, L1RN, P = 0.4620, GK, P =0.1461, JAK2, P = 0.0048). e UMAP projection of single-cell RNA-seq data of 28,983 cells from two mouse datasets. Seven clusters at UMAP of integrated dataset. f Expression of Hbegf shown in dot plot. g GSEA using top 200 expressed gene markers of each cell cluster in human synovium as gene sets and ranked gene lists from Hbegf+ fibroblasts in mouse synovium. Normalized enrichment scores of HBEGF+ fibroblasts, NES = 1.6723, P.adj = 0.0001. h Percentage of Hbegf+ fibroblasts in synovium under different states (health, STA, and STA after treatment)