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.
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).
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.
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.
Rheumatoid arthritis (RA) is a chronic and systemic autoimmune disease. It is characterized by chronic synovitis and progressive articular damage [1,2,3]. Accurate assessment of RA process has potential to provide optimal treatment strategies as current clinical biomarkers are unable to monitory the disease activity and a part of patients with negative inflammatory tests still have active disease [4,5,6]. Moreover, classification of RA synovium is possible to offer sensitive predictors of RA progress . A series of studies have demonstrated that increasing numbers of mononuclear phagocytes, synovial fibroblasts, B cells, and T cells participated in RA progression and caused destruction of articular cartilage and bone [8,9,10,11,12,13,14]. However, few researchers paid attention to the roles of synovial cells in RA remission. Recently, with the help of single-cell transcriptome sequencing technology, tissue resident macrophages that contributed to the remission of RA were found in synovium. In human, the population of MerTKposCD206pos synovial tissue macrophages was higher in RA remission synovium compared to active RA synovium. Moreover, MerTKposCD206pos synovial tissue macrophages showed the biofunction including mediator profile resolution and repair responses . In mice, locally renewing Cx3cr1+ tissue-resident macrophages form a tight-junction barrier to protect intra-articular structures from inflammatory reaction . These two studies strongly indicated that there possibly may be subpopulations of synovium-resident cells with similar function in synovium and reinstating synovial homeostasis.
Fibroblasts are important structural cells in the synovium and have key roles in maintaining the homeostasis. Previous literatures have reported that fibroblasts had an important role in RA pathological processes . Recent scRNA-seq analyses have also described fibroblast heterogeneity in the synovium  and identified specific cluster associated with RA such as THY1+ fibroblasts, which interacted with macrophages and endothelial cells and led to severe inflammatory arthritis [19, 20]. However, most of these researches concentrated on fibroblast subclasses in promoting inflammation in arthritis. Fibroblast subsets that communicate with MerTKposCD206pos synovial tissue macrophages and engage in tissue repair have not been reported in detail . Furthermore, because of the usefulness of a detailed description of fibroblast subclasses in addressing the mechanisms of RA remission and monitoring RA activity, we hypothesized that there is a distinct population of fibroblasts producing an active role in RA remission. In this study, we integrated single-cell RNA sequencing datasets and bulk RNA-seq datasets and found HBEGF+ fibroblasts had an important role in RA remission.
The total of 13 datasets, including 7 single-cell transcriptomics datasets and 6 bulk RNA-seq datasets, were collected from public datasets Gene Expression Omnibus (GEO) and NIH IMMPOR (Table 1). The single-cell transcriptomics datasets contained 1 dataset of peripheral blood mononuclear cells from RA patients, 1 dataset of synovial cells from RA patients , 1 dataset of CD45− synovial cells from RA patients , 1 dataset of chondrocytes from osteoarthritis (OA) patients , 1 dataset of meniscus cells from OA patients , 1 dataset of CD45+ synovial cells mice with different states of arthritis (health and K/BxN serum transfer arthritis (STA)) , and 1 data of synovial cells from mice with different states of arthritis (health, STA, and STA in remission) . The bulk RNA-seq datasets contained 1 dataset of fibroblasts  from human synovial tissue, 1 dataset of macrophages  from human synovial tissue, 5 multicenter datasets of synovial tissue from 102 people with different states of arthritis (health, RA, OA) [24, 25], and 1 dataset of synovial tissue from RA patients before and after drug treatment . All the datasets were processed in R (V.4.0.0), and the results were showed using ggplot2 R package (V.3.3.5) except where mentioned.
Human single-cell RNA sequencing analysis
The human single-cell transcriptomics datasets, composed of 3 datasets from RA patients and 2 datasets from OA patients, were analyzed based on the states of arthritis.
Integrating scRNA-seq datasets of cells from RA patients
According to the Seurat single-cell analysis standard workflow [27, 28], firstly, each dataset was used to create Seurat object. Specifically, cells with <500 measured genes and >5% mitochondrial contamination were defined as low-quality cells and cells with >4500 measured genes were identified as potential doublets. They were filtered out from each dataset. After being filtered, total 29,382 cells were selected for following processes. All RA Seurat objects were merged into a different RA state object. The merged object was normalized (function NormalizeData, method = “LogNormalize,” scale. factor = 10,000), the 3000 most variable genes were identified, and the expression levels of these genes were scaled before performing PCA in variable gene space. Next, batch effect was corrected and merged object was integrated by running Harmony (version 1.0) . The top 25 harmony dimensions were provided as an input for UMAP and visualized the first two UMAP dimensions at a clustering resolution of 0.2. All steps were performed using functions implemented in the Harmony package and Seurat package (NormalizeData, FindVariableFeatures, ScaleData, RunPCA, FindNeighbours, FindClusters, RunUMAP) with default parameters, except where mentioned.
Next, distinct cell types were labeled by canonical marker genes such as PRG4, PDPN (fibroblasts), THY1, MCAM (mural cells), CD34, VWF (endothelial cells), CD2, CD4 (CD4+ T cells), CD8A, GNLY, GZMB (CD8+ T cells), LTB, CD3D (natural killer cells), VSIG4, CD163 (synovium macrophages), CD68, LYZ (periperal blood macrophages), XBP1, CD27 (plasma cells), CD79A, CD37 (B cells), and STMN1 (STMN1+ cells). Gene expression of each cluster was visualized using Dotplot.
Integrating scRNA-seq datasets of cells from OA patients
The OA scRNA-seq datasets were integrated following the steps mentioned above with the same parameters. Total 6708 cells were involved for analysis after filtering and classified into distinct cell types referring to the source of cells. Expression of EGFR was showed using function FeaturePlot.
Cell-cell communication in RA synovium
After identifying cell types in RA synovium, cell-cell communication was analyzed by implementing the CellChat (V.1.1.3) pipeline . A new CellChat object was created from the merged Seurat object. The paracrine/autocrine signaling interaction dataset of CellChatDB was set as referencing database. Next, the communication probability was computed using a truncated mean of 20% (function computeCommunProb, type = "truncatedMean", trim = 0.2). After that, the cell-cell communication was inferred and the cell-cell communication network was aggregated with default parameters. The number of interactions was visualized to show the aggregated cell-cell communication network and signaling sent from each cell cluster. EGF signaling pathway network was showed in heatmap and ligand such as HBEGF, AREG, BTC, EGF, EREG, TGFA, and receptor EGFR, which involved in EGF signaling pathway were showed using function Featureplot based on the merged RA Seurat object.
Bioinformatics analysis of HEBGF+ fibroblasts
When changing clustering resolution from 0.2 to 0.5 for visualizing the first two UMAP dimensions, fibroblasts that highly expressed HBEGF were divided into one group. Therefore, this group of cells was defined as HBEGF+ fibroblasts (average expression of HBEGF higher than 1.5) while fibroblasts lowly expressed HBEGF were defined as HBEGF− fibroblasts. Macrophages that highly expressed HBEGF were also defined as HBEGF+ macrophages and macrophages that lowly expressed HBEGF were defined as HBEGF− macrophages. Differential gene expression that HBEGF+ fibroblasts compared to HBEGF− fibroblasts was calculated by the function FindMarkers (Seurat R package) and showed in volcano plot. Gene expression of each cell type was computed using function FindAllMarkers (Seurat R package) and top 10 of which were showed in heatmap using function Doheatmap (Seurat R package). Expression of CLIC5, PDGFD, BDH2, ENPP1, SLAMF8, GK, L1RN, and JAK2 were displayed using function Vlnplot (Seurat R package). Differential expression gene markers of HBEGF+ fibroblasts and HBEGF− fibroblasts with the log2 fold change >1 were selected to generate a gene list for GO enrichment analysis using the clusterProfiler R package (V.4.2.2) [31, 32].
Human bulk RNA sequencing analysis
The bulk RNA-seq datasets were consisted of 1 dataset of fibroblasts from human synovial tissue, 1 dataset of macrophages from human synovial tissue, 5 multicenter datasets of synovial tissue from 53 RA joints, 33 OA joints and 26 healthy joints, and 1 dataset of synovial tissue derived from 12 RA patients before and after classical treatment of combination of triple disease-modifying anti-rheumatic drugs (triple DMARD treatment, methotrexate, sulfasalazine, and hydroxychloroquine). Datasets from cells and synovial tissue were integrated for analysis separately. Dataset from RA joint before and after drug treatment was analyzed individually.
Datasets of fibroblasts and macrophages were merged and normalized before correcting batch effect. Sva R package (V. 3.42.0) was used to correct batch effect by running function ComBat. Then point plot was used to show the HBEGF expression in fibroblasts and macrophages. The HBEGF expression between fibroblasts and macrophages was tested by Student’s t test with a significance threshold of P < 0.05. Next, samples of fibroblasts were divided into three groups, high HBEGF group (fibroblasts with count of HBEGF higher than 1000), middle HBEGF group (fibroblasts with count of HBEGF between 100 and 1000), and low HBEGF group (fibroblasts with count of HBEGF between lower than 100). Differential gene expression analysis was performed using DESeq2 R package (V 1.36.0) between high and low HBEGF expression group and differential expression genes (log2 fold change >1 and p < 0.05) were showed in volcano plot.
Five multicenter bulk RNA-seq datasets were merged, normalized, and corrected following the steps mentioned above. Then the expression of HBEGF, AREG, BTC, EGF, EREG, and TGFA was showed in box plot. HBEGF and gene markers such as CLIC5, PDGFD, BDH2, ENPP1, GK, IL1B, L1RN, and SLAMF8 were tested by the Pearson correlation test and displayed in dot plot separately. Expression of HBEGF, CLIC5, PDGFD, BDH2, ENPP1, GK, IL1B, L1RN, and SLAMF8 between healthy joints and RA joints was tested by Student’s t test with a significance threshold of P < 0.05. Expression of HBEGF, CLIC5, PDGFD, BDH2, ENPP1, GK, IL1B, L1RN, and SLAMF8 before and after treatment was tested by paired Student’s t test and showed in box plot with a significance threshold of P < 0.05.
Mouse single-cell RNA sequencing analysis
The mouse single-cell transcriptomics datasets were composed of 1 dataset of CD45+CD11b+Ly6G- synovial cells from mice with healthy joints and STA and 1 dataset of synovial cells from mice with healthy joints, STA, and STA in remission.
Two mouse scRNA-seq datasets were integrated following the steps of integrating scRNA-seq datasets of cells from RA patients with the same parameters. After filtering, total 28983 cells were involved for analysis. Next, distinct cell types were identified by marker genes such as Prg4, Pdpn, Hbegf (Hbegf+ fibroblasts and Hbegf− fibroblasts), Thy1, Mcam (mural cells), Cd34, Vwf (endothelial cells), Cd2, Cd4 (T cells), Vsig4, Cd163 (synovium macrophages), Cd79A, Cd37 (B cells), and Stmn1 (Stmn1 + cells). Expression of Hbegf was showed using the function FeaturePlot. Differential gene expression was calculated for each cell type using the function FindAllMarkers. Next, differential expression gene markers of Hbegf+ fibroblasts were ordered by the log2 fold change to generate a gene list, which was then used as an input for GSEA analysis. The top 200 differential expression gene markers of each cell type from human were used as the gene sets when running fgsea (V.1.2.0). Finally, the percentage of Hbegf+ fibroblasts was calculated under different arthritis conditions.
Expression of HBEGF was downregulated in RA synovium but upregulated in RA remission synovium
A total of 29,382 cells were included in the downstream analysis after three scRNA-seq datasets had been integrated and corrected for batch effects by Harmony. Eleven clusters in RA synovium were identified and showed in uniform manifold approximation and projection (UAMP) (Fig. 1a) and defined as fibroblasts (PRG4 and PDPN), mural cells (THY1 and MCAM), endothelial cells (CD34 and VWF), CD4+ T cells (CD2 and CD4), CD8+ T cells (CD8A, GNLY, and GZMB), natural killer cells (LTB and CD3D), synovium macrophages (CD68 and LYZ), periperal blood macrophages (CD68 and LYZ), plasma cells, B cells (CD79A and CD37), and STMN1+ cells (STMN1) referring to the canonical marker genes (Fig. 1b).
Next, CellChat was employed to analyze the cell-cell communication between each cell type in the RA synovial microenvironment (Fig. 1c). We found that fibroblasts had an extensive communication network with other clusters and participated in various kinds of paracrine or autocrine signaling interactions.
Among all the signaling pathways that fibroblasts were involved, we discovered the EGF signaling pathway that was the ligand-receptor interaction between fibroblasts and synovium macrophages (Fig. 2a). HBEGF, one of the ligands in EGF signaling pathway, was mainly expressed in a part of fibroblasts and synovial macrophages (Fig. 2b). Other ligands interacting with EGFR, such as AREG, BTC, EGF, EREG, and TGFA, were lowly expressed by all clusters (Fig. 2d). EGFR, corresponding receptor to HBEGF, was only expressed by fibroblasts (Fig. 2c) in synovium. Then we noticed other tissues such as cartilage and meniscus in the joint cavity and found meniscal cells highly expressed EGFR, whereas chondrocytes did not (Fig. 2e, f). Next, we compared bulk RNA-seq data harvested from synovial cells under heathy, OA, and RA condition. HBEGF expression was significantly decreased in the RA state compared to the healthy and OA states (P = 0.0003, Fig. 2g). Moreover, comparing to the expression of HBEGF, expression of other EGFR ligands such as AREG, BTC, EGF, EREG, and TGFA were extremely low in all cases compared to the expression of HBEGF. After a 6-month triple DMARD treatment, the expression of HBEGF increased in RA synovium (P =0.05433, Fig. 2h).
HBEGF+ fibroblast was a subset with distinct biofunction in synovium
As can be seen from Fig. 2b, the fibroblasts which highly expressed HBEGF (average expression of HBEGF higher than 1.5) were distributed into a small group. Then we defined them as HBEGF+ fibroblasts and fibroblasts with low expression of HBEGF were defined as HBEGF− fibroblasts. Referring to the previous studies, macrophages with high or low expression of HBEGF were defined as HBEGF+ macrophages and HBEGF− macrophages separately (Fig. 3a).
By calculating the differential gene expression between HBEGF+ fibroblasts and HBEGF− fibroblasts, we detected 162 genes with distinct expression patterns between HBEGF+ fibroblasts and HBEGF− fibroblasts. Forty-four gene markers were highly expressed by HBEGF+ fibroblasts and 118 gene markers by HBEGF− fibroblasts. The main differential expressed gene of HBEGF+ fibroblasts was HBEGF whereas the differentially expressed gene in HBEGF− fibroblasts is MDK, which is related to the pathogenesis of RA (Fig. 3b). Next, we figured out genes mainly expressing by each cluster and found HEBGF+ fibroblasts were heterogeneous from other subpopulations with high expression of CRTAC1, ITGB8, SCARA5, THBS4, and ITGBL1 while HBEGF− fibroblasts highly expressed CXCL12 and MMP2 (Fig. 3c). In bulk RNA sequencing datasets of synovial fibroblasts and macrophages, we also found subsets of fibroblasts and macrophages highly expressing HBEGF in synovium (Fig. 3d). Comparing to low HBEGF group (fibroblasts with count of HBEGF lower than 100), we noted that the differential expression genes (DEG) of high HBEGF group (fibroblast with count of HBEGF higher than 1000) were more abundant in HBEGF+ fibroblasts.
Finally, we used GO enrichment analysis to figure out the difference between HBEGF+ fibroblasts and HBEGF− fibroblasts in biofunction. We found that HBEGF+ fibroblasts involved in positive regulation of cell migration and motility, cellular component movement and cell growth whereas HBEGF− fibroblasts engaged in collagen metabolic and catabolic process and angiogenesis (Fig. 3f, g).
Population of HBEGF+ fibroblasts was decreased in RA synovium but increased in RA remission synovium
To evaluate the population changes of HBEGF+ fibroblasts in different states of synovium, we investigated gene markers that are specially expressed by HBEGF+ fibroblasts and HBEGF+ macrophages. We found that CLIC5, PDGFD, BDH2, and ENPP1 were mainly expressed by HBEGF+ fibroblasts, and SLAMF8, GK, L1RN, and JAK2 were mainly expressed by HBEGF+ macrophages (Fig. 4a).
Then, we selected these gene markers for the Pearson correlation test in the integrated bulk RNA sequencing dataset and found that the expression of CLIC5, PDGFD, BDH2, and ENPP1 was positively correlated with expression of HBEGF while the expression of SLAMF8, GK, L1RN, and JAK2 was negatively correlated with expression of HBEGF (Fig. 4b). Expression of CLIC5, PDGFD, BDH2, and ENPP1 was decreased, while expression of SLAMF8, GK, L1RN, and JAK2 was increased compared to synovial membranes in healthy and OA states (Fig. 4c). Moreover, after triple DMARD treatment, the expression of CLIC5, PDGFD, BDH2, and ENPP1 had an increase in most of patients. Conversely, the expression of SLAMF8, GK, L1RN, and JAK2 had a decrease after treatment in most of patients (Fig. 4d).
Next, we analyzed single-cell RNA sequencing datasets from mouse with different states of arthritis (healthy, STA, and STA after treatment). Total 28,983 cells included fibroblasts, mural cell, endothelial cells, T cells, synovium macrophages, B cells, and STMN1+ cells (Fig. 4e). Hbegf+ fibroblasts also existed in mouse synovium (Fig. 4f). However, few macrophages with high expression of Hbegf could be found in mouse synovium. To prove that Hbegf+ fibroblasts from mice were similar to HBEGF+ fibroblasts from human, we ran GSEA analysis and found the phenotype of Hbegf+ fibroblasts from mice aligned closer to HBEGF+ fibroblasts than the other clusters from human (Fig. 4g). Finally, we found that the population change of Hbegf+ fibroblasts had a similar pattern as the expression change of HBEGF under different conditions. The proportion of Hbegf+ fibroblasts decreased in RA synovium comparing to healthy tissue. However, after treatment, the proportion of Hbegf+ fibroblasts returned to a healthy level (Fig. 4h).
In rheumatoid arthritis (RA), synovial fibroblasts have been considered as the key roles in regulation of joint homeostasis [10, 12, 14, 17,18,19,20]. Corresponding with previous researches, our study confirmed that fibroblasts had an extensive communication network with other clusters in RA synovium. Among the cell-cell communications between fibroblasts and other clusters, we found EGF singling pathway and HBEGF+ fibroblasts.
This study subdivided fibroblasts into 2 clusters based on the expression of HBEGF. Fibroblasts that highly expressed HBEGF (average expression of HBEGF higher than 1.5) were defined as HBEGF+ fibroblasts while the population with low expression of HBEGF were regarded as HBEGF− fibroblast. HBEGF, heparin-binding EGF-like growth factor, is one of the ligands for the ErbB family of epidermal growth factor receptors (including EGFR) . It stimulates the migration, differentiation, and proliferation of cells to fill the damaged area and repair tissue lesion. Recent literatures have reported the protective function of HBEGF in TNF-driven chronic intestinal inflammation  and cartilage degeneration diseases . And this study demonstrated HBEGF expression was downregulated in RA synovium and increased after classic therapeutic strategy—triple DMARD treatment. A similar pattern could be seen in the population of Hbegf+ fibroblasts in mice with different states of arthritis. The amount of Hbegf+ fibroblasts decreased in RA joints but increased in RA remission joints. And GO enrichment analysis showed that HBEGF+ fibroblasts played a role in cell growth and positive regulation of cell migration and motility cellular component movement while HBEGF− fibroblasts showed opposite biofunction and involved in collagen metabolic process and angiogenesis which were proven to promote inflammation in arthritis [36,37,38]. Therefore, we believe that HBEGF+ fibroblasts played essential roles in the remission of RA. And we also believe that further researches on HBEGF+ fibroblasts could help to address the mechanism of RA remission and may identify novel biomarkers for the prediction of RA process.
In synovium, a fraction of synovial macrophages also highly expressed HBEGF. Previous research defined them as HBEGF+ macrophages and described their function in RA synovium . HBEGF+ macrophages activated synovial fibroblasts and subsequently induced invasiveness in synovium. In order to figure out the primary source of HBEGF, we used the Pearson correlation test and found that expression of HBEGF was more relative to HBEGF+ fibroblasts than HBEGF+ macrophages. In mice, we could find Hbegf+ fibroblasts exist in synovium but Hbegf+ macrophages did not. Furthermore, the population change of Hbegf+ fibroblasts showed the similar pattern as the expression change of HBEGF. The population of Hbegf+ fibroblasts decreased in RA state synovium and had an increase after RA remission. That meant HBEGF+ fibroblasts, instead of HBEGF+ macrophages, were the primary source of HBEGF in synovium and the population change of HBEGF+ fibroblasts was a reason for the expression change of HBEGF.
This study has proven that expression of HBEGF was downregulated in RA state synovium. However, pervious literatures had pointed out that the importance of the ErbB family pathway in chronic pain and activating EGF singling pathway resulted in deterioration of RA. For instance, injecting HBEGF into the paw of mice caused painful mechanical hypersensitivity and severe pain and blocking the ErbB receptor could alleviate RA pain and joint inflammation [39, 40]. In vivo, there are two different structural forms of HBEGF including proHBEGF (transmembrane protein) and sHBEGF (soluble protein) [41, 42]. proHBEGF is a precursor for sHBEGF and can be cleaved at the plasma membrane to yield sHBEGF. proHBEGF takes part in juxtacrine activity and sHBEGF engages in paracrine activity. From scRNA-seq data, we noticed that HBEGF+ fibroblasts accounted for a small part comparing to the huge amount of HBEGF− fibroblasts. And specially, the ErbB receptor, EGFR, was highly expressed by HBEGF− fibroblasts, which also highly expressed MDK (Midkine) , and CXCL12 (C-X-C Motif Chemokine Ligand 12) [44, 45] that engage in the pathophysiology of RA. So it is possible that the different spatial location between HEBGF+ fibroblasts and HEBGF− fibroblasts in some ways prevents proHBEGF from contacting with EGFR and activating HBEGF− fibroblasts. Moreover, small quantities of HBEGF+ fibroblasts limit the production of sHBEGF and in other ways inhibit the interaction of sHBEGF with HEBGF− fibroblasts. Oppositely, injecting exogenous HBEGF, which served as sHBEGF, mediated EGFR expressed by HBEGF− fibroblasts directly and subsequently activated HBEGF− fibroblasts to secrete increasingly cytokines and chemokines.
In summary, population of HBEGF+ fibroblasts and expression of HBEGF decreased in RA synovium and increased after treatment. To conclude, 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.
Availability of data and materials
All the data of this manuscript was collected from public datasets such as GEO and NIH IMMPOR.
Heparin-binding EGF-like growth factor
Single-cell RNA sequencing
GO enrichment analysis:
Gene ontology enrichment analysis
Uniform manifold approximation and projection
Epidermal growth factor receptor
Epidermal growth factor receptor
Gene Expression Omnibus
Differential expression genes
triple DMARD treatment:
Combination of triple disease-modifying anti-rheumatic drugs
Sokka T, Pincus T. Erythrocyte sedimentation rate, C-reactive protein, or rheumatoid factor are normal at presentation in 35%-45% of patients with rheumatoid arthritis seen between 1980 and 2004: analyses from Finland and the United States. J Rheumatol. 2009;36(7):1387–90. https://doi.org/10.3899/jrheum.080770 Epub 2009 May 1. PMID: 19411389.
Orange DE, Agius P, DiCarlo EF, et al. Histologic and transcriptional evidence of subclinical synovial inflammation in patients with rheumatoid arthritis in clinical remission. Arthritis Rheum. 2019;71(7):1034–41. https://doi.org/10.1002/art.40878 Epub 2019 Jun 5. PMID: 30835943; PMCID: PMC6594902.
Laurent L, Anquetil F, Clavel C, et al. IgM rheumatoid factor amplifies the inflammatory response of macrophages induced by the rheumatoid arthritis-specific immune complexes containing anticitrullinated protein antibodies. Ann Rheum Dis. 2015;74(7):1425–31. https://doi.org/10.1136/annrheumdis-2013-204543 Epub 2014 Mar 11. PMID: 24618262.
Culemann S, Grüneboom A, Nicolás-Ávila JÁ, et al. Locally renewing resident synovial macrophages provide a protective barrier for the joint. Nature. 2019;572(7771):670–5. https://doi.org/10.1038/s41586-019-1471-1 Epub 2019 Aug 7. PMID: 31391580; PMCID: PMC6805223.
Ji Q, Zheng Y, Zhang G, et al. Single-cell RNA-seq analysis reveals the progression of human osteoarthritis. Ann Rheum Dis. 2019;78(1):100–10. https://doi.org/10.1136/annrheumdis-2017-212863 Epub 2018 Jul 19. PMID: 30026257; PMCID: PMC6317448.
Sun H, Wen X, Li H, et al. Single-cell RNA-seq analysis identifies meniscus progenitors and reveals the progression of meniscus degeneration. Ann Rheum Dis. 2020;79(3):408–17. https://doi.org/10.1136/annrheumdis-2019-215926 Epub 2019 Dec 23. PMID: 31871141; PMCID: PMC7034356.
Woetzel D, Huber R, Kupfer P, et al. Identification of rheumatoid arthritis and osteoarthritis patients by transcriptome-based rule set generation. Arthritis Res Ther. 2014;16(2):R84. https://doi.org/10.1186/ar4526 PMID: 24690414; PMCID: PMC4060460.
Broeren MG, de Vries M, Bennink MB, et al. Disease-regulated gene therapy with anti-inflammatory interleukin-10 under the control of the CXCL10 promoter for the treatment of rheumatoid arthritis. Hum Gene Ther. 2016;27(3):244–54. https://doi.org/10.1089/hum.2015.127 Epub 2016 Feb 16. PMID: 26711533.
Hao Y, Hao S, Andersen-Nissen E, et al. Integrated analysis of multimodal single-cell data. Cell. 2021;184(13):3573–3587.e29. https://doi.org/10.1016/j.cell.2021.04.048 Epub 2021 May 31. PMID: 34062119; PMCID: PMC8238499.
Korsunsky I, Millard N, Fan J, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods. 2019;16(12):1289–96. https://doi.org/10.1038/s41592-019-0619-0 Epub 2019 Nov 18. PMID: 31740819; PMCID: PMC6884693.
Jin S, Guerrero-Juarez CF, Zhang L, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun. 2021;12(1):1088. https://doi.org/10.1038/s41467-021-21246-9 PMID: 33597522; PMCID: PMC7889871.
Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7. https://doi.org/10.1089/omi.2011.0118 Epub 2012 Mar 28. PMID: 22455463; PMCID: PMC3339379.
Zhou L, Zhou W, Joseph AM, et al. Group 3 innate lymphoid cells produce the growth factor HB-EGF to protect the intestine from TNF-mediated inflammation. Nat Immunol. 2022;23(2):251–61. https://doi.org/10.1038/s41590-021-01110-0 Epub 2022 Jan 31. PMID: 35102343; PMCID: PMC8842850.
Gudmann NS, Junker P, Juhl P, et al. Type IV collagen metabolism is associated with disease activity, radiographic progression and response to tocilizumab in rheumatoid arthritis. Clin Exp Rheumatol. 2018;36(5):829–35 Epub 2018 May 8. PMID: 29745884.
Elshabrawy HA, Chen Z, Volin MV, et al. The pathogenic role of angiogenesis in rheumatoid arthritis. Angiogenesis. 2015;18(4):433–48. https://doi.org/10.1007/s10456-015-9477-2. Epub 2015 Jul 22. PMID: 26198292; PMCID: PMC4879881.
Tsai CH, Chen CJ, Gong CL, et al. CXCL13/CXCR5 axis facilitates endothelial progenitor cell homing and angiogenesis during rheumatoid arthritis progression. Cell Death Dis. 2021;12(9):846. https://doi.org/10.1038/s41419-021-04136-2 PMID: 34518512; PMCID: PMC8437941.
Wangzhou A, Paige C, Neerukonda SV, et al. A ligand-receptor interactome platform for discovery of pain mechanisms and therapeutic targets. Sci Signal. 2021;14(674):eabe1648. https://doi.org/10.1126/scisignal.abe1648 PMID: 33727337; PMCID: PMC8097872.
Hieda M, Isokane M, Koizumi M, et al. Membrane-anchored growth factor, HB-EGF, on the cell surface targeted to the inner nuclear membrane. J Cell Biol. 2008;180(4):763–9. https://doi.org/10.1083/jcb.200710022 PMID: 18299347; PMCID: PMC2373455.
Shindo E, Nanki T, Kusunoki N, et al. The growth factor midkine may play a pathophysiological role in rheumatoid arthritis. Mod Rheumatol. 2017;27(1):54–9. https://doi.org/10.1080/14397595.2016.1179860 Epub 2016 May 10. PMID: 27163292.
NC and BF analyzed the data and wrote the main manuscript text. They contributed equally to this work and should be considered co-first authors. ZH, XY, and JW prepared all the figures. All authors reviewed and approved the final manuscript.
Nachun Chen, resident doctor. Baoying Fan, attending physician. Zhiyong He, head of the First Department of Orthopaedics; Xinping Yu, chief physician. Jinjun Wang, associate chief physician. All authors work in First Department of Orthopaedics, Zhongshan City People’s Hospital Affiliated to Sun Yat-sen University, Zhongshan City, Guangdong Province, China.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
Chen, N., Fan, B., He, Z. et al. Identification of HBEGF+ fibroblasts in the remission of rheumatoid arthritis by integrating single-cell RNA sequencing datasets and bulk RNA sequencing datasets.
Arthritis Res Ther24, 215 (2022). https://doi.org/10.1186/s13075-022-02902-x