This study involved 815 patients with moderate-to-severe plaque psoriasis, diagnosed by certified dermatologists. Blood samples were collected from peripheral veins, and serum was extracted for further analysis. All participants had no history of autoimmune or systemic inflammatory conditions. Prior to participation, written informed consent was secured from each individual. The study adhered to the Declaration of Helsinki principles and received ethical approval from the Ethics Committee of Huashan Hospital, Fudan University under protocol number 2020-1202.

Researchers collected single-cell transcriptomic datasets encompassing six common skin conditions from public databases. They individually analyzed nine single-cell RNA sequencing datasets, which covered a variety of skin diseases. To avoid altered gene expression from tissue dissociation and ensure data comparability, single-cell datasets generated via fluorescence-activated cell sorting (FACS) were excluded.

During data processing, to ensure consistency, the research team excluded cells with a library size of less than 500, fewer than 200 detected genes, or a mitochondrial gene expression percentage exceeding 10%. They also utilized the DoubletFinder tool (https://github.com/chris-mcginnis-ucsf/DoubletFinder) to identify and remove potential doublets. Additionally, the research team converted Ensembl gene IDs to corresponding gene symbols using annotations from Biomart via the org.Hs.eg.db package (https://bioconductor.org/packages/org.Hs.eg.db/).

After initial quality control, the researchers independently processed the raw data for each skin condition. They adopted the standard workflow of Seurat (version 5.1.0, Satija Lab, New York, USA), utilizing a series of functions including NormalizeData, RunHarmony, RunUMAP, FindVariableFeatures, ScaleData, RunPCA, FindNeighbors, and FindClusters, with a clustering resolution of 0.5. To accurately identify NK cells and exclude the interference of other immune cell subpopulations, researchers defined cells co-expressing PTPRC (CD45), KLRB1 (CD161), and IL7R (CD127) as ILCs and cells co-expressing KLRD1 (CD94) as NK cells [15]. To further ensure the accuracy of identification, they also employed the automated annotation tool SingleR (https://bioconductor.org/packages/SingleR/) for independent validation. Through this process, the researchers identified NK cells and retained them for subsequent reintegration analysis.

Prior to dataset reintegration, only genes shared across all skin conditions were retained to ensure comparability. Harmony (v1.2.0) (https://github.com/immunogenomics/harmony) was applied to correct for batch effects among different donors. Based on parameters described in the previously mentioned Seurat workflow, subclusters consisting of fewer than 1000 cells and located distantly from main clusters were identified as potential doublets and removed. The remaining cells were designated as the final NK cell population.

The Harmony function was applied to the sorted NK cells for batch correction, and a clustering resolution of 0.5 was set in the FindClusters function. To annotate each subset, differentially expressed genes (DEGs) were identified using the Wilcox algorithm via the FindAllMarkers function.

Enrichment analyses for Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) were carried out using the enrichGO and enrichKEGG R packages, with a threshold of adjusted p-value 0.05. Canonical pathways from h.all.v6.2.symbols.gmt, sourced from the GSEA website (https://www.gsea-msigdb.org/gsea/), were analyzed using the AUCell R package (v1.12.0, https://github.com/aertslab/AUCell) with default settings to assess pathway activities.

To infer cellular developmental trajectories, pseudotime ordering and trajectory construction were performed using the core algorithms of Monocle2. Visualization was conducted with the R package Slingshot (v2.12.0, https://bioconductor.org/packages/Slingshot/).

Visualization of the single-cell datasets was conducted in the R environment using several key packages, including ggplot2, scRNAtoolVis, pheatmap, ComplexHeatmap, and ggvenn.

The tissue sections were first treated with eco-friendly dewaxing solutions I, II, and III (G1128, Servicebio, Wuhan, Hubei, China) for 10 minutes each. They were then dehydrated in three changes of absolute ethanol (100092183, Sinopharm Chemical Reagent, Shanghai, China), spending 5 minutes in each, before being rinsed with distilled water. Antigen retrieval was conducted by heating the slides in retrieval buffer for 30 minutes, ensuring minimal buffer evaporation and avoiding tissue dryness. After heating, the slides were cooled naturally to room temperature. Finally, the sections were washed in PBS (G0002, Servicebio, Wuhan, Hubei, China) (pH 7.4) on a shaking platform, with three 5-minute cycles.

After gently removing excess liquid, a hydrophobic barrier was drawn around the tissue using a PAP pen. Blocking was performed for 30 minutes with bovine serum albumin (BSA): for goat-derived primary antibodies, 10% donkey serum was used; for all other sources, 3% BSA was applied. Slides were incubated overnight at 4 °C in a humidified chamber with the following primary antibodies (dilution ratio: 1:500): IL13 (ABS120028, Absin, Shanghai, China), GZMB (GB15092, Servicebio, Wuhan, Hubei, China), and NCAM1 (GB12041, Servicebio, Wuhan, Hubei, China). The next day, slides were washed in PBS (pH 7.4) on a shaker three times, 5 minutes each.

Following this, the sections were incubated with the following secondary antibodies (dilution ratio: 1:300) in the dark at room temperature for 50 minutes: CY3-conjugated goat anti-rabbit IgG (GB21303, Servicebio, Wuhan, Hubei, China) and CY3-conjugated goat anti-mouse IgG (GB21301, Servicebio, Wuhan, Hubei, China). After incubation, the slides were washed three times with PBS (pH 7.4), for 5 minutes each. Nuclear staining was performed using DAPI, which was applied and allowed to incubate under the same dark conditions for 10 minutes, after which another series of PBS washes (3 $\times$ 5 minutes) was carried out. Autofluorescence was reduced by treating the sections with quencher solution B for 5 minutes, followed by a 10-minute rinse under running water. Finally, the samples were coverslipped using an anti-fade mounting medium.

This kit utilizes a double-antibody sandwich assay. An anti-Human IL18 antibody (A1010A0118, Biotnt, Shanghai, China) coats the enzyme-linked immunosorbent plate. Human IL18 in standards and samples binds to this monoclonal antibody. Biotinylated anti-Human IL18 is added to form a plate-bound immune complex. Horseradish peroxidase-conjugated streptavidin attaches to biotin. After adding the substrate solution, a blue color forms. The reaction stops upon adding a termination solution. The Human IL18 concentration correlates directly with the OD value at 450 nm. A standard curve allows determination of Human IL18 levels in samples. The experiment used 80 human peripheral blood samples.

Introduce 100 µL of relevant standards and samples to the wells, seal, mix well, and incubate at 37 °C for 60 minutes. Thoroughly rinse the plate four times with 1$\times$ wash buffer and dry on absorbent paper. Add 100 µL of 1$\times$ diluted Biotinylated Antibody to each well, reseal, and incubate for 45 minutes at 37 °C. Rinse again four times, dry, then add 100 µL of 1$\times$ Streptavidin-HRP to each well, seal once more, and incubate for 15 minutes at 37 °C. After another four rinses and drying, add 100 µL of chromogen solution to each well, leave the plate in the dark at room temperature for 5–30 minutes (monitor color change and stop promptly if color deepens excessively). Add 100 µL of stop solution, mix well. Configure the microplate reader to 450 nm (primary) and 630 nm (secondary) for absorbance measurement, and read within 30 minutes after stopping the reaction.

The single-cell datasets used in this study are available in the Gene Expression Omnibus and the Genome Sequence Archive (Table 1, Ref. [15, 16, 17, 18, 19, 20, 21, 22, 23, 24]). The disease data we selected were all collected before any treatment was administered. The code used for bioinformatics analysis can be found at https://github.com/issac0124/1/issues.

Statistical analyses were performed using the “stats” package in R software (version 3.6.2, R Core Team, Vienna, Austria) and GraphPad Prism (version 8.0, GraphPad Software, LLC, La Jolla, CA, USA). A p-value less than 0.05 was considered statistically significant.

To investigate the heterogeneity of skin-resident natural killer cells, we collected publicly available single-cell transcriptomic datasets comprising 40 skin samples across six skin conditions, including healthy skin [16] (n = 2), psoriasis [17, 18, 19] (n = 13), atopic dermatitis [15, 20] (n = 9), lupus [21, 22] (n = 3), melanoma [24] (n = 2), and vitiligo [16, 23] (n = 11). A total of 12,028 NK cells were identified from the 40 integrated samples for subsequent analysis (Fig. 1A). Subclustering analysis revealed that NK cells from healthy individuals and patients with various skin diseases were grouped into nine distinct subsets, all exhibiting high expression of NK cell lineage–associated genes (Fig. 1B,C). Based on the analysis of differentially expressed genes (DEGs), these subsets were identified and annotated using established canonical markers, including GZMB⁺, GZMB⁺ IL13⁺, GZMB⁺ Regulator of G-protein signaling 1 (RGS1)⁺, GZMB⁺ Keratin (KRT)⁺, RGS1⁺AHR⁺, Transcription factor 7 (TCF7)⁺, TCF7⁺RGS1⁺, RGS1⁺ and Crystallin Alpha B (CRYAB)⁺ NK cells (Fig. 1D,E).

Fig. 1.

Single-cell RNA sequencing (scRNA-Seq) analysis reveals the heterogeneity of NK cells in human skin diseases. (A) Schematic workflow of single-cell RNA sequencing analysis. (B) UMAP embedding shows nine NK cell clusters derived from 12,028 cells across six skin conditions. (C) Dot plot displaying the expression of NK cell–related genes and lineage-specific markers in the skin. (D) Density plot illustrating representative markers for each NK cell subset. (E) Heatmap depicting distinct marker genes for each NK cell cluster. NK, natural killer; UMAP, Uniform Manifold Approximation and Projection.

We next examined the distribution of these NK cell subsets in various skin diseases. Their abundance, visualized through 2D embeddings and bar plots, varied significantly across conditions (Supplementary Fig. 1). Notably, GZMB⁺ NK cells were present in PSO and also exhibited a certain level of expression in AD, which is consistent with previous reports [25]. GZMB⁺IL13⁺ NK cells were predominantly found in AD. In melanoma, all NK cells expressed RGS1, which is a marker associated with the prognosis of melanoma [26, 27]. High expression of RGS1 indicates a poor prognosis for melanoma [28]. TCF7⁺ NK cells are primarily expressed in lupus. Previous studies have identified TCF7 as a susceptibility locus for lupus [29]. Furthermore, RGS1⁺ NK cells were found to be unique to vitiligo. These results collectively highlight the dynamic and diverse phenotypic changes of NK cells in the progression of various skin disorders.

To assess the functional characteristics of different NK cell subsets, we identified differentially expressed genes (DEGs) for each subset (Fig. 2A). CRYAB⁺, GZMB⁺, GZMB⁺RGS1⁺, GZMB⁺IL13⁺, TCF7⁺RGS1⁺, and TCF7⁺ NK cells show upregulation of genes related to immune response, including immune cell chemotaxis (CCL2, CCL19, CXCL2, CXCL13, CCL3, CCR8, SELL, CCR7), lymphocyte development, differentiation, and activation (TCF7, PRDM1, NR4A1, FCER1G, FCMR), and cytotoxicity (FGFBP2, GZMH, PRF1, FCGR3A, CRTAM, EFHD2). GZMB⁺IL13⁺ NK cells also express a variety of cytokines (IL13, IL22), which are associated with the pathogenesis of atopic dermatitis. GZMB⁺KRT⁺ NK cells express genes associated with epidermal cells (CALML5, KRT2, S100A7, KRTDAP). RGS1⁺ and RGS1⁺AHR⁺ NK cells are associated with stress and metabolism (MT1G, MT1E, MT1X, ATP5E, ATP5L, ATP5I) (Fig. 2A,B). These findings suggest that CRYAB⁺, GZMB⁺, GZMB⁺RGS1⁺, GZMB⁺IL13⁺, TCF7⁺RGS1⁺, and TCF7⁺ NK cells may exhibit pro-inflammatory properties, thereby promoting skin inflammation. In addition, the GO and KEGG enrichment analysis of the selected DEGs indicated that GZMB⁺KRT⁺ NK cells are associated with the IL17 signaling pathway, while GZMB⁺RGS1⁺ NK cells are related to Th17 cell differentiation, which are respectively immune and differentiation pathways related to the pathogenesis of psoriasis. Meanwhile, GZMB⁺IL13⁺ NK cells are associated with the pathway of Bacterial invasion of epithelial cells, which may be related to the impaired skin barrier in atopic dermatitis. CRYAB⁺ NK cells are associated with cell-to-cell junctions, while RGS1⁺AHR⁺ NK cells are related to antigen processing and presentation. RGS1⁺ is associated with various autoimmune diseases, TCF7⁺ is mainly related to multiple signaling pathways, while TCF7⁺RGS1⁺ is primarily associated with various pathogen infections. In addition, GZMB⁺ NK cells are associated with chemotaxis, which may be related to the inflammatory response in the skin of psoriasis (Fig. 2C).

Fig. 2.

Transcriptional programs of the NK cell subsets in the human skin NK cell atlas. (A) Key genes for each NK cell subset across six skin conditions. (B) Top Gene Ontology (GO) pathways enriched in the top 100 upregulated genes in each subset, based on average log₂ fold change versus other skin NK cell subsets. (C) Top Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched in the top 100 upregulated genes in each subset. (D) Violin plots showing AUCell scores for subset-specific gene sets in each NK cell cluster.

To confirm the immune roles of NK cells, we evaluated their cytokine response profiles using area under the curve (AUC) analysis (Fig. 2D). Our results indicate that RGS1⁺ NK cells exhibit elevated oxidative phosphorylation activity, while GZMB⁺ and TCF7⁺ NK cells are associated with the complement system (Fig. 2D).

We further examined whether these cells adopt distinct transcriptional programs to promote inflammation in each disease. In psoriasis, GZMB⁺ NK cells show upregulation of genes associated with lymphocyte migration (e.g., CX3CR1, SPON2, CHST2) (Fig. 3A). The GO-enriched pathways are associated with the activation and effector functions of immune cells (Fig. 3B). KEGG enrichment analysis indicated that in psoriasis, GZMB⁺ NK cells retained the NK cell-mediated cytotoxicity signature and were involved in the infection process of pathogens and the migration of leukocytes across the endothelium (Fig. 3C). In atopic dermatitis, GZMB⁺ NK cells express genes associated with cell metabolism (e.g., AKR1C3, ODC1, GCLM) and cell proliferation and differentiation (e.g., ZEB2, LYAR) (Fig. 3D). The GO-enriched pathways are associated with immune response, negative regulation of lymphocyte activation, and cell killing (Fig. 3E). KEGG enrichment analysis indicated involvement in the Chemokine signaling pathway, Viral protein interaction with cytokine and cytokine receptor, and Efferocytosis (Fig. 3F).

Fig. 3.

Transcriptional changes of GZMB⁺ NK cells shared between psoriasis and atopic dermatitis. (A) Volcano plot showing DEGs between GZMB⁺ NK cells and other cells in psoriasis. (B) Bar chart showing GO-enriched pathways of DEGs between GZMB⁺ and other cells in psoriasis. (C) Dot plot showing KEGG pathways of DEGs between GZMB⁺ and other cells in psoriasis. (D) Volcano plot illustrating DEGs between GZMB⁺ NK cells and other cells in atopic dermatitis. (E) Bar chart showing GO-enriched pathways of DEGs between GZMB⁺ and other cells in atopic dermatitis. (F) Dot plot showing KEGG pathways of DEGs between GZMB⁺ and other cells in atopic dermatitis.

In psoriasis, NK cells expressing GZMB show upregulation of genes associated with cell motility (e.g., pFN1, CFL1) and chemotaxis (e.g., CCL3, CCL4, CCL5) (Fig. 4A). GO enrichment analysis revealed pathways involved in various immune signaling and immune activation (Fig. 4B). KEGG analysis based on the selected DEGs showed that in psoriasis, NK cells expressing GZMB are functionally enriched in Th17 cell differentiation and Lipid and atherosclerosis (Fig. 4C). GSEA enrichment analysis revealed that in psoriasis, NK cells not expressing GZMB show upregulation of the Cornified envelope formation, Coronavirus disease-COVID-19, and Ribosome pathways (Fig. 4D). GZMB⁺IL13⁺ NK cells show upregulation of genes associated with cell metabolism (e.g., ALDOA, EEF1G) and mitochondrial-related functions (e.g., MICOS10, MICOS13) (Fig. 4E). In atopic dermatitis, GZMB⁺IL13⁺ NK cells are associated with GO pathways related to T-cell activation and lymphocyte differentiation (Fig. 4F). KEGG analysis based on the selected DEGs showed that in atopic dermatitis, IL13⁺ NK cells are functionally enriched in the Bacterial invasion of epithelial cells pathway and various receptor pathways (Fig. 4G). Gene Set Enrichment Analysis (GSEA) revealed that in atopic dermatitis, GZMB⁺IL13^– NK cells are enriched in the Toll-like receptor signaling pathway, Human cytomegalovirus infection, and Efferocytosis pathways (Fig. 4H).

Fig. 4.

Disease-related functional differences of Granzyme B (GZMB)⁺ NK cells between psoriasis and atopic dermatitis. (A) Differentially expressed genes (DEGs) between psoriasis NK cells expressing GZMB and those not expressing GZMB. (B) Bar chart showing GO-enriched pathways of DEGs between NK cells expressing GZMB and those not expressing GZMB in psoriasis. (C) Dot plot showing KEGG pathways of DEGs between NK cells expressing GZMB and those not expressing GZMB in psoriasis. (D) GSEA enrichment plot showing representative pathways enriched in psoriasis NK cells expressing GZMB and those not expressing GZMB. (E) Volcano plot illustrating DEGs between GZMB⁺IL13⁺ NK and GZMB⁺IL13^– NK cells in atopic dermatitis. (F) Bar chart showing GO-enriched pathways of DEGs between GZMB⁺IL13⁺ NK and GZMB⁺IL13^– NK cells in atopic dermatitis. (G) Dot plot showing KEGG pathways of DEGs between GZMB⁺IL13⁺ NK and GZMB⁺IL13^– NK cells in atopic dermatitis. (H) Gene Set Enrichment Analysis (GSEA) enrichment plot showing representative pathways enriched in GZMB⁺IL13⁺ NK and GZMB⁺IL13^– NK cells in atopic dermatitis.

Psoriasis and atopic dermatitis exhibit high complexity; in atopic dermatitis, cytokines such as IL13, IL2, and IL22 are significantly elevated (Fig. 5A,B). Psoriasis is mainly enriched in pathways related to cell metabolism (Fig. 5C). In addition, KEGG analysis based on the selected DEGs showed that psoriasis is associated with the IL17 signaling pathway (Fig. 5D). Among the changes in the proportions of NK cell subpopulations in the two diseases, the most significant differences were observed in GZMB⁺ NK cells and GZMB⁺IL13⁺ NK cells (Fig. 5E). Pseudotime analysis indicated that the GZMB⁺ NK cell population could exhibit the expression of RGS1 during the differentiation process, but ultimately would show the expression of IL13 (Fig. 5F). During the simulated differentiation process, genes associated with the inflammatory response in atopic dermatitis (such as the S100 protein family) and genes related to immune regulation (TNFRSF18) were upregulated, which may contribute to the development of paradoxical reactions during psoriasis treatment (Fig. 5G). Taken together, GZMB⁺ NK cells in psoriasis patients may gradually differentiate to express the cytokine IL13, thereby promoting the development of paradoxical reactions.

Fig. 5.

Dynamic transcriptional features of GZMB⁺ NK cells in psoriasis. (A) Volcano plot showing differentially expressed genes (DEGs) in NK cells between psoriasis and atopic dermatitis. (B) Heatmap illustrating gene expression profiles of NK cells in psoriasis and atopic dermatitis. (C) Bar plot showing Gene Ontology (GO) enrichment pathways of NK cell DEGs in psoriasis. (D) Dot plot presenting KEGG enrichment pathways of NK cell DEGs in psoriasis. (E) Pie chart illustrating the proportion of NK cell subsets in psoriasis and atopic dermatitis. (F) Pseudotime trajectory analysis of NK cells across six skin conditions using Monocle2; lines represent inferred cellular trajectories. (G) Heatmap showing dynamic changes in gene expression during the differentiation of GZMB⁺ NK cells.

To verify our hypothesis, we conducted a retrospective study on samples from patients with moderate-to-severe plaque psoriasis who were treated with IL17A inhibitors. We identified 48 patients whose symptoms improved significantly and 32 who experienced paradoxical reactions (Fig. 6A,B). Our findings revealed substantially higher levels of GZMB and IL13 in the peripheral blood of patients with paradoxical reactions than in the control group (Fig. 6C). Immunofluorescence studies also validated the increased expression of GZMB in psoriasis and IL13 in atopic dermatitis (AD) (Fig. 6D). Moreover, both GZMB and IL13 signals were detected in the dermis and at the dermoepidermal junction in both conditions, implying a possible transformation between the two NK cell phenotypes. By analyzing the public single-cell transcriptomic data of psoriasis patients treated with biologics, such as IL17 and IL23 inhibitors, we observed an upregulation of GZMB expression post-treatment [17, 30] (Fig. 6E), which is consistent with our finding from the analysis of public databases that the proportion of GZMB is higher in atopic dermatitis than in psoriasis. In short, as psoriasis evolves into paradoxical reactions, GZMB⁺ NK cells might begin to express IL13, thereby playing a role in the development of these reactions.

Fig. 6.

Characteristics of natural killer cells in patients with paradoxical reactions. (A) Schematic diagram illustrating the workflow of clinical sample collection. (B) Baseline characteristics of patients with significant remission versus those who developed paradoxical reactions. (C) ELISA results showing the differential levels of GZMB and IL13 in patients with remission and those with paradoxical reactions. The blue axis is significant improvement (SI), and the red axis is paradoxical skin reaction (PSR). (D) Dot plots displaying changes in NK cell gene expression following treatment with IL23 inhibitors and IL17 inhibitors in patients with paradoxical reactions. (E) Immunofluorescence staining showing the abundance of GZMB⁺ and IL13⁺ cells in lesional tissues of psoriasis and atopic dermatitis. *, p $\le$ 0.05; **, p $\le$ 0.01; scale bar = 100 µm.