DEGs of IgG4-RD relative to HCs
In this study, we mainly explored the biological characteristics of IgG4-RD based on transcriptome data (Figure 1). The patients were categorized into IgG4-RD and HC groups. Data cleaning operations, such as normalization and annotation probes, were performed on the expression profiles of the two datasets. We performed principal component analysis to explore the effect of gene expression values on the IgG4-RD group relative to the HC group, which revealed differences in the principal components between the two groups (Figure 2A–B). The DEGs of the two datasets were obtained via differential analysis using the limma package. In the GSE40568 dataset, 626 DEGs were obtained, including 550 upregulated and 76 downregulated genes (Figure 2C, E). In the GSE66465 dataset, 158 DEGs were obtained, including 107 upregulated and 51 downregulated genes (Figure 2D, F). Volcano plots and heat maps were subsequently plotted for the DEGs, which could accurately distinguish the IgG4-RD group from HC group (Figure 2 C–F).
Identification of IgG4-RD-related gene modules using WGCNA
First, we performed WGCNA on LSG tissues of IgG4-RD and set a threshold for genes with the CV of >5%; accordingly, 2,207 disease group-related genes were included in this analysis. Samples were categorized into disease and control groups through systematic tree clustering analysis (Figure 3A). Based on the “pickSoftThreshold” function of the “WGCNA” package, we finally selected a soft threshold of 20 (Figure 3B). Consequently, we obtained seven gene set modules using topological overlap matrix and hierarchical average linkage clustering method analysis (Figure 3C–D). The turquoise module showed the strongest correlation with IgG4-RD (r = 0.89, p = 0.003), and 344 genes associated with IgG4-RD were obtained in LSG using WGCNA (Figure 3D). The DEGs obtained through GSE40568 dataset analysis of LSG samples were intersected with the genes in the IgG4-RD module, and 157 intersecting genes were obtained (Figure 3E). There were 839 genes with CV of >5% in PBMCs, and 10 modules were obtained using WGCNA (Figure 4A-D). Among them, turquoise was found to be associated with IgG4-RD, and 249 disease group-related genes were finally obtained. The DEGs obtained through GSE66456 dataset analysis of PBMC samples were intersected with the genes in the IgG4-RD module, providing 121 intersecting genes (Figure 4E).
Biological processes and signaling pathways enriched in IgG4-RD-related DEGs
We performed GO functional enrichment analysis on key DEGs to analyze the relationship between IgG4-RD-related DEGs and BP, MF, CC, biological pathways, and diseases. The results revealed that GO items coenriched with IgG4-RD-related gene modules of LSG and PBMC samples were mainly enriched in BP, such as cGMP metabolism, coagulation regulation, glucocorticoid response, corticosteroid response, platelet degranulation, fatty acid metabolism, inflammatory response regulation, lymphocyte activation regulation, and steroid hormone response (Figure 5A, D–E). Further, they were enriched in CC, such as platelet alpha granule, plasma membrane raft, collagen-containing extracellular matrix (ECM), platelet alpha granule membrane, and caveola (Figure 5A, D–E), as well as MF, such as ECM binding, cytokine binding, heparin-binding, peptidase regulator activity, and glycosaminoglycan binding (Figure 5A, D–E).
Next, KEGG pathway enrichment analysis was performed on the DEGs, which revealed that DEGs in LSG were enriched in biological pathways, such as malaria, peroxisome proliferator-activated receptor signaling pathway, ECM–receptor interaction, focal adhesion, and proteoglycans in cancer (Figure 5B), and the most significantly enriched part of each pathway was also displayed (Figure 5F). The KEGG pathway enrichment results of the gene modules showed the highest association with IgG4-RD in PBMCs, such as malaria, ECM–receptor interaction, viral protein interaction with cytokine and cytokine receptor, platelet activation, and chemokine signaling pathway (Figure 5C, G).
Then, we analyzed GSEA based on PBMC and LSG samples and noted a possible association with IgG4-RD and NABA secreted factors, BIOCARTA eryth pathway, REACTOME formation of the fibrin clot clotting cascade, WP human complement system, REACTOME response to elevated platelet cytosolic CA2, WP adipogenesis, REACTOME L1CAM interactions, REACTOME platelet activation signaling and aggregation, WP focal adhesion PI3K AKT MTOR signaling pathway, REACTOME signaling by RHO GTPASES, and REACTOME signaling by receptor tyrosine kinases (Figure 6A–D).
PPI construction and hub gene identification
We constructed a PPI network of key IgG4-RD-related DEGs using STRING data to analyze the PPI of DEGs and visualized the interaction. Cytoscape was used to evaluate the obtained data and visualize the PPI network. The larger the circle, the higher the fold-increase in differential expression. Blue represents downregulated genes, whereas red represents upregulated genes (Figure 7A, C). Further, we used the MCC algorithm of the cytoHubba plugin to select the top 15 genes as key genes (Figure 7B, D). miRNAs that might bind to the hub gene were predicted through the website, and a regulatory network was constructed. In total, 4 hub genes and 71 miRNAs were included in the network (Figure 7E). Moreover, the DEGs and WGCNA analysis revealed 5 common significant hub genes in both LSG and PBMC samples of IgG4-RD patients (Figure 8A). These 5 hub genes (GNG11, PDE5A, PLK2, PROS1, and THBS1) were upregulated in both LSG and PBMCs (Figure 8C, D). Surprisingly, THBS1 was the only key gene in both GSE66456 dataset from PBMC samples and GSE40568 dataset from LSG samples, indicating the important role of THBS1 in IgG4-RD (Figure 8B).
Plasma THBS1 levels are significantly higher in IgG4-RD
We determined the presence of THBS1 in the plasma samples of patients with IgG4-RD using ELISA to verify whether THBS1 is upregulated in IgG4-RD. This study included 48 patients with IgG4-RD who did not receive treatment and 24 HCs who were matched by age and sex. Detailed clinical characteristics of the participants are summarized in Table S1. Plasma THBS1 levels were significantly elevated in patients with IgG4-RD (42.59 μg/mL; interquartile range [IQR], 32.13–61.44) compared with those in HCs (17.63 μg/mL; IQR, 8.22–34.80; p < 0.0001) (Figure 8E). Additionally, the receiver operating characteristic curve analysis was used to determine the cutoff level of THBS1 for identifying patients with IgG4-RD based on the Youden index, which provided a value of 27.94 μg/mL and area under the curve of 0.8351 (95% confidence interval: 0.739–0.931, p < 0.0001) (Figure 8F).
Changes in the infiltration components of immune cells in IgG4-RD
The gene expression matrix data were analyzed via CIBERSORTx, and an immune cell infiltration matrix was obtained, which showed the distribution of various types of immune cell infiltration in the samples of each group (Figure 9A). After filtering out cells that were not infiltrated in any sample, a matrix of 21 immune cell infiltrates was obtained. The differences in immune cell infiltration between the control and IgG4-RD groups were also analyzed, where M2 macrophages, neutrophils, and follicular helper T cells were significantly different between the two groups (Figure 9B). Concurrently, the correlation between various types of immune cell infiltrations and THBS1 expression in IgG4-RD samples was analyzed. A negative correlation was observed between the expression of THBS1 and number of M0 macrophages, CD4-naive T cells, gamma/delta T cells, and memory B cells (Figure 10A–C, G). Eosinophils, resting natural killer cells, and CD8 T cell were positively correlated with THBS1 expression (Figure 10D–F). The composition of various types of immune cells was also significantly correlated in IgG4-RD (Figure 10H).
Differential THBS1 Expression and Enhanced Classical Monocyte Proportions in IgG4-RD Patients
After stringent quality control steps, we harvested a total of 52,176 cells from the PBMCs of four patients diagnosed with IgG4-RD and three healthy controls (Figure 11A). The identified cells comprised a range of cell types including CD8+ NKT-like cells, Classical Monocytes, Effector CD4+ T cells, Effector CD8+ T cells, Memory B cells, Memory CD4+ T cells, Myeloid Dendritic cells, Naive CD4+ T cells, Naive CD8+ T cells, Natural killer cells, Non-classical monocytes, Plasma B cells, Plasmacytoid Dendritic cells, Pre-B cells, and Pro-B cells (Figure 11A-B). Further analysis of the proportion of each cell type in both the disease and control groups yielded informative results. We observed that certain cell types, specifically Non-classical monocytes, Plasma B cells, Effector CD4+ T cells, and Classical Monocytes, constituted a higher proportion of the total cell population in the PBMCs of IgG4-RD patients compared to other cell types (Figure 11C-D). Our analysis was extended to explore gene expression patterns within these identified cell types. We discovered that the gene THBS1, implicated in various cellular processes, is primarily expressed in Classical Monocytes. Notably, this cell type was found to be more abundant in the PBMCs of IgG4-RD patients. (Figure 11E) These findings underscore the potential role of these cell types, and specifically the THBS1-expressing Classical Monocytes, in the pathogenesis of IgG4-RD, paving the way for more targeted therapeutic strategies.