EBV reactivation priming of the peripheral immune system in multiple sclerosis relapse

Abstract Despite decades of research, the cellular and molecular events preceding multiple sclerosis (MS) relapse remain incompletely understood. Here, in this observational study of longitudinal blood samples from patients with relapsing-remitting MS, we used single-cell RNA sequencing, bulk transcriptomics, multiparameter flow cytometry and targeted viral reverse transcription quantitative polymerase chain reaction (RT−qPCR) to construct a time-resolved atlas of immune perturbations surrounding relapse. A reproducible pre-relapse signature in monocytes and B cells, emerging up to 3 months before clinical onset, was enriched for host genes responsive to Epstein−Barr virus (EBV) lytic reactivation factors. RT−qPCR confirmed elevated EBV LMP-1 transcripts in pre-relapse B cells, and flow cytometry demonstrated expansion of CD11c+ atypical B cell populations displaying EBV surface protein gp350. Pre-relapse transcriptional modules overlapped with MS genome-wide association study (GWAS) risk loci and EBNA-2-bound enhancers, suggesting that inherited MS susceptibility and EBV-responsive programs operate through shared regulatory elements. How this peripheral activation relates to central nervous system lesion formation remains to be established. These findings nonetheless suggest that EBV reactivation, when occurring within a genetically predisposed peripheral immune environment, is a proximal precursor of MS relapse. Main Multiple sclerosis (MS) is a chronic inflammatory disease of the central nervous system (CNS) and a leading cause of neurologic disability in young adults, affecting approximately 2.9 million people worldwide1. The disease course, particularly in early stages, is characterized by symptomatic flare-ups separated by periods of remission whose timing is highly variable and difficult to predict. A central challenge in identifying the triggers of relapse is uncertainty in the timing of molecular events that precede clinical symptoms. Multiple inflammatory processes participate in relapse pathophysiology, many of which are transient. In viral meningitis, for example, neutrophilic pleocytosis in the cerebrospinal fluid may be detectable for only hours after infection, whereas lymphocytic dominance persists weeks later2. In MS, the inability to sample immune activity during such short-lived presymptomatic windows has limited efforts to capture the molecular events that initiate relapses. A growing body of evidence supports serum neurofilament light chain (NfL) as a useful biomarker of disease activity in MS3,4,5. Moreover, proteomic studies have identified elevated cytokines that accompany relapse6,7. However, these markers reflect tissue damage that has already occurred and an inflammatory process already underway, rather than the upstream factors that instigate relapse. In parallel, B cell dysregulation has emerged as a central feature of MS relapse biology. Since B cell depletion was shown to reduce relapse rates8, experimental and human studies have demonstrated that B cells function not only as antibody producers but also as antigen-presenting cells and sources of inflammatory cytokines that activate pathogenic T cell programs9,10. Among B cell subsets, CD11c+T-bet+ atypical memory B cells (ABCs) expand in autoimmune diseases and chronic viral infections, including Epstein–Barr virus (EBV), HIV and hepatitis C11,12,13. In MS, ABCs accumulate in blood and CNS-adjacent compartments, display elevated expression of interferon-stimulated genes and antigen presentation machinery14,15,16 and are less effectively depleted by anti-CD20 therapies17. Whether ABCs expand before relapse and whether they contribute to the inflammatory cascade remain undetermined. MS is further shaped by host genetics. Genome-wide association studies (GWASs) have identified more than 200 independent risk loci, many mapping to regulators of immune activation, B cell differentiation and antigen presentation18. HLA-DRB1*15:01 confers the largest individual risk effect19, and additional non-HLA loci are concentrated in B cell and innate immunoregulatory pathways20,21,22,23. Notably, the EBV nuclear antigen EBNA-2 binds preferentially to enhancers at MS risk loci and remodels host transcriptional landscapes24,25,26, providing a direct molecular link between viral biology and host genetic susceptibility. These studies, however, address lifetime disease risk; whether the same alleles shape the timing and instigation of relapses is unknown. Here we present a time-resolved, single-cell atlas of the peripheral immune system surrounding MS relapse. We identify reproducible immune perturbations emerging weeks to months before clinical relapse, centered on monocytes and CD11c+ atypical B cell populations, enriched for host genes regulated by early-lytic EBV factors and overlapping with MS GWAS risk loci and EBNA-2-bound regulatory elements. Although the observational design precludes definitive causal inference, our findings implicate viral reactivation and host genetic susceptibility in the immune cascade that precedes relapse. Results Study design To investigate immune mechanisms preceding MS relapse, we leveraged the Comprehensive Longitudinal Investigation in MS at Brigham and Women’s Hospital (CLIMB) study, an ongoing prospective cohort that has followed thousands of patients with MS for more than 25 years through annual magnetic resonance imaging (MRI), standardized neurological examinations and systematic biospecimen collection27. As part of routine CLIMB enrollment, participants donate blood samples at annual or scheduled clinic visits, and peripheral blood mononuclear cells (PBMCs) are cryopreserved alongside detailed clinical metadata. Pre-relapse samples are exceedingly rare, identifiable only retrospectively when a clinically stable patient happens to suffer a relapse after a routine blood draw. Decades of systematic biobanking were required to accumulate sufficient specimens for molecular investigation of the presymptomatic immune state. We screened 85 participants with available biospecimens and identified 15 individuals for our initial single-cell discovery profiling based on the availability of remission, pre-relapse and/or acute relapse samples captured within defined peri-attack windows. An additional 99 patients were evaluated using one or more orthogonal methodologies: bulk RNA sequencing (RNA-seq) of sorted CD19+ B cells (n = 70 patients, one of whom also contributed to the discovery cohort), targeted viral transcript detection by reverse transcription quantitative polymerase chain reaction (RT−qPCR) (n = 30) or multiparameter flow cytometry (n = 23) (Supplementary Fig. 1 and Supplementary Table 1). EBNA-1 IgG ELISA was performed on 59 participants with available serum; all 38 patients with MS were seropositive, as were 19 of 21 healthy controls (Supplementary Table 2). These 21 healthy controls served as a comparative reference population for single-cell sequencing (Fig. 1a–e). All samples were collected at least 30 days after corticosteroid exposure, and disease-modifying therapy (DMT) was unchanged between paired timepoints in all but two patients (Supplementary Table 3). Age distributions were comparable between the MS and healthy cohorts, with a higher proportion of female participants in the MS cohort, consistent with MS epidemiology (Fig. 1c). Relapses were confirmed independently by two neurologists and required a contrast-enhancing lesion on MRI. Remission required flanking clinical and radiographic stability (Fig. 1d). A time-resolved single-cell atlas identifies prioritized immune states PBMCs from the discovery cohort were subjected to 10x Genomics 3′ single-cell RNA sequencing (scRNA-seq), yielding 281,799 high-quality cells after quality control filtering and doublet exclusion, including 154,709 cells from patients with MS across remission, pre-relapse and relapse timepoints and 127,090 cells from healthy donors (Fig. 2a and Extended Data Fig. 1). Datasets were integrated using Harmony batch correction while preserving disease-associated transcriptional structure28. Reference-based annotation complemented by careful manual assignment resolved major circulating lineages, including classical and non-classical monocytes, naive and memory B cells, plasmablasts, T cell subsets, dendritic cells and natural killer cells (Fig. 2a,b and Extended Data Fig. 2). Key patient and technical metadata, including donor, DMT, sex, age and experimental batch, were inspected on the integrated uniform manifold approximation and projection (UMAP) (Extended Data Fig. 3), and donor was modeled as a random effect in downstream analyses. To quantify cell-state-specific transcriptional shifts while accounting for donor-level heterogeneity, we applied the scDist mixed-effects model29. Negative control contrasts (healthy versus healthy, remission versus remission) produced negligible perturbations, confirming robustness (Extended Data Fig. 4). A complementary sensitivity analysis restricted to donors with paired pre-relapse and remission samples recapitulated the primary perturbation pattern (Extended Data Fig. 5). In striking contrast, the pre-relapse versus remission comparison revealed widespread transcriptional divergence across nearly all immune populations (Fig. 2c), with the most pronounced perturbations in FGFBP2+ cytotoxic γδ T cells, monocytes and activated type 2 conventional dendritic cells (cDC2s) (Fig. 2d). Multiple B cell subsets also ranked among the top third of perturbed populations in the pre-relapse contrast. Notably, perturbation magnitudes were generally larger in the pre-relapse state than during overt relapse, possibly indicating that peripheral immune dysregulation peaks during the presymptomatic phase before attenuating once clinical symptoms emerge. To integrate genetic risk information at single-cell resolution, we performed patient-aware single-cell disease relevance scoring (scDRS) using MS GWAS summary statistics18 (Fig. 2e,f). Across PBMC states, pre-relapse perturbation magnitude was not correlated with the pre-relapse change in MS genetic risk enrichment (Spearmanʼs ρ = −0.14, P = 0.53). B cell subsets were the exception: intermediate, switched memory and TCL1A+ naive B cells showed both substantial pre-relapse perturbation (scDist 5.79−6.47) and the three largest pre-relapse increases in scDRS score of any population (Δ = 0.44−0.54), rendering B cells a population that we dissect at high resolution below (Fig. 2g). An atypical B cell axis defines the pre-relapse immune state We reclustered approximately 30,000 high-confidence B cells to resolve within-lineage heterogeneity. This analysis delineated 14 discrete populations, including naive, switched memory, activated switched memory, activated memory, atypical/ABC-like memory, ABC, CD11c++ activated memory, plasmablast and cycling plasmablast states, among others (Fig. 3a–c and Extended Data Fig. 6). Imputed surface protein markers were consistent with subset organization and distinguished naive, memory, activated and plasmablast compartments (Fig. 3c). Subset-level scDist analysis revealed that naive B cells exhibited the largest transcriptional perturbations in the pre-relapse condition, followed by atypical/ABC-like memory, ABC, switched memory and CD11c++ activated memory B cells (Fig. 3d). Differential expression revealed enrichment for type I interferon signaling, antigen processing and viral defense pathways, with canonical interferon-stimulated genes (IFITM2, IRF7) and major histocompatibility complex class II (MHC-II) components (HLA-DQA1, HLA-DQA2) upregulated in pre-relapse naive B cells and HLA-DQA1 additionally upregulated in atypical/ABC-like memory, switched memory, anergic/naive-leaning and ABC subsets (Extended Data Fig. 7). Viral mRNA translation was the top Reactome pathway upregulated in pre-relapse B cells (Benjamini–Hochberg-adjusted P = 6.3 × 10−90). These enriched pathways implicate early activation of innate antiviral programs as a defining feature of the pre-relapse immune state. Notably, these enrichments persisted after regression of cell-cycle-associated transcripts, excluding generalized proliferation as the basis and supporting a bona fide activation phenotype rather than non-specific clonal expansion. Patient-aware scDRS analysis showed that ABCs carried the strongest enrichment (estimated marginal mean (EMM) = 1.16, 95% confidence interval (CI): 1.06−1.25), matched by cycling plasmablasts (1.16, 0.81−1.50, although estimated from only 70 cells), followed by CD11c++ activated memory B, activated switched memory B, plasmablast and atypical/ABC-like memory B cells (Fig. 3e). Plotting pre-relapse perturbation against the within-patient change in genetic risk enrichment revealed that ABC, atypical/ABC-like memory B, activated switched memory B, switched memory B and transitional/immature-like B states occupied the upper-right quadrant of the scDist−scDRS plot, where pre-relapse transcriptional perturbation and inherited MS susceptibility co-occur (Fig. 3f). Atypical and activated antigen-experienced B cell states are thus enriched for both MS genetic risk programs and prerelapse transcriptional changes. ABC-like and EBV gp350+ B cells expand before relapse To validate transcriptomic findings at the protein level, we performed multiparameter flow cytometry on paired pre-relapse and remission PBMCs from a validation cohort using a deep B cell phenotyping panel (Fig. 4a and Extended Data Fig. 8a–c). After gating on live singlet CD19+ B cells, FlowSOM clustering identified metaclusters corresponding to known B cell subsets, including naive, memory, activated memory, ABC-like and plasmablast populations (Fig. 4b,c). Core clustering markers included CD20, CD11c, CD21, CD27, CD18, CD38, CD172a, CD5 and CXCR4, with additional overlay markers used for phenotypic characterization. We label the CD11c+ metaclusters (13 and 15, analyzed together) ‘ABC-like’ given the surface marker phenotype but note that without intracellular staining for T-bet and ZEB2, a canonical ABC identity is not confirmed. Within-patient analysis of 23 paired pre-relapse and remission samples demonstrated a significant increase in ABC-like B cell abundance during the pre-relapse window (21 of 23 patients increased; median 1.6-fold; two-sided Wilcoxon signed-rank P = 8.8 × 10−5) (Fig. 4d). Using fluorescence-minus-one (FMO) thresholding, we assessed expression of the EBV late-lytic glycoprotein gp350 across metaclusters. gp350+ cells localized preferentially with activated memory and ABC-like metaclusters (Fig. 4e,f). Critically, the frequency of gp350+ B cells among CD19+ cells was significantly increased in pre-relapse relative to remission (binomial mixed model with donor random intercepts and experiment day as a fixed effect: odds ratio = 1.41, 95% CI: 1.11−1.79, P = 0.0043; model-predicted frequency 0.263% versus 0.186%), with frequencies peaking in the 0−90-day pre-relapse window (odds ratio = 1.47, 95% CI: 1.12−1.94, Benjamini–Hochberg-adjusted P = 0.011) (Fig. 4g and Extended Data Fig. 8d,e). These data provide orthogonal validation that ABC-like B cells expand before relapse, accompanied by an increased frequency of EBV gp350+ B cells. EBV-associated viral and host transcriptional programs are elevated before relapse To directly assess EBV transcript expression during disease activity windows, we performed RT−qPCR on FACS-sorted CD19+ B cells from 30 patients with MS with paired pre-relapse and remission samples (25 passed quality control), targeting LMP-1, BZLF1, BLLF1 (gp350), EBER2 and additional EBV transcripts (Fig. 5a). Principal component analysis (PCA) of normalized ΔCt values revealed separation along principal component 1 of pre-relapse from remission samples based on EBV transcript profiles alone (Fig. 5b). Among individual transcripts, LMP-1 was significantly elevated in pre-relapse samples compared to paired remission specimens (odds ratio = 1.50, 95% CI: 1.26−1.93, Benjamini−Hochberg-adjusted P = 0.0022) (Fig. 5c). BZLF2 (encoding gp42) transcripts showed a similar pattern that did not reach significance (P = 0.058). Prior in vitro studies have shown that BZLF1 transcripts peak 8−12 hours after lytic induction, and the protein product Zta persists for several days30, implying a narrow sampling window for capturing lytic transcripts in clinical specimens. LMP-1, by contrast, is expressed during both pro-growth latency II/III programs and the lytic cycle31,32,33, potentially affording a wider detection window. To systematically associate transcriptional responses observed during the pre-relapse window to defined phases of EBV activity, we integrated our data with a curated EBV−human co-expression atlas derived from 201 lymphoblastoid cell line transcriptomes spanning latent and lytic viral states34. In this reference, host genes are grouped into modules based on their correlation with specific EBV factors. Fold over-representation analysis of genes driving the pre-relapse ABC scDist signature revealed enrichment for host gene modules correlated with both latent growth program and early-lytic EBV factors. The three strongest were the modules for BSLF2/BMLF1, LMP-1 and BZLF1. Of the 79 modules tested, 13 were significantly enriched, comprising nine of 32 early lytic, three of 10 latent and one of six unassigned modules, whereas none of the 31 true late or leaky-late modules reached significance (Fig. 5d and Extended Data Fig. 9). Bulk RNA-seq validation in independently sorted CD19+ B cells from the validation cohort confirmed enrichment of an LMP-1 host response signature in pre-relapse (0−90 days before onset) relative to remission samples (linear mixed model with donor random intercepts: difference in singscore enrichment = 0.0094, 95% CI: 0.0056−0.0131, P = 8.7 × 10⁻6) (Fig. 5e). Within the 0−90-day pre-relapse window, enrichment of the LMP-1 response peaked approximately 30−90 days before relapse onset and declined during peri-relapse intervals, mirroring temporal dynamics observed in the discovery cohort. These patterns of immune response to early-lytic reactivation without true late-lytic programs, which strictly depend on viral DNA replication for their expression35,36, are consistent with predominantly abortive EBV reactivation (Extended Data Fig. 10). Surface gp350 detected by flow cytometry suggests that a minority of cells nonetheless complete the cycle or retain residual late-lytic protein, so abortive and completed outcomes may coexist (or, alternatively, late-lytic responses may be transient and resolved by the time of sampling). In any case, EBV reactivation might not need to be productive to exert pathogenic effects on B cells. Because EBV-correlated host gene modules overlap with B cell proliferation signatures34, we tested whether the pre-relapse ABC signature could instead reflect general B cell proliferation. It remained separable from host genes whose promoters are bound by the latency-associated nuclear antigen (LANA) of Kaposiʼs sarcoma-associated herpesvirus (KSHV), a related gammaherpesvirus that also drives B cell proliferation (Fig. 6a), and persisted after regressing out cell cycle variation, together supporting an EBV-specific rather than a proliferation-driven origin for these transcriptional changes. EBNA-2-bound regulatory elements and MS genetic risk define pre-relapse B cell states To determine if EBV-associated transcriptional programs and pre-relapse perturbation co-occur in specific B cell states, we performed pre-ranked gene set enrichment analysis (fgsea) using external EBNA-2 chromatin immunoprecipitation followed by sequencing (ChIP−seq) promoter occupancy gene sets, together with EBNA-2-associated activation and repression signatures derived from intersected RNA-seq studies, across all reclustered B cell subsets (Fig. 6a). Promoter occupancy sets were either shared between EBV types 1 and 2 or specific to each type. KSHV LANA promoter occupancy gene sets served as non-EBV herpesvirus controls. ABC, switched memory B, transitional/immature-like B and plasmablast subsets showed significant EBNA-2-associated enrichment across three gene sets each at false discovery rate (FDR)-corrected thresholds, and four further subsets (activated switched memory B, activated memory B, anergic/naive-leaning B and TCL1A+ naive B) across two, whereas KSHV LANA enrichment was not significant in any subset, supporting viral specificity of the observed transcriptional programs. Considered together, pre-relapse perturbation (scDist), EBNA-2 enrichment (the mean supportive normalized enrichment score across prioritized gene sets, driven particularly by the EBV type 1 set) and MS genetic prioritization (scDRS) placed ABC states at the intersection of all three, followed by activated switched memory B cells (Fig. 6b). Figure 6c presents a working model in which EBV-associated activity in the reservoir memory B cell compartment promotes preferential survival of MS risk-enriched ABC-like and CD11c++ B cell states, with associated peripheral immune activation that precedes CNS inflammation. Discussion Decades of prospective biobanking have now made it possible to capture critical windows of disease activity in MS, revealing that EBV activity and associated immune activation precede relapse by weeks to months. This observation helps reconcile the epidemiologic requirement of EBV infection for MS onset37 with the intermittent, unpredictable nature of relapses and is consistent with the clinical efficacy of B cell depletion therapies8, possibly by reducing the latent viral reservoir or intercepting a transient pathogenic state. LMP-1 transcripts were selectively upregulated in pre-relapse B cells. LMP-1 encodes LMP-1, which has historically been understood to alter B cell behavior by mimicking constitutive CD40 signaling to drive NF-κB-dependent survival, antigen presentation and cytokine responses. However, recent work has demonstrated that LMP-1ʼs TES1 domain activates non-canonical NF-κB signaling through proteasomal degradation of cIAP1, cIAP2 and TRAF2, a mechanism fundamentally distinct from B cell receptor-mediated canonical NF-κB activation38. Non-canonical NF-κB signaling was among the most enriched transcriptomic signatures in pre-relapse CD19+ B cells. Beyond LMP-1, most assayed EBV genes showed positive pre-relapse effect sizes (Fig. 5c), a pattern consistent with coordinated viral transcriptional activity captured at variable sampling intervals before relapse. CD11c+ ABC-like B cells that expand before relapse represent a plausible intermediary between the EBV reservoir and subsequent infiltration of the CNS by peripheral immune cells. ABC-like cells exhibit a CXCR4-high/CXCR5-low chemokine receptor profile that shifts responsiveness from lymphoid follicle homing toward CXCL12 (refs. 39,40), which, in active MS lesions, redistributes from the abluminal to the luminal endothelial surface41 via CXCR7-mediated scavenging42. This potentially establishes a chemotactic gradient favoring CNS-directed B cell recruitment. The combination of pre-relapse transcriptional perturbation, MS genetic risk enrichment and EBNA-2-associated regulation within ABC and activated memory B cell states (Figs. 3f and 6b) suggests that these subsets are probable pathogenic intermediates rather than bystanders of systemic inflammation. In absolute terms, gp350+ B cells remained rare, constituting less than 0.2% of CD19+ cells during remission before increasing by more than 40% pre-relapse (Fig. 4g). ABC-like cells similarly represented a minority of the total B cell compartment, yet their expansion was observed in most paired samples (Fig. 4d). Instigator cells may not need to be abundant to exert outsized effects when occupying privileged antigen-presenting positions. Among the strongest transcriptional perturbations were observed in CD14+ and CD16+ monocytes. These were dominated by interferon-stimulated and antigen presentation genes including SP140 and IRF8, both MS risk loci. Because monocytes are not canonical EBV reservoirs, the enrichment of EBV-associated host genes in these cells is suggestive of paracrine signaling: type I interferons and inflammatory cytokines produced by virally activated B cells propagate to bystander myeloid populations43. The tight coupling of monocyte and B cell activation in the pre-relapse window supports feedforward amplification in which virus-driven activity in a minority of B cells elicits disproportionate systemic consequences, consistent with earlier reports that relapse is preceded by interferon-associated signatures44. Host genetics likely modulate the threshold at which this sequence becomes pathogenic. EBNA-2 binds to regulatory elements near MS susceptibility loci and remodels enhancer landscapes24,25,26, and our finding that EBNA-2-bound genes are enriched in pre-relapse ABC transcriptional gene sets implies that susceptibility alleles may amplify cellular responses to episodic EBV reactivation cues. EBNA-2 transcripts were not significantly elevated in pre-relapse B cells; this may reflect transient EBNA-2 expression kinetics relative to durable downstream effects on host enhancer landscapes or indicate that LMP-1-driven signaling operates partly independently of EBNA-2 in this context. Under this logic, the same viral event that resolves innocuously in most individuals may drive pathogenic immune escalation in those with permissive chromatin architectures, effectively lowering the threshold for CNS-directed inflammation. Among EBNA-2 gene sets, type 1-derived sets showed more consistent enrichment than type 2-specific sets (Fig. 6a,b), possibly reflecting predominance of type 1 EBV in Western populations. Whether this reflects genuine allelic differences in enhancer remodeling at MS risk loci requires further study. These observations collectively support a two-component model of relapse initiation: a time-limited perturbation, such as systemic infection, physiological stress or a stochastic event that provokes EBV reactivation within latently infected memory B cells, and a genetically and epigenetically defined susceptibility state that governs the magnitude and persistence of the ensuing immune response. When both components co-occur, interferon and cytokine cascades may surpass the level required for CNS infiltration to initiate relapse. The selective, coordinated activation within B cell and monocyte compartments is consistent with host responses observed in settings of herpesvirus activity45,46. A reproducible pre-relapse signature carries translational implications. Flow cytometric measurements of ABC-like frequency and gp350 positivity, together with molecular detection of LMP-1 and related EBV-responsive transcripts, represent candidate biomarkers that could complement MRI and serum NfL, although effect sizes are modest, and formal evaluation of sensitivity, specificity and added value over existing markers will require prospective validation with predefined cutoffs. Such markers could be informative in two settings: patients not on B cell depletion, where circulating B cell subsets remain assessable, and the months following depletion, where tracking B cell repopulation could reveal whether the ABC signature re-emerges. Because ABCs are less effectively depleted by anti-CD20 therapy17, they may reconstitute earlier than the broader B cell pool. In this setting, improved monitoring could help personalize care by triaging which patients require more frequent re-dosing. More broadly, the data prioritize nodes between viral signaling and host amplification, particularly LMP-1-driven non-canonical NF-κB signaling and downstream B cell activation pathways, as candidate intervention points. Long-term immunotherapies targeting EBV could reduce the latent reservoir, to potentially shift MS management from suppression toward prevention. Several limitations warrant consideration. Although the temporal precedence and multimodal consistency of these findings support a directional relationship, the observational design cannot establish causality, and at least three alternative models remain consistent with the data. First, subclinical CNS injury could precede and trigger peripheral immune activation and secondary EBV reactivation, although this is less likely given the radiographic stability at pre-relapse sampling. Second, within the B cell pool, naive B cells exhibited the largest pre-relapse transcriptional perturbations despite not being canonical EBV latency reservoirs, consistent with a model in which an upstream trigger acts on the naive compartment first and drives differentiation into activated and ABC-like states, with EBV reactivation emerging as an amplifying rather than an initiating event. Third, generic B cell activation could permissively upregulate EBV transcription as a bystander effect, a possibility that the current analyses cannot fully exclude. Viral RNA abundance was low and required pre-amplification for sensitive RT−qPCR detection, and viral transcripts were not reliably captured in scRNA-seq, likely reflecting the rarity of EBV+ instigator cells and known limitations of standard 3′ workflows47. The discordance between gp350 protein detection by flow cytometry and the absence of significant BLLF1 transcript elevation by RT−qPCR may reflect the longer half-life of surface glycoproteins relative to their encoding mRNAs or could indicate antibody cross-reactivity with a structurally related epitope; orthogonal confirmation would strengthen this finding. Moreover, enrichment of gp350+ events within ABC-like metaclusters represents population-level co-occurrence, not single-cell co-localization; the gp350 assay lacked Fc receptor blocking and EBV− cell line controls to exclude non-specific binding at the observed low positivity rates (approximately 0.2%), and the absence of activation markers in the panel precluded adjustment for generalized B cell activation. Although KSHV LANA served as a non-EBV herpesvirus control, additional controls using host gene sets from other B-lymphotropic herpesviruses would further substantiate EBV specificity. B cell receptor and T cell receptor clonotype analysis was not performed and remains an important direction for establishing clonal relationships between EBV-reactive and autoreactive lymphocytes. Future studies incorporating functional perturbation of the LMP-1 axis (for example, NIK/IKKα blockade) will be needed to determine whether LMP-1-driven signaling is required for ABC expansion or is a correlate of it, and inclusion of non-MS autoimmune comparators such as systemic lupus erythematosus (SLE), in which ABC and EBV dysregulation are also prominent, will be essential to establish disease specificity of the pre-relapse EBV signature. Finally, this study is confined to peripheral blood. Whether analogous virus-driven activation states operate within the CNS, and how peripheral immune priming couples to lesion formation, remain open questions. In summary, these data support a model in which episodic EBV-associated signaling in B cells, manifesting as LMP-1-linked activation and detectable viral antigen, precedes relapse by weeks to months and is accompanied by expansion of ABC-like B cells and systemic interferon-dominant immune priming that, within a permissive genetic background, increases the likelihood of relapse. By resolving these events in the pre-symptomatic window, our findings extend the role of EBV in MS beyond the requirement of initial infection for disease onset to episodic reactivation as a candidate driver of relapse biology. Translation will require higher-sensitivity viral detection, prospective sampling and functional perturbation of the implicated networks; that viral and genetic signals intersect in a defined B cell state provides a tractable next target. Methods Cohort and sample processing Participants were drawn from the CLIMB longitudinal MS cohort at Brigham and Women’s Hospital, in which enrolled patients donate blood at annual or routine clinic visits and undergo standardized MRI27. The Mass General Brigham institutional review board provided ethical approval of this work. All included participants signed an informed consent form. Inclusion criteria required a confirmed diagnosis of MS, absence of corticosteroid use for at least 1 month before sampling and stable DMT exposure across pre-relapse and remission samples. Any exceptions are noted for individual assays. Pre-relapse samples were identified retrospectively. Specifically, these were specimens drawn at routine visits from patients who were clinically stable at the time of collection but who subsequently experienced a relapse (within 90 days for scRNA-seq and bulk RNA-seq samples and within 150 days for flow cytometry and RT−qPCR samples). Relapse samples were collected at or after clinical onset; remission samples were drawn at least 100 days after any prior relapse and at least 180 days before any subsequent relapse, with complete clinical and radiographic stability. PBMCs were isolated by Ficoll gradient, cryopreserved in 10% DMSO/90% FBS and stored in vapor-phase liquid nitrogen. Matched healthy controls were included under identical processing conditions. Single-cell RNA-seq Cryopreserved PBMC samples were thawed and prepared as single-cell suspensions; viability exceeded 85% for all samples used for scRNA-seq profiling. Suspensions were combined into pools containing samples from up to three donors (some pools included samples from unrelated studies, which were identified by genotype during demultiplexing and excluded). Approximately 10,000 viable cells were targeted per sample, and each pool was loaded into one well of a Chromium Next GEM Single Cell 3′ v3.1 chip (10x Genomics). In total, 52 PBMC vials representing 42 unique patient timepoints were multiplexed across 28 microfluidics wells and processed in eight batches (Supplementary Table 4). After GEM recovery and library construction according to the manufacturer’s protocol, libraries were sequenced on the Illumina NovaSeq 6000 platform to a median of 34,930 mean reads per cell per library (range, 21,902−89,643). The resulting FASTQ files were aligned to the GRCh38-2020-A reference transcriptome and processed into count matrices with Cell Ranger (version 7.2.0, 12 pools; version 8.0.1, 16 pools). BAM files from each multiplexed pool were demultiplexed to donors with Vireo48 using cellSNP49 allele counts and donor genotypes (gVCFs) from whole-genome sequencing (WGS); the five single-donor pools were assigned directly. Downstream analyses were performed in Seurat version 5 (refs. 50,51). Cells were filtered using fixed and per-library adaptive quality control thresholds (scater); doublets were removed as the union of genotype-based (Vireo) and expression-based (scDblFinder52) calls; and residual low-quality clusters were removed by iterative clustering (Extended Data Fig. 1). Final Harmony-integrated clusters were labeled with their majority Azimuth51 annotation from the cellular indexing of transcriptomes and epitopes by sequencing (CITE−seq) PBMC reference when the top label comprised more than 80% of cells in the cluster; remaining clusters were annotated manually after marker inspection (Extended Data Fig. 2). Surface protein abundances were imputed from the reference antibody-derived tag (ADT) panel during Azimuth reference mapping. Differential and perturbation analysis Cell-state-specific transcriptional shifts between conditions were quantified with scDist29. For each cell state, scDist was fit on log-normalized expression with donor as a random effect. The method estimates per-principal-component condition effects (20 principal components) with linear mixed models, applies empirical Bayes shrinkage to the estimated effects and summarizes the posterior distribution of the between-condition distance D. We report the posterior median with 95% credible intervals (2.5th−97.5th posterior percentiles) and P values from scDist’s simulation-based test of the null hypothesis D = 0, Benjamini−Hochberg adjusted across cell states. Cell states were retained if they comprised ≥50 cells (PBMC atlas) or ≥20 cells (B cell subsets) and ≥2 donors per condition. Specificity was assessed by split-null analyses in which donors within a single condition were randomly partitioned into pseudogroups and reanalyzed with scDist (100 iterations; Extended Data Fig. 4). Gene-level scDist effect sizes were obtained by projecting the per-principal-component condition effects into gene space and were used to rank genes for enrichment analyses. Differential expression within each B cell subset used NEBULA53 (NBGMM; raw counts were modeled as a function of condition with a donor random intercept and library size offset; genes expressed in ≥10% of cells in either condition), with Benjamini−Hochberg correction within subset (FDR < 0.05). Gene set enrichment employed pre-ranked fgsea54 and over-representation analyses using custom EBV factor−host gene modules34 and EBNA-2 ChIP−seq-derived gene sets. Molecular Signatures Database (MSigDB)55 version 2023 hallmark sets were used to construct host activation control sets. Patient-aware scDRS, which quantifies enrichment of MS genetic risk-associated expression in individual cells, used MS GWAS summary statistics from the International Multiple Sclerosis Genetics Consortium18. MAGMA version 1.10 gene-level analysis defined the disease gene set (top 1,000 genes), and scDRS computed per-cell normalized disease relevance scores with total counts, detected genes and mitochondrial percentage as covariates. Scores were averaged per sample and cell state for the PBMC atlas and per donor and subset for the B cell object (combinations with ≥5 cells retained), and cell state estimated marginal means (EMMs) were obtained from a cell-number-weighted linear mixed model with donor random intercepts (two-sided t-tests, Benjamini−Hochberg correction across cell states). B cell reclustering and annotation Cells annotated to B lineage states in the PBMC atlas (TCL1A+ naive B, naive B, intermediate B, switched memory B and plasmablast; 33,847 cells) were extracted and reclustered. Within the B cell compartment, normalization, variable feature selection, scaling and PCA were recomputed, and cells were clustered on the first 30 principal components. Seven first-pass contaminant clusters (T/natural killer, monocyte, red blood cell, dendritic cell and platelet) were removed and the workflow was repeated, followed by Harmony integration across batches, clustering and UMAP embedding; one residual T-cell-like cluster was removed at this stage. The remaining clusters were manually annotated into 14 subsets (transitional/immature-like, TCL1A+ naive, naive, anergic/naive-leaning, interferon-stimulated naive, switched memory, activated switched memory, activated memory, early activated memory, atypical/ABC-like memory, CD11c++ activated memory, ABC, plasmablast and cycling plasmablast) on the basis of canonical RNA markers (Fig. 3b and Extended Data Fig. 6b), with Azimuth-imputed surface protein abundances providing orthogonal support for the subset assignments. Viral module construction EBV factor−host gene modules were derived from the EBV transcriptome atlas of Arvey et al.34, in which EBV and human gene expression were profiled in lymphoblastoid cell lines. For each of the 79 annotated EBV genes, Spearmanʼs rank correlation coefficients were computed against each of 22,467 human genes across the 201 lymphoblastoid sample atlas. Host genes were then ranked in descending order of correlation with each EBV gene, and a fixed number of the most positively correlated host genes were retained as that EBV factorʼs module. No correlation magnitude or significance threshold was applied, so all 79 modules are of equal size, and enrichment testing is rank based rather than threshold based. Modules were annotated post hoc by the lifecycle stage of the EBV gene from which they derive (latent (n = 10), early lytic (n = 32), leaky late (n = 15), true late (n = 16) and unassigned (n = 6)), and this annotation was used for grouping and display of enrichment results. The resulting gene set collection was used for over-representation and GSEA testing of B cell differential expression rankings and, for bulk RNA-seq, as input to module scoring. Multiparameter flow cytometry Paired pre-relapse and remission PBMCs (n = 23 patients, 46 paired samples) were thawed and stained with a viability dye and a deep B cell phenotyping panel that included CD19, CD20, CD11c, CD21, CD27, CD38, CD5, CD18, CD172a, CXCR4 and EBV gp350. Cells were acquired on an LSRFortessa flow cytometer (BD Biosciences). Gating was performed sequentially on lymphocytes, singlets, viable cells and CD3−CD19+ B cells. Compensated marker intensities were arcsinh transformed (cofactor 150). A 10 × 10 self-organizing map was trained with FlowSOM56 on core marker intensities (robust median/median absolute deviation (MAD) scaled; 1,000 cells sampled per file) and grouped into 16 metaclusters by consensus metaclustering and then projected onto all cells. One metacluster of likely non-B-cell debris was removed. Core markers defining metaclusters included CD20, CD11c, CD21, CD27, CD18, CD38, CD172a, CD5 and CXCR4; additional overlay markers were used for phenotypic characterization; and gp350 positivity thresholds were defined per acquisition day from FMO controls at a 1% false-positive rate. Because the panel did not include intracellular T-bet or Zeb2, the CD11c+CD21lo metacluster is designated ‘ABC-like’ to reflect surface-phenotypic rather than transcription-factor-confirmed identification. Within-patient changes in metacluster frequencies were analyzed using two-sided paired Wilcoxon signed-rank tests. gp350+ frequencies were modeled with a binomial generalized linear mixed model (GLMM) (gp350+ counts among CD19+ B cells were modeled as a function of condition and acquisition day, with patient random intercepts). Marginal condition estimates were obtained with EMMs. Temporal dynamics were assessed by re-fitting the model with pre-relapse samples binned by time before relapse onset (0−90 days and 90−150 days) versus remission. Supplementary Table 5 lists the flow cytometry assay details, and Supplementary Table 6 lists the antibodies used for both multiparameter flow cytometry assay and bulk CD19+ cell isolation by FACS. Bulk RNA-seq RNA-seq libraries from FACS-sorted CD19+ B cells were sequenced on NovaSeq 6000 and NovaSeq X Plus instruments (Illumina; paired-end 150-bp reads). Reads were quantified against a decoy-aware Salmon index built with k = 31 from transcript sequences extracted with gffread from the same GRCh38 reference used for scRNA-seq (GRCh38-2020-A), using the full genome as decoy sequence. Quantification used salmon quant (Salmon version 1.12.0) with automatic library type detection and GC-bias and sequence-bias correction. Transcript-level estimates were summarized to genes with tximport (countsFromAbundance = ‘lengthScaledTPM’) using a transcript-to-gene map derived from the same annotation. Immunoglobulin and T cell receptor V, D, J and C segment genes were removed before normalization to prevent clonal transcript levels from distorting per-sample scaling. Genes were filtered with edgeR::filterByExpr, and libraries were normalized by the trimmed mean of M-values (TMM) method; log2 counts per million were computed with a prior count of 1. For the LMP-1 host response analysis, genes were ranked within each sample (singscore::rankGenes), and each sample was scored against the LMP-1 host module of the EBV factor−host co-expression atlas by single-sample rank-based enrichment (singscore::simpleScore, unidirectional up-set). Specificity of the score was assessed against 500 random gene sets matched to the LMP-1 module on mean expression in 5% quantile bins; the observed across-sample standard deviation and the donor intra-class correlation were compared with the empirical null distributions from these matched gene sets. Supplementary Table 7 lists the bulk RNA-seq samples. RT−qPCR FACS-sorted CD19+ B cells from 30 patients with MS with paired pre-relapse and remission samples and radiographic evidence of disease activity (at least one new or enlarging gadolinium-enhancing lesion) were subjected to RNA extraction (Qiagen, RNeasy Mini Kit), reverse transcription and targeted pre-amplification using a pool of TaqMan assays to enhance sensitivity for low-abundance viral transcripts. RT−qPCR was performed across four 384-well plates targeting EBV factors LMP-1, BZLF1, BLLF1 (gp350), BZLF2 (gp42), EBER1, EBER2, EBNA-1, EBNA-2 (type 1 and type 2), EBNA-3A, LMP-2B, BALF1, BGLF4, BGLF5, BNLF2a and BRLF1, with B2M as an endogenous control. No-template control wells were excluded before analysis. Five donors (MS030, MS047, MS050, MS085 and MS086) were excluded for documented assay quality control failure as was any donor for whom more than 50% of wells were unreliable. Targets for which more than 80% of wells were unreliable in either condition were excluded (LMP-2B). Only donors retaining both a pre-relapse and a remission sample after these exclusions were analyzed, leaving 25 patients and 50 samples. Technical replicate wells were averaged within each donor, condition and target; residual missing values were imputed with the condition mean for that target. Expression was measured relative to B2M as ΔCt = Ct(target) − Ct(B2M) and plotted on the 40 − ΔCt scale, on which higher values indicate greater abundance. Unsupervised structure was assessed by PCA of the per-sample ΔCt matrix (prcomp, centered, not scaled). The association between each transcript and pre-relapse status was assessed separately for each target by univariable logistic regression of pre-relapse status on 40 − ΔCt (glm, binomial family), with profile-likelihood 95% CIs. Two-sided Wald P values are reported unadjusted alongside Benjamini−Hochberg-adjusted P values. Supplementary Table 8 lists the qPCR sample details, and Supplementary Table 9 lists the primer−probe panel. EBNA-2 and KSHV integration EBNA-2 promoter occupancy gene sets were derived from external EBNA-2 ChIP−seq datasets for EBV types 1 and 2. EBNA-2 activation and repression signatures were defined from intersected external RNA-seq studies. KSHV LANA promoter occupancy gene sets served as non-EBV viral controls. Enrichment was assessed using pre-ranked fgsea on scDist-derived gene rankings for each B cell subset. Statistical analyses The unit of study for all statistical comparisons was the participant blood draw; one sample per participant per disease state was analyzed, so n samples correspond to n biologically independent participants per group unless stated otherwise. No technical replicates were treated as independent observations; where technical replicates existed, they were combined before analysis: specifically, replicate scRNA-seq vials of the same sample were pooled at the cell level, and technical qPCR wells were averaged within donor, condition and target. All P values are two-tailed unless stated otherwise; scDistʼs simulation-based test of D = 0 rejects for large values of the test statistic and is one-sided by construction. Group comparisons employed Wilcoxon tests, univariable logistic regression or (generalized) linear mixed models with donor random intercepts. Effect sizes are expressed as log2 fold changes for expression or odds ratios for logistic and binomial mixed models. Analyses were performed in R version 4.5.2, MAGMA version 1.10 was run as a standalone binary and scDRS was run in Python 3.12.3. Salmon quantification was run under WSL Ubuntu. Reporting summary Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article. Data availability Raw sequencing reads (FASTQ) and processed, deidentified count matrices for the single-cell and bulk RNA sequencing data generated in this study, together with the sample-level metadata required to reproduce the analyses, are deposited in the Gene Expression Omnibus (GEO) (GSE344578) and made publicly available without restriction at the time of publication. Flow cytometry FCS files, qPCR amplification data, genotype-based demultiplexing assignments (Vireo donor_ids.tsv) and derived analysis objects are deposited in Zenodo (DOI pending). WGS reads and individual-level clinical data from the CLIMB cohort are potentially identifying and are governed by participant informed consent and Mass General Brigham institutional review board policies; these data are, therefore, subject to controlled access. WGS data were used solely for genotype-based sample demultiplexing, and the resulting donor assignments are publicly deposited, so all results are reproducible without access to the genomes. Deidentified individual-level data will be made available to qualified investigators for non-commercial academic research from the date of publication and for a minimum of 5 years thereafter, contingent on a data use agreement and any institutional review board approvals required by the host and requesting institutions. Requests should be directed to the corresponding author (T.C.; tchitnis@bwh.harvard.edu), who will coordinate review with Mass General Brigham data governance officials; requesters can expect an initial response within 4 weeks. Publicly available datasets used in this study include MS GWAS summary statistics from the International Multiple Sclerosis Genetics Consortium, the EBV−human expression atlas of Arvey et al.34, EBNA-2 and KSHV LANA ChIP−seq data (GEO: GSE246060, GSE29498 and GSE56144) and the Azimuth human PBMC reference. Code availability Custom code used to process the sequencing data and to generate the analyses and figures reported in this study is available through a public GitHub repository (https://github.com/TNRC-MGB/RRMS-Biomarkers) and archived with a citable DOI via Zenodo (pending), without restriction, under an open source (MIT) license. References Walton, C. et al. Rising prevalence of multiple sclerosis worldwide: insights from the Atlas of MS, third edition. Mult. Scler. 26, 1816–1821 (2020). Jaijakul, S., Salazar, L., Wootton, S. H., Aguilera, E. & Hasbun, R. The clinical significance of neutrophilic pleocytosis in cerebrospinal fluid in patients with viral central nervous system infections. Int. J. Infect. Dis. 59, 77–81 (2017). Thebault, S. et al. Serum neurofilament light chain predicts long term clinical outcomes in multiple sclerosis. Sci. Rep. 10, 10381 (2020). Canto, E. et al. Association between serum neurofilament light chain levels and long-term disease course among patients with multiple sclerosis followed up for 12 years. JAMA Neurol. 76, 1359–1366 (2019). Chitnis, T. et al. Neurofilament light chain serum levels correlate with 10-year MRI outcomes in multiple sclerosis. Ann. Clin. Transl. Neurol. 5, 1478–1491 (2018). Akesson, J. et al. Proteomics reveal biomarkers for diagnosis, disease activity and long-term disability outcomes in multiple sclerosis. Nat. Commun. 14, 6903 (2023). Chitnis, T. et al. Inflammatory and neurodegenerative serum protein biomarkers increase sensitivity to detect clinical and radiographic disease activity in multiple sclerosis. Nat. Commun. 15, 4297 (2024). Hauser, S. L. et al. B-cell depletion with rituximab in relapsing-remitting multiple sclerosis. N. Engl. J. Med. 358, 676–688 (2008). Ramesh, A. et al. A pathogenic and clonally expanded B cell transcriptome in active multiple sclerosis. Proc. Natl Acad. Sci. USA 117, 22932–22943 (2020). Lanz, T. V. et al. Clonally expanded B cells in multiple sclerosis bind EBV EBNA1 and GlialCAM. Nature 603, 321–327 (2022). Jenks, S. A. et al. Distinct effector B cells induced by unregulated Toll-like receptor 7 contribute to pathogenic responses in systemic lupus erythematosus. Immunity 49, 725–739 (2018). Portugal, S., Obeng-Adjei, N., Moir, S., Crompton, P. D. & Pierce, S. K. Atypical memory B cells in human chronic infectious diseases: an interim report. Cell Immunol. 321, 18–25 (2017). SoRelle, E. D., Reinoso-Vizcaino, N. M., Horn, G. Q. & Luftig, M. A. Epstein-Barr virus perpetuates B cell germinal center dynamics and generation of autoimmune-associated phenotypes in vitro. Front. Immunol. 13, 1001145 (2022). Jelcic, I. et al. T-bet+ CXCR3+ B cells drive hyperreactive B-T cell interactions in multiple sclerosis. Cell Rep. Med. 6, 102027 (2025). Claes, N. et al. Age-associated B cells with proinflammatory characteristics are expanded in a proportion of multiple sclerosis patients. J. Immunol. 197, 4576–4583 (2016). SoRelle, E. D. et al. Early multiple sclerosis activity associated with TBX21+CD21loCXCR3+ B cell expansion resembling EBV-induced phenotypes. JCI Insight 10, e188543 (2025). El Mahdaoui, S. et al. CD11c+ B cells in relapsing-remitting multiple sclerosis and effects of anti-CD20 therapy. Ann. Clin. Transl. Neurol. 11, 926–937 (2024). International Multiple Sclerosis Genetics Consortium. Multiple sclerosis genomic map implicates peripheral immune cells and microglia in susceptibility. Science 365, eaav7188 (2019). Hollenbach, J. A. & Oksenberg, J. R. The immunogenetics of multiple sclerosis: a comprehensive review. J. Autoimmun. 64, 13–25 (2015). Smets, I. et al. Multiple sclerosis risk variants alter expression of co-stimulatory genes in B cells. Brain 141, 786–796 (2018). Cenit, M. C. et al. STAT3 locus in inflammatory bowel disease and multiple sclerosis susceptibility. Genes Immun. 11, 264–268 (2010). Parnell, G. P. & Booth, D. R. The multiple sclerosis (MS) genetic risk factors indicate both acquired and innate immune cell subsets contribute to MS pathogenesis and identify novel therapeutic opportunities. Front. Immunol. 8, 425 (2017). Skarlis, C., Papadopoulos, V., Raftopoulou, S., Mavragani, C. P. & Evangelopoulos, M. E. B-cell activating factor gene variants in multiple sclerosis: possible associations with disease susceptibility among females. Clin. Immunol. 257, 109847 (2023). Keane, J. T. et al. The interaction of Epstein-Barr virus encoded transcription factor EBNA2 with multiple sclerosis risk loci is dependent on the risk genotype. EBioMedicine 71, 103572 (2021). Mechelli, R. et al. Epstein-Barr virus genetic variants are associated with multiple sclerosis. Neurology 84, 1362–1368 (2015). Mechelli, R. et al. A disease-specific convergence of host and Epstein−Barr virus genetics in multiple sclerosis. Proc. Natl Acad. Sci. USA 122, e2418783122 (2025). Gauthier, S. A., Glanz, B. I., Mandel, M. & Weiner, H. L. A model for the comprehensive investigation of a chronic autoimmune disease: the multiple sclerosis CLIMB study. Autoimmun. Rev. 5, 532–536 (2006). Korsunsky, I. et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat. Methods 16, 1289–1296 (2019). Nicol, P. B. et al. Robust identification of perturbed cell types in single-cell RNA-seq data. Nat. Commun. 15, 7610 (2024). Ersing, I. et al. A temporal proteomic map of Epstein-Barr virus lytic replication in B cells. Cell Rep. 19, 1479–1493 (2017). Chang, Y. et al. Induction of Epstein-Barr virus latent membrane protein 1 by a lytic transactivator Rta. J. Virol. 78, 13028–13036 (2004). Carter, K. L., Cahir-McFarland, E. & Kieff, E. Epstein-barr virus-induced changes in B-lymphocyte gene expression. J. Virol. 76, 10427–10436 (2002). Ahsan, N., Kanda, T., Nagashima, K. & Takada, K. Epstein-Barr virus transforming protein LMP1 plays a critical role in virus production. J. Virol. 79, 4415–4424 (2005). Arvey, A. et al. An atlas of the Epstein-Barr virus transcriptome and epigenome reveals host-virus regulatory interactions. Cell Host Microbe 12, 233–245 (2012). Buschle, A. & Hammerschmidt, W. Epigenetic lifestyle of Epstein-Barr virus. Semin. Immunopathol. 42, 131–142 (2020). Price, A. M. & Luftig, M. A. To be or not IIb: a multi-step process for Epstein-Barr virus latency establishment and consequences for B cell tumorigenesis. PLoS Pathog. 11, e1004656 (2015). Bjornevik, K. et al. Longitudinal analysis reveals high prevalence of Epstein-Barr virus associated with multiple sclerosis. Science 375, 296–301 (2022). Sun, Y. et al. Epstein-barr virus latent membrane protein 1 targets cIAP1, cIAP2 and TRAF2 for proteasomal degradation to activate the non-canonical NF-κB pathway. PLoS Pathog. 22, e1013898 (2026). Fujiwara, M., Kondo, R., Sugiyama, Y., Maruyama, M. & Nishikimi, A. Increased Fascin1 and Pak1 expressions enhance age-associated B-cell actin cytoskeleton remodeling and motility. Cell Biochem. Funct. 43, e70090 (2025). Vidal-Pedrola, G. et al. Characterization of age-associated B cells in early drug-naive rheumatoid arthritis patients. Immunology 168, 640–653 (2023). McCandless, E. E. et al. Pathological expression of CXCL12 at the blood-brain barrier correlates with severity of multiple sclerosis. Am. J. Pathol. 172, 799–808 (2008). Cruz-Orengo, L. et al. CXCR7 influences leukocyte entry into the CNS parenchyma by controlling abluminal CXCL12 abundance during autoimmunity. J. Exp. Med. 208, 327–339 (2011). Thorley-Lawson, D. A. & Gross, A. Persistence of the Epstein−Barr virus and the origins of associated lymphomas. N. Engl. J. Med. 350, 1328–1337 (2004). Comabella, M. et al. A type I interferon signature in monocytes is associated with poor response to interferon-β in multiple sclerosis. Brain 132, 3353–3365 (2009). SoRelle, E. D. et al. Epstein-Barr virus reactivation induces divergent abortive, reprogrammed, and host shutoff states by lytic progression. PLoS Pathog. 20, e1012341 (2024). Munz, C. Latency and lytic replication in Epstein−Barr virus-associated oncogenesis. Nat. Rev. Microbiol. 17, 691–700 (2019). Younis, S. et al. Epstein-Barr virus reprograms autoreactive B cells as antigen-presenting cells in systemic lupus erythematosus. Sci. Transl. Med. 17, eady0210 (2025). Huang, Y., McCarthy, D. J. & Stegle, O. Vireo: Bayesian demultiplexing of pooled single-cell RNA-seq data without genotype reference. Genome Biol. 20, 273 (2019). Huang, X. & Huang, Y. Cellsnp-lite: an efficient tool for genotyping single cells. Bioinformatics 37, 4569–4571 (2021). Hao, Y. et al. Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat. Biotechnol. 42, 293–304 (2024). Hao, Y. et al. Integrated analysis of multimodal single-cell data. Cell 184, 3573–3587 (2021). Germain, P. L., Lun, A., Garcia Meixide, C., Macnair, W. & Robinson, M. D. Doublet identification in single-cell sequencing data using scDblFinder. F1000Res. 10, 979 (2021). He, L. et al. NEBULA is a fast negative binomial mixed model for differential or co-expression analysis of large-scale multi-subject single-cell data. Commun. Biol. 4, 629 (2021). Korotkevich, G. et al. Fast gene set enrichment analysis. Preprint at bioRxiv https://doi.org/10.1101/060012 (2021). Liberzon, A. et al. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 1, 417–425 (2015). Van Gassen, S. et al. FlowSOM: using self-organizing maps for visualization and interpretation of cytometry data. Cytometry A 87, 636–645 (2015). Acknowledgements We thank the participants and staff of the CLIMB study, M. Polgar-Turcsanyi for data management and colleagues in the Translational Neuroimmunology Research Center for critical discussions. Funding This work was supported by the US Department of Defense (to T.C.), the Water Cove Charitable Foundation (to T.C.), the Cindy Larsen-Chugg Chair (to T.C.) and a Broad Institute Founders Grant (to B.E.G.) and by generous support from George and Sandra K. Schussel to B.E.G. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript. Author information Authors and Affiliations Contributions T.C. and D.A.K. conceived and designed the work. All authors made contributions to the acquisition and/or analysis of the data. T.C., D.A.K., D.C., L.E.S., S.S. and B.E.G. contributed to the interpretation of data. D.A.K. and T.C. drafted the work or substantively revised it. All authors have approved the submitted version. Corresponding author Ethics declarations Competing interests T.C. has received research support from the US Department of Defense, the National Institutes of Health, the National Multiple Sclerosis Society, Genentech Novartis, Octave Biosciences, Sanofi and Tiziana Life Sciences. She has served as a consultant/advisory board member for Genentech, Lilly, Novartis, Octave Biosciences and Sanofi. H.L.W. has received research support from the Cure Alzheimer’s Fund, the National Institutes of Health and the National Multiple Sclerosis Society. He has served on the advisory board for Sanofi and for Tiziana Life Sciences. Peer review Peer review information Nature Medicine thanks Christian Münz and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editor: Jerome Staal, in collaboration with the Nature Medicine team. Additional information Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Extended data Extended Data Fig. 1 PBMC Single-cell RNA-seq quality control. PBMC Single-cell RNA-seq quality control, doublet removal, and batch correction. a, PBMC scRNA-seq and whole-genome sequencing (WGS) data were processed through raw count generation and sample demultiplexing (Cell Ranger and Vireo, respectively), per-library quality control (QC) filtering (scater), expression-based doublet detection (scDblFinder), and combined genotype- and expression-based doublet removal with exclusion of unassigned or mis-assigned cells, followed by first-pass clustering. b, Dot plot showing initial cluster evaluation based on expression of canonical lineage markers for T cells (TRAC, CD3D, CD3E, IL7R, CCR7), NK/cytotoxic cells (NKG7, GNLY, PRF1, GZMB, FCER1G), B cells (CD19, MS4A1, CD79A, BANK1), plasma cells (MZB1, JCHAIN, XBP1, IGKC), monocytes (CD14, S100A8, S100A9, FCGR3A, LYZ), RBC/platelets (HBB, HBA1, HBA2, PPBP, PF4, ITGA2B), vascular cells (COL14A1, ANGPT1, EGFL7, IGFBP7), potential ambient signatures (MALAT1, JUN), and remaining doublet score (scDblFinder_score). Dot size represents the percentage of cells expressing each gene and color intensity indicates scaled average expression (z-scored across clusters). Accompanying panels show cell totals per cluster, total UMI counts (nCount), number of detected genes (nFeature), and percentage of mitochondrial reads (pct_mt). Orange-highlighted rows denote clusters flagged for removal. c, UMAP embedding of first-pass clustering results with numbered cluster identities. d, UMAP highlighting the 2.68% of cells removed from cluster-level QC, with annotations indicating non-PBMC contaminants (endothelial/vascular and erythroid), monocyte/T cell ambient RNA clusters, and platelet-containing doublets (platelet + monocyte, platelet + T cell, and platelet + B cell). e, Heatmap displaying the percentage of cells removed per sample, stratified by disease condition (Healthy, PreRelapse, Relapse, and Remission). f, UMAP embedding of second-pass clustering after removal of prior-flagged clusters. g, UMAP colored by sample to assess batch effects prior to correction. h, UMAP embedding after Harmony batch correction, with renumbered cluster identities. The final atlas comprised 281,799 cells (127,090 from 21 healthy donors and 154,709 from 15 donors with MS), derived from 52 individual PBMC vials. Extended Data Fig. 2 Annotation of PBMC clusters to known cell types. a, Stacked bar chart showing the fraction of cells in each Harmony-integrated cluster (y-axis) assigned to their top predicted label after reference mapping to a published PBMC dataset (Azimuth51; Hao et al.50, Cell, 2021). Bars are colored by predicted cell type identity. The dashed vertical line at 0.8 denotes the threshold above which clusters were directly assigned their top Azimuth label; clusters falling below this threshold were manually annotated based on canonical marker gene expression. b, Confusion matrix comparing the final manual PBMC annotations used in this study (y-axis) with Azimuth-predicted labels (x-axis). Color intensity indicates the percentage of cells in each manually defined category assigned to a given Azimuth label, with darker blue denoting higher concordance. c, UMAP embedding of the final annotated PBMC atlas with cell type labels. Extended Data Fig. 3 Sample composition and covariate distribution across the PBMC atlas. a, Bar chart showing the total number of cells recovered per sample (y-axis, log10 scale), faceted by disease condition (Healthy, PreRelapse, Relapse, and Remission). The horizontal dashed line indicates the minimum target of 1,500 cells per library at microfluidics loading. b, Stacked bar chart displaying the cell type composition (cell fraction, y-axis) of each sample (x-axis), colored by cell annotation. c-f, UMAP embeddings of the final PBMC atlas colored by experimental day (eight batches, A1-A2 and B1-B6, each corresponding to a single 10x Chromium chip) (c), donor identity (d), disease-modifying therapy (Avonex, Copaxone, Rebif, or Untreated; note, both healthy donors and untreated MS donors are represented in the Untreated category) (e), and age (continuous scale, years) (f), demonstrating that cells intermingle across covariates without evident confounding structure. Extended Data Fig. 4 Split-null validation of scDist transcriptomic distances. a,b, Representative null iterations from split-null analyses in which donors within a single condition were randomly assigned to two pseudogroups and compared using scDist with donor as a random effect. a, Healthy versus healthy split-null. b, MS remission versus MS remission split-null. Each subplot displays one iteration, with cell types ranked along the y-axis by scDist distance (x-axis); horizontal bars indicate 95% credible intervals. Dot color represents −log10(Benjamini-Hochberg-adjusted P value), with values capped at 3 (adjusted P ≤ 0.001). Across 100 null iterations, most cell types yield near-zero distances, consistent with the absence of a true transcriptomic perturbation between pseudogroups. c, Histogram of the maximum scDist distance across all PBMC cell types for each healthy split-null iteration (gray bars). The red vertical line indicates the observed maximum pre-relapse scDist distance, which falls far outside the null distribution. Extended Data Fig. 5 Paired-sample sensitivity analysis of scDist transcriptomic distances. a,b, scDist analysis restricted to MS donors with paired relapse-associated and remission samples, comparing relapse versus remission (a) and pre-relapse versus remission (b). Cell types are ranked along the y-axis by scDist distance (x-axis); horizontal bars indicate 95% credible intervals. Dot color represents −log10(Benjamini-Hochberg-adjusted P value). Results recapitulate the primary analysis, with the largest transcriptomic perturbations observed in myeloid, innate lymphoid, and B cell populations. Extended Data Fig. 6 Quality control and granular analysis of B cell subsets. a, UMAP embedding of B cells (TCL1A+ Naive B, Naive B, Intermediate B, Switched Memory B, and Plasmablast) isolated from the PBMC atlas after first-pass re-clustering. b, Dot plot evaluating B cell cluster identity using the canonical lineage marker panel described in Extended Data Fig. 1b. Dot size represents the percentage of cells expressing each gene and color intensity indicates scaled average expression (z-scored across clusters). Accompanying panels show cell totals per cluster, total UMI counts (nCount), number of detected genes (nFeature), and percentage of mitochondrial reads (pct_mt). Orange-highlighted rows denote clusters flagged for removal, with annotations indicating the basis for exclusion (low quality, monocyte doublet, T/NK contamination, RBC contamination, DC contamination, or CD4 T contamination). c, UMAP embedding of second-pass clustering after removal of contaminating clusters. d, UMAP colored by experimental batch (A1-A2, B1-B6) post harmony batch correction. e, Final B cell cluster identities (0-14). Extended Data Fig. 7 Differential gene expression across B cell clusters between pre-relapse and remission. Volcano plots showing differential gene expression (pre-relapse versus remission) for each of the B cell subsets identified after re-clustering and manual annotation (see Fig. 3a and Extended Data Fig. 6). Differential expression was computed using NEBULA with a negative binomial gamma mixed model (NBGMM) and donor as a random effect. The x-axis represents log fold change (logFC) and the y-axis represents −log10(FDR), where FDR was calculated by Benjamini-Hochberg correction within each subset. The horizontal dashed line indicates FDR = 0.05. Genes are colored by direction of effect: red, upregulated in pre-relapse; green, upregulated in remission; blue, not significant. Labeled genes (positive logFC only) were selected from the concordant intersection of NEBULA and scDist results within each subset, thus highlighting genes that were both significantly differentially expressed and among the strongest contributors to the cluster-level transcriptomic distance identified by scDist. Genes expressed in ≥10% of cells in either condition were tested; library size was included as an offset. Extended Data Fig. 8 FACS gating strategy and FlowSOM clustering of CD19 + B cells. a, Representative gating strategy for flow cytometry to extract CD19 + B cells for FlowSOM analysis. Sequential gates were applied to isolate lymphocytes (FSC-A versus SSC-A), singlets (FSC-A versus FSC-H), live cells (viability dye versus SSC-A), CD3− cells (CD3 versus SSC-A), and CD19 + B cells (CD19 versus SSC-A). CD19+ events (with compensation applied) were imported into R for downstream analysis. b, UMAP projections of CD19 + B cells coloured by arcsinh-transformed (cofactor = 150) expression of ten surface markers: CD5, CD11C, CD18, CD19, CD20, CD21, CD27, CXCR4, CD38, and CD172 (color scale clamped at 0-7). c, UMAP colored by FlowSOM metacluster identity. Cells were downsampled to 1,000 per sample, robust z-scored (median/MAD), and clustered using a 10 Å~ 10 self-organizing map with consensus metaclustering into 16 clusters; one cluster (cluster 6) was removed as likely non B cell debris. The remaining 15 metaclusters were manually annotated as shown. d, Modelpredicted percentage of gp350+ B cells (y-axis) by condition (pre-relapse vs. remission) stratified by experimental day. Predictions and 95% confidence intervals were derived from estimated marginal means of a binomial generalised linear mixed model (GLMM) with condition and experimental day as fixed effects and donor as a random intercept (n = 46 samples from 23 patients with paired pre-relapse and remission draws). gp350 positivity thresholds were determined per experimental day at a 1% false-positive rate. e, Forest plot of odds ratios (log scale; x-axis) from a binomial GLMM in which pre-relapse samples were further stratified by proximity to relapse onset (0-90 days and 90-150 days before relapse), with experimental day as a fixed effect and donor as a random intercept. Error bars indicate 95% Wald confidence intervals; dashed line marks an odds ratio of 1. Extended Data Fig. 9 EBV host response signatures are significantly enriched within the ABC-developmental trajectory. a, Schematic of the analytical workflow. Human host gene sets associated with individual EBV factors were extracted from a previously published EBV-human expression atlas (Arvey et al.34), retaining the top 500 genes per factor. From this study, genes were ranked by NEBULA Wald statistic (log fold change / standard error; pre-relapse versus remission) within each B cell subset. Gene set enrichment analysis (GSEA) was then performed using fgseaMultilevel (Korotkevich et al.54) to test for enrichment of each EBV factor-derived host gene set in the pre-relapse transcriptomics data. b, Schematic of the EBV-human expression atlas. RNA-seq was performed on 201 lymphoblastoid cell lines (LCLs) with matched EBV and human transcriptomes, yielding host gene signatures associated with latent (Host-EBVlatent) and lytic (Host-EBVlytic) EBV gene expression. c, Reference key mapping individual EBV genes to their lifecycle stage: other (orange), latent (red), early (blue), leaky late (green), and true late (purple). d, GSEA dot plots showing normalized enrichment scores (NES; x-axis) for each EBV factor-derived host gene set (rows, ordered by lifecycle stage as in c) across the 13 B cell subsets. Dot color indicates the -log10(adjusted P value), after Benjamini-Hochberg adjustment. Enrichment was tested in the positive direction only (genes upregulated pre-relapse relative to remission). B cell subsets along the ABC-developmental trajectory show the most prominent enrichment of EBV host response signatures. Extended Data Fig. 10 EBV-host responses implicate early lytic reactivation in pre-relapse B cells. Violin plots showing the distribution of normalized enrichment scores (NES; y-axis) for EBV factor-derived host gene sets grouped by EBV lifecycle stage (x-axis: latent, early, leaky late, true late, other) across B cell subsets. Each dot represents one EBV factor-derived gene set with positive NES (top 500 host genes per factor from the Arvey et al.34, expression atlas; same data as Extended Data Fig. 9d, filtered to NES > 0). Dot color indicates -log10(adjusted P value) on a grey-to-red scale (capped at 5), where P values were adjusted using the Benjamini-Hochberg method. Positive NES indicates upregulation in pre-relapse relative to remission. GSEA was performed using fgseaMultilevel (Korotkevich et al.54) with genes ranked by NEBULA Wald statistic, as described in Extended Data Fig. 9. The highest NES values and most significant enrichments are concentrated among early and leaky late lytic EBV factor signatures rather than true late lytic or latent signatures. Supplementary information Supplementary Information (download PDF ) Tables 1−9 and Fig. 1. Rights and permissions Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/. About this article Cite this article King, D.A., Saxena, S., Caefer, D. et al. EBV reactivation priming of the peripheral immune system in multiple sclerosis relapse. Nat Med (2026). https://doi.org/10.1038/s41591-026-04665-3 Received: Accepted: Published: Version of record: DOI: https://doi.org/10.1038/s41591-026-04665-3