In depth transcriptomic profiling defines a landscape of dysfunctional immune responses in patients with VEXAS syndrome

Abstract

VEXAS (Vacuoles, E1 enzyme, X-linked, Autoinflammatory, Somatic) syndrome is caused by inactivating somatic mutations in the UBA1 gene. Here, we characterize the immunological landscape of VEXAS syndrome by performing multi-omics single-cell RNA analysis, cytokine multiplex assays, and in vitro functional assays on patients' peripheral blood. Our data reveals a broad immune system activation with upregulation of multiple inflammatory response pathways and proinflammatory cytokines. Unexpectedly, we find that monocytes have dysfunctional features irrespective of UBA1 mutation status, exhibiting impaired efferocytosis and blunted cytokine production in vitro. In contrast, UBA1-mutated NK cells show an upregulation of the inflammation pathways and enhanced cytotoxicity. Within the lymphocyte subsets, predominantly UBA1 wild-type, we identify clonal expansion of effector memory CD8⁺ T cells and skewed B cell differentiation with loss of transitional B cells and expansion of plasmablasts. Thus, our analysis indicates that VEXAS syndrome is characterized by profound alterations in both adaptive and innate immune systems, accounting for the complex pathophysiology of the disease, and provides a basis to understand the marked clinical heterogeneity and variable disease course.

Fig. 1. Multi-omics single-cell sequencing analysis of PBMCs from VEXAS patients and healthy donors.

a Overview of the experimental workflow. A figure was created with BioRender. Mizumaki, H. (2025) [https://BioRender.com/n37z398]. b Hemoglobin levels (HGB), platelet counts (PLT), white blood cell counts (WBC), neutrophil counts, monocyte counts, B cell counts, NK cell counts, and T cell counts from VEXAS patients (n = 9). Background shading shows a normal reference range for each parameter. Data are presented as mean and SD. c A Uniform Manifold Approximation and Projection (UMAP) plot of 178,168 cells from all subjects (n = 14, left). UMAP plots of 125,806 peripheral blood mononuclear cells (PBMCs) derived from VEXAS patients (n = 9, upper right) and 52,362 PBMCs derived from healthy donors (n = 5, lower right). Leiden clusters based on 5’ gene expression are shown and colored by major cell types. d A violin plot showing expression distributions of selected canonical marker genes in 10 major cell types. Rows and columns represent selected marker genes and cell types, respectively. e A UMAP plot of cells with projected mutation status assignment for wild-type UBA1 (wtUBA1; n = 18,140 cells) and mutated UBA1 (mtUBA1; n = 1402 cells). f Normalized frequency of mtUBA1 cells in major cell types. Cell types with more than 100 cells were analyzed. gdT gamma delta T cells, MAIT mucosa-associated invariant T cell. Source data are provided as a Source Data file.

Fig. 2. Compositional changes and activation of the inflammatory pathways in immune cells from VEXAS patients.

a A neighborhood graph of PBMCs using Milo differential abundance testing. Nodes represent neighborhoods from the PBMC population. A color scale indicates log2-fold change (FC) differences between VEXAS patients and healthy donors. Significant changes are colored in blue and red. Nondifferential abundance Nhoods (a false discovery rate [FDR] ≥ 0.10) are indicated in white. b Beeswarm and box plots showing the distribution of log2FC differences between VEXAS (n = 9) and healthy donors (n = 5) in neighborhoods in different cell type clusters. Colors are represented similarly to (a). A box plot shows median and interquartile ranges (IQR); lower and upper hinges correspond to the first and third quartiles, respectively. An upper whisker extends from the hinge to the largest value no further than 1.5*IQR from the hinge. A lower whisker extends from the hinge to the smallest value at most 1.5*IQR from the hinge. c UMAP plots overlaid with projections of the IFN-α, IFN-γ, and TNF response pathway gene module scores. d Gene module scores of IFN-α, IFN-γ, and TNF responses in PBMCs. A heatmap depicting differences between average scores of VEXAS patients and those of healthy donors in each of 10 cell types. Gene module scores of cells of VEXAS patients were compared to those of healthy donors in each cell type using the two-sided Welch’s t-test and shown as *P < 1.0 × 10^–10, **P < 1.0 × 10^–100. e Predicted ligand-receptor interactions of significantly upregulated genes in VEXAS patients as compared to healthy donors. f A volcano plot displaying normalized protein expression of cytokines in plasma measured by the Olink Target 96 Inflammation panel immunoassay. The plot shows 92 plasma cytokines differentially expressed between VEXAS patients and healthy donors. x and y axes represent a magnitude of a cytokine’s log2FC and a significance scale by the -log10 (p adjusted [padj]), respectively. A black dotted line indicates padj = 0.05. Significantly dysregulated cytokines are highlighted in red. gdT gamma delta T cells, MAIT mucosa-associated invariant T cell. Source data are provided as a Source Data file.

Fig. 3. Dysfunctional monocytes in VEXAS.

a A UMAP plot of 24,887 monocytes and 704 dendritic cells (DC). b A neighborhood graph of monocytes/DCs, which was generated similarly as in Fig. 2a. c Beeswarm and box plots of monocytes/DCs for VEXAS (n = 9) and healthy donors (n = 5), which were generated similarly as in Fig. 2b. d A volcano plot of differentially expressed genes of monocytes in VEXAS as compared to healthy donor: upregulated and downregulated genes in red and blue, respectively. A horizontal dotted line represents padj = 0.05, vertical dotted lines indicate absolute log2FC values = +–0.25. P values were calculated with the two-sided Wilcoxon rank-sum test and Bonferroni correction for multiple comparisons. e A dot plot showing gene set enrichment scores of top 10 upregulated hallmark pathways across subtypes. Dot sizes indicate mean normalized enrichment score (NES) differences between VEXAS and healthy donors, and the color scale indicates FDR values; non-significant pathways (FDR ≥ 0.10) in gray. f Gene module scores of HLA class II genes and sepsis-associated monocyte signatures in monocytes across VEXAS (n = 9), healthy donors (n = 5), and other multisystem autoimmune diseases (Systemic lupus erythematosus (SLE), n = 33; microscopic polyangiitis (MPA), n = 8; Behcet’s disease (BD), n = 4). Data are presented as mean with SD. P values were calculated with the two-sided unpaired Mann–Whitney U test. g A UMAP plot of cells with projected UBA1 mutation status as wtUBA1 (981 cells) and mtUBA1 monocytes (651 cells). h Dynamic changes of mtUBA1 and wtUBA1 monocytes ratios to all monocytes in VEXAS patients (log scale on y axis) along differentiation (pseudotime ordering from classical to non-classical monocytes on x axis). Data are presented as mean with 1.96*SE. i Frequency of CSFE⁺ bait cells among CytoTell Blue⁺CD14⁺ monocytes (VEXAS, n = 6; healthy donors, n = 6). Data are presented as mean with SD. P values were calculated using the two-sided unpaired Mann–Whitney U test. j Cytokine detection of IFN-γ, IL-6, and IL-1β in culture supernatants of purified CD14⁺ monocytes (VEXAS, n = 10; healthy donors, n = 10). Data are presented as mean with SD. P values were calculated using the two-sided unpaired Mann–Whitney U test. cMono classical monocytes, intMono intermediate monocytes, ncMono non-classical monocytes, cDC conventional DCs, pDC plasmacytoid DCs. Source data are provided as a Source Data file.

Fig. 4. Compositional and transcriptomic landscapes of mtUBA1 NK cells.

a A UMAP plot of 11,846 NK cells. b A neighborhood (Nhood) graph of NK cells,  which was generated similarly as in Fig. 2a. c Beeswarm and box plots of NK cells for VEXAS (n = 9) and healthy donors (n = 5), which were generated similarly as in Fig. 2b. d A UMAP plot of cells with projected UBA1 mutation status assignment for wtUBA1 (574 cells) and mtUBA1 NK cells (125 cells). e Dynamic changes of mtUBA1 and wtUBA1 NK cell ratios to all NK cells in VEXAS patients along differentiation. x axis, pseudotime ordering from CD56^bright NK cells to adaptive-like NK cells estimated by slingshot; y axis, ratios of cell numbers of mtUBA1 or wtUBA1 NK cells to all VEXAS NK cells on a log scale. Data are presented as mean with 1.96*SE. f A dot plot showing gene set enrichment scores of top 10 upregulated hallmark pathways across subtypes. Dot sizes indicate mean normalized enrichment score (NES) differences between VEXAS patients and healthy donors, and the color scale indicates FDR values. Non-significant pathways (FDR ≥ 0.10) are in gray. g A volcano plot of differentially expressed genes between mtUBA1 NK cells and wtUBA1 NK cells. Genes upregulated and downregulated in VEXAS patients are highlighted in red and blue, respectively. A horizontal dotted line, a padj value = 0.05; vertical dotted lines, absolute log2FC values = −0.25 and 0.25. h Gene module scores of the IFN-γ response, TNF response, leukocyte mediated cytotoxicity, and apoptosis pathways in mtUBA1 and wtUBA1NK cells. P values were calculated with the two-sided unpaired Mann–Whitney U test. Source data are provided as a Source Data file.

Fig. 5. Clonal expansion of effector CD8⁺ T cells with enhanced inflammation and cytotoxicity in VEXAS.

a A UMAP plot of 128,236 T cells. b A neighborhood (Nhood) graph of T cells, which was generated similarly as in Fig. 2a. c Beeswarm and box plots of T cells for VEXAS (n = 9) and healthy donors (n = 5), which were generated similarly as in Fig. 2b. d Distributions of T cell clone states in each subtype between VEXAS patients and healthy donors. e UMAP embedding of T cells from VEXAS patients (right) and healthy donors (left) colored by clonal expansion sizes. f Gini index and Shannon diversity index of TCR clonality in CD4⁺ T cells, CD8⁺ T cells, and total T cells from VEXAS (n = 9) and healthy donors (n = 5). Data are presented as mean with SD. P values were calculated with the two-sided unpaired Mann–Whitney U test. g Dynamic changes of gene module scores of the IFN-γ response, T-cell mediated cytotoxicity, and exhaustion pathways in CD8⁺ T cells from VEXAS patients and healthy donors along differentiation. x axis, pseudotime ordering from naïve CD8⁺ T cells to effector memory CD8⁺ T cells estimated by slingshot; y axis, gene module scores for each pathway. Data are presented as mean with 1.96*SE. gdT gamma delta T cells, MAIT mucosa-associated invariant T cell, proT proliferative T, Treg regulatory T cells, TCR T cell receptor. Source data are provided as a Source Data file.

Fig. 6. Skewed differentiation of B cells and clonal expansion of plasmablasts in VEXAS.

a A UMAP plot of 11,335 B cells. b A neighborhood (Nhood) graph of B cells, which was generated similarly as in Fig. 2a. c Beeswarm and box plots of each B cell subtype for VEXAS (n = 9) and healthy donors (n = 5), which were generated similarly as in Fig. 2b. d Proportions of B cell subtypes relative to the total number of B cells in VEXAS (n = 11, red dots) and healthy donors (n = 8, blue dots) by flow cytometry. Data are presented as mean with SD. P values were calculated with the two-sided unpaired Mann–Whitney U test. e UMAP plots of B cells from VEXAS patients (right) and healthy donors (left) colored by clonal expansion sizes. f Distributions of B cell clone states in each subtype in VEXAS patients and healthy donors. g Isotype distributions of immunoglobulins expressed by each transcriptionally defined B cell subpopulation within all subjects. h A bar chart showing gene set enrichment scores of top 10 upregulated hallmark pathways in plasmablasts from VEXAS patients. A horizontal axis displays normalized enrichment score (NES) differences between VEXAS patients and healthy donors. A color scale indicates FDR values. BCR B-cell receptor. Source data are provided as a Source Data file.