Contribution of germline and somatic mutations to risk of neuromyelitis optica spectrum disorder

Abstract

Neuromyelitis optica spectrum disorder (NMOSD) is a rare autoimmune disease characterized by optic neuritis and transverse myelitis, with an unclear genetic background. A genome-wide meta-analysis of NMOSD in Japanese individuals (240 patients and 50,578 controls) identified significant associations with the major histocompatibility complex region and a common variant close to CCR6 (rs12193698; p = 1.8 × 10^-8, odds ratio [OR] = 1.73). In single-cell RNA sequencing (scRNA-seq) analysis (25 patients and 101 controls), the CCR6 risk variant showed disease-specific expression quantitative trait loci effects in CD4⁺ T (CD4T) cell subsets. Furthermore, we detected somatic mosaic chromosomal alterations (mCAs) in various autoimmune diseases and found that mCAs increase the risk of NMOSD (OR = 3.37 for copy number alteration). In scRNA-seq data, CD4T cells with 21q loss, a recurrently observed somatic event in NMOSD, showed dysregulation of type I interferon-related genes. Our integrated study identified novel germline and somatic mutations associated with NMOSD pathogenesis.

Keywords: clonal hematopoiesis; eQTL; genome-wide association study; human leukocyte antigen; mosaic chromosomal alterations; neuromyelitis optica spectrum disorder; scRNA-seq; single-cell RNA sequencing; somatic mutation.

Graphical abstract

Figure 1

NMOSD GWAS meta-analysis of Japanese individuals (A) Manhattan plot of the NMOSD GWAS meta-analysis. The red horizontal line indicates the genome-wide significance threshold of p = 5.0 × 10⁻⁸. (B) Regional association plot around the lead variant, rs12193698, at CCR6. The purple diamond indicates the lead variant, rs12193698. Other dots represent SNPs and are colored according to linkage disequilibrium (r²) with the lead variant.

Figure 2

Regional association plots of HLA variants with NMOSD risk NMOSD GWAS meta-analysis of Japanese individuals (A) Regional association plots in nominal results (upper) and those conditioned on HLA-DRB1, -DQA1, and -DQB1 (lower). In each plot, diamonds represent the HLA variants, including SNPs, classical alleles, and amino acid polymorphisms of the tested HLA genes. The dashed horizontal lines represent the genome-wide significance threshold of p = 5.0 × 10⁻⁸. (B) Association plots in HLA-DRB1, -DQA1, and -DQB1 genes. Diamonds represent the amino acid polymorphisms (purple) and classical alleles (blue) for the tested HLA genes. For amino acid polymorphisms, the smallest p values among those for results of binary and omnibus tests at each position are shown. An allele or position of the smallest p values at each step is displayed in a red circle. The dashed horizontal lines represent the genome-wide significance threshold of p = 5.0 × 10⁻⁸.

Figure 3

Cell-type- and disease-specific expression profiles of CCR6 Regional association plots of HLA variants with NMOSD risk NMOSD GWAS meta-analysis of Japanese individuals. (A) Uniform manifold approximation and projection (UMAP) visualization of 1,004,361 PBMCs from the NMOSD cases (n = 25) and healthy controls (n = 101). (B) Cell-type-specific expression of CCR6. (C) Percentage of CCR6-expressing cells and CCR6 expression levels in each cell type. (D) The disease-specific eQTL effects of the CCR6 risk variant in five cell types. Boxes denote the interquartile range, and horizontal lines denote the median. Whiskers extend to 1.5 times the interquartile range. †FDR <0.1. (E) Correlations of the eQTL effect sizes of rs12193698 and 95% confidence intervals between NMOSD and healthy controls.

Figure 4

Identification of Th17 cells among CD4T cells and eQTL analysis using Th17 cells (A–C) Uniform manifold approximation and projection (UMAP) visualization of 372,313 CD4T cells projected onto the reference dataset: manually annotated (the same as Figure 3A) in (A), annotated based on the reference data in (B), and annotated to indicate whether Th17 or other cell types in (C). (D) Tile plot showing the proportion of the manually annotated four cell types in the reference-based 18 cell types. (E) Percentage of CCR6-expressing cells and CCR6 expression levels in the 18 reference-based cell types. (F) Beeswarm plot showing the distribution of adjusted log₂ fold change in abundance between NMOSD and controls in Nhoods according to Th17 cells and other cell types excluding Th17 cells. Boxes denote the interquartile range, and vertical lines denote the median. Whiskers extend to 1.5 times the interquartile range. (G) Correlations of the eQTL effect sizes of rs12193698 and 95% confidence intervals between NMOSD and healthy controls in CD4 memory cells. The plot includes all CD4 memory cells (black), Th17 cells in CD4 memory cells (red), and CD4 memory cells excluding Th17 cells (light blue).

Figure 5

Integrated analysis into the association between somatic mutations and NMOSD (A) Forest plot of the associations between disease risks and mCAs (CNA and CN-LOH). The odds ratios and 95% confidence intervals were calculated using a multivariate logistic regression adjusted for age and sex. (B) Schematic overview of our integrated analysis of genomics and single-cell transcriptomics. (C) Uniform manifold approximation and projection (UMAP) embeddings of mNMOSD1 scRNA-seq data colored by cell type (left) and clone (right). (D) UMAP embeddings of mNMOSD1 scRNA-seq data colored by the density of the 21q loss cells (left) and normal cells (right). (E) Heatmap showing in-sample odds ratios of each cell type containing mutated cells. Dendritic cells (DCs) in mNMOSD4 are masked in gray due to a limited cell count (<5). (F) DEG analysis between the 21q loss and normal CD4T cells in mNMOSD1 and mNMOSD2. Significant DEGs satisfying FDR (adjusted p values via the Benjamini-Hochberg method) < 0.1 are colored in light blue or pink, and DEGs on 21q are colored in navy. The upward triangle denotes a data point beyond the axis limit. The genes included in the “response to virus” pathway were annotated. (G) Top 10 enriched biological pathways of the downregulated DEGs in CD4T cells with 21q loss. The dot color indicates the adjusted enrichment p values via the Benjamini-Hochberg method, and the dot size represents the gene count annotated to each term.