Front. Immunol.Frontiers in ImmunologyFront. Immunol.1664-3224Frontiers Media S.A.10.3389/fimmu.2025.1665086ImmunologyOriginal ResearchB cell subpopulations and their role in the pathogenesis of primary Sjögren’s syndrome: insights from single-cell RNA sequencingZhangXiaoyu1†HanYuanwei2†ZhangJing3SongYu4ZhouYan5WangJia6TianJinhai6XinXinran1LiuBin1*FangLi78*ZhuHong5*1Department of Rheumatology, Affiliated Hospital of Qingdao University, Qingdao, Shangdong, China2The First Clinical Medical College of Ningxia Medical University, Yinchuan, China3Department of Endodontics, Stomatology Hospital, General Hospital of Ningxia Medical University, Yinchuan, China4Department of Oral and Maxillofacial Surgery, Stomatology Hospital, General Hospital of Ningxia Medical University, Yinchuan, China5Department of Rheumatology, General Hospital of Ningxia Medical University, Yinchuan, China6Biochip Center, General Hospital of Ningxia Medical University, Yinchuan, China7School of Health Management, Ningxia Vocational and Technical College for Minorities, Wuzhong, China8School of Nursing, Ningxia Medical University, Yinchuan, China
Edited by: Ting Wang, University of Washington, United States
Reviewed by: Yanmei Li, Tianjin Medical University General Hospital, China
Xiujia Yang, Southern Medical University, China
*Correspondence: Bin Liu, binliu72314@163.com; Li Fang, 1301705392@qq.com; Hong Zhu, nxzhuh@126.com
†These authors have contributed equally to this work
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Objective
This study aims to investigate the potential role of B cells in the pathogenesis of Primary Sjögren’s Syndrome (pSS) by analyzing cell types, differentially expressed genes, and associated signaling pathways using single-cell RNA sequencing.
Methods
Peripheral blood mononuclear cells (PBMCs) from 3 pSS patients and 3 healthy controls (HCs) were collected. Single-cell transcriptomic analysis was performed, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, transcription factor analysis, pseudotime analysis, cell communication analysis, and B cell receptor (BCR) repertoire analysis. Genes and pathways potentially involved in the pathogenesis of pSS were identified, and key genes were validated by qRT-PCR. Statistical significance was assessed using T-tests and the Wilcoxon rank-sum test, with a p-value < 0.05 considered statistically significant.
Results
Single-cell RNA sequencing of peripheral blood B cells from three patients with primary Sjögren’s syndrome (pSS) and three healthy controls (HCs) identified three subpopulations: memory B (Bmem), naïve B (NaiveB), and plasma cells (PlasmaCells). In pSS, differentially expressed genes were enriched in Type I interferon signaling, antigen processing/presentation, and MHC class II binding. Transcription factors related to interferon responses, including NR2F6, IRF5, STAT2, and IRF9, were upregulated. Cell–cell communication analysis highlighted frequent interactions via TNFSF10–TNFRSF10C and TGFB1–TGFBR3. Pseudotime analysis indicated accelerated NaiveB differentiation along the effector branch. B cell receptor repertoire analysis revealed increased IGHV4-34 usage and higher IGHJ4/IGHJ6 usage in PlasmaCells, with reduced IGHV1-3, IGHV1-69D, and IGHV2-7D usage. qRT-PCR validation in 22 pSS patients and 22 HCs confirmed significant ISG15 upregulation (p < 0.0001).
Conclusion
B cells contribute to the pathogenesis of pSS through the Type I IFN signaling pathway mediated by genes such as ISG15, alterations in BCR clonality, IGHV-J gene rearrangements, and abnormal gene usage.
pSS is a chronic inflammatory autoimmune disease characterized by lymphocyte proliferation and progressive damage to exocrine glands. Patients with pSS have multiple autoantibodies in their serum. In addition to dysfunction of the salivary and lacrimal glands, multi-organ and multi-system involvement may occur (1). The etiology of pSS remains unclear, but studies suggest that disruption of the innate immune barrier, through mechanisms involving IFN pathways, plays a key role in the pathogenesis of SS, particularly in the early stages of the disease (2). Research indicates that B cells play a central role in the pathogenesis of pSS (3), although there are still differing views on the specific mechanisms by which B cells contribute to pSS. Clarifying the pathogenesis is crucial.
Single-cell RNA sequencing (scRNA-seq) is a new technology for high-throughput RNA sequencing and analysis at the single-cell level. scRNA-seq provides insights by analyzing differences between cells and subpopulations of cells. Studies have used scRNA-seq to analyze differential expression in synovial cells of rheumatoid arthritis patients, offering new insights into the pathology and heterogeneity of rheumatoid arthritis, and providing information for novel targeted therapies (4). However, research on scRNA-seq in pSS is limited. This study performed scRNA-seq on peripheral blood cells from pSS patients and healthy individuals, and analyzed the B cell subpopulations in depth, revealing differences in expression across B cell subgroups. This provides new insights into the involvement of B cells in the pathogenesis of pSS.
Materials and methodsStudy subjects
A total of 25 patients diagnosed with pSS for the first time at the Rheumatology and Immunology Department of Ningxia Medical University General Hospital were enrolled as the pSS group. Additionally, 25 healthy individuals from the same hospital who underwent physical examinations were included as the HC group. Three subjects from each group were randomly selected for single-cell sequencing, while the remaining 22 subjects in each group were used for qRT-PCR analysis. The pSS diagnostic criteria followed the 2002 and 2016 classification standards established by the American College of Rheumatology (ACR) or the European League Against Rheumatism (EULAR) (5). Exclusion criteria: patients with other systemic autoimmune diseases (including systemic lupus erythematosus, dermatomyositis, rheumatoid arthritis, adult-onset Still’s disease, etc.), vascular diseases, infectious diseases, hematologic diseases, tumors, neuropsychiatric disorders, pregnancy, and other conditions. The study protocol was approved by the Ethics Committee of Ningxia Medical University General Hospital, and all participants provided informed consent (KYLL-2024-0327).
Research methodsPreparation of single-cell suspension
PBMCs were isolated using density gradient centrifugation with lymphocyte separation medium (Ficoll-Paque Plus, GE Healthcare) and washed with PBS without calcium and magnesium. To remove red blood cells, 2 mL of GEXSCOPE® red blood cell lysis buffer (RCLB, Singleron) was added at 25 °C for 10 minutes. The solution was then centrifuged at 500×g for 5 minutes and resuspended in PBS. The blood sample was centrifuged at 400g for 5 minutes at 4 °C, and the supernatant was discarded. After removal of red blood cells, the PBMCs were separated by centrifugation at 400g for 10 minutes at 4 °C. The supernatant was discarded, and the PBMCs were resuspended in PBS to obtain a single-cell suspension. The cell viability was assessed using Trypan Blue staining, with cell viability greater than 90% based on microscopy counting.
RT & Amplification & Library Construction
For single-cell sequencing, 3 pSS patient samples and 3 healthy control samples (2×105 cells/mL, 100 μL) were loaded onto a SCOPE-chip™ microfluidic chip. Libraries were constructed following the protocol of sCircle® Single Cell Full Length Immuno_BCR Library Kit (Biotechnologies). Specifically, poly(A) tails were captured by magnetic beads with molecular markers. Cells and mRNA were labeled after the cells were lysed. The magnetic beads in the chip were collected, and mRNAs were reverse-transcribed into complementary DNA (cDNA) and amplified. After local cDNAs were fragmented and spliced, transcriptome sequencing libraries suitable for the Illumina sequencing platform were constructed. The remaining cDNA was enriched to the full length immune receptor (BCR) by three rounds. Then the enriched products were fragmented and spliced to construct the BCR sequencing libraries, suitable for the Illumina sequencing platform. Finally, sequencing of the libraries was performed on Illumina Nova 6000, with a pair-end length of 150bp.
Quality control, dimensionality reduction, and clustering
The raw sequencing reads were processed using CeleScope v1.14.0 (Singleron Biotechnology), with default parameters. Briefly, Barcodes and UMIs were extracted from R1 reads and corrected. Adapter sequences and poly A tails were trimmed from R2 reads and the trimmed R2 reads were aligned against the GRCh38 (hg38) transcriptome using STAR(v2.6.1b). Uniquely mapped reads were then assigned to exons with FeatureCounts(v2.0.1) (6). Successfully Assigned Reads with the same cell barcode, UMI and gene were grouped together to generate the gene expression matrix for further analysis. For each dataset, quality control, dimensionality reduction, and clustering analysis were conducted using Scanpy v1.8.2 in the Python 3.7 environment (7). The following filtering criteria were applied: Exclude cells with gene counts below 200 or in the top 2% of gene counts; Exclude cells with unique molecular identifier (UMI) counts in the top 2%; Exclude cells with mitochondrial gene content exceeding 30%; Exclude genes expressed in fewer than 5 cells.
After filtering, 65670 cells were retained for the downstream analyses. The raw count matrix was normalized by total counts per cell and logarithmically transformed into normalized data matrix. Top 2000 variable genes were selected by setting flavor = ‘seurat_v3’. Principle Component Analysis (PCA) was performed on the scaled variable gene matrix, and top20 principle components were used for clustering and dimensional reduction. Batch effect between samples was removed by Harmony v1.0 (8). Cells were separated by using Louvain algorithm and setting resolution parameter at 1.2. Cell clusters were visualized by using Uniform Manifold Approximation and Projection (UMAP) (9).
Differential gene identification and cell type annotation
To identify differentially expressed genes (DEGs), we used the scanpy.tl.rank_genes_groups() function based on Wilcoxon rank sum test with default parameters, and selected the genes expressed in more than 10% of the cells in either of the compared groups of cells and with an average log(Fold Change) value greater than 1 as DEGs. Adjusted p value was calculated by benjamini-hochberg correction and the value 0.05 was used as the criterion to evaluate the statistical significance.
The cell type identification of each cluster was determined according to the expression of canonical markers from the reference database SynEcoSysTM (Singleron Biotechnology). SynEcoSysTM contains collections of canonical cell type markers for single-cell seq data, from CellMakerDB, PanglaoDB and recently published literatures. Cell doublets were estimated based on the expression pattern of canonical cell markers. Any clusters enriched with multiple cell type-specific markers were excluded for downstream analysis.
Pathway enrichment analysis
To investigate the potential functions , Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were used with the “clusterProfiler” R package v 3.16.1 (10, 11). Pathways with p_adj value less than 0.05 were considered as significantly enriched. Selected significant pathways were plotted as bar plots.
Transcription factor regulatory network analysis
Transcription factor network was constructed by pyscenic (v0.11.0) using scRNA expression matrix and transcription factors in AnimalTFDB. First, GRNBoost2 predicted a regulatory network based on the co-expression of regulators and targets. CisTarget was then applied to exclude indirect targets and to search transcription factor binding motifs. After that, AUCell was used for regulon activity quantification for every cell. Cluster-specific TF regulons were identified according to Regulon Specificity Scores (RSS) and the activity of these TF regulons were visualized in heatmaps (12).
Cell-cell communication network analysis
Cell-cell interaction (CCI) were predicted based on known ligand–receptor pairs by Cellphone DB (v2.1.0) version (13). Permutation number for calculating the null distribution of average ligand-receptor pair expression in randomized cell identities was set to 1000. Individual ligand or receptor expression was thresholded by a cutoff based on the average log gene expression distribution for all genes across each cell type. Predicted interaction pairs with p value <0.05 and of average log expression > 0.1 were considered as significant and visualized by heatmap_plot and dot_plot in CellphoneDB.
Pseudotime analysis
Cell differentiation trajectory was reconstructed with the Monocle2 v 2.10.0. For constructing the trajectory, top 2000 highly variable genes were selected by FindVairableFeatures, and dimension-reduction was performed by DDRTree. The trajectory was visualized by plot_cell_trajectory function in Monocle2 (14, 15).
VDJ analysis
ScBCR clonotype assignment were performed using CeleScope vdj pipeline v1.14.0 (Singleron Biotechnology), with GRCh38 as reference. In brief, a BCR diversity metric, containing clonotype frequency and barcode information, was obtained. For the BCR, only cells with one productive IGH chain and one productive IGK/IGL chain were kept for further analysis. Each unique IGH-IGK/IGL pair was defined as a clonotype.If one clonotype was present in at least two cells, cells harboring this clonotype were considered to be clonal and the number of cells with such pairs indicated the degree of clonality of the clonotype.
B cell isolation by Ficoll gradient and CD19-positive selection
PBMCs were obtained from fresh EDTA-anticoagulated peripheral blood of pSS patients and healthy controls by Ficoll-Paque density gradient centrifugation (Cat#LTS1077, TBD). PBMCs were then incubated with CD19 MicroBeads (Cat#17954, STEMCELL Technologies) and positively selected according to the manufacturer’s protocol. Purity of the isolated B cells was determined by flow cytometry using an anti-human CD19-APC antibody (Cat#309512, BioLegend) and consistently exceeded 95%. Total RNA was extracted from the purified B cells using TRIzol reagent (Cat#15596018CN, Invitrogen), quantified and assessed for purity, and subsequently subjected to quantitative PCR (qRT-PCR) analysis.
qRT-PCR
For qRT-PCR analysis, 22 pSS patient samples and 22 healthy control blood samples (≥5 mL) were collected. Total RNA was extracted using the Biotech RNA extraction kit RP4002, and RNA concentration and purity were measured using a NanoDrop 2000 spectrophotometer (Thermo Scientific). Samples with an A260/A280 ratio of 1.8–2.0 and a concentration >50 ng/μL were stored at -80 °C. Quantitative reverse transcription was performed using the PrimeScript RT kit (Takara Bio). 200 ng RNA, 2 μL PrimeScript RT Master Mix (Takara Bio), and RNase-free water were mixed to a final volume of 10 μL to complete the reverse transcription and obtain cDNA. qPCR was performed using 2 μL cDNA, 10 μL TB-Green Premix Ex Taq (Takara Bio), 6.4 μL RNase-free water, and 0.8 μL primers (forward and reverse). The reaction was performed on a Lightcycler 480 System with the following conditions: 95 °C for 3 s, followed by 40 cycles of denaturation at 95 °C for 5 s and extension at 60 °C for 30 s. Each sample was tested in triplicate, and the difference in Ct values (fluorescence threshold cycle number) was less than 0.5. The relative expression level of ISG15 was analyzed using the 2-ΔΔCt method. The primer sequences used were as follows: ISG15 forward primer: 5’-CGCAGATCACCCAGAAGATCG-3’, reverse primer: 5’-TTCGTCGCATTTGTCCACCA-3 ‘; GAPDH forward primer: 5’-CCACGGCTGCTTCCAGCTCC-3’, reverse primer: 5’-GGACTCCATGCCCAGGAAGGAA-3 ‘.
Statistical analysis
Cell distribution comparisons between two groups were performed using the unpaired two-tailed Wilcoxon rank-sum test. Comparisons of gene expression or gene features between two groups were performed using the unpaired two-tailed Student’s t-test. Non-normally distributed quantitative data were presented as median and percentiles [M(P25, P75)], and group comparisons were performed using the Mann-Whitney U test. Statistical analyses and data presentation were performed using R and Python. A p-value of <0.05 was considered statistically significant.
ResultsAnalysis of cell types and differential genes in peripheral blood single cells of pSS and HC groups
A total of 65,670 cells were captured, with 34,221 from the pSS group and 31,449 from the HC group (Supplement table1). Unsupervised clustering of all cells identified 20 clusters (Supplementary Figure S1A), which were subsequently annotated as eight major cell types based on canonical marker gene expression (Supplementary Figure S5B): B cells, T cells, natural killer (NK) cells, neutrophils, basophils, mononuclear phagocytes (MPs), plasmacytoid dendritic cells (pDCs), and erythrocytes (Figures 1A, B). Among the 2,938 captured B cells, unsupervised of all cells identified 4 clusters (Supplementary Figure S1C). Based on canonical marker gene expression,the results showed (Supplementary Figure S1D): Cluster 1 and Cluster 3 were identified as NaiveB (IGHD,IL4R,IGHM),Cluster 2 was classified as Bmem(CD27,TNFRSF13B,ANXA2), and Cluster 4 was determined to be PlasmaCells(JCHAIN,MZB1, IGHG1) (Figure 1C). In comparison to the HC group, the pSS group exhibited trends of increased median proportions in B cells (4.40% vs 3.59%), MPs (11.53% vs 2.17%), NK cells (7.37% vs 4.29%), basophils (0.91% vs 0.49%), and pDCs (0.07% vs 0.06%). Conversely, trends of decreased median proportions were observed in T cells (24.87% vs 31.92%) and neutrophils (49.98% vs 58.47%) (Figure 1D). Compared to the HC group, the pSS group exhibited elevated median proportions of NaiveB (71.70% vs 67.65%) and PlasmaCells (2.70% vs 0.67%), while the median proportion of Bmem was reduced (25.61% vs 32.35%) (Figure 1E).The top 10 differentially expressed genes in B cells were IGHM, IGLC2, MS4A1, IGHA1, IGLC3, IGHD, CD74, CD79A, IGHG1, and BANK1 (Figure 1F). A heatmap of the top 10 upregulated genes in each B cell subtype revealed that NaiveB expressed high levels of TCL1A, IGHD, and FCER2, PlasmaCells showed high expression of IGHA1, JCHAIN, and MZB1, while Bmem cells had elevated levels of AIM2, LINC01781, and ITGB1 (Figure 1G).
PBMC Cluster Types and Differential Genes in the pSS Group and HC Group. (A) UMAP plot of cell dimensionality reduction for the HC group; (B) UMAP plot of cell dimensionality reduction for the pSS group; (C) UMAP plot of B cell dimensionality reduction; (D) Proportion of PBMC cell composition; (E) Proportion of B cell composition; (F) Heatmap of top 10 differentially expressed genes in each PBMC cell cluster: yellow indicates high expression, purple indicates low expression; (G) Heatmap of top 10 differentially expressed genes in B cell subsets.
UMAP plots A and B depict clusters of immune cells, including B cells, NK cells, T cells, neutrophils, basophils, MPs, and pDCs. Plot C shows clusters of Bmen, NaiveB, and PlasmaCells. Bar chart D compares cell proportions in healthy controls (HC) and primary Sjögren's syndrome (pSS), while bar chart E highlights differences in Bmen, NaiveB, and PlasmaCells. Heatmap F shows gene expression across various cell types, and heatmap G details gene expression in Bmen, NaiveB, and PlasmaCells.Differential gene comparison and functional enrichment analysis of B cell subtypes in the pSS group
In this study, differentially expressed genes (DEGs) in the three B cell subtypes (Bmem, NaiveB, and PlasmaCells) from patients with pSS were analyzed using GO and KEGG pathway enrichment. For the Bmem subtype, upregulated genes such as MX1, IFI44L, ISG15, and STAT1 were primarily enriched in the Type I interferon signaling pathway and immune response-related processes. Downregulated genes, including H1-10, FOSB, and JUN, were mainly involved in ribosomal functions and mRNA binding (Figure 2A). KEGG analysis showed that upregulated genes were enriched in pathways such as "Proteasome," "Antigen processing and presentation," and "Epstein-Barr virus infection" (Figure 2B), while downregulated genes were associated with pathways like "Ribosome," "COVID-19," and "IL-17 signaling pathway" (Figure 2C).GO analysis revealed that upregulated genes were enriched in "Type I interferon signaling pathway," "Regulation of innate immune response," and "Viral response" at the biological process (BP) level (Figure 2D), while downregulated genes were enriched in "Ribosome" and "RNA binding" at the cellular component (CC) and molecular function (MF) levels (Figure 2E).
Inter-group Comparison and Functional Analysis of Differential Genes in the Bmem Subset. (A) Heatmap of top 20 differential gene expression comparisons in the Bmem subset; (B) KEGG pathway analysis of upregulated genes in the Bmem subset of the pSS group; (C) KEGG pathway analysis of downregulated genes in the Bmem subset of the pSS group; (D) GO term analysis of upregulated genes in the Bmem subset of the pSS group; (E) GO term analysis of downregulated genes in the Bmem subset of the pSS group.
Five-panel scientific figure detailing gene expression and pathway analysis: A. Heatmap showing gene expression levels across samples, with a gradient from blue to yellow indicating low to high expression. B. Bar chart listing pathways like Proteasome and Influenza A with corresponding significance values. C. Bar chart highlighting pathways significant for Coronavirus and IL-17 signaling. D. Dot plot of biological processes, cellular components, and molecular functions with color-coded ontology categories. E. Dot plot focusing on protein targeting and translation processes, categorized by ontology.
In the NaiveB subtype, upregulated genes such as IFI44L, ISG15, and STAT1 were associated with Type I interferon responses and viral immune responses. The downregulated genes, including RASAL2, FOSB, and KLF6, were enriched in translation initiation and ribosomal processes (Supplementary Figure S2A). KEGG analysis highlighted the enrichment of upregulated genes in "Prion disease" "Oxidative phosphorylation" and "Antigen processing and presentation" (Supplementary Figure S2B), while downregulated genes were primarily enriched in "Ribosome" "COVID-19" and "Th17 cell differentiation" pathways (Supplementary Figure S2C). GO analysis indicated that upregulated genes in NaiveB cells were mainly enriched in "Response to type I interferon" and "Response to virus" at the BP level, "Phagosome" and "Mitochondrial inner membrane" at the CC level, and "MHC class II protein complex binding" at the MF level (Supplementary Figure S2D). Downregulated genes were enriched in "Ribosome" and "RNA binding" at the CC and MF levels (Supplementary Figure S2E).
For PlasmaCells, the upregulated genes such as ISG15, IFI6, and STAT1 were linked to antigen processing, MHC class II binding, and peroxidase activity. Downregulated genes like FRAT2, SLC25A27, and DIP2B were associated with ribosomal functions and translation initiation (Supplementary Figure S3A). KEGG pathway analysis revealed that upregulated genes in PlasmaCells were enriched in "Protein processing in the endoplasmic reticulum" "Antigen processing and presentation" and "Parkinson's disease" (Supplementary Figure S3B), while downregulated genes were primarily associated with "Osteoclast differentiation" and "Endocytosis" (Supplementary Figure S3C) .GO analysis showed that upregulated genes were enriched in "Processing and presentation of exogenous peptide antigens" at the BP level, "Endoplasmic reticulum protein complex" at the CC level, and "Molecular carrier activity" and "MHC class II protein complex binding" at the MF level (Supplementary Figure S3D). Downregulated genes were mainly involved in "Ribosome" and "RNA binding" processes (Supplementary Figure S3E).
Transcription factor analysis of B cell subtypes in the pSS group
Compared to the HC group, the transcription factors NR2F6 and IRF5 were significantly upregulated in Bmem cells in the pSS group, while KLF4 and JUNB were significantly downregulated. In NaiveB cells, the transcription factors STAT2, IRF9, THRB, and BATF3 were significantly upregulated, whereas TGIF2, KLF4, FOS, BOCH2, and FOXO1 were significantly downregulated. In PlasmaCells, the transcription factors AHR, ZFP64, ARID3A, and XBP1 were significantly upregulated, while ELF2, FOXO1, and POLR3G were significantly downregulated (Figure 3A).
Transcriptional Regulation and Cellular Communication in B Cell Subsets. (A) Heatmap of average transcription factor expression comparisons between groups in B cell subsets; (B) Interaction pair count plot between two cell types: The outer ring represents cell types; the inner ring shows red for ligands and blue for receptors; line clarity is positively correlated with the number of interaction pairs between the two cell types, with clearer lines indicating more interactions; (C) Heatmap of ligand-receptor pair counts between two cell types in the pSS group: Darker colors represent a higher number of interaction pairs between the two cell types; (D) Heatmap of ligand-receptor pair counts between two cell types in the HC group.
Four-panel image showing heatmaps and a chord diagram related to cell types. Panel A is a heatmap displaying gene expression across different conditions with a color scale from blue to red. Panel B is a circular chord diagram illustrating interactions between various cell types with colored connections. Panels C and D are heatmaps depicting interactions among specific cell types, with color intensity indicating interaction levels. Each panel includes color bars and legends for reference.Cell communication network analysis
The interactions between B cells and other cell types in the pSS group were analyzed, including 30 ranked interaction pairs, chemokine interaction pairs, immune checkpoint interaction pairs, growth factor interaction pairs, and cytokine interaction pairs. The results showed that the interaction pairs between B cells and MPs were the most numerous (Figures 3B–D). As ligands, B cells primarily communicate through the following pairs: LGALS9-HAVCR2, CXCL8-CXCR2, TNFSF10-TNFRSF10C, TGFB1-TGFBR3, HLA-E-KLRC1, FAM3C-HLA-C, among others (Supplementary Figure S4A-E). When B cells act as receptors, the main communication occurs through pairs such as CD28-CD86, CXCL8-CXCR2, TNFSF138-TNFRSF13C, TGFB1-TGFBR3, ARP-CD74, COPA-CD74, etc. (Supplementary Figure S5A-E).
Pseudo-time analysis of B cell subtypes in the pSS and HC groups
Pseudotime trajectory analysis of B cell subtypes in the pSS group (Figure 4A) revealed potential developmental relationships among these subsets. The analysis identified 5 key cellular state nodes and two main developmental branch paths. Naive B cells were predominantly located at the starting point of the trajectory and distributed along the initial segments of both branches. The first branch path was primarily enriched with Naive B cells, while the second branch path indicated a progression toward effector B cell differentiation, with Bmem cells enriched in intermediate states along this branch (Figure 4B). Cell states along the pseudotime trajectory were classified into 11 distinct states. State 9 was in the early left branch, followed by State 8 before the branching point. Transitional States 3–7, 10, and 11 clustered around the branch root. The long right branch was mainly State 1 at late pseudotime, and State 2 appeared rarely at the end of the left branch (Figure 4C). Comparison with the HC group showed that in the pSS group, the number of Naive B cells was relatively reduced at the trajectory starting point but significantly increased at the initial segment of the second branch, which functions as the effector branch. PlasmaCells were localized at the terminus of the second branch in both the HC and pSS groups (Figure 4D). In the HC group, naive B cells were mainly at the origin and decreased over time. In pSS, fewer naive B cells were present at the origin, showing a rise then decline. Memory B cells localized to intermediate nodes in both groups, but were more abundant at terminal nodes in pSS. Plasma cells were terminal in both, with higher proportion in pSS. Integrated all_group data showed dynamic changes: naive B cells increased then decreased, memory B cells were mid-distributed with late presence, and plasma cells were terminal (Figure 4E). Further analysis of marker gene expression patterns along pseudotime classified them into 6 clusters (Figure 4F): Early-stage clusters (1–2) included JUN, NFKBIA, DUSP1, FOS (stress-responsive/MAPK, NF-κB) and FCER2, IGHD, TCL1A (B cell activation/maturation), showing progressive downregulation. Late-stage clusters (4–6) comprised immunoglobulin/plasma cell genes (IGKV4-1, IGHG, IGHA, IGLC, JCHAIN, XBP1, MZB1), B2M (MHC I), and ITGB1, HBB, HBA2 (adhesion/oxygen transport), all upregulated toward pseudotime end. Cluster 3 (AC099560.1, RPS3A) remained low, indicating suppression of ribosomal/metabolic programs. Overall, expression shifted from early stress and immune activation to late antibody production, antigen presentation, and structural remodeling.
Pseudo-temporal trajectory plots of B cell subpopulations. (A) Pseudo-temporal trajectory plot of B cell evolution; (B) Distribution characteristics of each B cell subpopulation in the pseudo-temporal differentiation trajectory; (C) Distribution of different cell subpopulations in the pseudo-temporal differentiation trajectory; (D) Distribution of pseudo-temporal trajectories of each B cell subpopulation across different groups; (E) Distribution of major time points of B cell clusters in different groups (HC, pSS, and all_group combined) along pseudo-time; the y-axis represents pseudo-time (from bottom to top), the x-axis shows the proportion of cell types at different time points, and different colors represent different cell types; The width of each colored band reflects the proportion of the corresponding cell type at a given pseudotime position; (F) Gene expression changes along pseudo-time; the x-axis represents pseudo-time, and the y-axis shows the gene expression levels.
Panel A shows a scatter plot of pseudotime values for components one and two, illustrating a progression over time. Panel B depicts clusters labeled Bmem, NaiveB, and PlasmaCells with distinct color codes. Panel C presents states one through eleven in various colors using similar components. Panel D compares HC and pSS clusters within the same component space. Panel E uses violin plots for HC, pSS, and all groups, highlighting distributions of Bmem, NaiveB, and PlasmaCells. Panel F displays a heatmap of gene expressions, including JUN and IGHG1, with clusters one to seven indicated by color.BCR analysis
BCR data analysis of 6 samples from the pSS and HC groups revealed that the clonal frequencies of all B-cell subsets in both groups were predominantly single-clonotype-dominated . The proportion of large-scale clones (clonotype frequency > 10) was 0% across all subsets. The proportion of medium-scale clones (clonotype frequency > 1 and ≤ 10) was low in both groups (Figure 5A). Both NaiveB and Bmem cells exhibited clonal expansion in samples from both groups; however, neither their single-clonotype frequencies nor medium-scale clonotype frequencies showed statistically significant differences between the two groups. In contrast, PlasmaCells displayed clonal expansion in all pSS samples (predominantly single-clonotype), whereas in the HC group, clonality was detected only in HC03 (Supplement table 2). Further analysis of BCR clonal diversity using the D50 diversity index demonstrated no statistically significant differences between the pSS and HC groups (p = 0.7) (Figure 5B).
Comparison of BCR clonotype distribution and IGHV/IHJ gene usage in B cell subpopulations. (A) The number and proportion of different BCR clonotypes between groups: Circle size represents the total number of clonotypes (Single represents a single clonotype; Medium represents clonotypes with a frequency greater than 1 but less than or equal to 10; Large represents clonotypes with a frequency greater than 10). The colors within the circle represent the proportion of clonotypes with different frequencies. (B) Inter-group comparison of BCR clonal diversity: Analyzed using D50 diversity index. The D50 diversity index represents the number of clonotypes required to account for 50% of the total BCR repertoire; lower values suggest reduced diversity; (C) Comparison of IGHV gene usage frequency in B cell clones between the two groups. (D) Comparison of IGHJ gene usage frequency in B cell clones between the two groups. (E) IGHV gene usage frequency in B cell subpopulations of the HC group. (F) IGHV gene usage frequency in B cell subpopulations of the pSS group.
Panel A shows pie charts of clone sizes in pSS and HC groups across different cell types, categorized as single, medium, or large. Panel B is a bar graph comparing the D50 diversity index between HC and pSS. Panel C presents a bar graph showing percentages of different genes in HC and pSS. Panel D compares IGHJ gene percentages between HC and pSS. Panel E and F are heatmaps illustrating gene usage across different cell types, with color intensity indicating frequency.
Further analysis of the IGHV and IGHJ genes of B cells from the pSS and HC groups showed that in the pSS group, the top 5 IGHV gene segments with the highest clonotype usage frequencies were: IGHV3-23, IGHV4-34, IGHV4-39, IGHV3-33, and IGHV4-59. In the HC group, the top 5 IGHV gene segments were: IGHV3-33, IGHV3-23, IGHV4-39, IGHV4-59, and IGHV1-69D. Compared to the HC group, the frequency of IGHV4-34 usage was higher in the pSS group (median = 9.3% vs 3.9%), although the difference did not reach statistical significance (p = 0.100, Mann-Whitney U = 0.000) (Figure 5C).
In terms of IGHJ gene segments, both pSS and HC groups showed relatively high usage of IGHJ3, IGHJ4, IGHJ5, and IGHJ6 (Figure 5D). Further analysis of the usage frequency of heavy chain variable region (V) genes between the two groups in each B cell subtype showed that in the PlasmaCells subtype, the frequency of IGHJ4 and IGHJ6 usage was higher in the pSS group compared to the HC group, while the frequency of IGHJ1 and IGHJ5 usage was lower (Figure 5E). The usage frequency of IGHV gene segments in the IGHV1-3, IGHV1-69D, and IGHV2-7D subtypes was relatively decreased in the pSS group. No significant changes in heavy chain V region gene usage frequencies were observed in the Bmem and NaiveB subtypes (Figure 5F).
qRT-PCR analysis
From the above data, it was found that in the B cell subtypes of pSS patients, type I IFN-related genes such as ISG15, IFI44L, and IFI44 were significantly upregulated. Both functional enrichment and cell communication analyses suggested that the IFN pathway may be involved in the pathogenesis of pSS. To further validate the reliability of these results, the ISG15 gene was selected for qRT-PCR verification. ISG15 was signifificantly overexpressed in the peripheral blood of pSS patients (p < 0.0001) and in purifified B cells (p = 0.0076) (Figure 6).
Comparison of ISG15 expression levels between the HC group and the pSS group.
Two dot plots comparing the relative expression of ISG15 in peripheral blood mononuclear cells (PBMC) and B cells. The left plot shows a significant increase in the pSS group compared to HC (P<0.0001). The right plot also shows higher expression for pSS versus HC (P=0.0076).Discussion
pSS is a complex and heterogeneous disease, and its pathogenesis remains unclear. RNA-seq technology retains the transcriptional differences of different cells, providing a significant advantage in identifying cell subtype characteristics and cell interactions. In this study, we performed comprehensive bioinformatics analysis of the transcriptome to study cell clustering, gene differences, and related pathways in pSS. First, we clustered a total of 65,670 cells from 6 samples of peripheral blood from the pSS and HC groups. Compared to the HC group, the pSS group exhibited an elevated median proportion of B cells, suggesting that B cell activation might be related to the onset of pSS. Subsequently, we annotated the B cell populations and identified three B cell subtypes: Bmem, NaiveB, and PlasmaCells. Additionally, the median proportions of both NaiveB and PlasmaCells were higher in the pSS group. A comparison of differential genes among the three B cell subtypes revealed that, compared to the HC group, the upregulated genes in both Bmem and NaiveB cells in the pSS group included ISG15, IFI44L, IFI44, STAT1, IFI6, IFI3, IFIT3, and MX1. The upregulated genes common to Bmem, NaiveB, and PlasmaCells in the pSS group included ISG15, IFI6, and STAT1. These genes are mainly involved in functions such as “encoding interferon-induced proteins,” “encoding ISG15 ubiquitin-like modifiers,” and “encoding STAT1 signal transduction and transcription activation factor 1.” Differential gene comparison showed that the upregulated genes in all three B cell populations of the pSS group were involved in the type I interferon response pathway. Luo S et al. found that type I IFN stimulates monocyte differentiation and induces immature dendritic cells to express chemokines and costimulatory molecules, promoting the onset of SLE (16). Cui Y et al. used single-cell RNA-seq to explore the common molecular mechanisms between SLE and primary pSS and found that IFN response and ITGB2 signaling pathways play crucial roles in both diseases (17). It is speculated that the type I interferon signaling pathway mediates B cell involvement in the pathogenesis of pSS.
To further investigate the role of the type I IFN response in B cells in pSS, GO and KEGG pathway enrichment analyses of differential genes in B cell subpopulations showed that, compared to the HC group, the differentially upregulated genes in NaiveB and Bmem cells in the pSS group were involved in the type I interferon signaling pathway at the BP level. In this study, we observed significant upregulation of interferon signaling pathway–related genes (such as ISG15, IFI44L, and IFI44) in B cells from pSS patients, along with marked increases in upstream transcription factors (STAT2, IRF9, and IRF5). Previous studies have shown that IFI44L is a type I IFN-stimulated gene, and its upregulation in pSS patients has been observed (18). ISG15 is a ubiquitin-like protein that, after activation by IFN-α and IFN-β, conjugates with intracellular target proteins. Cui Y et al. found that ISG15 expression levels in the saliva and serum of pSS patients were higher than those in controls and identified IFI44L and ISG15 as common hub genes in both pSS and SLE (17). Given that ISG15 is a key marker of IFN-I activation and can representatively reflect pathway activity, we selected ISG15 for subsequent validation. To further validate this, we performed qRT-PCR analysis of ISG15 expression in the peripheral blood of 22 pSS patients and 22 healthy controls, showing that ISG15 expression was significantly higher in the peripheral blood of pSS patients compared to the HC group. In addition, we examined sorted B cells and found that ISG15 expression was also markedly elevated in B cells from pSS patients. This suggests that B cells may be involved in the pathogenesis of pSS through IFN-related genes such as ISG15.
Our study further performed transcription factor analysis of B cell subpopulations. In Bmem cells, NR2F6 and IRF5 were significantly upregulated. The nuclear receptor transcription factor NR2F6 is a member of orphan nuclear receptors. Natascha Hermann-Kleiter et al. have demonstrated that NR2F6 antagonizes the ability of Th0 and Th17 CD4(+) T cells to induce IL-2 and IL-17 expression, suggesting that NR2F6 may be involved in the pathogenesis of pSS (19). IRF5 is a regulator of type I IFN and IFN-stimulated genes (ISGs). IFN induces STAT activation, which in turn triggers ISG expression. STAT1, STAT2, and IRF9 amplify the JAK-STAT signaling pathway to enhance the IFN response, and the JAK-STAT pathway transduces intracellular signals for various cytokines, which is crucial for the pathogenesis of autoimmune diseases (20, 21). Ivashkiv LB et al. found that a common feature of SLE patients is the increased expression of type I IFN, and IRF5 expression was significantly elevated in the peripheral blood mononuclear cells of SLE patients (22). Su Song et al. found that inhibiting IRF5 expression can prevent the onset and severity of SLE. Our study showed that IRF5 is highly expressed in B cells, suggesting that the increased expression of IRF5 in Bmem cells in the pSS group may be related to the disease's development (23). KLF4 was significantly downregulated in Bmem cells. KLF4 is an evolutionarily conserved zinc finger transcription factor that regulates various cellular processes, such as cell growth, proliferation, and differentiation. Tao H et al. found that KLF4 promotes dentinogenesis and odontoblast differentiation through the regulation of TGF-β signaling and interaction with histone acetylation (24). Some pSS patients exhibit extensive tooth loss and multiple dental caries, suggesting that KLF4 might be involved in the mechanisms of dental damage in pSS patients. Our study showed that KLF4 is significantly downregulated in pSS, while it is upregulated in the normal group. Thus, we speculate that the downregulation of KLF4 in pSS may suppress dentin growth, leading to tooth loss and caries in these patients.
In NaiveB cells, STAT2, IRF9, THRB, and BATF3 were significantly upregulated. STAT is an essential transcription factor in the type I IFN-mediated signaling pathway. STAT2 is defined as an auxiliary factor that participates exclusively in type I IFN (IFN-α, -β, -τ, -ω) signaling transduction; IFN signals through the JAK-STAT pathway, activating the transcription of ISGs. STAT1-2 heterodimers bind to IRF9 to form the activated transcription complex ISGF3. In PlasmaCells, AHR, ZFP64, ARID3A, and XBP1 were significantly upregulated. These findings suggest that multiple transcription factors in B cells participate in the pathogenesis of pSS through the type I IFN-mediated signaling pathway.
Pseudotime trajectory analysis of B cell populations revealed that Naive B cells were predominantly localized at the starting point of the developmental trajectory and distributed along the initial segments of both branches, suggesting their role as the origin of B cell differentiation. The first branch was primarily enriched with Naive B cells, potentially representing a relatively quiescent or self-maintaining state. The second branch delineated the classical differentiation path from Naive B through Bmem to Plasma Cells. Compared to the HC group, the pSS group exhibited a reduced proportion of Naive B cells at the trajectory origin but a significant increase at the initial segment of the second (effector) branch. This uneven distribution implies that more Naïve B cells in pSS may be primed to prematurely enter the effector differentiation path. Collectively, the pseudotime analysis constructed a dynamic map of B cell evolution in pSS, initiating from Naive B states and diverging into two paths: (1) Naive B maintenance and (2) effector. The pSS group was characterized by accelerated and skewed entry of Naïve B cells into the effector path, highlighting a dysregulated B cell maturation trajectory in pSS.
Notably, analysis of marker gene expression patterns along pseudotime revealed dynamic transcriptional changes along the trajectory. The early stage was enriched for stress-response and immediate-early genes (JUN, FOS, NFKBIA) as well as B cell activation–related genes (FCER2, IGHD, TCL1A), whose expression gradually declined over pseudotime. In the intermediate stage, ribosomal and metabolic genes (AC099560.1, RPS3A) remained consistently low. The late stage showed marked upregulation of immunoglobulin and plasma cell differentiation–related genes (IGKV4-1, IGHG, IGHA, IGLC, JCHAIN, XBP1, MZB1), antigen presentation–related genes (B2M), and adhesion and oxygen transport–related genes (ITGB1, HBB, HBA2), indicating a transcriptional profile associated with antibody production and immune effector functions at the trajectory end. These results suggest that cells in the early pseudotime phase are dominated by stress and immune activation–related transcriptional programs, which subsequently shift toward an immune effector–related expression profile.
Cell communication analysis suggested that the number of interactions between B cells and MPs cells was the highest. When B cells acted as ligands in the pSS group, the cell communication with MPs primarily occurred via TGF-β and TNF. TGF-β ligands bind various TGF-β receptors, leading to the recruitment and activation of SMAD family transcription factors that regulate gene expression. These proteins can modulate the expression and activation of interferon γ and tumor necrosis factor α.
BCR analysis of peripheral blood B cells from the pSS and HC groups revealed that both groups exhibited predominantly monoclonal expansion across B cell subsets. While clonal expansion was observed in Naive B and Bmem in both groups, no statistically significant differences were detected between pSS and HC. Notably, plasma cells showed clonal expansion in all pSS samples but were detected in only one HC sample. The lack of significant intergroup differences in BCR clonal diversity may be attributed to the limited sample size (n=3 per group) and substantial individual heterogeneity.
Further analysis of the IGHV and IGHJ genes showed that the use of these genes had changed in the pSS group compared to the HC group. Compared to the HC group, the median usage frequency of IGHV4-34 was higher in pSS (median = 9.3% vs. 3.9%). Arbuckle, Odendahl, and others found that the use of the IGHV4 family increased in peripheral blood B cells in SLE, particularly IGHV4-34 (25, 26). Doorenspleet and others also found an increased frequency of IGHV4-34 use in synovial B cells in early RA (27). In the peripheral blood plasmaCells subset of pSS patients, the frequencies of IGHJ4 and IGHJ6 use were higher, while frequencies of IGHV1-3, IGHV1-69D, and IGHV2-7D subtypes were relatively lower. No significant changes were observed in the heavy chain V region gene usage in NaiveB and Bmem subpopulations. Primary immune thrombocytopenia (ITP) is a disease caused by IgG antibodies against platelets. Studies analyzing BCR libraries in ITP found B cell clones carrying IGHV4-28/IGHJ4 in all ITP patients. There is a close immunological connection between pSS and lymphoma, with the risk of lymphoma significantly increased in long-term pSS patients (28). Xuemin Xue and others found that in B cell lymphoma, rearrangements of IGHV and IGHJ genes were present (29). These studies indicate that abnormal use and rearrangements of IGHV and IGHJ gene segments may contribute to autoimmune diseases and lymphoma. We can also speculate that the rearrangement and abnormal usage of IGHV and IGHJ genes in B cells in pSS patients may contribute to the disease's onset and the development of pSS to lymphoma.
However, this study has certain limitations. The scRNA-Seq data were derived from PBMCs rather than purified B cells, which may limit the resolution of transcriptional features within B cell subpopulations. In addition, experimental validation focused on downstream interferon-effector genes such as ISG15, but not on upstream regulators. Future studies with purified B cell subsets and functional assays will be needed to clarify the role of IFN-I signaling in aberrant B cell activation in pSS.
In summary, this study used scRNA-Seq and chip data analysis to reveal the key genes and related signaling pathways in B cells involved in the pathogenesis of pSS. B cells participate in the pathogenesis of pSS not only through the type I IFN signaling pathway mediated by genes like ISG15, IFI44L, IFI44, STAT1, and IFI6, but also through changes in BCR clonotypes, rearrangements, and abnormal usage of IGHV-J genes.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found here: https://www.jianguoyun.com/p/DTSjHPUQyZLgDRjKnYoGIAA.
Ethics statement
This study was approved by the Medical Ethics Committee of Ningxia Medical University General Hospital (Ethics No. KYLL-2021-585). The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation in this study was provided by the participants’ legal guardians/next of kin. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.
Author contributions
XZ: Funding acquisition, Writing – original draft. YH: Writing – original draft. JZ: Data curation, Writing – original draft. YS: Data curation, Writing – original draft. YZ: Investigation, Writing – original draft. JW: Software, Writing – original draft. XX: Methodology, Writing – review & editing. BL: Writing – review & editing, Methodology. JT: Writing – original draft, Supervision. LF: Writing – review & editing. HZ: Funding acquisition, Investigation, Writing – original draft.
Funding
The author(s) declare financial support was received for the research and/or publication of this article. The study was supported by nature science foundation of Ningxia province (2025AAC020097), Shandong Provincial Natural Science Foundation (ZR2023MHO66), Qingdao Medical and health scientific research project (2024-WJKY160), ‘Clinical Medicine +X’ research project of Affiliated Hospital of Qingdao University (QDFY+X2024210), Shandong Postdoctoral Science Foundation (SDBX2024043) and the Postdoctoral Fellowship Program of CPSF under Grant Number GZC20251551.
Acknowledgments
Xiaoyu Zhang and Yuanwei Han edited the manuscript meticulously. Jing Zhang and Yu Song and Xinran Xin was responsible for figure preparation and data interpretation. Yan Zhou, Jia Wang and Jinhai Tian collected and analyzed data. Hong Zhu, Li Fang and Bin Liu signed and instituted critical revisions. Unanimously, all contributors endorsed and consented to the manuscript's submission.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2025.1665086/full#supplementary-material
ReferencesMarietteXCriswellLA. Primary Sjögren′s Syndrome [J. N Engl J Med. (2018) 378:931–9. doi: 10.1056/NEJMcp1702514, PMID: 29514034Brito-ZerónPBaldiniCBootsmaHBowmanSJJonssonRMarietteX. Sjogren syndrome [J. Nat Rev Dis Prime. (2016) 2:16047. doi: 10.1038/nrdp.2016.47, PMID: 27383445NocturneGMarietteX. B cells in the pathogenesis of primary Sjögren syndrome [J. Nat Rev Rheumatol. (2018) 4:133–45. doi: 10.1038/nrrheum.2018.1, PMID: 29416129ZhangFJonssonAHNathanAMillardNCurtisMXiaoQ. Deconstruction of rheumatoid arthritis synovium defines inflammatory subtypes [J. Nature. 623:616–24. doi: 10.1038/s41586-023-06708-y, PMID: 37938773ShiboskiCHShiboskiSCSerorRCriswellLALabetoulleMLietmanTM. American College of Rheumatology/European League Against Rheumatism Classification Criteria for Primary Sjögren's Syndrome: A Consensus and Data-Driven Methodology Involving Three International Patient Cohorts [J. Arthritis Rheumatol. (20162017) 69:35–45. doi: 10.1136/annrheumdis-2016-210571, PMID: 27789466DobinADavisCASchlesingerFDrenkowJZaleskiCJhaS. STAR: ultrafast universal RNA-seq aligner[J. Bioinformatics. 291:15–21. doi: 10.1093/bioinformatics/bts635, PMID: 23104886WolfFAAngererPTheisFJ. SCANPY: large-scale single-cell gene expression data analysis [J. Genome Biol. (2018) 19:15. doi: 10.1186/s13059-017-1382-0, PMID: 29409532KorsunskyIMillardNFanJSlowikowskiKZhangFWeiK. Fast, sensitive and accurate integration of single-cell data with Harmony[J. Nat Methods. (2019) 16:1289–96. doi: 10.1038/s41592-019-0619-0, PMID: 31740819BechtEMcInnesLHealyJDutertreC-AImmanuelWKwokH. Dimensionality reduction for visualizing single-cell data using UMAP [J. Nat Biotechnol. (2018), 3. doi: 10.1038/nbt.4314, PMID: 30531897YuGWangLGHanYHeQY. clusterProfiler: an R package for comparing biological themes among gene clusters [J. OMICS. (2012) 16:284–7. doi: 10.1089/omi.2011.0118, PMID: 22455463KanehisaMGotoS. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. (2000) 28:27–30. doi: 10.1093/nar/28.1.27, PMID: 10592173Van de SandeBFlerinCDavieKDe WaegeneerMHulselmansGAibarS. A scalable SCENIC workflow for single-cell gene regulatory network analysis [J. Nat Protoc. (2020) 15:2247–76. doi: 10.1038/s41596-020-0336-2, PMID: 32561888EfremovaMVento-TormoMTeichmannSAVento-TormoR. CellPhoneDB: inferring cell–cell communication from combined expression of multi-subunit ligand–receptor complexes [J. Nat Protoc. (2020) 15:1484–506. doi: 10.1038/s41596-020-0292-x, PMID: 32103204QiuXHillAPackerJLinDMaYATrapnellC. Single-cell mRNA quantification and differential analysis with Census [J. Nat Methods. (2017) 14:309–15. doi: 10.1038/nmeth.4150, PMID: 28114287QiuXMaoQTangYWangLChawlaRPlinerHA. Reverse graph embedding resolves complex single-cell developmental trajectories [J. Nat Methods. (2017) 14:979–82. doi: 10.1038/nmeth.4402, PMID: 28825705LuoSWuRLiQZhangG. Epigenetic regulation of IFI44L expression in monocytes affects the functions of monocyte-derived dendritic cells in systemic lupus erythematosus [J. J Immunol Res. (2022) 2022:4053038. doi: 10.1155/2022/4053038, PMID: 35592687CuiYZhangHWangZGongBAl-WardHDengY. Exploring the shared molecular mechanisms between systemic lupus erythematosus and primary Sjögren's syndrome based on integrated bioinformatics and single-cell RNA-seq analysis [J. Front Immunol. (2023) 14:1212330. doi: 10.3389/fimmu.2023.1212330, PMID: 37614232JaraDCarvajalPCastroIBarreraMJAguileraSGonzálezS. Type I interferon dependent hsa-miR-145-5p downregulation modulates MUC1 and TLR4 overexpression in salivary glands from Sjogren’s syndrome patients [J. Front Immunol. (2021) 12:685837. doi: 10.3389/fimmu.2021.685837, PMID: 34149728Hermann-KleiterNGruberTLutz-NicoladoniCThuilleNFresserFLabiV. The nuclear orphan receptor NR2F6 suppresses lymphocyte activation and T helper 17-dependent autoimmunity [J. Immunity. (2008) 29:205–16. doi: 10.1016/j.immuni.2008.06.008, PMID: 18701084WackATerczyńska-DylaEHartmannR. Guarding the frontiers: the biology of type III interferons [J. Nat Immunol. (2015) 16:802–9. doi: 10.1038/ni.3212, PMID: 26194286IvashkivLBDonlinLT. Regulation of type I interferon responses [J. Nat Rev Immunol. (2014) 14:36–49. doi: 10.1038/nri3581, PMID: 24362405FengDStoneRCElorantaMLSangster-GuityNNordmarkGSigurdssonS. Genetic variants and disease-associated factors contribute to enhanced interferon regulatory factor 5 expression in blood cells of patients with systemic lupus erythematosus [J. Arthritis Rheumatol. (2010) 62:562–73. doi: 10.1002/art.27223, PMID: 20112383SongSDeSNelsonVChopraSLaPanMKamptaK. Inhibition of IRF5 hyperactivation protects from lupus onset and severity [J. J Clin Invest. (2020) 130:6700–17. doi: 10.1172/JCI120288, PMID: 32897883TaoHLinHSunZPeiFZhangJChenS. Klf4 Promotes Dentinogenesis and Odontoblastic Differentiation via Modulation of TGF-β Signaling Pathway and Interaction With Histone Acetylation [J. J Bone Miner Res. (2019) 34:1502–16. doi: 10.1002/jbmr.3716, PMID: 31112333ArbuckleMRMcClainMTRubertoneMVScofieldRHDennisGJJamesJA. Development of autoantibodies before the clinical onset of systemic lupus erythematosus. [J]. N Engl J Med. (2003) 349:1526–33. doi: 10.1056/NEJMoa021933, PMID: 14561795OdendahlMJacobiAHansenAFeistEHiepeFBurmesterGR. Disturbed peripheral B lymphocyte homeostasis in systemic lupus erythematosus [J. J Immunol. (2000) 165:5970–9. doi: 10.4049/jimmunol.165.10.5970, PMID: 11067960DoorenspleetMEKlarenbeekPLde HairMJvan SchaikBDEsveldtREvan KampenAH. Rheumatoid arthritis synovial tissue harbours dominant B-cell and plasma-cell clones associated with autoreactivity [J. Ann Rheum Dis. (2014) 73:756–62. doi: 10.1136/annrheumdis-2012-202861, PMID: 23606709Brito-ZerónPKostovBFraileGCaravia-DuránDMaureBRascónFJ. Characterization and risk estimate of cancer in patients with primary Sjögren syndrome [J. J Hematol Oncol. (2017) 10:90. doi: 10.1186/s13045-017-0464-5, PMID: 28416003XueXFuLQiuTCaoZWangXRaoW. Losing CD45 and various B-cell markers in a case of MYC-driven pediatric high-grade B-cell lymphoma, not otherwise specified that transformed from Burkitt's lymphoma during rituximab-containing treatments: a case report. [J]. Virchows Arch. (2023) 483:111–6. doi: 10.1007/s00428-022-03433-1, PMID: 36383247