Multimodal analyses of vitiligo skin identify tissue characteristics of stable disease

Abstract

Vitiligo is an autoimmune skin disease characterized by the destruction of melanocytes by autoreactive CD8+ T cells. Melanocyte destruction in active vitiligo is mediated by CD8+ T cells, but the persistence of white patches in stable disease is poorly understood. The interaction between immune cells, melanocytes, and keratinocytes in situ in human skin has been difficult to study due to the lack of proper tools. We combine noninvasive multiphoton microscopy (MPM) imaging and single-cell RNA-Seq (scRNA-Seq) to identify subpopulations of keratinocytes in stable vitiligo patients. We show that, compared with nonlesional skin, some keratinocyte subpopulations are enriched in lesional vitiligo skin and shift their energy utilization toward oxidative phosphorylation. Systematic investigation of cell-to-cell communication networks show that this small population of keratinocyte secrete CXCL9 and CXCL10 to potentially drive vitiligo persistence. Pseudotemporal dynamics analyses predict an alternative differentiation trajectory that generates this new population of keratinocytes in vitiligo skin. Further MPM imaging of patients undergoing punch grafting treatment showed that keratinocytes favoring oxidative phosphorylation persist in nonresponders but normalize in responders. In summary, we couple advanced imaging with transcriptomics and bioinformatics to discover cell-to-cell communication networks and keratinocyte cell states that can perpetuate inflammation and prevent repigmentation.

Keywords: Autoimmunity; Dermatology; Diagnostic imaging; Expression profiling; Skin.

Figure 1. In vivo MPM images of vitiligo lesional and nonlesional skin showing metabolic changes with depth independent of sun exposure.

(A) Representative en-face MPM images from the stratum granulosum in nonlesional (A1) and lesional skin (A2) and from the basal layer in nonlesional (B1) and lesional skin (B2) of 1 vitiligo patient. Average mitochondrial clustering values (β values) based on Z stacks from all vitiligo patients (n = 12) as a function of depth for nonlesional (top right) and lesional (bottom right) skin are shown as spline fits. Data are shown as mean ± SD of the measurements for the images in all the Z stacks at each area. The labels A1, A2, B1, and B2 within the mitochondrial clustering panels represent the mitochondrial clustering values extracted from the panel’s respective labeled images. Scale bars: 20 μm. (B) Representative en-face MPM images from the stratum granulosum in sun-exposed (A1) and non–sun-exposed skin (A2) and from the basal layer in sun-exposed (B1) and non–sun-exposed skin (B2) of 5 healthy volunteers. (C) Distribution of the median β values (left) and β variability values (right) in nonlesional and lesional skin of 12 vitiligo patients; each value corresponds to a Z stack of images acquired in nonlesional and lesional skin. *P < 0.05 by 2-tailed t test.(D) Distribution of the median β values (left) and β variability values (right) in sun-exposed and non–sun-exposed skin of 5 healthy volunteers; each value corresponds to a Z stack of images acquired in non–sun-exposed and sun-exposed areas.

Figure 2. Single-cell isolation of nonlesional and lesional skin of vitiligo patients for scRNA-Seq.

(A) Schematic diagram of single-cell isolation and scRNA-Seq data analyses. (B) UMAP plot of the cells from all patients in both nonlesional (left) and lesional skin (right).

Figure 3. scRNA-Seq analyses of lesional and nonlesional skin reveal unique keratinocyte cell populations in vitiligo patients.

(A) Feature plots showing expression of the selected markers in the UMAP space of all cells. (B) High-density plot showing relative gene expression of key marker genes in different cell subpopulations. Each density plot is composed by bar charts, and bar height corresponds to the relative expression level of the gene in cells that is ordered from low to high. (C) Percentages of cell subpopulations across patients, lesional skin, and nonlesional skin (left). Comparison of the percentages of each cell subpopulation across lesional and nonlesional skin (middle). Comparison of the percentages of major cell types including keratinocytes, stress keratinocytes, melanocytes, and immune cells across lesional and nonlesional skin (right). The bar plot shows that the percentages of keratinocytes and melanocytes decrease, while the percentages of stress keratinocytes and immune cells increase in lesional skin compared with nonlesional skin.

Figure 4. Stress keratinocytes have a unique gene signature and are the main source of CXCL9 and CXCL10.

(A) Heatmap of scaled expression levels of top 10 differentially expressed genes between nonlesional and lesional keratinocytes and enriched Hallmark pathways of the highly expressed genes in lesional keratinocytes. (B) Dot plots of signature scores of WNT signaling and OxPhos pathway between nonlesional and lesional skin. The size represents the percentage of expressing cells, and colors indicate the scaled signature scores. (C) Heatmap of scaled expression levels of differentially expressed genes between stress keratinocytes and other keratinocytes. (D) Enriched Hallmark pathways of highly expressed genes in stress keratinocytes and other keratinocytes, respectively. (E) The composition of stress keratinocytes and other keratinocytes in nonlesional and lesional skin. (F) Dot plot of stress associated markers in nonlesional skin, lesional skin, and stress keratinocytes. The size represents the percentage of expressing cells, and colors indicate the scaled expression.

Figure 5. Stress keratinocytes have altered energy utilization and shift toward oxidative phosphorylation.

(A) Violin plots of signature scores of OxPhos, glycolysis, WNT signaling, IFN-γ, IFN-α, and inflammatory response across nonlesional skin, lesional skin, and stress keratinocytes. The 2-sided Wilcoxon rank-sum test was used to evaluate whether there are significant differences in the computed signature scores. (B) Enrichment analysis of 21 metabolic pathways in stress keratinocytes versus other keratinocytes. Each dot represents 1 pathway. The x axis represents the differential gene signature scores of each metabolic pathway between stress keratinocytes and other keratinocytes. The y axis represents the Pearson’s correlation of gene signature scores between each metabolic pathway and stress response. Gene signature scores of stress response were computed based on the marker genes of stress keratinocytes. Colors represent the P values from 2-sided Wilcoxon rank-sum tests of each gene signature score between stress keratinocyte and other keratinocytes. (C) The number of differentially expressed OxPhos- and glycolysis-related genes in stress keratinocytes versus other keratinocytes. (D) Heatmap showing the average expression of top 18 differentially expressed OxPhos-related genes between stress keratinocytes and other keratinocytes. The top green bars represent the difference in the proportion of expressed cells between stress keratinocytes and other keratinocytes. (E) RNAscope demonstrating Krt6A and Krt10 in situ hybridization in patient-matched lesional and nonlesional punch grafting tissue. DAPI (cyan) was used to stain nuclei and second harmonic generation (blue) demonstrating location of collagen.

Figure 6. Cell-to-cell communication analysis reveals major signaling changes between nonlesional and lesional vitiligo skin.

(A) Number of inferred interactions among all cell subpopulations between nonlesional (NL) and lesional (LS) skin. (B) The relative information flow of all significant signaling pathways within the inferred networks between nonlesional and lesional skin. Signaling pathways labeled in green represent pathways enriched in nonlesional skin, the middle ones colored by black are equally enriched in both nonlesional and lesional skin, and the ones colored by purple are more enriched in lesional skin. (C) Visualization of outgoing and incoming interaction strength of each cell subpopulation in the inferred cell-to-cell communication network of nonlesional (top) and lesional skin (bottom). The dot sizes are proportional to the number of total interactions associated with each cell subpopulation. Dashed circle indicates the most altered cell subpopulations when comparing the outgoing and incoming interaction strength between nonlesional and lesional skin. (D) The signaling changes associated with the 3 most altered cell subpopulations.

Figure 7. Keratinocyte-melanocyte and keratinocyte–immune cell communication is altered in lesional vitiligo skin compared with nonlesional skin.

(A) Bubble plot in left panel shows the decreased signaling from keratinocyte and immune subpopulations to melanocytes (nonlesional versus lesional skin). Bubble plot in right panel shows all significant signaling from stress keratinocyte to melanocytes and immune subpopulations. (B) Inferred cell-to-cell communication networks of WNT and CXCL signaling in nonlesional and lesional skin, respectively. The edge width is proportional to the inferred communication probabilities. The relative contribution of each ligand-receptor pair to the overall signaling pathways.

Figure 8. Pseudotemporal dynamics reveal transition dynamics of stress keratinocytes.

(A) Projection of keratinocytes onto the PHATE space revealed the potential lineage relationships between different keratinocyte subpopulations in nonlesional (NL, left panel) and lesional (LS, right panel) skin. Cells were colored by the annotated cell identity. (B) The inferred pseudotemporal trajectories of all cells using Monocle 3. Cells were colored by the inferred pseudotime. Pseudotemporal trajectory analysis revealed 2 potential transitional paths, as indicated by Path 1 and Path 2. (C) Pseudotemporal dynamics of all pseudotime-dependent genes along the Path 1 and Path 2. Each row (i.e., gene) is normalized to its peak value along the pseudotime. These genes were clustered into 5 groups with the average expression patterns (middle) and representative genes (right). Solid and dashed lines indicate the average expression of a particular gene group in Path 1 and Path 2, respectively. The number of genes in each gene group is indicated in parenthesis. (D) Enriched biological processes of the 5 gene groups in C.

Figure 9. Upregulation of stress response and OxPhos are seen in the reconstructed pseudotemporal dynamics of stress keratinocytes.

(A) The reconstructed pseudotemporal dynamics of selected marker genes along the inferred pseudotime in Path 1 and Path 2, respectively. Black lines represent the average temporal patterns that were obtained by fitting a cubic spline. Cells were colored by the inferred pseudotime. (B) Pseudotemporal dynamics of the pseudotime-dependent genes related with the stress response and along the inferred pseudotime in Path 1 and Path 2. (C) Pseudotemporal dynamics of the pseudotime-dependent genes related with OxPhos along the inferred pseudotime in Path 1 and Path 2.

Figure 10. Keratinocyte energy utilization normalize in vitiligo patients who respond to punch grafting treatment but persist in nonresponders.

(A) Representative clinical images from vitiligo patients undergoing punch grafting treatment. Clinical responder on top and nonresponder on the bottom. (B) Average mitochondrial clustering values (β values) based on Z stacks from 6 vitiligo patients as a function of depth for responders and nonresponders at baseline are shown as spline fits. Patients were followed and imaged again after 10 weeks. Average mitochondrial clustering values (β values) for clinical responders (n = 3) and nonresponders (n = 3) are shown. (C) Distribution of β variability values (right) in punch grafting responders and nonresponders (n = 6); each value corresponds to a Z stack of images acquired. *P < 0.05 by t test.