Мультиэнштейная геномная ассоциация и мультиомикальный анализ позволяют выявить специфичные особенности многоклерозной генетики.

Main

Multiple sclerosis (MS), the most prevalent chronic inflammatory disease of the central nervous system (CNS), is characterized by demyelination disseminated in space and time, resulting in a range of neurological symptoms 1 . Globally, an estimated 2.8 million individuals are affected by MS, with prevalence varying geographically; for example, the prevalence among people of European ancestry (EUR) is over 10 times higher than among those of East Asian ancestry 2 , 3 . These regional differences suggest that both genetic and environmental factors contribute to the pathogenesis of MS 4 . However, the genetic background of MS and differences between ancestral populations remain incompletely understood 5 .

To date, genome-wide association studies (GWAS) have detected more than 200 associations with MS genetic susceptibility 6 – 10 . Nevertheless, MS GWAS for non-EUR populations have been limited to a few biobanks 11 , 12 , and there have been none among East Asian individuals, partially owing to the low prevalence in this population. Thus, at present, the only clearly identified risk variants for MS in East Asian individuals are within the major histocompatibility complex (MHC) region 13 , 14 . Notably, recent studies have shown that performing GWAS meta-analysis across populations and biobanks can increase statistical power and reveal heterogeneity 15 . Thus, large-scale cross-population GWAS meta-analyses for MS may help to advance our understanding of its pathogenesis.

Previous MS GWAS results have prioritized cell types in which MS risk variants exhibit functional roles and have highlighted the importance of peripheral immune cells, including T cells, B cells and microglia (that is, immune cells specifically found in the CNS) 10 , 16 . However, these studies mainly utilized bulk transcriptome and epigenome datasets, and the roles of nonimmune cell types in susceptibility to MS have remained insufficiently explored. Single-cell and single-nucleus RNA sequencing (sc/snRNA-seq) is useful to address these issues but has not been fully utilized to interpret MS GWAS results. Furthermore, recent advances in spatial transcriptomics (ST) technologies have shown considerable potential for capture of cellular pathogenesis through consideration of the two- or three-dimensional structures of tissues 17 and seem to be well-suited for study of MS, which is characterized by spatial heterogeneity.

Here we performed a non-EUR MS GWAS recruiting a total of 688 cases and 205,199 controls from the Japanese population and identified a novel Japanese-specific risk locus located upstream of ETS1 . To maximize statistical power, we used a total of 29,374 cases and 1,843,563 controls from 4 ancestral populations (Japanese ancestry, EUR, African ancestry (AFR) and American ancestry (AMR)), performed single-population GWAS meta-analyses for the EUR and AFR populations as well as cross-population GWAS meta-analysis and identified a total of 22 novel loci. Next, we performed fine-mapping of human leukocyte antigen (HLA) variants in the Japanese population. Furthermore, by integrating GWAS and sc/snRNA-seq data, we detected cell types associated with MS among peripheral immune cells and brain cells. In addition, spatial and temporal heterogeneity in associations with MS genetic background was observed in ST datasets for MS subcortical lesions.

Results

Meta-analysis of Japanese MS GWAS

We conducted a meta-analysis of two Japanese MS GWAS datasets. For the first GWAS (dataset 1), we enrolled MS cases mainly from the Japan MS/NMOSD Biobank 14 , 18 and the University of Osaka Hospital, and included controls from the BioBank Japan (BBJ) second cohort 19 . For the second GWAS (dataset 2), we enrolled cases and controls from the BBJ first cohort 19 . As the BBJ is a multi-institutional hospital-based biobank in Japan 19 , we retrieved MS cases based on ICD-10 codes from the BBJ. In total, 688 cases and 205,199 controls were analyzed, involving 8,858,013 variants that fulfilled the postimputation criteria (minor allele frequency (MAF) > 0.5% and imputation quality R 2 > 0.7). The quantile–quantile plot and genomic inflation factor (of 1.026) indicated little genomic inflation (Extended Data Fig. 1a ).

We identified two loci that satisfied the genome-wide significance threshold of P < 5.0 × 10 −8 in the Japanese MS GWAS meta-analysis (Fig. 1a and Table 1 ). The lead variant of the most significant locus was identified within the MHC region ( rs80079921 , P = 4.9 × 10 −29 , odds ratio (OR) = 2.71). The other was a risk variant at 11q24 ( rs2199759 , P = 2.1 × 10 −8 , OR = 1.58) that had not been previously reported in any populations and was located within an intergenic region (670 kb upstream of KIRREL3 and 914 kb downstream of ETS1 ; Extended Data Fig. 1b ). ETS1 , the nearest coding gene downstream of rs2199759 , encodes a transcription factor that has an essential role in differentiation of lymphoid lineages 20 and has been genetically and biologically associated with several immune-related traits and autoimmune diseases, including MS 10 , 21 . In addition, downregulation of ETS1 in brain endothelial cells has been reported to cause blood–brain barrier (BBB) dysfunction 22 , a key pathological feature of MS. These facts suggest that this locus may act as a distal enhancer-like regulator of ETS1 and contribute to MS pathogenesis through immune dysregulation and BBB impairment.

EUR, AFR and cross-population MS GWAS meta-analyses

To maximize statistical power, we systematically collected and processed GWAS summary statistics from our Japanese population study, other studies 9 , 10 , the UK Biobank 11 , MVP 12 , FinnGen 23 , All of Us 24 and Genomics England 25 and conducted a cross-population meta-analysis involving four ancestral populations: Japanese ancestry, EUR, AFR and AMR (Supplementary Table 1 ). We also performed single-population meta-analyses for which multiple datasets were available (EUR and AFR). The sample sizes were 27,572 cases and 1,436,801 controls for EUR; 819 cases and 155,904 controls for AFR; and 29,374 cases and 1,843,563 controls for the cross-population meta-analysis.

In the EUR and AFR GWAS meta-analyses, we identified 153 and 2 significant loci, respectively, including 17 and 2 novel loci (Fig. 1b,c , Table 2 and Supplementary Tables 2 and 3 ). In the cross-population GWAS meta-analysis, we identified 156 significant loci, including three novel loci that were not significant in the single-population meta-analyses (Fig. 1d , Table 2 and Supplementary Table 4 ). Several novel variants were located near genes with experimentally established associations with MS, such as rs41287361 ( SEMA4D ) 26 , 27 and rs142239370 ( PRF1 ) 28 . SEMA4D encodes semaphorin 4D, which acts as an axonal guidance factor during neuronal development and also has a crucial role in immune responses, and Sema4d knockout mice are resistant to development of experimental autoimmune encephalomyelitis 26 , 27 . Notably, PRF1 was reported as a male-specific MS-associated gene 28 , and rs142239370 was nominally significant ( P < 0.05) only in the MVP EUR and AFR cohorts, of which 90% of participants were males 12 .

Next, to investigate sex differences in the genetic architecture of MS, we performed sex-stratified GWAS in cohorts with available individual-level genotype data (Japanese datasets 1 and 2, UK Biobank, All of Us and Genomics England) and conducted sex-stratified GWAS meta-analyses in the Japanese and EUR populations. For rs142239370, which was nominally significant only in the MVP cohorts, the male-specific effect did not reach significance ( P male = 0.768, P female = 0.652), suggesting instead a potential environment-specific association. We systematically evaluated the sex-specificity of lead variants identified in the Japanese and EUR populations using Cochran’s Q test (Extended Data Fig. 2a and Supplementary Table 5 ). Among 145 lead variants identified in the Japanese or EUR meta-analysis that were subsequently tested in the sex-stratified GWAS, 11 were nominally significant ( P < 0.05), but none remained significant after Bonferroni correction ( P < 0.05/145).

The largest EUR MS GWAS conducted by the International Multiple Sclerosis Genetics Consortium (IMSGC) 10 to date included 47,429 cases and 68,374 controls, of which 14,802 cases and 26,703 controls in the discovery cohort underwent genome-wide analysis (we hereafter refer to this discovery cohort as IMSGC2019). To investigate factors contributing to detection of associations, we compared the effect sizes of novel variants identified in our meta-analyses with those reported by IMSGC2019 (Extended Data Fig. 2b and Supplementary Table 6 ). The effect sizes for the 17 novel variants in our EUR meta-analysis were concordant with IMSGC2019, indicating that the detection of these variants relied primarily on our larger sample size. In our AFR meta-analysis, one novel variant (rs191128996) was AFR-specific (allele frequency (AF) in the 1000 Genomes Project (AF 1KG) AFR = 0.019, AF 1KG EUR = 0.001), whereas the other ( rs11627878 ) had sufficient AF in EUR (AF AFR = 0.036, AF EUR = 0.339) but was not nominally significant in IMSGC2019 ( P IMSGC2019 = 0.772). In our cross-population meta-analysis, one novel variant ( rs11771510 ) showed a concordant effect, one (rs142239370) was nominally significant only in the MVP cohorts but not in IMSGC2019 as mentioned above, and one ( rs4848520 ) was not included in IMSGC2019.

Finally, we examined the differences in effect sizes between the GWAS results for the Japanese ancestry and other populations. Of the 125 loci significant in the EUR population that were tested in the Japanese ancestry GWAS, 24 were nominally significant ( P < 0.05) in the Japanese population, and 11 showed nominally significant heterogeneity (Extended Data Fig. 2c and Supplementary Table 7 ). The effect sizes of these loci showed a significant positive correlation between the Japanese and EUR populations (Spearman’s ρ = 0.67, P = 6.3 × 10 −18 ). By contrast, the association with rs2199759 , a novel variant detected in the Japanese GWAS, was not observed in the other populations. This signal was distinct from the MS-associated locus located within the ETS1 gene that was identified in previous studies and in our meta-analyses for other populations (Extended Data Fig. 3 ). These findings indicate that there are both shared and population-specific factors in the genetic background of MS susceptibility.

In summary, our results indicate that some associations are difficult to detect in a single cohort for several reasons: (1) limited sample size, (2) population-specific variants (that is, low AF in other populations), (3) population-specific causality (for example, gene–environment interactions and gene–gene interactions), (4) sex-specific effects or (5) variant filtering issues (for example, imputation quality). By leveraging multiple populations and biobank datasets, we were able to overcome these limitations in our meta-analyses and identify additional associations.

HLA fine-mapping analysis emphasized the associations between MS and class II HLA genes

To fine-map risk variants within MHC for MS susceptibility, we performed HLA imputation on our Japanese genotype data using an HLA reference panel of the Japanese population 29 . After imputation and quality control (QC) (MAF > 0.5%, r 2 by DEEP*HLA > 0.7 and classical HLA genes), we obtained 53 two-digit and 87 four-digit HLA alleles and 603 HLA amino acid polymorphisms. The most significant association was for phenylalanine at HLA-DQβ1 amino acid position 9 ( P = 4.9 × 10 −37 ; Fig. 2a , Table 3 and Supplementary Table 8 ). By stepwise conditional analysis, we successively detected HLA-DQA1*03:01, HLA-DRB1*15:01 and glycine at HLA-DQβ1 position 125 as independent associations (Fig. 2b–d and Table 3 ). After conditioning on these four associations, no HLA variants satisfied the genome-wide significant threshold, whereas the single nucleotide polymorphism (SNP) associations remained significant ( rs9266359 , P = 2.6 × 10 −8 ; Fig. 2e ). To address the potentially insufficient conditioning process owing to deviations from additivity, we performed conditional analyses incorporating nonadditive effects of the conditional terms; these identified the same four associations (Extended Data Fig. 4 and Supplementary Table 9 ). Furthermore, we performed nonadditive analyses for the four independent associations identified in our HLA fine-mapping and found that glycine at HLA-DQβ1 position 125 exhibited a nominally significant nonadditive effect ( P = 0.0016; Supplementary Table 10 ).

In addition, to compare the relative contributions of HLA alleles and amino acid polymorphisms, we analyzed them separately (HLA only and amino acid only). In the HLA only analyses, HLA-DRB1*15:01, HLA-DRB1*04 and HLA-DRB1*08 were identified as independent associations (Extended Data Fig. 5a –d and Supplementary Table 11 ). In the amino acid only analyses, phenylalanine at HLA-DQβ1 position 9, alanine at HLA-DQβ1 position 57 and aspartic acid at HLA-DQβ1 position −21 were identified as independent associations (Extended Data Fig. 5e –h and Supplementary Table 12 ). HLA-DQA1*03:01 and alanine at HLA-DQβ1 amino acid position 57 were in high linkage disequilibrium (LD) ( r 2 = 0.877) (Extended Data Fig. 6 ), whereas the other pairs showed at most moderate LD ( r 2 < 0.6). Therefore, the four associations detected in the original analysis were not fully conditioned by other variants in the HLA only and amino acid only analyses, suggesting that they probably represent independent associations. We evaluated model fitting for each model; however, it was difficult to statistically determine the best-fitting model (Supplementary Fig. 1 ).

Our HLA fine-mapping for MS susceptibility used the largest sample sizes and the most comprehensive HLA variants in the Japanese population to date. We identified four HLA risk variants that were localized to class II HLA genes. These findings were consistent with those of previous studies 30 , 31 , confirming the crucial role of class II HLA genes in MS. HLA-DRB1*15:01 is considered to be among the most common MS risk variant in the MHC region across ancestries 14 , 30 , 31 , whereas phenylalanine at HLA-DQβ1 position 9 has been reported only by Ogawa et al. 13 (Supplementary Table 13 ). Therefore, our confirmation of this association in a larger Japanese cohort will be valuable for future investigation of MS risk in the Japanese population. Furthermore, although HLA-DRB1*04:05 has been reported as an East Asian-specific MS risk variant 13 , 14 , although this allele was initially significant ( P = 2.5 × 10 −8 ) in our main HLA analysis, its association did not reach genome-wide significance after conditioning on HLA-DQβ1 amino acid position 9 ( P = 0.0012). Instead, HLA-DRB1*04 was detected only in the model that included HLA alleles alone. These results suggest that the association with DRB1*04:05 could be largely explained by HLA-DQβ1 amino acid position 9. In summary, our analysis highlights the importance of high-resolution HLA fine-mapping including amino acid polymorphisms.

Cell type specificity analysis using sc/snRNA-seq data from PBMCs and MS brain lesions

The majority of studies exploring functional annotation of MS-associated variants have used bulk transcriptome and epigenome information 9 , 10 , 16 . Recently, single-cell technologies have enabled interpretation of disease pathology with higher resolution 18 , 32 . As immune cells have a key biological role in the pathogenesis of MS as well as that of other autoimmune diseases 33 , we performed scRNA-seq of peripheral blood mononuclear cells (PBMCs) from 20 patients with MS and obtained 163,541 cells. We manually annotated these cells based on expression of known marker genes and defined 8 major cell types and 29 finer cell types, according to Azimuth-based annotations 34 (Extended Data Fig. 7 ).

To detect cell types associated with MS susceptibility at single-cell resolution, we applied scDRS 35 to the cross-population GWAS meta-analysis mentioned above and our PBMC scRNA-seq data (Fig. 3a–c ). scDRS identifies disease-associated cell types by detecting cells that exhibit excess expression of disease-associated genes, as defined by MAGMA 36 , compared with background cells. Among the eight main cell types, CD4 + T cells showed significant association with MS ( P = 0.0020). Further investigations of the finer CD4 + T cell types detected significant associations of naive CD4 + T cells ( P = 0.015), central memory T cells ( P = 0.0020), effector memory T cells ( P = 0.010) and memory regulatory T cells (memory T reg cells; P = 0.0040). As the CD4 + T helper cell lineage has a central role in MS pathogenesis 37 , it was reasonable that naive, central memory and effector memory CD4 + T cells would show strong associations. Although the role of T reg cells in the pathogenesis of MS has remained unclear, some studies have suggested heterogeneous functions of these cells in MS 38 . Our results, which showed a stronger association of memory T reg cells than naive T reg cells, support this possibility. The significance of CD4 + T cells and their subtypes was replicated in a large-scale PBMC scRNA-seq dataset from the Asian Immune Diversity Atlas (AIDA) 39 , suggesting that the associations between MS susceptibility and these cell types are robust (Supplementary Fig. 2 ).

As well as immune cells, the brain, which is a critical organ in MS, is important for elucidation of the pathogenesis of MS. To detect MS-associated cell types in the brain, we applied scDRS to snRNA-seq data for MS subcortical lesions (74,686 cells) 40 . T cells ( P = 0.0010), myeloid cells ( P = 0.0010), endothelial cells ( P = 0.0040) and B cells ( P = 0.0070) showed significant associations (Fig. 3c,d ). After we had excluded T cells, B cells and myeloid cells from the snRNA-seq data to enable us to focus on nonimmune brain cell types, endothelial cells remained significant ( P = 0.0010), and astrocytes and stromal cells also showed significant associations ( P = 0.0090 and P = 0.013, respectively) (Supplementary Fig. 3 ). Endothelial cells are the main component of the BBB, the disruption of which is a contributor to the development of MS 41 , 42 . Astrocytes have recently been identified as early and highly active players in lesion formation 43 , 44 . Thus, it is unsurprising that these cell types would show strong associations with MS susceptibility. Overall, our results suggest that in addition to traditionally implicated cell types (that is, immune cells and microglia), endothelial cells and astrocytes may have substantial roles in determining susceptibility to MS.

Integration of GWAS with ST highlights spatial and temporal heterogeneity in MS brain lesions

We detected brain cell types associated with MS genetic risk at single-cell resolution. snRNA-seq cannot capture spatial information, but the recent development of ST methods has enabled profiling of cell-specific gene expression levels at specific spatial positions 17 . As MS lesions in the CNS are characterized by spatial spread 1 , brain ST analysis could help us to understand MS pathogenesis. Furthermore, as MS lesions change over time into chronic active lesions and chronic inactive lesions (hereafter referred to as ‘active’ and ‘inactive’) 45 , it is important to compare ST data at different time points.

To investigate the association between the MS genetic architecture and the spatial regions of MS brain lesions, we applied gsMap 46 to ST data of eight active and four inactive MS subcortical lesions 40 . gsMap can integrate ST data with GWAS results for spatially resolved mapping of trait-associated cells or bins. We uniformly analyzed all ST slides using the joint analysis mode in gsMap, which is applicable when multiple ST slides are available and allows fair comparisons across them. As gsMap was designed for single-population GWAS, we applied it to the EUR meta-analysis GWAS results. We used the six niches defined in the original study by gene expression and deconvoluted cell type proportions as spatial annotations; these niches comprised lesion rim (LR), lesion core (LC), vascular infiltrating (VI), periplaque white matter, ependyma and gray matter cells.

In the active lesion, the VI, LR, LC and periplaque white matter niches showed significant associations ( P = 2.9 × 10 −12 , 5.0 × 10 −11 , 5.2 × 10 −7 and 3.9 × 10 −6 , respectively; Fig. 3e,f and Extended Data Fig. 8a ). Similarly, in the inactive lesion, these niches showed significant associations ( P = 2.0 × 10 −9 , 4.5 × 10 −5 , 0.0010 and 0.019, respectively), but these associations were weaker than in the active lesion. The active lesion-specific pattern was particularly pronounced for VI and LR (Fig. 3f,g ). The VI niche is primarily located within the LC one and is characterized by endothelial cells, stromal cells and immune cells (Extended Data Fig. 8b ). It thus resembles lesion-associated perivascular spaces, and it is therefore reasonable that MS genetic risk would be enriched in the VI niche, as such cell types reflect a fundamental pathological process of MS, namely immune cell infiltration into the CNS. The LR was characterized by myeloid cells, consistent with mediation of demyelination by macrophages and monocytes in the development of MS.

To evaluate whether these results were specific to MS, we conducted gsMap analyses using GWAS data for other immune-mediated traits (rheumatoid arthritis 47 , systemic lupus erythematosus 48 , ulcerative colitis 49 , Crohn’s disease 49 and allergy 50 ) and control traits (Alzheimer’s disease 51 , depression 52 and low-density lipoprotein 53 ). Associations with VI and LR in active lesions were notably significant for MS, indicating that gsMap could capture MS-specific spatial features in MS lesions (Fig. 3g and Extended Data Fig. 8c ). We further compared the strength of associations detected by gsMap with deconvoluted cell type proportions. We observed positive correlations for some cell types, including VI–stromal cells ( ρ = 0.43, P = 2.9 × 10 −12 ; Extended Data Fig. 9a ) and LR–myeloid cells ( ρ = 0.69, P < 1.0 × 10 −300 ; Extended Data Fig. 9b ); however, these associations were not driven solely by any single cell type (Extended Data Fig. 9 ). These findings suggest that the associations identified by gsMap are determined by complex combinations of cell type proportions per spot and cell–cell spatial interactions, rather than by single-cell expression profiles alone.

Collectively, our findings demonstrate the spatial and temporal heterogeneity of MS brain lesions through an approach integrating genetics and ST.

Discussion

In this study, we performed a GWAS of MS in the Japanese population and identified a population-specific risk variant. ETS1 , the nearest coding gene downstream of this variant, encodes a transcription factor that plays an essential part in the differentiation of lymphoid lineages 20 . Next, we curated GWAS data from existing studies and global biobanks, conducted EUR, AFR and cross-population GWAS meta-analyses, and identified a total of 22 novel loci, including loci that were difficult to detect in a single cohort. Our results suggest that these associations may not have been detected in previous studies owing to smaller sample sizes, population-specific variants, population-specific effects, sex-specific effects and variant filtering. We propose that GWAS meta-analyses leveraging multiple populations and biobank datasets could help to overcome these limitations. As MS susceptibility is influenced by environmental exposures such as Epstein–Barr virus infection and latitude 54 , gene–environment interaction may represent a critical factor contributing to differences in MS genetic architecture between populations or cohorts. In the HLA fine-mapping, we identified four independent class II HLA risk alleles in the Japanese population. The three risk alleles other than HLA-DRB1*15:01, the most common MS risk HLA allele, have not been identified in EUR populations 10 , 30 , 31 .

We constructed a large PBMC scRNA-seq dataset from patients with MS. By integrating our GWAS results with sc/snRNA-seq datasets for PBMCs and subcortical lesions, we detected enrichment of genome-wide host genetic risk in specific cell types. In addition to the CD4 + T helper cell lineage, one of the traditionally implicated cell types 37 , memory T reg cells, endothelial cells, astrocytes and stromal cells showed significant associations with MS genetic risk. Although the role of T reg cells in the pathogenesis of MS has remained unclear, some studies have suggested heterogeneous functions of T reg cells in MS 38 . Moreover, the enrichment of MS genetic risk in endothelial cells, a key component of the BBB 41 , 42 , and that of astrocytes, which are early and highly active players in lesion formation 43 , 44 , are plausible findings.

In addition, we integrated GWAS and MS subcortical lesion ST data and described spatial and temporal heterogeneity of MS brain lesions from a genetics perspective. Chronic inflammation within the CNS is strongly associated with increased disability in MS. This condition, known as ‘smoldering’, is thought to be relatively resistant to disease-modifying treatments 55 . In our analysis, genetic risk factors for MS were strongly enriched in smoldering lesions. This suggests that genetic risk factors for MS may also have an impact on disease progression in the brain. This analytical approach to understanding the spatiocellular features of disease genetics at single-cell spatial resolution by integrating MS GWAS with ST datasets of MS brain lesions is expected to become an increasingly important element of GWAS.

Our study had several potential limitations. First, our GWAS meta-analyses primarily relied on biobank data, and part of the sample collection was based on ICD-10 codes, which may have introduced uncertainty in the medical diagnoses. A framework for accurate collection of data from patients with MS on a biobank scale needs to be developed. Second, we included MS cases regardless of disease stage or clinical subtype to increase sample size. To better understand genetic contributions to MS pathogenesis and progression, more clinical data from patients with MS should be collected for use in within-case GWAS. Third, ST data for brain are scarce, and analysis methods have not been standardized. Future expansion of ST data will resolve these problems. Fourth, the detection of cell types and regions related to MS using scDRS and gsMap used a polygenic approach; identification of the specific genes and pathways that drive the observed associations remains a challenge for future studies.

In summary, we conducted large-scale MS GWAS meta-analyses using multiple populations and integrative analyses with sc/snRNA-seq and ST data. Our results highlight the importance of studying disease genetic background across diverse populations and omics modalities.

Methods

Study participants and genotyping in the Japanese datasets

For the Japanese MS GWAS, we constructed two case–control genotype datasets (hereafter referred to as ‘dataset 1’ and ‘dataset 2’). Dataset 1 consisted of samples obtained from the Japan MS/NMOSD Biobank, the University of Osaka Hospital and the BBJ second cohort 19 . Dataset 2 comprised samples exclusively from the BBJ first cohort 19 . MS cases in the Japan MS/NMOSD Biobank and the University of Osaka Hospital were diagnosed using the 2010 Revised McDonald criteria 56 , whereas those in the BBJ were identified based on ICD-10 code (G35). The remaining BBJ samples were used as controls. All participants provided written informed consent, and approval was given by the ethics committees of the University of Osaka and the University of Tokyo.

The samples in dataset 1 were genotyped using an Illumina Asian Screening Array chip, whereas those in dataset 2 were genotyped with an Illumina HumanOmniExpressExome BeadChip or a combination of the Illumina HumanOmniExpress and HumanExome BeadChips.

QC and imputation in the Japanese datasets

We applied the same stringent QC to both Japanese datasets. We excluded samples with a genotyping call rate <0.98. We included only samples of estimated East Asian ancestry, based on principal component analysis (PCA) with the HapMap project samples 57 . Duplicate samples were removed based on genotype information (KING kinship coefficient 58 , as calculated by PLINK2 59 , greater than 0.354). We excluded controls with autoimmune or allergic diseases. After QC filtering of samples, dataset 1 included 627 cases and 52,019 controls, and dataset 2 included 61 cases and 153,180 controls.

We applied standard QC filters to variants, excluding those with SNP call rate < 0.99, minor allele count < 5, P < 1.0 × 10 −6 (1.0 × 10 −10 for dataset 2) for Hardy–Weinberg equilibrium (HWE), or >5% AF difference compared to the Tohoku Medical Megabank Project dataset 60 and in-house reference panel described below.

We conducted genome-wide genotype imputation to estimate untyped variants. We used an in-house reference panel consisting of Japanese whole-genome sequencing data ( n = 11,754) as a haplotype reference for genotype imputation and performed haplotype estimation using SHAPEIT4 61 , followed by genotype imputation with Minimac4 62 .

GWAS of the Japanese datasets

We conducted GWAS using a generalized linear mixed model (GLMM) with random effects implemented in SAIGE (v.0.45) 63 for the Japanese datasets. We included age, sex and the top ten principal components as covariates. We used variants with MAF > 0.5% and imputation quality R 2 > 0.7. We set the genome-wide significance threshold to P < 5.0 × 10 −8 .

UK Biobank

UK Biobank is a large-scale biomedical database and research resource containing genetic, lifestyle and health information from half a million UK participants. Detailed characteristics of the cohort and genotype–phenotype data have been described elsewhere 11 . Briefly, we used genomic data based on genotyping by either the Applied Biosystems UK BiLEVE Axiom Array or Applied Biosystems UK Biobank Axiom Array and imputation using a Genomics England haplotype reference panel 64 . We identified MS cases based on ICD-10 code (G35). We excluded controls with autoimmune or allergic diseases. Individuals were included based on the following criteria: (1) concordance with sex information, (2) genotyping call rate > 0.95, (3) heterozygosity rate < +3 s.d. and (4) the white British subset previously reported 10 . After QC, 1,950 cases and 376,163 controls remained for association analysis. Variant QC was conducted based on the following criteria: (1) missing data < 0.05, (2) MAF > 0.5%, (3) P > 1 × 10 −15 for HWE and (4) imputation quality R 2 > 0.7. We performed genome-wide association tests using a GLMM as implemented in REGENIE (v.3.3) 65 . We included age, sex and the top ten principal components as covariates.

All of Us

The All of Us Research Program was established in 2015 to enroll a diverse group of at least one million people across the USA for biomedical research and human health improvement 24 . We identified MS cases based on ICD-10 codes (G35). We excluded controls with autoimmune or allergic diseases. Individuals were included based on the following criteria: (1) genotyping call rate > 0.98, (2) heterozygosity rate < +3 s.d. and (3) EUR, AFR and AMR individuals identified by PCA-based criteria. After QC, 1,260 cases and 99,143 controls for EUR, 295 cases and 45,659 controls for AFR, and 182 cases and 34,839 controls for AMR remained for association analysis. Variant QC was conducted based on the following criteria: (1) missing data < 0.03, (2) MAF > 1% and (3) P > 1 × 10 −15 for HWE. We performed genome-wide association tests using a GLMM as implemented in REGENIE (v.3.2.4) 65 . We applied microarray based genotyped data for step 1 and whole-genome sequencing data for step 2 of REGENIE. Age, sex and the top ten principal components were included as covariates.

Genomics England

Genomics England was established by the UK Department of Health and Social Care to conduct the 100,000 Genomes Project, focusing on patients with rare diseases and cancer 25 . We used whole-genome sequencing data of aggregated variant calls ( https://re-docs.genomicsengland.co.uk/site_qc/ ), including 78,195 post-QC germline genomes. We identified MS cases based on ICD-10 codes. Controls included patients with cancer or unaffected relatives of rare diseases without autoimmune or allergic diseases. Additional variant QC was conducted based on the following criteria: (1) Missing data < 0.02, (2) MAF > 0.5% and (3) P > 1 × 10 −6 for HWE. We included EUR individuals, identified by PCA-based criteria. Then, 244 cases and 28,846 controls remained for association analysis. We performed genome-wide association tests using a GLMM as implemented in SAIGE (v.1.0.7) 63 . We included age, sex and the top ten PCs as covariates.

GWAS meta-analysis

For the GWAS meta-analyses, we also collected GWAS summary statistics from previous studies and biobank datasets (Supplementary Table 1 ). To align the human genome reference build across cohorts with compatibility with the previous GWAS, we converted GWAS summary statistics in hg38 to hg19 using UCSC LiftOver 66 . We ensured that all summary statistics were aligned to the same strand of the reference genome by calculating the consensus rate across genome-wide nonpalindromic variants, following the procedure implemented in gwas-sumstats-harmoniser in the GWAS catalog 67 . We included variants with MAF > 0.5% and imputation quality R 2 > 0.7 in each dataset and used only those shared by at least two datasets. We applied a standard fixed-effects approach implemented in METAL (version released on 2011-03-25) 68 for GWAS meta-analyses. For studies with overlapping samples, we employed the Lin–Sullivan method 69 in RE2C (v.1.0.6) 70 . We used FUMA (v.1.6.1) 71 to define independent associated loci with default parameters. We considered a locus to be novel if the genomic region within ±250 kb of the lead variant did not include variants that had been previously reported.

Sex-stratified GWAS meta-analysis

Sex-stratified GWAS meta-analyses were performed for the Japanese ancestry and EUR populations. For the EUR analysis, cohorts with available individual-level genotype data (UK Biobank, All of Us and Genomics England) were included. QC filtering was conducted as described above for each cohort. Within each cohort, GWAS was performed separately for male and female individuals using the same covariates as in the overall GWAS (excluding sex), and the results were subsequently meta-analyzed within each sex. In the Japanese analyses, 201 male cases and 112,285 male controls, and 487 female cases and 91,294 female controls were included. In the EUR analyses, 886 male cases and 231,229 male controls, and 2,568 female cases and 272,652 female controls were included. The heterogeneity of effect sizes for variants identified in the overall Japanese ancestry and EUR meta-analyses was evaluated using Cochran’s Q test.

HLA imputation and association study

We performed HLA imputation with the HLA reference panel of the Japanese population ( n = 1,118) 29 . We applied DEEP∗HLA 72 , a method for deep learning-based HLA imputation, to the genotyped SNPs in the MHC region and imputed the two- and four-digit HLA alleles and amino acid variants of class I and II HLA genes ( HLA-A , HLA-C , HLA-B , HLA-DRB1 , HLA-DQA1 , HLA-DQB1 , HLA-DPA1 and HLA-DPB1 ). The variants imputed with an imputation quality score of r 2 > 0.7 in 10-fold cross-validation and MAF > 0.5% were used for further analyses.

For HLA fine-mapping, we evaluated associations of HLA variants with risk of MS using a logistic regression model with a glm() function implemented in R (v.4.3.1). Accounting for the employment of a logistic regression model, we excluded third-degree or more closely related individuals (KING kinship coefficient 58 , as calculated by PLINK2 59 , greater than 0.0884) from dataset 1 for the Japanese GWAS. We estimated the significance of multiallelic amino acid polymorphisms using an omnibus test for each amino acid position through log-likelihood ratio test, comparing the likelihood of the fitted model with the null model. The significance of the improvement in model fitting was evaluated based on the deviance, which follows a chi-squared distribution with m − 1 degrees of freedom for an amino acid position with m polymorphic residues. To evaluate independent risk among and across the HLA genes, we conducted a forward-type stepwise conditional regression analysis that also included the associated variant genotypes as covariates. When conditioning on a specific HLA amino acid position, we included all the amino acid polymorphisms belonging to that position. We otherwise used the same covariates as those used in the GWAS. We adopted a genome-wide association significance threshold of P = 5.0 × 10 −8 .

We also performed conditional analyses incorporating nonadditive effects of the conditional terms. We defined the nonadditive term δ i , j as 1 − |1 − x i , j |, where x i , j is the dosage of allele j in individual i . We included both additive and nonadditive effects as conditioning terms, while retaining additive effects as the main effect terms. All other covariates and testing methods were kept the same as in our original analysis. Further, we performed nonadditive analyses for the four independent associations identified in our HLA fine-mapping. At each conditioning step of our original analysis, we compared the likelihoods of a model with only the additive effect versus one with both additive and nonadditive effects for the main effects.

We evaluated model fit for three models (combined HLA and amino acid, HLA only and amino acid only) using two approaches. First, we compared the Akaike information criterion and Bayesian information criterion across conditional models. Second, we performed ten-fold cross-validation and compared predictive performance using the area under the receiver operating characteristic curve (AUC). As the individuals involved in MS cases were younger than the controls, the nominal model including covariates (age, sex and 10 genetic principal components) achieved a high AUC (>0.9). To avoid this covariate-driven inflation, we also evaluated AUC for models excluding these covariates.

Study participants and specimen collection of PBMCs for scRNA-seq

Peripheral blood samples were obtained from patients with MS ( n = 20) recruited at The University of Osaka Hospital. We collected peripheral blood into heparin tubes and isolated PBMCs using Leucosep (Greiner Bio-One) density gradient centrifugation according to the manufacturer’s instructions. We processed blood samples within 3 h of collection and stored them at −80 °C until use.

Droplet-based single-cell sequencing

We processed single-cell suspensions through a 10x Genomics Chromium Controller. Oil droplets of encapsulated single cells and barcoded beads (GEMs) were subsequently reverse transcribed in a Veriti Thermal Cycler, resulting in complementary DNA (cDNA) tagged with a cell barcode and unique molecular index. cDNA was amplified to generate single-cell libraries according to the manufacturer’s protocol. Quantification was performed using a Bioanalyzer High Sensitivity DNA assay (High Sensitivity DNA Kit, Agilent, 5067–4626). Amplified cDNA was enzymatically fragmented, end-repaired and poly-A tagged. Cleanup and size selection were performed on amplified cDNA using SPRIselect magnetic beads (Beckman-Coulter, SPRIselect, B23317 ). Next, Illumina sequencing adapters were ligated to the fragments and cleaned up using SPRIselect magnetic beads. Then, sample indices were selected and amplified, followed by double-sided size selection using SPRIselect magnetic beads. We assessed the final library quality using a Bioanalyzer High Sensitivity DNA assay and sequenced samples on an Illumina NovaSeq 6000 in paired-end mode with 90-bp read2 and 20,000 read depth. Four samples for detection of chromosomal alteration were sequenced with 270-bp read2 to achieve 80,000 reads per cell for gene expression.

Alignment, quantification and QC of PBMC scRNA-seq data

We processed droplet libraries using Cell Ranger v.7 (10x Genomics). Sequencing reads were aligned with STAR using the GRCh38 human reference genome. We created raw count matrices with the cellranger_count function and processed them with the remove-background function of CellBender (v.0.1.0) to remove cell-free ambient RNAs in the droplets and call cells. We removed cells that had fewer than 500 unique molecular indices, fewer than 200 feature counts in each sample, more than 10% of reads from mitochondrial genes or more than 10% of reads from hemoglobin genes. In addition, putative doublets were removed using Scrublet (v.0.2.3) 73 and scds (v.1.18.0) 74 for each sample.

Manual annotation of PBMC scRNA-seq data

Python package Scanpy (v.1.9.8) 75 was used for data scaling, transformation, clustering, dimensionality reduction and most of the visualization. Expression values were normalized to a sum of 1 × 10 4 per cell and log-transformed with an added pseudocount of 1. For integration, we identified 3,000 shared highly variable genes using the scanpy.pp.highly_variable_genes() function. PCA was run on gene expression, followed by batch correction using the scanpy.pp.harmony_integrate() function 76 . Uniform manifold approximation and projection (UMAP) dimension reduction was performed based on the first 30 harmony-adjusted principal components 77 . We manually annotated cell types, accounting for Azimuth-based annotations 34 and previous studies of PBMCs (Extended Data Fig. 7 ).

Detection of cell types associated with MS

Using scDRS (v1.0.2) 35 , we assessed the associations between MS and cell types in the sc/snRNA-seq datasets for PBMCs and MS brain subcortical lesions. The PBMC scRNA-seq dataset comprised 163,541 cells from 20 patients with MS. The MS subcortical lesion snRNA-seq dataset comprised 74,686 cells from seven lesions of six patients with MS, with 66,402 cells remaining after removal of immune cells (T cells, B cells and myeloid cells). First, we used MAGMA 36 to calculate the gene P values and z -scores from the summary statistics of the cross-population MS GWAS meta-analysis, selecting the top 1,000 genes outside the MHC region (chr. 6: 25–34 Mb) as a set of putative disease genes following default settings. Second, we used scDRS to calculate a disease score for each cell by aggregating the expression of the putative disease genes and computed a set of 1,000 Monte Carlo control scores using random gene sets. Third, we used scDRS to normalize these scores and calculate P values for individual cells. We used the compute_score function with default parameters and included the number of genes per cell, age and sex as covariates. Finally, we performed analysis at the cell type level with scDRS to associate putative gene sets with cell types in the scRNA-seq datasets using the compute_downstream function with default settings. The false discovery rate (FDR) was calculated via the Benjamini–Hochberg method across all cell types in each dataset, and associations with FDR < 0.05 were considered to be statistically significant.

To replicate our findings in PBMC, we used the AIDA phase 1 Data Freeze v.2 PBMC scRNA-seq dataset 39 . This publicly available dataset includes four hierarchical annotation levels (levels 1 to 4). We used 1,219,548 cells annotated across these levels. To most closely match the annotations in our PBMC dataset, we created major cell types from levels 1 and 2 and used level 3 as finer cell types. We applied scDRS following the procedure as described above, except that the number of genes per cell and sex were included as covariates. We conducted a sensitivity analysis by varying the number of MAGMA genes used for scDRS (500, 750, 1,000, 1,250 and 1,500) and confirmed that the cell type specificity observed in PBMC and subcortical lesions remained robust (Supplementary Fig. 4 ).

Spatially resolved mapping of regions associated with MS

To detect regions associated with MS in subcortical lesion ST data generated by the 10x Genomics Visium ST platform, we applied gsMap 46 to the EUR GWAS meta-analysis results. We used a total of 12 MS subcortical lesion ST slides (8 active lesions and 4 inactive lesions) and considered the six niches defined in the original study 40 for annotations in gsMap. We implemented the gsMap joint analysis mode to produce unified results for all slides. First we generated the slice mean across 12 slides using the cleate_slice_mean command. Second, latent representations for spots were learned. Third, we identified homogeneous spots for each spot based on their latent representations and generated gene specificity scores for each spot using the slice mean. Fourth, we assigned gene specificity scores to SNPs using transcription start site information and computed the stratified LD score. Fifth, spatial LD score regression was performed to assess the associations with MS in individual spots. Sixth, we aggregated P values of individual spots across slides for each disease condition by the Cauchy combination method and obtained a P value for each niche annotation. All steps were conducted using default parameters.

We compared the strength of associations detected by gsMap with cell type proportions deconvoluted in the original publication of the ST dataset 40 . Furthermore, to identify characteristic cell types for each niche, we averaged cell type proportions per slide across all ST slides, including six control slides to match the conditions in the original publication, and compared them with those in other niches using Wilcoxon rank-sum test. P values were adjusted using the Benjamini–Hochberg method, and an FDR < 0.05 was considered to indicate statistical significance.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Online content

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41588-026-02741-5 .

Supplementary information

Supplementary Figs. 1–4.

Reporting Summary

Supplementary Tables 1–13.

Acknowledgements

We thank all the members and participants of the Japan MS/NMOSD Biobank, BBJ, IMSGC, UK Biobank, FinnGen and MVP. We also thank the All of Us participants for their contributions, without which this research would not have been possible, and the All of Us Research Program of the National Institutes of Health for making available the participant data examined in this study. This research was also made possible through access to data in the National Genomic Research Library, which is managed by Genomics England Limited (a wholly owned company of the Department of Health and Social Care). The National Genomic Research Library holds data provided by patients and collected by the NHS as part of their care and data collected as part of their participation in research and is funded by the National Institute for Health Research and NHS England. The Wellcome Trust, Cancer Research UK and the Medical Research Council have also funded research infrastructure.

Extended data

Author contributions

R.F., K.O., T.O., N.I. and Y. Okada designed the study and wrote the paper. R.F., K.O., S.N., Y. Ogawa, R.E., K.S., S.T., T.Y., Y. Shirai, Y.Y., G.S., C.K., T.N. and Y. Okada conducted data analysis. K.O., R.E., T.Y., Y. Shirai and K. Matsuda performed experiments. M.W., A.H., M.Y., F.S., M.K., M.M., M. Nakamori, Y. Shimizu, I.K., K. Miyamoto, M. Niino, A.K., Y.N., T.M., H.M. and J.-i.K. contributed to sample and data collection. T.O. N.I. and Y. Okada supervised the study. All authors read and approved the final paper.

Peer review

Peer review information

Nature Genetics thanks Harm-Jan Westra and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Funding

K.O. was supported by JSPS KAKENHI (JP24K18685), AMED (JP23tm0524002), Takeda Science Foundation, Kanzawa Medical Research Foundation and Kobayashi Foundation. M.W. was supported by JSPS KAKENHI (JP25K10629). T.O. was supported by JSPS KAKENHI (JP23K27517, JP25H01863) and AMED (JP22gm1810003). N.I. was supported by JSPS KAKENHI (24K02371), AMED (JP25wm0625504s0102) and the Health and Labour Sciences Research Grant on Intractable Diseases (Neuroimmunological Diseases) from the Ministry of Health, Labour and Welfare of Japan (23FC1009 and 26FC0201). Y. Okada was supported by JSPS KAKENHI (25H01057); AMED (JP223fa627001, JP223fa627010, JP223fa627011, JP22zf0127008, JP23tm0524002, JP24wm0625504, JP24gm1810011, JP25kk0305022, JP26tm0135237, JP26ek0410158); JST Moonshot R&D (JPMJMS2021, JPMJMS2024); Ono Pharmaceutical Foundation for Oncology, Immunology, and Neurology; a Baelz Research Grant; and the RIKEN TRIP initiative (AGIS).

Data availability

GWAS summary statistics and single-cell sequencing data for PBMCs are available at the National Bioscience Database Center (NBDC) Human Database with accession ID hum0197 . GWAS genotype data from the BBJ are available at the NBDC Human Database (research ID: hum0014 and hum0311 ). The IMSGC GWAS summary statistics were obtained from https://imsgc.net/ , and those of ref. 9 and for other traits were obtained from the GWAS Catalog ( https://www.ebi.ac.uk/gwas/ ). The UK Biobank analysis was conducted via application number 47821 ( https://www.ukbiobank.ac.uk/ ); individual-level phenotypic and genetic data from UK Biobank can be requested from the UK Biobank Access Management System ( https://www.ukbiobank.ac.uk/enable-your-research/apply-for-access ). Individual-level phenotypic and genetic data from the All of Us Research Program can be accessed following application to use the All of Us Researcher Workbench, a cloud-based computing platform ( https://www.researchallofus.org/register/ ), and data from the 100,000 Genomes Project can be accessed following application to join the Genomics England Research Network ( https://www.genomicsengland.co.uk/research/academic/join-research-network ). The GWAS summary statistics for FinnGen Data Freeze 11 were obtained from https://www.finngen.fi/en , and those for MVP were obtained through dbGaP (accession number phs002453 ). The AIDA phase 1 Data Freeze v2 PBMC scRNA-seq data were obtained from CZ CELLxGENE ( https://cellxgene.cziscience.com/collections/ced320a1-29f3-47c1-a735-513c7084d508 ), and the MS subcortical lesion snRNA-seq and ST data were obtained from the UCSC Cell Browser ( https://cells-test.gi.ucsc.edu/?ds=ms-subcortical-lesions ).

Code availability

We used publicly available software for the analyses, as described in Methods .

Competing interests

The authors declare no competing interests.

Footnotes

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

These authors contributed equally: Rintaro Fujimoto, Kotaro Ogawa.

These authors jointly supervised this work: Tatsusada Okuno, Noriko Isobe, Yukinori Okada.

Lists of authors and their affiliations appear at the end of the paper.

Contributor Information

Tatsusada Okuno, Email: okuno@neurol.med.osaka-u.ac.jp.

Noriko Isobe, Email: isobe.noriko.342@m.kyushu-u.ac.jp.

Yukinori Okada, Email: yuki-okada@m.u-tokyo.ac.jp.

Japan MS/NMOSD Biobank:

Mitsuru Watanabe , Fumitaka Shimizu , Yuko Shimizu , Izumi Kawachi , Katsuichi Miyamoto , Masaaki Niino , Yuji Nakatsuji , Takuya Matsushita , Jun-ichi Kira , Tatsusada Okuno , and Noriko Isobe

the BioBank Japan Project:

Shinichi Namba , Koichi Matsuda , and Yukinori Okada

Extended data

is available for this paper at https://doi.org/10.1038/s41588-026-02741-5 .

Supplementary information

The online version contains supplementary material available at https://doi.org/10.1038/s41588-026-02741-5 .

References

1. Woo, M. S., Engler, J. B. & Friese, M. A. The neuropathobiology of multiple sclerosis. Nat. Rev. Neurosci. 25 , 493–513 (2024).

2. Walton, C. et al. Rising prevalence of multiple sclerosis worldwide: insights from the Atlas of MS, third edition. Mult. Scler. 26 , 1816–1821 (2020).

3. Eskandarieh, S., Heydarpour, P., Minagar, A., Pourmand, S. & Sahraian, M. A. Multiple sclerosis epidemiology in East Asia, South East Asia and South Asia: a systematic review. Neuroepidemiology 46 , 209–221 (2016).

4. Goodin, D. S., Khankhanian, P., Gourraud, P. A. & Vince, N. The nature of genetic and environmental susceptibility to multiple sclerosis. PLoS ONE 16 , e0246157 (2021).

5. Jacobs, B. M. et al. Towards a global view of multiple sclerosis genetics. Nat. Rev. Neurol. 18 , 613–623 (2022).

6. Sawcer, S. et al. Genetic risk and a primary role for cell-mediated immune mechanisms in multiple sclerosis. Nature 476 , 214–219 (2011).

7. Beecham, A. H. et al. Analysis of immune-related loci identifies 48 new susceptibility variants for multiple sclerosis. Nat. Genet. 45 , 1353–1360 (2013).

8. Isobe, N. et al. An ImmunoChip study of multiple sclerosis risk in African Americans. Brain 138 , 1518–1530 (2015).

9. Andlauer, T. F. M. et al. Novel multiple sclerosis susceptibility loci implicated in epigenetic regulation. Sci. Adv. 2 , 54 (2016).

10. Patsopoulos, N. A. et al. Multiple sclerosis genomic map implicates peripheral immune cells and microglia in susceptibility. Science 365 , eaav7188 (2019).

11. Bycroft, C. et al. The UK Biobank resource with deep phenotyping and genomic data. Nature 562 , 203–209 (2018).

12. Verma, A. et al. Diversity and scale: genetic architecture of 2068 traits in the VA Million Veteran Program. Science 385 , eadj1182 (2024).

13. Ogawa, K. et al. Next-generation sequencing identifies contribution of both class I and II HLA genes on susceptibility of multiple sclerosis in Japanese. J. Neuroinflammation 16 , 162 (2019).

14. Watanabe, M. et al. HLA genotype-clinical phenotype correlations in multiple sclerosis and neuromyelitis optica spectrum disorders based on Japan MS/NMOSD Biobank data. Sci. Rep. 11 , 607 (2021).

15. Zhou, W. et al. Global Biobank Meta-analysis Initiative: powering genetic discovery across human disease. Cell Genomics 2 , 100192 (2022).

16. Harroud, A. et al. Locus for severity implicates CNS resilience in progression of multiple sclerosis. Nature 619 , 323–331 (2023).

17. Cheng, M. et al. Spatially resolved transcriptomics: a comprehensive review of their technological advances, applications, and challenges. J. Genet. Genomics 50 , 625–640 (2023).

18. Yata, T. et al. Contribution of germline and somatic mutations to risk of neuromyelitis optica spectrum disorder. Cell Genomics 5 , 100776 (2025).

19. Nagai, A. et al. Overview of the BioBank Japan Project: study design and profile. J. Epidemiol. 27 , S2–S8 (2017).

20. Garrett-Sinha, L. A. Review of Ets1 structure, function, and roles in immunity. Cell. Mol. Life Sci. 70 , 3375 (2013).

21. Ann Garrett-Sinha, L. & Garrett-Sinha, A. An update on the roles of transcription factor Ets1 in autoimmune diseases. WIREs Mech. Dis. 15 , e1627 (2023).

22. Luo, Y. et al. Endothelial ETS1 inhibition exacerbate blood–brain barrier dysfunction in multiple sclerosis through inducing endothelial-to-mesenchymal transition. Cell Death Dis. 13 , 462 (2022).

23. Kurki, M. I. et al. FinnGen provides genetic insights from a well-phenotyped isolated population. Nature 613 , 508–518 (2023).

24. Bick, A. G. et al. Genomic data in the All of Us Research Program. Nature 627 , 340–346 (2024).

25. Caulfield, M. et al. National Genomic Research Library. figshare 10.6084/m9.figshare.4530893.v7 (2020).

26. Kumanogoh, A. et al. Requirement for the lymphocyte semaphorin, CD100, in the induction of antigen-specific T cells and the maturation of dendritic cells. J. Immunol. 169 , 1175–1181 (2002).

27. Okuno, T. et al. Roles of Sema4D–Plexin-B1 interactions in the central nervous system for pathogenesis of experimental autoimmune encephalomyelitis. J. Immunol. 184 , 1499–1506 (2010).

28. Camiña-Tato, M. et al. Gender-associated differences of perforin polymorphisms in the susceptibility to multiple sclerosis. J. Immunol. 185 , 5392–5404 (2010).

29. Hirata, J. et al. Genetic and phenotypic landscape of the major histocompatibilty complex region in the Japanese population. Nat. Genet. 51 , 470–480 (2019).

30. Moutsianas, L. et al. Class II HLA interactions modulate genetic risk for multiple sclerosis. Nat. Genet. 47 , 1107–1113 (2015).

31. Hollenbach, J. A. & Oksenberg, J. R. The immunogenetics of multiple sclerosis: a comprehensive review. J. Autoimmun. 64 , 13–25 (2015).

32. Namkoong, H. et al. DOCK2 is involved in the host genetics and biology of severe COVID-19. Nature 609 , 754–760 (2022).

33. Attfield, K. E., Jensen, L. T., Kaufmann, M., Friese, M. A. & Fugger, L. The immunology of multiple sclerosis. Nat. Rev. Immunol. 22 , 734–750 (2022).

34. Stuart, T. et al. Comprehensive integration of single-cell data. Cell 177 , 1888–1902.e21 (2019).

35. Zhang, M. J. et al. Polygenic enrichment distinguishes disease associations of individual cells in single-cell RNA-seq data. Nat. Genet. 54 , 1572–1580 (2022).

36. Leeuw, C. A. de, Mooij, J. M., Heskes, T. & Posthuma, D. MAGMA: generalized gene-set analysis of GWAS data. PLOS Comput. Biol. 11 , e1004219 (2015).

37. Kunkl, M., Frascolla, S., Amormino, C., Volpe, E. & Tuosto, L. T helper cells: the modulators of inflammation in multiple sclerosis. Cells 9 , 482 (2020).

38. Verma, N. D. et al. Multiple sclerosis patients have reduced resting and increased activated CD4 + CD25 + FOXP3 + T regulatory cells. Sci. Rep. 11 , 10476 (2021).

39. Kock, K. H. et al. Asian diversity in human immune cells. Cell 188 , 2288–2306 (2025).

40. Lerma-Martin, C. et al. Cell type mapping reveals tissue niches and interactions in subcortical multiple sclerosis lesions. Nat. Neurosci. 27 , 2354–2365 (2024).

41. Abbott, N. J., Patabendige, A. A. K., Dolman, D. E. M., Yusof, S. R. & Begley, D. J. Structure and function of the blood–brain barrier. Neurobiol. Dis. 37 , 13–25 (2010).

42. Ortiz, G. G. et al. Role of the blood–brain barrier in multiple sclerosis. Arch. Med. Res. 45 , 687–697 (2014).

43. Brosnan, C. F. & Raine, C. S. The astrocyte in multiple sclerosis revisited. Glia 61 , 453–465 (2013).

44. Ponath, G., Park, C. & Pitt, D. The role of astrocytes in multiple sclerosis. Front. Immunol. 9 , 337102 (2018).

45. Dal-Bianco, A. et al. Slow expansion of multiple sclerosis iron rim lesions: pathology and 7 T magnetic resonance imaging. Acta Neuropathol. 133 , 25–42 (2017).

46. Song, L., Chen, W., Hou, J., Guo, M. & Yang, J. Spatially resolved mapping of cells associated with human complex traits. Nature 641 , 932–941 (2025).

47. Ishigaki, K. et al. Multi-ancestry genome-wide association analyses identify novel genetic mechanisms in rheumatoid arthritis. Nat. Genet. 54 , 1640–1651 (2022).

48. Bentham, J. et al. Genetic association analyses implicate aberrant regulation of innate and adaptive immunity genes in the pathogenesis of systemic lupus erythematosus. Nat. Genet. 47 , 1457–1464 (2015).

49. de Lange, K. M. et al. Genome-wide association study implicates immune activation of multiple integrin genes in inflammatory bowel disease. Nat. Genet. 49 , 256–261 (2017).

50. Ferreira, M. A. et al. Shared genetic origin of asthma, hay fever and eczema elucidates allergic disease biology. Nat. Genet. 49 , 1752–1757 (2017).

51. Bellenguez, C. et al. New insights into the genetic etiology of Alzheimer’s disease and related dementias. Nat. Genet. 54 , 412–436 (2022).

52. Als, T. D. et al. Depression pathophysiology, risk prediction of recurrence and comorbid psychiatric disorders using genome-wide analyses. Nat. Med. 29 , 1832–1844 (2023).

53. Graham, S. E. et al. The power of genetic diversity in genome-wide association studies of lipids. Nature 600 , 675–679 (2021).

54. Alfredsson, L. & Olsson, T. Lifestyle and environmental factors in multiple sclerosis. Cold Spring Harb. Perspect. Med. 9 , a028944 (2019).

55. Giovannoni, G. et al. Smouldering multiple sclerosis: the ‘real MS’. Ther. Adv. Neurol. Disord. 15 , 17562864211066751 (2022).

56. Polman, C. H. et al. Diagnostic criteria for multiple sclerosis: 2010 revisions to the McDonald Criteria. Ann. Neurol. 69 , 292–302 (2011).

57. International HapMap Consortium A haplotype map of the human genome. Nature 437 , 1299–1320 (2005).

58. Manichaikul, A. et al. Robust relationship inference in genome-wide association studies. Bioinformatics 26 , 2867–2873 (2010).

59. Chang, C. C. et al. Second-generation PLINK: rising to the challenge of larger and richer datasets. GigaScience 4 , 7 (2015).

60. Tadaka, S. et al. 3.5KJPNv2: an allele frequency panel of 3552 Japanese individuals including the X chromosome. Hum. Genome Var. 6 , 28 (2019).

61. Delaneau, O., Zagury, J.-F., Robinson, M. R., Marchini, J. L. & Dermitzakis, E. T. Accurate, scalable and integrative haplotype estimation. Nat. Commun. 10 , 5436 (2019).

62. Das, S. et al. Next-generation genotype imputation service and methods. Nat. Genet. 48 , 1284–1287 (2016).

63. Zhou, W. et al. Efficiently controlling for case-control imbalance and sample relatedness in large-scale genetic association studies. Nat. Genet. 50 , 1335–1341 (2018).

64. Shi, S. et al. A Genomics England haplotype reference panel and imputation of UK Biobank. Nat. Genet. 56 , 1800–1803 (2024).

65. Mbatchou, J. et al. Computationally efficient whole-genome regression for quantitative and binary traits. Nat. Genet. 53 , 1097–1103 (2021).

66. Hinrichs, A. S. et al. The UCSC Genome Browser Database: update 2006. Nucleic Acids Res. 34 , D590–D598 (2006).

67. Cerezo, M. et al. The NHGRI-EBI GWAS Catalog: standards for reusability, sustainability and diversity. Nucleic Acids Res. 53 , D998–D1005 (2025).

68. Willer, C. J., Li, Y. & Abecasis, G. R. METAL: fast and efficient meta-analysis of genomewide association scans. Bioinformatics 26 , 2190–2191 (2010).

69. Lin, D.-Y. & Sullivan, P. F. Meta-analysis of genome-wide association studies with overlapping subjects. Am. J. Hum. Genet. 85 , 862–872 (2009).

70. Lee, C. H., Eskin, E. & Han, B. Increasing the power of meta-analysis of genome-wide association studies to detect heterogeneous effects. Bioinformatics 33 , i379–i388 (2017).

71. Watanabe, K., Taskesen, E., van Bochoven, A. & Posthuma, D. Functional mapping and annotation of genetic associations with FUMA. Nat. Commun. 8 , 1826 (2017).

72. Naito, T. et al. A deep learning method for HLA imputation and trans-ethnic MHC fine-mapping of type 1 diabetes. Nat. Commun. 12 , 1639 (2021).

73. Wolock, S. L., Lopez, R. & Klein, A. M. Scrublet: computational identification of cell doublets in single-cell transcriptomic data. Cell Syst. 8 , 281–291 (2019).

74. Bais, A. S. & Kostka, D. scds: computational annotation of doublets in single-cell RNA sequencing data. Bioinformatics 36 , 1150–1158 (2020).

75. Wolf, F. A., Angerer, P. & Theis, F. J. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 19 , 15 (2018).

76. Korsunsky, I. et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat. Methods 16 , 1289–1296 (2019).

77. McInnes, L., Healy, J. & Melville, J. UMAP: uniform manifold approximation and projection for dimension reduction. Preprint at https://arxiv.org/abs/1802.03426 (2020).

Associated Data

Supplementary Materials

Supplementary Figs. 1–4.

Reporting Summary

Supplementary Tables 1–13.

Data Availability Statement

GWAS summary statistics and single-cell sequencing data for PBMCs are available at the National Bioscience Database Center (NBDC) Human Database with accession ID hum0197 . GWAS genotype data from the BBJ are available at the NBDC Human Database (research ID: hum0014 and hum0311 ). The IMSGC GWAS summary statistics were obtained from https://imsgc.net/ , and those of ref. 9 and for other traits were obtained from the GWAS Catalog ( https://www.ebi.ac.uk/gwas/ ). The UK Biobank analysis was conducted via application number 47821 ( https://www.ukbiobank.ac.uk/ ); individual-level phenotypic and genetic data from UK Biobank can be requested from the UK Biobank Access Management System ( https://www.ukbiobank.ac.uk/enable-your-research/apply-for-access ). Individual-level phenotypic and genetic data from the All of Us Research Program can be accessed following application to use the All of Us Researcher Workbench, a cloud-based computing platform ( https://www.researchallofus.org/register/ ), and data from the 100,000 Genomes Project can be accessed following application to join the Genomics England Research Network ( https://www.genomicsengland.co.uk/research/academic/join-research-network ). The GWAS summary statistics for FinnGen Data Freeze 11 were obtained from https://www.finngen.fi/en , and those for MVP were obtained through dbGaP (accession number phs002453 ). The AIDA phase 1 Data Freeze v2 PBMC scRNA-seq data were obtained from CZ CELLxGENE ( https://cellxgene.cziscience.com/collections/ced320a1-29f3-47c1-a735-513c7084d508 ), and the MS subcortical lesion snRNA-seq and ST data were obtained from the UCSC Cell Browser ( https://cells-test.gi.ucsc.edu/?ds=ms-subcortical-lesions ).

We used publicly available software for the analyses, as described in Methods .