Jialing Xiao1†
Xiaojuan Guo2†Lingxue Gong3†Qinhui Li4Kai Dong1Xiaoxin Guo1Huan Li1Ting Wang1Liang Wang1Weijia Wu1Chengzi Gan1Mingze Li1
Bo Gong1Yixuan Jiang5Yixiao Wu5Yifan Hu5Liang Han1*Jianing Yang2*
Yujie Mao2*- 1Department of Human Disease Genes Key Laboratory of Sichuan Province and Institute of Laboratory Medicine, Sichuan Academy of Medical Sciences and Sichuan Provincial People’s Hospital, School of Medicine, University of Electronic Science and Technology of China, Chengdu, Sichuan, China
- 2Department of Dermatology, Sichuan Provincial People’s Hospital, School of Medicine, University of Electronic Science and Technology of China, Chengdu, China
- 3College of Clinical Medicine, Guizhou Medical University, Guiyang, China
- 4Medical Services Department, Sichuan Provincial People’s Hospital, School of Medicine, University of Electronic Science and Technology of China, Chengdu, Sichuan, China
- 5University of Electronic Science and Technology of China, School of Medicine, Chengdu, Sichuan, China
Background: Non-segmental vitiligo (NSV) is an autoimmune disorder characterized by irregular depigmented skin patches due to melanocyte loss, which causes considerable psychosocial burden. Although localized mechanisms underlying vitiligo pathogenesis have been studied extensively, investigations into peripheral blood mononuclear cells (PBMCs), key mediators of autoimmune diseases, remain limited.
Methods: To address this gap, we performed single-cell RNA sequencing (scRNA-seq) on peripheral blood samples from 3 untreated patients with generalized, progressive non-segmental vitiligo (GP-NSV) and 3 healthy controls. Findings were validated using flow cytometry in an additional cohort of 7 GP-NSV patients and 30 controls. Computational analyses, including pseudotime trajectory reconstruction and pathway enrichment, were employed to characterize immune cell subsets and their functional states.
Results: Vitiligo patients exhibited striking heterogeneity in PBMC subsets. KLRC2+ NK cells were markedly reduced and enriched in tumor necrosis factor (TNF) and apoptotic signaling pathways, a finding further confirmed by flow cytometry. Pseudotime analysis indicated that NK cells underwent negative regulation of DNA metabolic processes alongside activation of granzyme-mediated programmed cell death. In addition, the frequency of FCGR3A+ Cytotoxic CD8+T cell was reduced, with enrichment in T cell activation and differentiation signatures. STAM+ regulatory T cells (Tregs) were increased, whereas EGR1+ B cells were decreased, both subsets showing enrichment in pathways linked to osteoclast differentiation and calcium ion metabolism, suggesting a potential role of calcium homeostasis dysregulation in disease pathogenesis.
Conclusions: This study provides the single-cell atlas of PBMCs in GP-NSV, uncovering profound transcriptional and compositional alterations across multiple immune cell subsets in active vitiligo. These findings offer novel insights into systemic immune dysregulation in GP-NSV and pave the way for novel targeted therapeutic strategies.
1 Introduction
Vitiligo is an autoimmune skin disease characterized by patches of white skin that can develop anywhere on the body and is caused mainly by CD8+ T cells attacking the epidermal melanocytes (1). On the basis of the rate of lesion expansion and its distribution characteristics, vitiligo can be categorized into two types: segmental vitiligo is characterized by unilateral and segmental lesion distribution, early onset, and rapid stabilization, whereas Non-segmental vitiligo (NSV) typically presents as bilateral, acral, or generalized symmetric lesions that evolve gradually over time (2). Notably, non-segmental vitiligo accounts for 84–95% of all vitiligo cases (3). It affects approximately 1% of the world’s population and have a huge psychosocial impact on people, such as shame, depression and anxiety (4) and are often associated with various systemic diseases including thyroid disease, rheumatoid arthritis, and alopecia areata (5). These research findings suggest that vitiligo may be a systemic immune disease. PBMCs serve as central sentinels of systemic immunity, acting both as frontline effectors of innate immunity and as bridges to adaptive immune responses. Through cytokine secretion, antigen presentation, and cytotoxic activity, they coordinate systemic immune reactions and can serve as sensitive indicators of disease progression and immune dysregulation in autoimmune conditions such as vitiligo. Nevertheless, their roles in vitiligo remain poorly understood.
Both lesional and peripheral immune dysregulation have been implicated in vitiligo. Cytotoxic CD8+ T cells, particularly those recognizing melanocyte-specific antigens, are considered key effectors mediating melanocyte destruction. These cells produce interferon-γ (IFN-γ), which drives chemokine production (CXCL9 and CXCL10) and recruits additional autoreactive T cells to the skin (6). NK cells and Innate Lymphoid Cell type 1(ILC1) also contribute to melanocyte injury through the IFN-γ–CXCR3B axis (7). Studies of peripheral blood using conventional immunophenotyping have revealed that vitiligo patients exhibit altered CD4+/CD8+ T cell ratios, expansion of activated T cell subsets, and changes in NK cell frequencies, reflecting systemic immune dysregulation associated with the disease (8). In addition, multiple studies have reported a reduction in both the frequency and suppressive function of Tregs in vitiligo patients, accompanied by decreased expression of FOXP3, IL−10, and TGF−β in the skin and peripheral blood, particularly during active disease, further underscoring the presence of systemic immune dysregulation (9). Emerging evidence suggests that peripheral immune subsets like B cells, γδ T cells, and tissue-resident-like T cells may play a role in immunity, but their contribution to vitiligo remains largely unexplored.
ScRNA-seq has revolutionized the ability to explore cellular heterogeneity at the individual cell level and uncover previously unrecognized immune cell subsets. It allows unbiased characterization of cell type composition and functional states, at single-cell resolution. Most scRNA-seq studies in vitiligo to date have focused on skin lesions, revealing tissue-resident memory T cells, melanocyte stress responses, and inflammatory signaling networks. For examples, regulatory T cells exert regulatory control over the development of vitiligo by modulating the interaction between the CCL5-CCR5 signaling pathways and CD8+ T cells, resulting in a trendic cycle (10). Simultaneously, Xu et al. (11) highlighted the distinctive role of epidermal fibroblasts in this process, as they secrete chemotactic factors to recruit and activate CD8+ cytotoxic T cells. Despite the current understanding that some cells play a role in the recruitment or activation of CD8+ cells, there remains a need for a more in-depth exploration of the functional characterization of other cells.
However, comprehensive profiling of the peripheral immune landscape in vitiligo remains limited, and the relationship between systemic immune states and cutaneous pathology is poorly defined.
In this study, PBMCs were isolated from three GP-NSV patients and three healthy controls for scRNA-seq to construct a high-resolution map of immune cell populations. We identified distinct immune cell subsets showing significant alterations in the patient group, including KLRC2+ NK cell and FCGR3A+ CD8+ Tem cells, explored their transcriptional programs and functional characteristics, and partially validated in subset by flow cytometry. Our findings provide new insights into systemic immune dysregulation and inform the development of targeted immunomodulatory therapies in vitiligo.
2 Methods and Materials
2.1 Ethics statement
This study was conducted in accordance with the Declaration of Helsinki. Ethical approval was obtained from the Institutional Review Board of Sichuan Provincial People’s Hospital (No. 2022731). Written informed consent was obtained from all participants prior to their inclusion in the study. The main research subjects were human peripheral blood samples. All personal data were anonymized to ensure participant confidentiality.
2.2 Subjects and clinical sample collection
In order to obtain data that better reflect the systemic immune status of vitiligo, we recruited patients with GP-NSV (12), who had depigmented areas covering more than 5% of body surface area and had the appearance of new lesions or extension of existing lesions within the last 3 months. Key exclusion criteria were previous use of any systemic treatment within 6 months. Healthy controls were recruited based on the absence of personal history of vitiligo, psoriasis, or any other autoimmune diseases. Additionally, all controls were in good general health with no significant comorbidities, including no history of systemic diseases or malignancy. The demographic characteristics of patients and healthy controls are detailed in Supplementary Table 1.
For single-cell sequencing, peripheral blood samples (10 mL/donor) were collected from three GP-NSV patients and three healthy controls using EDTA anticoagulant tubes (BD Vacutainer). For PBMC isolation, whole blood was diluted 1:1 with phosphate-buffered saline (PBS, Biosharp, Cat. No. BL302A) then layered over lymphocyte separation solution (Ficoll-Paque PREMIUM, Cat. No. 17544203) in 50 mL centrifuge tubes. Centrifugation was performed at 2000 rpm for 20 minutes at 20°C with brake disabled in a horizontal rotor. The intermediate white layer (peripheral blood mononuclear cells, PBMCs) was carefully aspirated and washed twice with PBS (300 × g, 10min), followed by resuspension in RPMI 1640 medium containing 0.04% BSA (Corning, Cat. No. 10-040-CVR). Cell concentration and viability were assessed using a Luna cell counter or Trypan blue staining.
For flow cytometry validation, peripheral blood samples from a larger cohort of 30 healthy controls and 7 GP-NSV patients were collected. Each sample consisted of 5 mL of peripheral blood. To maintain sample integrity, the collected blood samples were kept at room temperature and analyzed within 24 hours.
2.3 10x Genomics single-cell sample processing and cDNA library preparation
The freshly prepared single-cell suspension was adjusted to a concentration of 700–1200 cells/μl. According to the 10× Genomics Chromium Next GEM Single Cell 3′ Reagent Kit v3.1 User Guide (Cat. No. 1000268), the suspension was processed for loading and library construction. The constructed libraries were subjected to high-throughput sequencing on the Illumina NovaSeq 6000 platform using the PE150 platform. The raw data generated from high-throughput sequencing were in fastq format. These data were processed using the official 10x Genomics software, CellRanger, for quality statistics and alignment to the reference genome. The software quantifies high-throughput single-cell transcriptome data by recognizing cell-specific barcodes and unique molecular identifiers (UMIs) for different mRNA molecules within each cell. Data on cell viability, capture efficiency, and sequencing depth for each sample are provided in Supplementary Table S1.
2.4 Single-cell data preprocessing, gene expression quantification and cell-type determination
Based on the initial quality control from CellRanger, further quality control was performed using the Scrublet software package (13). Based on the UMI count matrix obtained from CellRanger, doublets were first identified and excluded using the Scrublet software. The matrix was then imported into Seurat (version 3.2.3) for subsequent quality control. The following filtering criteria were applied: (1) cells with a mitochondrial gene ratio exceeding 20%; (2) cells with fewer than 500 detected genes or a total UMI count below 1,000; and (3) after the above steps, cells with total UMI counts deviating beyond two standard deviations (mean ± 2 SD) from the overall mean were further removed to eliminate extreme outliers. To correct for batch effects across samples, canonical correlation analysis (CCA) in Seurat was employed for data integration. Prior to integration, the expression matrix of each sample was normalized using the NormalizeData function with the parameter normalization.method = “LogNormalize” and scale.factor = 10000. The top 2,000 highly variable genes were selected via the FindVariableFeatures function for downstream analysis. Subsequently, based on the first 10 principal components, cell clustering was performed using the FindNeighbors and FindClusters functions with a resolution of 0.8. Dimensionality reduction and visualization were achieved using Uniform Manifold Approximation and Projection (UMAP). Finally, the FindAllMarkers function was utilized to identify significantly upregulated marker genes in each cell cluster, and automated cell type annotation was performed using the SingleR package with reference datasets to minimize subjective bias. Differential expression analysis between groups was performed using the FindMarkers function in Seurat, with the parameters set to min.pct = 0.1 and test.use = “wilcox”. The identified marker genes were visualized using the VlnPlot and FeaturePlot functions. Using the SingleR package (14), cell type identification was performed on the basis of a single-cell reference expression quantitative public dataset. The expression profile of the cells to be identified was correlated with the reference dataset, and the cell type with the highest correlation in the reference dataset was assigned to the cells being identified. This approach reduces human subjective interference to some extent. The identification principle involves calculating the Spearman correlation between the expression profile of each cell in the sample and each annotated cell in the reference dataset. The cell type with the highest correlation to the sample cell’s expression in the reference dataset is selected as the final identified cell type.
2.5 Pseudotemporal trajectory analysis
Pseudotemporal trajectory analysis was performed on the target cell subpopulations using Monocle2 (15). Briefly, the UMI expression matrix and corresponding metadata from the Seurat object were converted into a CellDataSet object, with the expression family specified as the negative binomial distribution (negbinomial.size()). Size factors and dispersion values were then estimated using the estimateSizeFactors and estimateDispersions functions, respectively. Highly variable genes across cells were selected via the dispersionTable function for downstream ordering analysis. The data were subsequently reduced to two dimensions using the reduceDimension function, and cells were ordered along the inferred trajectory with the orderCells function. Finally, trajectory topology, dynamic gene expression patterns, and pseudotemporal clusters were visualized using the plot_cell_trajectory, plot_genes_in_pseudotime, and plot_pseudotime_heatmap functions, respectively.
2.6 KEGG and GO pathway analysis
Gene Set Preparation: Background gene set files were downloaded and organized from the KEGG database (https://www.kegg.jp/) and the GO database using the GSEABase package (v1.44.0). Pathway Activity Scoring: The GSVA package (v1.30.0) was used to calculate pathway activity scores for single cells. GSVA (Gene Set Variation Analysis) is an unsupervised gene set enrichment analysis method that calculates gene set enrichment scores for individual samples on the basis of gene set variation analysis. Differential Pathway Analysis: The LIMMA (16) package (v3.38.3) was used to calculate differences in pathway activity between different groups. LIMMA (Linear Models for Microarray Data) is a linear model method used to analyze differential expression in high-throughput gene expression data.
2.7 RNA extraction and library construction
Total RNA was extracted using the TRIzol reagent according to the manufacturer’s instructions. RNA purity and quantification were assessed using a NanoDrop 2000 spectrophotometer (Thermo Scientific, USA), and RNA integrity was evaluated with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Transcriptome libraries were constructed using the VAHTS Universal V6 RNA-seq Library Prep Kit following the manufacturer’s protocol.
2.8 RNA sequencing and differential gene expression analysis
Library sequencing was performed on the Illumina NovaSeq 6000 platform, generating 150 bp paired-end reads. Each sample yielded approximately 50 million raw reads. Raw reads in FASTQ format were processed using fastp (17) to remove low-quality reads, obtaining clean reads for subsequent analysis. Read alignment to the reference genome was conducted using HISAT2 (18), followed by gene expression quantification in FPKM (19). The read counts for each gene were obtained using HTSeq-count (20).Principal component analysis (PCA) and visualization were performed using R (v3.2.0) to evaluate the biological reproducibility of the samples.
Differential gene expression analysis was conducted using DESeq2 (21), where genes with a q-value<0.05 and a fold change > 2 or<0.5 were defined as differentially expressed genes (DEGs). Hierarchical clustering analysis of DEGs was performed using R (v3.2.0) to illustrate gene expression patterns across different groups and samples. Additionally, radar plots were generated for the top 30 DEGs using the R package ggradar to visualize the expression changes of upregulated and downregulated genes.
2.9 Flow cytometry validation of KLRC2+ subset
Fresh heparin-anticoagulated whole blood samples were inverted 10 times at room temperature to mix thoroughly. Then, 5 μL of antibody and 25 μL aliquots of blood were added to a flow cytometry tube, followed by the addition of 1 μL KLCR2 antibody (BD, Cat# 748169) and 4 μL of a multi-color peripheral blood flow cytometry detection antibody kit (CD3/CD16 + 56/CD45/CD19, labeled with FITC/PE/PerCP/APC; Mindray, REF#105-004448-00). After vortex mixing, the samples were incubated in the dark for 15 minutes. Subsequently, 450 μL of red blood cell lysis buffer was added, vortexed, and incubated at room temperature for 5 minutes. Then, 2 mL of physiological saline was added, vortexed, and the samples were centrifuged at 450 g for 5 minutes. The supernatant was discarded, and 500 μL of physiological saline was added and stored. Finally, the samples were analyzed using the BD FACSCanto™ II flow cytometer.
2.10 Statistical analysis
FlowJo (version 10) was used to gate different cell subsets expressing KLRC2 in PBMCs, and the percentages of positive cells within each subset were exported. The data were then imported into SPSS (version 27) for statistical analysis. Based on the data distribution, the independent-samples Mann-Whitney U test was employed for comparative analysis, with a two-tailed P-value<0.05 considered statistically significant. GraphPad Prism (version 8.0.2) was used to generate bar graphs.
3 Results
3.1 scRNA-seq reveals five cell types in PMBC
The workflow for the single-cell sequencing study in this section encompasses clinical sample collection, isolation and preparation of PBMCs, library construction, sequencing, and subsequent data visualization and interpretation of biological significance (Figure 1A). We initially isolated 59,700 cells from PBMCs obtained from six samples for single-cell sequencing. These samples comprised three cases of patients with generalized, progressive non-segmental vitiligo (28,799 cells) and three healthy control cases (30,921 cells). Following cell collection, rigorous quality control procedures were applied, including the removal of doublets, empty droplets, and low-quality cells, before proceeding to subsequent analysis. Unsupervised clustering using Seurat revealed a total of 19 clusters (Figure 1B). The cell groupings for each subject are shown in (Supplementary Figure 1A). Cell cycle analysis revealed that the cell cluster in the lower-right corner was predominantly in the G1 phase (Figure 1C). To impart biological significance to these subpopulations, we first identified four major cell types on the basis of classical cell markers: T cells (CD3E, CD3D, TRAC), B cells (MS4A1, CD79A, and CD9B), NK cells (NKG7, GZMB, and IFNG), and monocytes (LYZ, S100AB, and S100A9) (Figures 1D, E). Notably, Clusters 6 and 13 exhibited expression profiles encompassing both T cell markers (CD3D, CD3E, and TRAC) and NK cell markers (NKG7 and GZMB) (Figure 1E). Considering this dual expression pattern, we classified these two cell subpopulations as NKT cells. In the initial identification of 19 cellular subtypes, we observed a decrease in Cluster 11 at the NK cell (P<0.05) and an increase in Cluster 13 at the NKT cell (P<0.05) within the vitiligo group (Figure 1F), suggesting a potential association with the onset and progression of vitiligo.
Figure 1. Single-cell transcriptomic analysis of peripheral blood mononuclear cells (PBMCs) from healthy controls and non-segmental vitiligo patients. (A) Schematic of the single-cell transcriptome sequencing workflow for the healthy control and patient groups. (B) UMAP plot showing the distribution of all cells, divided into 19 subclusters by clustering. After quality control, a total of 59,700 cells were obtained. (C) Cell cycle analysis of the total cell population. (D) Identification of cell types in the UMAP plot using Seurat and cell markers, including monocytes, T cells, B cells, and NK cells. (E) In each UMAP plot, classical marker genes were used to identify cell types. (F) Comparison of cell frequencies across different cell subclusters between the healthy control and patient groups. Data are presented as mean ± SEM, and paired t-test was used. *P <0.05.
3.2 Heterogeneity in NK and NKT cells in vitiligo patients
To delve deeper into the roles of NK and NKT cellular subtypes in the occurrence and development of vitiligo, we conducted differential gene expression analysis. The NK and NKT cell groupings for each subject are shown in (Supplementary Figure 1B). The results revealed that NK Cluster exhibited elevated expression of genes such as MYOM2, IFITM1, CX3CR1, GIMAP7, XIST, TXNIP, and SLFN5 (log2FC>1; P<0.05), whereas genes such as EGR1, KDM6B, DUSP2, IFNG, HLA-B, BTG2, FOSB, DUSP1, NR4A2, ZFP36, CXCR4, MYADM, JUND, FOS, JUN, TNFAIP3, CD69, and NFKBIA were downregulated (log2FC<-1; P<0.05) (Figure 2A). Further KEGG and GO analyses revealed enrichment in apoptotic signaling pathways in Cluster 11, explaining the observed reduction in NK cell numbers in vitiligo patients compared with those in the control group. Additionally, this cellular subtype was enriched in inflammatory pathways such as TNF signaling, IL-17 signaling, MAPK signaling, antigen processing and presentation, and Th17 cell differentiation (specific pathways to be investigated). GO analysis further revealed significant enrichment in T cell activation for this cellular subtype, suggesting a potential relevance to the occurrence of vitiligo (Figures 2B, C).
Figure 2. Gene expression and functional enrichment differences analysis of primary cell subsets with spectral differences in vitiligo and health control. (A) and (D) represent volcano plots of differentially expressed genes in Cluster 11 and Cluster 13 between normal controls and vitiligo patients. (B) and (E) show the KEGG functional enrichment analysis for Cluster 11 and Cluster 13, with the x-axis representing gene counts and the y-axis representing the relative functional categories of gene enrichment. (C) and (F) represent GO analysis for Cluster 11 and Cluster 13, with the red box indicating biological processes (BP), blue representing cellular components (CC), and green representing molecular functions (MF).
In contrast, NKT Cluster exhibited a significant increase in cell number compared with the normal control group. Differential gene expression analysis identified upregulated genes, including B2M, IFITM1, MT-ATP8, MTCO3P12, and MYOM2 (log2FC>1; P<0.05), whereas genes such as HLA-B, NFKBIA, TNFAIP3, CXCR4, FOS, DUSP2, AC020916.1, LINC-PINT, ZFP36, RGCC, JUND, CD69, NR4A2, and MYADM were downregulated (log2FC<-1; P<0.05) (Figure 2D). Through combined KEGG and GO analyses, we found that Cluster 13 was associated with biological processes such as the IL-17 signaling pathway, TNF signaling pathway, B cell receptor signaling pathway, T cell receptor signaling pathway, and hematopoietic regulation (Figures 2E, F). This may explain the observed increase in cell numbers and suggest a regulatory role in the immune system.
3.3 The reduction of the NK cell subset containing KLRC2
To identify the specific NK and NKT cell subpopulations potentially involved in vitiligo, we extracted 9,257 cells from all NT and NKT cluster and re-clustered them into 15 distinct subtypes (Figure 3A). Notably, we observed a decrease in the prevalence of several cellular subtypes, namely Clusters 3, 6, and 8, in vitiligo patients (Figure 3B). We analyzed the expressions of CD3D, CD3E, CD7, and GZMB to determine whether these cells exhibit characteristics of NK cells and NKT cells (Figure 3C). In the quest for common markers among these cellular subtypes, we identified the gene marker KRLC2 through violin plots (Figures 3D, E), which was expressed consistently across Clusters 3, 6, and 8. Although the specificity of this cellular marker is not high, it currently serves as a relatively robust marker for labeling these subtypes. Consequently, we defined this cellular subtype as KLRC2+ NK cells. This cell groupings for each subject are shown in (Supplementary Figure 2A). In this clusters, certain genes, such as MYOM2, XIST, CX3CR1, and IFITM1, were totally significantly upregulated (log2FC>1; P<0.05), whereas others, including BTG2, FOS, CXCR4, CD69, TNFAIP3, and NFKBIA, were significantly downregulated (log2FC<-1; P<0.05). Further KEGG and GO analyses revealed significant enrichment of T cell activation, IL-17 signaling, apoptosis, and oxidative stress response pathways in Clusters 3 and 6 (Figures 3F, G).
Figure 3. Transcriptional profiling of NK cell subpopulations in non-segmental vitiligo and healthy controls. (A) UMAP plot showing the reclustering of NK cell subpopulations into 19 distinct clusters, with each cluster color-coded and numbered. (B) Cell frequency analysis of NK cell subpopulations in healthy controls and non-segmental vitiligo patients, highlighting differences in clusters 3, 6, and 8 between the groups. (C) Further classification of NK cell subpopulations using classical cell markers to identify NK cells. (D) Violin plots identifying specific differentially expressed genes within NK cell subpopulations. (E) UMAP plot illustrating the specific expression of KLRC2 in clusters 3, 6, and 8. (F) KEGG pathway enrichment analysis of KLRC2+ NK cells. The x-axis represents the enrichment score, and the y-axis lists the pathways. Pathways are ordered by significance, with the most significantly enriched pathways at the top. (G) GO enrichment analysis of DEGs in KLRC2+ NK cells, highlighting key categories: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). The x-axis represents enrichment scores or gene counts, and the y-axis lists the GO terms. Terms are ranked by significance, with the most significant terms at the top. Data are presented as mean ± SEM, and paired t-test was used. *P<0.05.
To further validate the reduction of KLRC2+ NK cell subsets, we collected samples from 7 vitiligo patients of different age groups and 30 healthy controls (Supplementary Table 2). The frequency of KLRC2+ NK cells was assessed using flow cytometry, with the results shown in (Supplementary Figures 3A, B). Interestingly, flow cytometry data revealed a significant reduction in KLRC2 expression on the surface of T cells, B cells, and monocytes in the patient cohort, whereas the frequency of KLRC2+ NK cells showed no significant difference compared with healthy controls. Consistent with this, we performed bulk RNA sequencing on samples corresponding to those used in single-cell sequencing, and the results showed a significant decrease in KLRC2 expression (Supplementary Figures 3C, D). Consistent with vitiligo, a reduction in the CD56+ NK cell subpopulation was also observed in other autoimmune diseases, including hyperthyroidism and alopecia areata (22, 23).
3.4 NK cell subsets monocle analysis to delve its lifespan progression in vitiligo
The observed decrease in inhibitory receptors (CD158a+) and increase in activating receptors on NK cells have been implicated in the pathogenesis of vitiligo, suggesting a potential role for NK cell activation in disease development and progression (24). Similarly, we found that the NK cell subset KLRC2+ is reduced in non-segmental vitiligo patients, suggesting that different NK cell subsets may play distinct roles in disease development. Therefore, we used pseudotime analysis to investigate the branching of NK cell subsets to explore the impact of each subset on vitiligo. Pseudotime scores and sample distributions were used to determine the direction of NK cell differentiation (Figures 4A–D). Furthermore, we identified pseudotime-dependent genes across four distinct modules (modules 1–4), each characterized by unique expression patterns (Figure 4E). In the late stage of disease in NK cells, Module 1 and Module 3 exhibited high expression of genes such as GZMA, DYNLL1, and ATM, LCP1 which are associated with the negative regulation of DNA metabolic processes and the granzyme-mediated programmed cell death signaling pathway. In contrast, in the early stage of disease in NK cells, Module 2 and Module 4 showed high expression of genes such as HSPA5, RGCC, BIRC3, HSP90AA1, NFKBIZ, ZFP36L2, KLF6, and JMJD6, which are associated with programmed cell death and the positive regulation of macromolecule metabolic processes.
Figure 4. Pseudotime analysis reveals the developmental trajectories of NKT cell subpopulations. (A) Pseudotime trajectory plot of different NKT cell subpopulations showing developmental paths. Distinct branches represent different trajectories of cell development. (B) Pseudotime trajectory plot based on individual samples from healthy controls and non-segmental vitiligo patients, highlighting developmental paths and branching patterns. (C) Trajectory plot showing NKT cell developmental changes across different pseudo-time states. Various branches indicate different developmental trajectories. (D) Pseudo-time trajectory plot based on Pseudotime scoring, illustrating cell development along different paths. (E) Heatmap displaying the expression patterns of significantly varying genes in NKT cells on the basis of Pseudotime assessment, with emphasis on gene expression trends across different pseudo-time points.
Taken together, these findings indicate that NK cell functions undergo distinct shifts during the course of vitiligo, transitioning from early-stage profiles in NKT cells enriched for genes linked to programmed cell death and macromolecule metabolism, to late-stage NK cell profiles characterized by genes associated with the suppression of DNA metabolic processes and granzyme-mediated cytotoxicity, suggesting a potential role for these functional changes in disease progression.
3.5 PBMC scRNA-seq reveals functional and phenotypic changes in Treg and FCGR3A+ CD8+ T cells in vitiligo patients
T cell subsets play critical roles in the pathogenesis and progression of vitiligo. For example, CD8+ T cells primarily mediate melanocyte death (25). However, the functions of other T cell subsets remain unclear. Therefore, we reclustered the T cell subsets and identified 13 new subsets (Figure 5A). The T cell groupings for each subject are shown in (Supplementary Figure 1C). Using the SingleR package and classical cell markers (Figure 5B), we manually defined CD4+ T cell, Naïve CD4+ T, Naïve CD8+T, Cytotoxic CD8+T, MAIT T cell and Treg subsets (Figure 5C) (26). First, we conducted cell spectrum analysis on various subpopulations and found that the number of cells in Cluster 10 (Figure 5D), which contains Tregs and exclusively expresses FOXP3, had increased (P = 0.012) (Figure 5F). This finding is consistent with previous reports (27). We used a gene expression heatmap to identify some uniquely expressed marker genes in Treg cells, including RTKN2, IKZF2, TIGIT, IL2RA, CTLA4, and STAM. These genes may serve as markers for identifying Treg cells in vitiligo. The specific expression of STAM in Treg cells allowed us to define a novel STAM+ Treg subset in vitiligo, which may contribute to disease development. The STAM+Treg cell groupings for each subject are shown in (Supplementary Figure 2C). By examining the functions of these genes in Treg cells, we found that they are involved in multiple biological processes and pathways, including apoptosis, EBV infection, the TNF signaling pathway, and the IL-17 signaling pathway, as well as the negative regulation of protein phosphorylation, regulation of hematopoiesis, and modulation of the adaptive immune response, highlighting the immunoregulatory role of Treg cells (Figures 5I, J). Interestingly, Treg cells exhibit significant osteoclast differentiation activity, suggesting a potential link between vitiligo and osteoporosis (28). Similarly, by analyzing cell frequency differences and correlations, we found that the cell subpopulations in Cluster 5(P = 0.047) and Cluster 9(P = 0.034) were significantly reduced in patients (Figures 5D, E), and they identified as Cytotoxic CD8+T cells (Figure 5C). Violin plots indicated that FCGR3A may serve as a specific marker for this subpopulation (Figure 5F); accordingly, we defined FCGR3A+ Cytotoxic CD8+ T cells and observed a reduction of this subset in patients. This cell groupings for each subject are shown in (Supplementary Figure 2B). Functional enrichment analysis revealed that this cell subset is associated with biological processes related to mononuclear cell differentiation, T cell differentiation, and lymphocyte differentiation. These differentiation processes may contribute to the progression of vitiligo (Figures 5G, H).
Figure 5. Re-clustering transcriptomic analysis of T cell subpopulations in healthy controls and non-segmental vitiligo patients. (A) UMAP plot showing the re-clustering of T cell subpopulations into 13 new clusters. (B) Dot plot displaying the marker genes used to delineate major immune cell subsets. (C) UMAP visualization illustrating the identification of T cell subsets. (D) Bar chart depicting cell frequencies across different groups. (E) Heatmap illustrating the correlation between T cell subpopulations, with colors indicating the strength of the correlation from negative (blue) to positive (red). (F) Violin plot showing the expression of specific genes in different T cell subpopulations. (G) GO pathway enrichment analysis for FCGR3A+ Cytotoxic CD8 + T cells, with the x-axis representing enrichment scores and the y-axis listing pathways, ordered by significance. (H) KEGG analysis of differential gene expression in FCGR3A+ Cytotoxic CD8 + T cells, showing enrichment in biological processes (BP), cellular components (CC), and molecular functions (MF). The x-axis represents enrichment scores or gene counts, and the y-axis lists GO terms, ordered by significance. (I) GO pathway enrichment analysis for STAM+Treg cells, with the x-axis representing enrichment scores and the y-axis listing pathways, ordered by significance. (J) KEGG analysis of differential gene expression in STAM+Treg cells, showing enrichment in BP, CC, and MF. The x-axis represents enrichment scores or gene counts, and the y-axis lists GO terms, ordered by significance. Data are presented as mean ± SEM, and paired t-test was used. *P<0.05.
3.6 Single-cell RNA sequencing of PBMCs unravels altered functional states in the EGR1+ naïve B cell compartment in vitiligo
The role of B cells in the development of vitiligo is not yet clear. However, B cell activating factors may activate self-reactive B cells to produce auto-antibodies against melanocytes. These auto-antibodies may function as cellular adjuvants, activating CD4+ T cells and thereby enhancing their helper effect on the activation of CD8+ T cells (29). To investigate the role of B cells in vitiligo, we re-clustered B cells into 10 new subpopulations (Figure 6A) and used singleRs and musical marker (Figure 6B) to manually define them as memory B cells, naive B cells and plasma (Figure 6C). The B cell groupings for each subject are shown in (Supplementary Figure 1D). We calculated the frequencies of these cells and found that Clusters 6(P = 0.001) and 7(P = 0.028) were reduced in vitiligo patients (Figure 6D). Correlation analysis revealed that these subpopulations were highly correlated (Figure 6E). Violin plots revealed that Cluster 6 specifically expressed EGR1, which regulates CD44 transcription in B cells, promoting B cell homing and migration; thus, we defined these cells as EGR1+ B cells (Figure 6F). This cell groupings for each subject are shown in (Supplementary Figure 2D). Functional studies indicated that memory B cells might be related to calcium ion cellular responses (Figure 6G). This is related to H2O2-mediated oxidative stress in the epidermis, leading to altered calcium binding and a significant decrease in calcium ATPase activity (30). Additionally, GSVA analysis indicated that EGR1+ memory B cells are enriched in cytokine-cytokine receptor interactions and hematopoietic cell lineage pathways (Figure 6H).
Figure 6. Re-clustering transcriptomic analysis of B cell subpopulations in healthy controls and non-segmental vitiligo patients. (A) UMAP plot showing B cell re-clustering into 10 new subpopulations, with cells marked by color and number. (B) Dot plot displaying the marker genes used to delineate B cell subsets. (C) UMAP visualization illustrating the identification of B cell subsets. (D) Bar chart of B cell frequencies across groups. (E) Correlation heatmap of B cell subpopulations, with colors indicating the strength of the correlation from negative (blue) to positive (red). (F) Violin plot of specific gene expression in B cell subpopulations. (G) KEGG pathway enrichment analysis for naive B cells, with x-axis showing enrichment scores and y-axis listing pathways, ordered by significance. (H) Gene Set Enrichment Analysis (GSEA) of B cell subsets in vitiligo patients. The x-axis represents enriched pathways, and the y-axis represents cell subsets.
These findings suggest a functional role for these B cells in vitiligo; however, further research is needed to elucidate their exact mechanisms of action.
3.7 PBMC scRNA-seq identifies twelve monocyte subsets
Current research on the role of monocytes in the development of vitiligo is limited. Previous studies have indicated that monocytes accumulate a significant amount of hydrogen peroxide (31), leading to altered enzyme activity. These findings suggest that oxidative stress plays a crucial role in the pathogenesis of vitiligo. To further investigate the potential role of monocytes in the development and progression of vitiligo, we re-clustered the monocytes into 12 subpopulations (Figures 7A, B). Using classical cell markers and the SingleR method (14), we defined these subpopulations as monocytes, CD14+ monocytes, CD16+ monocytes (32–34), and myeloid dendritic cells (Figure 7C) (35). The monocyte cell groupings for each subject are shown in (Supplementary Figure 1E). We subsequently performed statistical analysis of the frequencies of the re-clustered cell subpopulations and found a significant difference between Cluster 1(P = 0.003) and Cluster 6 (P = 0.004) (Figure 7D). Correlation analysis further revealed a positive correlation between these two subpopulations (Figure 7E). We attempted to identify specific marker genes for these two subpopulations, but violin plot analysis did not reveal any significantly different marker genes within the clusters (Figure 7F). We then applied KEGG and GO analyses to explore their functional pathways (Figures 7G–J), which indicated that these two subpopulations are involved in pathways related to leukocyte cell-cell adhesion, the regulation of cell-cell adhesion, the regulation of leukocyte cell-cell adhesion, and the regulation of T cell activation. These pathways suggest that leukocyte cell-cell adhesion is a critical process involving multiple pathways and molecules, including integrins, selectins, ICAM, and VCAM, which play key roles in inflammation and immune system functions. Additionally, these processes are accompanied by T cell activation, indicating that monocytes may play a role in the pathogenesis of vitiligo. However, further studies are required to confirm this potential involvement.
Figure 7. Transcriptomic analysis of monocyte subpopulations re-clustering in healthy controls and non-segmental vitiligo patients. (A) UMAP plot illustrating the re-clustering of monocyte subpopulations into 12 distinct clusters, with individual cells denoted by color and numerical labels. (B) Comparison of monocyte frequencies across samples from healthy controls and non-segmental vitiligo patients. (C) UMAP plot showing the automatic annotation of monocyte subpopulations and other cell types using SingleR. (D) Bar chart representing the distribution of monocyte frequencies between the healthy controls and non-segmental vitiligo groups. (E) Heatmap depicting the correlations between monocyte subpopulations, with color intensity indicating the strength of the correlation, ranging from negative (blue) to positive (red). (F) Violin plot displaying the expression patterns of specific genes across different monocyte subpopulations. (G–J) KEGG and GO enrichment analyses of CD14+ monocytes and cluster 6 monocytes, highlighting the functional enrichment and differential gene expression within biological processes (BP), cellular components (CC), and molecular functions (MF). The x-axis represents enrichment scores or gene counts, and the y-axis lists pathway or GO term names, ranked by significance.
4 Discussion
Vitiligo is a depigmenting skin disorder caused by the loss of melanocytes, affecting approximately 1% of the global population. The NSV accounts for 84–95% of cases and imposes a substantial psychological burden on patients. Beyond its cutaneous manifestations, vitiligo is frequently accompanied by systemic autoimmune diseases, underscoring its value as a model for studying systemic immune dysregulation. Current models posit that disease spread is primarily driven by the IFN-γ to JAK/STAT axis, leading to elevated CXCL9/CXCL10 expression and subsequent recruitment of CXCR3+ cytotoxic CD8+ T cells to new skin sites (36). However, the “selectivity” underlying long-range lesion dissemination remains poorly understood. To address this gap, we performed single-cell RNA sequencing of PBMCs, a compartment reflecting integrated systemic immune activity, to gain a new immunological perspective on vitiligo pathogenesis. Our analysis revealed dynamic alterations and functional roles of KLRC2+ NK cells, Treg cells, FCGR3A+ Cytotoxic CD8 + T cells, and EGR1+ B cells, opening new avenues for systemic immunology-oriented research in vitiligo.
The balance of immune cells plays a crucial role in the progression of vitiligo (37–39). The inconsistency in reported NK cell abundance in vitiligo patients across different studies may be attributed to variations in disease stages and individual patient differences (24, 40). Cytokines derived from NK cells and their cytotoxic functions regulate immune responses and may contribute to various immune-mediated diseases (41), such as ankylosing spondylitis, Behçet’s disease, and multiple sclerosis.
Here, we report a reduction in NK cells accompanied by an increase in NKT cells in vitiligo patients. Notably, the diminished NK cell population was characterized by a distinctive marker profile—KLRC2—allowing us to define this subset as KLRC2+ NK cells. Functional enrichment analysis indicated that this subset is associated with apoptosis, T cell activation, TNF signaling, and lymphocyte differentiation. Consistent with our single-cell data, the enrichment of apoptosis-related genes aligns with the observed decrease in KLRC2+ NK cell abundance. Moreover, bulk RNA-seq analysis of PBMCs also revealed a reduction of KLRC2 expression, whereas the roles of this subset in other pathways remain to be further elucidated. Although flow cytometric validation of KLRC2 expression in bulk NK cell showed no statistical significance, this discrepancy may stem from two methodological considerations: (1) Structural homology between KLRC2 and KLRC1 (differing by only two amino acids) may compromise antibody specificity, potentially masking true expression differences; (2) Upon re-clustering the original differential spectrum cell subsets, Cluster 11 and Cluster 13, we further identified 15 clusters. Using canonical lineage markers (NK cells: CD7, GZMB; T cells: CD3D, CD3G), three intermediate clusters (Clusters 3, 6, 8) were identified at the NK-T cell interface, suggesting a transitional immunophenotype. Cross-compartment analysis demonstrated no significant KLRC2 modulation in sorted NK/NKT populations, but revealed marked downregulation in T cells from patients. This phenomenon may be due to Clusters 3, 6, and 8 being in a transitional state between NK and NKT cells, representing a T cell subset that was not distinctly classified at the single-cell level but became evident through flow cytometry analysis. This result aligns with single-cell sequencing data, further supporting the existence of a transitional T cell subset between NK and NKT cells. Moreover, the decreased KLRC2 expression in this T cell subset in patients suggests its potential as a diagnostic marker for vitiligo, highlighting its relevance in disease pathogenesis and biomarker development. Aberrant expression of KLRC2 in NK cells is critical for maintaining the microenvironment essential for tumor initiation and progression (42). Its absence also plays a pivotal role in immune regulation during transplantation (43)or infection (44). Upregulation of KLRC2 in vitiligo skin lesions has been reported, suggesting its potential involvement in local immune responses (45). However, the fate, function, and underlying mechanisms of the activated NK cell subset represented by KLRC2 in peripheral blood remain largely unexplored. Here, we demonstrate for the first time that KLRC2+ NK cells are markedly reduced in vitiligo patients and exhibit features enriched for apoptosis-related and immunoregulatory pathways. These findings not only fill a critical gap in the current literature but also provide new evidence and insights for understanding the pathogenesis of vitiligo from a systemic immunological perspective.
Disease progression is a dynamic process. To explore the transcriptional alterations in NK and NKT cells during this process, we investigated this issue using pseudotime analysis. In the early stage of the disease, NK cells exhibited high expression of genes such as GZMA (46), DYNLL1, ATM, and LCP1, which are involved in the negative regulation of DNA metabolism and granzyme-mediated programmed cell death. These findings are consistent with existing studies and suggest that NK cells in early-stage vitiligo may eliminate aberrant melanocytes through granzyme-mediated cytotoxicity while simultaneously limiting their own DNA metabolic activity to maintain cellular homeostasis. In the late stage of the disease, NK cells showed elevated expression of HSPA5, RGCC, BIRC3, HSP90AA1, NFKBIZ, ZFP36L2, KLF6, and JMJD6, genes associated with programmed cell death and the positive regulation of macromolecule metabolism. The expression patterns of these genes align with previously reported NK cell functional characteristics, indicating that in late-stage vitiligo, NK cells may adapt to chronic inflammatory conditions or sustained tissue damage by modulating metabolic and stress-response pathways. Future studies should focus on elucidating the activation mechanisms of NK cells, the molecular basis underlying their functional transitions, and their specific roles in disease progression, with the aim of informing strategies for early diagnosis and targeted therapy in vitiligo.
In vitiligo, the infiltration of CD8+ T cells is mediated primarily by the secretion of chemokines such as IFN-γ by dermal cells (11, 47, 48). These chemokines recruit and activate CD8+ T cells, which then infiltrate the dermal tissue and target melanocytes. Treg cells, through the CCL5-CCR5 axis, exert a protective effect on melanocytes by modulating the activity of CD8+ T cells (10). Existing evidence indicates that a subset of tissue-resident memory T (TRM) cells patrols the epidermis and papillary dermis, and are reactivated by homologous antigens, leading to relapse (49). Although the localized pathogenesis and recurrence of vitiligo have been extensively studied, there is still a lack of research on immune cell subsets in the context of vitiligo spread (50). Investigating the characteristics of immune cells in peripheral blood of progressive patients may be a crucial step in addressing the issue of vitiligo expansion. Here, we identified 13 cell subsets and found that, compared with those in healthy individuals, Treg cells in progressive vitiligo patients are more abundant and express genes such as RTKN2, IKZF2, IL2RA, and STAM. STAM is a cytosolic signal-transducing adaptor molecule that regulates JAK-STAT signaling through its role in cytokine receptor transport and serves as a critical negative regulator of the cGAS-STING pathway, thereby restraining hyperactivation of the innate immune response (51, 52).In addition to STAM, these genes are critical in autoimmune diseases and Treg activation (53–55). Notably, the STAM gene has not been previously reported in Treg cells, suggesting that it may play a significant role in Treg regulation in vitiligo. There, we defined this cell subset as the STAM+ Treg cell subset.
During the pathogenesis of vitiligo, regulatory T cells (Tregs) recognize the CCL5 ligand secreted by CD8+ T cells through their own surface CCR5 receptor, thereby exerting immunosuppressive functions and restraining the overactivation of CD8+ T cells (10). However, despite the increased abundance of Tregs in both the epidermis and peripheral blood of patients, they fail to effectively curb disease progression, suggesting that Tregs may reside in a “functionally impaired” state (56, 57), characterized by Th1-Treg polarization (58). Our single-cell transcriptomic analysis revealed a significant upregulation of the transcription factor STAT1 in patient-derived Tregs (56, 57), along with activation of the tumor necrosis factor (TNF)-related signaling pathway, further corroborating the shift of Tregs toward a pro-inflammatory Th1 phenotype (59). Notably, STAM, functioning as an intracellular signal-transducing adaptor molecule, may facilitate STAT1 activation via modulation of the JAK-STAT pathway. Collectively, these alterations lead to diminished immunosuppressive capacity and enhanced pro-inflammatory properties of Tregs, thereby impairing their regulatory control over CD8+ effector T cells and ultimately exacerbating the initiation and progression of vitiligo.
The secretion of cytokines may promote osteoporosis and is associated with skin diseases, suggesting an immune-mediated connection between them (28, 60). The role of CD8+ T cells as the primary cell population that targets melanocytes has not yet been fully elucidated in peripheral blood. Our data revealed a significant reduction in the FCGR3A+ Cytotoxic CD8 + T cell subset in patients. Functional enrichment analysis also indicated that this cell subset is associated with apoptosis-related signaling pathways, corroborating the reduction in cell abundance observed in patients. Additionally, the cells expressed pathways related to T cell differentiation and activation, suggesting an abnormal activation state of immune cells. FCGR3A, which is commonly associated with NK cells, mediates antibody-dependent cellular cytotoxicity (ADCC), enabling T cells to perform ADCC functions (61, 62). But, for the CD8+ T cell subset expressing FCGR3A, previous study (63)has shown that it is a terminally differentiated cell with the innate-like characteristics of CD8+ T cells in the SARS-CoV-2 infected. In the smoking population, the number of FCGR3A+ Cytotoxic CD8 + T cell subset was significantly elevated, suggesting immune activation and implying a potential association between smoking and the development of vitiligo (64). However, in our patients, this subset was markedly reduced, which may indicate immune dysregulation in individuals with vitiligo.
Upon stimulation with interleukin-4 (IL-4), B cells undergo proliferation and regulate immunoglobulin class switching (such as IgG1 and IgE), thereby facilitating T-cell development, which may contribute to the pathogenesis of vitiligo (65). However, no significant differences in the frequency of B cell subsets were observed in our data. Nevertheless, the role of B cell subpopulation heterogeneity in the pathogenesis and progression of vitiligo remains inadequately explored. In this study, we identified a subpopulation within memory B cells characterized by high EGR1 expression, which we defined as the EGR1+ memory B cell subpopulation. This subpopulation is enriched in pathways related to the cellular response to calcium ions and the regulation of hemopoiesis. GO analysis revealed significant enrichment of genes related to osteoclast differentiation and the IL-17 signaling pathway. Several genes that are not fully annotated, such as AC245014.3, SNORD3B-1, and AL021155.2, are also expressed in the EGR1+ memory B cell subpopulation. The expression of these genes may serve as critical factors in promoting B cell activation, thereby mediating the development of vitiligo. EGR1 mediates oxidative phosphorylation metabolic reprogramming in B-cell malignancies, thereby inducing resistance to ibrutinib (66). Studies in mouse models have shown that knockout of EGR1 inhibits the development of normal B cells and affects the differentiation of pro-B cells and immature B cells (67, 68). However, the reduction in EGR1 cell subpopulations and their functional characteristics in patients with vitiligo require further investigation and characterization.
Among the innate immune cell types, monocytes are less studied in vitiligo. A previous study reported that the proportion of CD80+ monocytes was significantly higher in the peripheral blood of vitiligo patients than in that of healthy controls, suggesting that monocytes may play a role in the development of vitiligo (24). In this study, we re-clustered monocytes and identified two subpopulations that were significantly reduced in vitiligo patients compared with normal controls. However, differential gene expression analysis did not distinguish specific subpopulations. Compared with the healthy controls, monocyte subpopulations expressed high levels of HLA-DRA, FOS, and the uncharacterized gene AC253572.2 in vitiligo patients. The literature also reports a significant increase in FOS expression in patients (69). Pathways related to leukocyte cell-cell adhesion and the regulation of cell-cell adhesion were notably active, as were pathways associated with viral infection.
This study has several limitations inherent to its design and implementation. First, the limited sample size of the discovery cohort for single-cell RNA-seq may compromise statistical robustness and limit the generalizability of our findings. Second, technical constraints precluded the complete validation of some transcriptionally defined cell clusters at the protein level, and the key molecular mechanisms we propose remain speculative without direct functional validation. Finally, the conclusion of impaired Treg cell function is primarily inferred from transcriptomic features suggestive of a Th1-like polarization, which requires confirmation through future in vitro functional assays.
In summary, our study delineates the intricate transcriptional landscape of immune cells in vitiligo, emphasizing key cell types and pathways that may play a pivotal role in its pathogenesis. Notably, we observed a reduction of KLRC2 expression in transitional T cells within NK cell and NKT cell subsets. These findings offer valuable insights into potential therapeutic targets for vitiligo and advance our understanding of its underlying immune mechanisms. However, further validation through larger sample sizes is required to substantiate our results.
Data availability statement
The original data presented in this study are openly available in the Genome Sequence Archive (GSA) for Human at the National Genomics Data Center (NGDC), China National Center for Bioinformation. The data can be accessed via the following URL: https://ngdc.cncb.ac.cn/gsa-human/ under the accession number HRA010691. All analysis code is provided in Data Sheet 2.
Ethics statement
The studies involving humans were approved by Ethics Committee of Sichuan Provincial People’s Hospital. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.
Author contributions
JX: Writing – original draft, Writing – review & editing, Resources. XJG: Writing – original draft, Resources. LG: Writing – original draft, Resources. QL: Writing – review & editing, Visualization, Funding acquisition. KD: Writing – review & editing, Visualization. XXG: Writing – review & editing. HL: Writing – review & editing, Visualization. TW: Writing – review & editing, Visualization. LW: Data curation, Writing – review & editing. WW: Writing – review & editing, Data curation. CG: Formal analysis, Writing – review & editing, Data curation. ML: Writing – review & editing, Data curation, Formal analysis. BG: Writing – review & editing, Supervision. YJ: Validation, Writing – review & editing. YW: Writing – review & editing, Validation. YH: Writing – review & editing, Software, Validation. LH: Software, Supervision, Writing – review & editing. JY: Resources, Writing – review & editing, Supervision. YM: Conceptualization, Funding acquisition, Project administration, Supervision, Writing – review & editing.
Funding
The author(s) declared that financial support was received for this work and/or its publication. This work was supported by the National Natural Science Foundation of China (Grant Numbers 82501268 and 82502798), the China Postdoctoral Science Foundation (Grant Number 2025M772611), the Key Research and Development Project of the Sichuan Provincial Science and Technology Department (Grant Number 2021YFS0201), and the Internal Research Fund of Sichuan Provincial People's Hospital (Grant Numbers 2023QN20, 24QNPY023 and 24QNPY030).
Acknowledgments
The authors would like to express their gratitude to the patients who visited the Department of Dermatology at Sichuan Provincial People’s Hospital. We are especially thankful for their valuable contribution of blood samples, which greatly supported this study. Furthermore, we thank OE Biotech Co., LTD (Shanghai, China) for providing single-cell sequencing for this study. We also acknowledge Sheng Pan, Suling Zhang and Qingzhen Zi for their assistance with bioinformatics analysis.
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 Generative AI was used in the creation of this manuscript. During the preparation of this work, the author used ChatGPT-4o only for the purpose of language polishing and refinement. After using this tool, the author thoroughly reviewed and edited the content as needed and takes full responsibility for the entirety of the published work, including its academic integrity and accuracy.
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.1698566/full#supplementary-material
Supplementary Figure 1 | UMAP visualization was employed to analyze and compare the overall cellular distributions, as well as the specific patterns of NK cells, NKT cells, B cells, monocytes, and T cells, between vitiligo patients and healthy controls. (A) UMAP visualization showing the overall distribution of total cells from each individual sample. (B) Distribution of NK and NKT cells across different individuals. (C) Distribution of T cell populations across different individuals. Each color represents a distinct cell cluster annotated based on transcriptional signatures. (D) Distribution of B cell populations across different individuals. (E) Distribution of monocyte populations across different individuals.
Supplementary Figure 2 | UMAP visualization depicting the expression distribution and inter-sample heterogeneity of the filtered marker genes (A-D) within their corresponding immune cell subsets. (A) Expression of KLRC2, a marker for cytotoxic NK cells, within the NK cell subset. (B) Expression of FCGR3A (encoding CD16) within the CD8+ T cell subset. (C) Expression of STAM within the regulatory T cell (Treg) subset. (D) Expression of the early growth response gene EGR1 within the naïve B cell subset. All panels illustrate the distribution of gene expression levels across different samples.
Supplementary Figure 3 | KLRC2 Expression in Peripheral Blood Immune Subsets and Total Leukocytes. (A) Gating strategy for flow cytometric analysis of KLRC2 expression in peripheral blood immune cell subsets, including monocytes, T cells, NK cells, B cells, and NKT cells. (B) Bar graph depicting KLRC2 expression in immune cell subsets of normal controls (HCs) and vitiligo patients, with the x-axis representing specific cell subsets and the y-axis indicating the proportion of KLRC2 expression in these subsets. (C) Heatmap of gene expression from bulk RNA sequencing of peripheral blood samples from vitiligo patients and normal controls. (D) Bar graph depicting the expression of KLRC2 at the bulk transcriptional level. Asterisks denote statistical significance, with * indicating P<0.05 and ** indicating P<0.01.
References
1. Cheuk S, Schlums H, Gallais Serezal I, Martini E, Chiang SC, Marquardt N, et al. CD49a expression defines tissue-resident CD8(+) T cells poised for cytotoxic function in human skin. Immunity. (2017) 46:287–300. doi: 10.1016/j.immuni.2017.01.009
2. Ezzedine K, Eleftheriadou V, Whitton M, and van Geel N. Vitiligo. Lancet. (2015) 386:74–84. doi: 10.1016/S0140-6736(14)60763-7
3. Silverberg NB. Update on childhood vitiligo. Curr Opin Pediatr. (2010) 22:445–52. doi: 10.1097/MOP.0b013e32833b6ac3
4. Zhang Y, Cai Y, Shi M, Jiang S, Cui S, Wu Y, et al. The prevalence of vitiligo: A meta-analysis. PloS One. (2016) 11:e0163806. doi: 10.1371/journal.pone.0163806
5. Lee JH, Ju HJ, Seo JM, Almurayshid A, Kim GM, Ezzedine K, et al. Comorbidities in patients with vitiligo: A systematic review and meta-analysis. J Invest Dermatol. (2023) 143:777–789.e6. doi: 10.1016/j.jid.2022.10.021
6. Riding RL and Harris JE. The role of memory CD8(+) T cells in vitiligo. J Immunol. (2019) 203:11–9. doi: 10.4049/jimmunol.1900027
7. Tulic MK, Cavazza E, Cheli Y, Jacquel A, Luci C, Cardot-Leccia N, et al. Innate lymphocyte-induced CXCR3B-mediated melanocyte apoptosis is a potential initiator of T-cell autoreactivity in vitiligo. Nat Commun. (2019) 10:2178. doi: 10.1038/s41467-019-09963-8
8. Hegab DS and Attia MA. Decreased circulating T regulatory cells in Egyptian patients with nonsegmental vitiligo: correlation with disease activity. Dermatol Res Pract. (2015) 2015:145409. doi: 10.1155/2015/145409
9. Giri PS, Mistry J, and Dwivedi M. Meta-analysis of alterations in regulatory T cells’ Frequency and suppressive capacity in patients with vitiligo. J Immunol Res. (2022) 2022:6952299. doi: 10.1155/2022/6952299
10. Gellatly KJ, Strassner JP, Essien K, Refat MA, Murphy RL, Coffin-Schmitt A, et al. scRNA-seq of human vitiligo reveals complex networks of subclinical immune activation and a role for CCR5 in T(reg) function. Sci Transl Med. (2021) 13:eabd8995. doi: 10.1126/scitranslmed.abd8995
11. Xu Z, Chen D, Hu Y, Jiang K, Huang H, Du Y, et al. Anatomically distinct fibroblast subsets determine skin autoimmune patterns. Nature. (2022) 601:118–24. doi: 10.1038/s41586-021-04221-8
12. Pigmentary Disorder Group, C.o.T. and W.M. Dermatology. Consensus on the diagnosis and treatment of vitiligo (2021 version). Chin J Dermatol. (2021) 54:105–9. doi: 10.35541/cjd.20200785
13. Wolock SL, Lopez R, and Klein AM. Scrublet: computational identification of cell doublets in single-cell transcriptomic data. Cell Syst. (2019) 8:281–291.e9. doi: 10.1016/j.cels.2018.11.005
14. Aran D, Looney AP, Liu L, Wu E, Fong V, Hsu A, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol. (2019) 20:163–72. doi: 10.1038/s41590-018-0276-y
15. Qiu X, Mao Q, Tang Y, Wang L, Chawla R, Pliner HA, et al. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods. (2017) 14:979–82. doi: 10.1038/nmeth.4402
16. Smyth GK. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. (2004) 3:Article3. doi: 10.2202/1544-6115.1027
17. Chen S, Zhou Y, Chen Y, and Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. (2018) 34:i884–90. doi: 10.1093/bioinformatics/bty560
18. Kim D, Langmead B, and Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. (2015) 12:357–60. doi: 10.1038/nmeth.3317
19. Roberts A, Trapnell C, Donaghey J, Rinn JL, and Pachter L. Improving RNA-Seq expression estimates by correcting for fragment bias. Genome Biol. (2011) 12:R22. doi: 10.1186/gb-2011-12-3-r22
20. Anders S, Pyl PT, and Huber W. HTSeq–a Python framework to work with high-throughput sequencing data. Bioinformatics. (2015) 31:166–9. doi: 10.1093/bioinformatics/btu638
21. Love MI, Huber W, and Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8
22. Wang Y, Zhang X, Ge J, Jin J, Zheng Z, Li J, et al. Single-cell landscape of peripheral blood mononuclear cells in patients with graves disease. Endocrinology. (2025) 166. doi: 10.1210/endocr/bqaf038
23. Gay-Mimbrera J, Gómez-Arias PJ, Álvarez-Heredia P, Batista-Duharte A, Rivera-Ruiz I, Aguilar-Luque M, et al. Integrated single-cell chromatin and transcriptomic analyses of peripheral immune cells in patients with alopecia areata. Front Immunol. (2025) 16:1565241. doi: 10.3389/fimmu.2025.1565241
24. Basak PY, Adiloglu AK, Koc IG, Tas T, Akkaya VB, et al. Evaluation of activatory and inhibitory natural killer cell receptors in non-segmental vitiligo: a flow cytometric study. J Eur Acad Dermatol Venereol. (2008) 22:970–6. doi: 10.1111/j.1468-3083.2008.02681.x
25. Mahmoud F, Abul H, al-Saleh Q, Haines D, Burleson J, and Morgan G. Peripheral T-cell activation in non-segmental vitiligo. J Dermatol. (1998) 25:637–40. doi: 10.1111/j.1346-8138.1998.tb02472.x
26. Ben Ahmed M, Zaraa I, Rekik R, Elbeldi-Ferchiou A, Kourda N, Belhadj Hmida N, et al. Functional defects of peripheral regulatory T lymphocytes in patients with progressive vitiligo. Pigment Cell Melanoma Res. (2012) 25:99–109. doi: 10.1111/j.1755-148X.2011.00920.x
27. Dwivedi M, Kemp EH, Laddha NC, Mansuri MS, Weetman AP, and Begum R. Regulatory T cells in vitiligo: Implications for pathogenesis and therapeutics. Autoimmun Rev. (2015) 14:49–56. doi: 10.1016/j.autrev.2014.10.002
28. Sirufo MM, De Pietro F, Bassino EM, Ginaldi L, and De Martinis M. Osteoporosis in skin diseases. Int J Mol Sci. (2020) 21. doi: 10.3390/ijms21134749
29. Lin X, Tian H, and Xianmin M. Possible roles of B lymphocyte activating factor of the tumour necrosis factor family in vitiligo autoimmunity. Med Hypotheses. (2011) 76:339–42. doi: 10.1016/j.mehy.2010.10.034
30. Schallreuter KU, Gibbons NC, Zothner C, Abou Elloof MM, and Wood JM. Hydrogen peroxide-mediated oxidative stress disrupts calcium binding on calmodulin: more evidence for oxidative stress in vitiligo. Biochem Biophys Res Commun. (2007) 360:70–5. doi: 10.1016/j.bbrc.2007.05.218
31. Hasse S, Gibbons NC, Rokos H, Marles LK, and Schallreuter KU. Perturbed 6-tetrahydrobiopterin recycling via decreased dihydropteridine reductase in vitiligo: more evidence for H2O2 stress. J Invest Dermatol. (2004) 122:307–13. doi: 10.1046/j.0022-202X.2004.22230.x
32. Eljaszewicz A, Kleina K, Grubczak K, Radzikowska U, Zembko P, Kaczmarczyk P, et al. Elevated numbers of circulating very small embryonic-like stem cells (VSELs) and intermediate CD14++CD16+ Monocytes in igA nephropathy. Stem Cell Rev Rep. (2018) 14:686–93. doi: 10.1007/s12015-018-9840-y
33. Moniuszko M, Bodzenta-Lukaszyk A, Kowal K, Lenczewska D, and Dabrowska M. Enhanced frequencies of CD14++CD16+, but not CD14+CD16+, peripheral blood monocytes in severe asthmatic patients. Clin Immunol. (2009) 130:338–46. doi: 10.1016/j.clim.2008.09.011
34. Saleh MN, Goldman SJ, LoBuglio AF, Beall AC, Sabio H, McCord MC, et al. CD16+ monocytes in patients with cancer: spontaneous elevation and pharmacologic induction by recombinant human macrophage colony-stimulating factor. Blood. (1995) 85:2910–7. doi: 10.1182/blood.V85.10.2910.bloodjournal85102910
35. Geissmann F, Manz MG, Jung S, Sieweke MH, Merad M, and Ley K. Development of monocytes, macrophages, and dendritic cells. Science. (2010) 327:656–61. doi: 10.1126/science.1178331
36. Liu H, Wang Y, Le Q, Tong J, and Wang H. The IFN-γ-CXCL9/CXCL10-CXCR3 axis in vitiligo: Pathological mechanism and treatment. Eur J Immunol. (2024) 54:e2250281. doi: 10.1002/eji.202250281
37. Chen Y, Griffiths CEM, and Bulfone-Paus S. Exploring mast cell-CD8 T cell interactions in inflammatory skin diseases. Int J Mol Sci. (2023) 24. doi: 10.3390/ijms24021564
38. Luo L, Zhu J, Guo Y, and Li C. Mitophagy and immune infiltration in vitiligo: evidence from bioinformatics analysis. Front Immunol. (2023) 14:1164124. doi: 10.3389/fimmu.2023.1164124
39. Zhang J, Yu S, Hu W, Wang M, Abudoureyimu D, Luo D, et al. Comprehensive analysis of cell population dynamics and related core genes during vitiligo development. Front Genet. (2021) 12:627092. doi: 10.3389/fgene.2021.627092
40. Yeh S, Li Z, Sen HN, Lim WK, Gill F, Perkins K, et al. Scleritis and multiple systemic autoimmune manifestations in chronic natural killer cell lymphocytosis associated with elevated TCRalpha/beta+CD3+CD4-CD8- double-negative T cells. Br J Ophthalmol. (2010) 94:748–52. doi: 10.1136/bjo.2009.171264
41. Kucuksezer UC, Aktas Cetin E, Esen F, Tahrali I, Akdeniz N, Gelmez MY, et al. The role of natural killer cells in autoimmune diseases. Front Immunol. (2021) 12:622306. doi: 10.3389/fimmu.2021.622306
42. Zhang Z, Ji W, Huang J, Zhang Y, Zhou Y, Zhang J, et al. Characterization of the tumour microenvironment phenotypes in Malignant tissues and pleural effusion from advanced osteoblastic osteosarcoma patients. Clin Transl Med. (2022) 12:e1072. doi: 10.1002/ctm2.1072
43. Vietzen H, Dohler B, Tran TH, Susal C, Halloran PF, Eskandary F, et al. Deletion of the natural killer cell receptor NKG2C encoding KLR2C gene and kidney transplant outcome. Front Immunol. (2022) 13:829228. doi: 10.3389/fimmu.2022.829228
44. Vietzen H, Zoufaly A, Traugott M, Aberle J, Aberle SW, and Puchhammer-Stockl E. Deletion of the NKG2C receptor encoding KLRC2 gene and HLA-E variants are risk factors for severe COVID-19. Genet Med. (2021) 23:963–7. doi: 10.1038/s41436-020-01077-7
45. Yu R, Broady R, Huang Y, Wang Y, Yu J, Gao M, et al. Transcriptome analysis reveals markers of aberrantly activated innate immunity in vitiligo lesional and non-lesional skin. PloS One. (2012) 7:e51040. doi: 10.1371/journal.pone.0051040
46. Arias M, Martínez-Lostao L, Santiago L, Ferrandez A, Granville DJ, and Pardo J. The untold story of granzymes in oncoimmunology: novel opportunities with old acquaintances. Trends Cancer. (2017) 3:407–22. doi: 10.1016/j.trecan.2017.04.001
47. Le Poole IC, van den Wijngaard RM, Westerhof W, and Das PK. Presence of T cells and macrophages in inflammatory vitiligo skin parallels melanocyte disappearance. Am J Pathol. (1996) 148:1219–28.
48. Harris JE, Harris TH, Weninger W, Wherry EJ, Hunter CA, and Turka LA. A mouse model of vitiligo with focused epidermal depigmentation requires IFN-gamma for autoreactive CD8(+) T-cell accumulation in the skin. J Invest Dermatol. (2012) 132:1869–76. doi: 10.1038/jid.2011.463
49. Dijkgraaf FE, Matos TR, Hoogenboezem M, Toebes M, Vredevoogd DW, Mertz M, et al. Tissue patrol by resident memory CD8(+) T cells in human skin. Nat Immunol. (2019) 20:756–64. doi: 10.1038/s41590-019-0404-3
50. Searle T, Al-Niaimi F, and Ali FR. Vitiligo: an update on systemic treatments. Clin Exp Dermatol. (2021) 46:248–58. doi: 10.1111/ced.14435
51. Takeshita T, Arita T, Higuchi M, Asao H, Endo K, Kuroda H, et al. STAM, signal transducing adaptor molecule, is associated with Janus kinases and involved in signaling for cell growth and c-myc induction. Immunity. (1997) 6:449–57. doi: 10.1016/S1074-7613(00)80288-5
52. Liang J and Yin H. STAM transports STING oligomers into extracellular vesicles, down-regulating the innate immune response. J Extracell Vesicles. (2023) 12:e12316. doi: 10.1002/jev2.12316
53. Choi SC, Park YP, Roach T, Jimenez D, Fisher A, Zadeh M, et al. Lupus susceptibility gene Pbx1 controls the development, stability, and function of regulatory T cells via Rtkn2 expression. Sci Adv. (2024) 10:eadi4310. doi: 10.1126/sciadv.adi4310
54. Butcher MJ, Park YP, Roach T, Jimenez D, Fisher A, Zadeh M, et al. Atherosclerosis-driven treg plasticity results in formation of a dysfunctional subset of plastic IFNgamma+ Th1/tregs. Circ Res. (2016) 119:1190–203. doi: 10.1161/CIRCRESAHA.116.309764
55. Ohkura N and Sakaguchi S. Transcriptional and epigenetic basis of Treg cell development and function: its genetic anomalies or variations in autoimmune diseases. Cell Res. (2020) 30:465–74. doi: 10.1038/s41422-020-0324-7
56. Feng T, Cao AT, Weaver CT, Elson CO, and Cong Y. Interleukin-12 converts Foxp3+ regulatory T cells to interferon-γ-producing Foxp3+ T cells that inhibit colitis. Gastroenterology. (2011) 140:2031–43. doi: 10.1053/j.gastro.2011.03.009
57. Liu X, Zhang W, Han Y, Cheng H, Liu Q, Ke S, et al. FOXP3(+) regulatory T cell perturbation mediated by the IFNγ-STAT1-IFITM3 feedback loop is essential for anti-tumor immunity. Nat Commun. (2024) 15:122. doi: 10.1038/s41467-023-44391-9
58. Okamoto M, Kuratani A, Okuzaki D, Kamiyama N, Kobayashi T, Sasai M, et al. IFN-γ-induced Th1-Treg polarization in inflamed brains limits exacerbation of experimental autoimmune encephalomyelitis. Proc Natl Acad Sci U.S.A. (2024) 121:e2401692121. doi: 10.1073/pnas.2401692121
59. Nie H, Zheng Y, Li R, Guo TB, He D, Fang L, et al. Phosphorylation of FOXP3 controls regulatory T cell function and is inhibited by TNF-α in rheumatoid arthritis. Nat Med. (2013) 19:322–8. doi: 10.1038/nm.3085
60. Liang W, Chen Q, Cheng S, Wei R, Li Y, Yao C, et al. Skin chronological aging drives age-related bone loss via secretion of cystatin-A. Nat Aging. (2022) 2:906–22. doi: 10.1038/s43587-022-00285-x
61. Clemenceau B, Vivien R, Berthome M, Robillard N, Garand R, Gallot G, et al. Effector memory alphabeta T lymphocytes can express FcgammaRIIIa and mediate antibody-dependent cellular cytotoxicity. J Immunol. (2008) 180:5327–34. doi: 10.4049/jimmunol.180.8.5327
62. Clemenceau B, Congy-Jolivet N, Gallot G, Vivien R, Gaschet J, Thibault G, et al. Antibody-dependent cellular cytotoxicity (ADCC) is mediated by genetically modified antigen-specific human T lymphocytes. Blood. (2006) 107:4669–77. doi: 10.1182/blood-2005-09-3775
63. Schreibing F, Hannani MT, Kim H, Nagai JS, Ticconi F, Fewings E, et al. Dissecting CD8+ T cell pathology of severe SARS-CoV-2 infection by single-cell immunoprofiling. Front Immunol. (2022) 13:1066176. doi: 10.3389/fimmu.2022.1066176
64. Martos SN, Campbell MR, Lozoya OA, Wang X, Bennett BD, Thompson IJB, et al. Single-cell analyses identify dysfunctional CD16(+) CD8 T cells in smokers. Cell Rep Med. (2020) 1. doi: 10.2139/ssrn.3517537
65. Imran M, Laddha NC, Dwivedi M, Mansuri MS, Singh J, Rani R, et al. Interleukin-4 genetic variants correlate with its transcript and protein levels in patients with vitiligo. Br J Dermatol. (2012) 167:314–23. doi: 10.1111/j.1365-2133.2012.11000.x
66. Liu Y, Kimpara S, Hoang NM, Daenthanasanmak A, Li Y, Lu L, et al. EGR1-mediated metabolic reprogramming to oxidative phosphorylation contributes to ibrutinib resistance in B-cell lymphoma. Blood. (2023) 142:1879–94. doi: 10.1182/blood.2023020142
67. Gururajan M, Simmons A, Dasu T, Spear BT, Calulot C, Robertson DA, et al. Early growth response genes regulate B cell development, proliferation, and immune response. J Immunol. (2008) 181:4590–602. doi: 10.4049/jimmunol.181.7.4590
68. Dinkel A, Warnatz K, Ledermann B, Rolink A, Zipfel PF, Burki K, et al. The transcription factor early growth response 1 (Egr-1) advances differentiation of pre-B and immature B cells. J Exp Med. (1998) 188:2215–24. doi: 10.1084/jem.188.12.2215
Keywords: GP-NSV, single-cell RNA sequencing, NK cell, KLRC2, Cytotoxic CD8⁺ T cell, FCGR3A, EGR1, B cell
Citation: Xiao J, Guo X, Gong L, Li Q, Dong K, Guo X, Li H, Wang T, Wang L, Wu W, Gan C, Li M, Gong B, Jiang Y, Wu Y, Hu Y, Han L, Yang J and Mao Y (2025) Single-cell transcriptomics reveals systemic immune dysregulation in non-segmental vitiligo. Front. Immunol. 16:1698566. doi: 10.3389/fimmu.2025.1698566
Received: 03 September 2025; Accepted: 19 November 2025; Revised: 15 November 2025;
Published: 08 December 2025.
Edited by:
Suneel Kumar, Rutgers, The State University of New Jersey, United StatesReviewed by:
Manoj Kumar Tembhre, All India Institute of Medical Sciences, IndiaMai Sharabi, Fayoum University, Egypt
Copyright © 2025 Xiao, Guo, Gong, Li, Dong, Guo, Li, Wang, Wang, Wu, Gan, Li, Gong, Jiang, Wu, Hu, Han, Yang and Mao. 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.
*Correspondence: Yujie Mao, eXVqaWVtYW9AeWVhaC5uZXQ=; Jianing Yang, ZHJ5YW5nam5AMTI2LmNvbQ==; Liang Han, NTQ2MjMxMTk0QHFxLmNvbQ==
†These authors have contributed equally to this work