Dimensionality reduction to disentangle the heterogeneity of MS neuropathology

To identify patterns in MS neuropathology, we established a computational workflow that consists of four stages (Fig. 1a). In stage 1, we collected previously generated quantitative and qualitative neuropathological data, predominantly related to the white matter [33]. Stage 2 consisted of processing this data, including transformation, imputation and an initial FAMD (factor analysis of mixed data). We selected the first three dimensions for further analysis and subsequently removed two donors with outlier values who contributed highly to a dimension. In stage 3, FAMD was performed on the definite dataset of 13 variables across 226 donors. Of the final FAMD dimensions, the first three together accounted for almost half of the total variance in the dataset, indicating that closer examination could provide insight into general neuropathological patterns observed in our MS cohort. In stage 4, we used orthogonal data types (clinical [33, 40], neuropathological [33, 58], immunological [13, 14], and genetics) to evaluate and interpret the dimensions.

Fig. 1

Overview of the study and the three independent dimensions. a Graphical outline of this study created with BioRender.com, consisting of 4 stages: (1) collecting previously generated input data, (2) data preparation, (3) exploratory FAMD, with scatter plots showing how donors score on the resulting dimensions, and (4) validation and exploration. b Dot plot displaying the relation between FAMD input variables and the first three dimensions. Dot size indicates the squared cosine, reflecting the proportion of variance in a variable explained by a dimension. Higher values correspond to a better quality of representation of the variable by the dimension. Dot color indicates whether the correlation between the dimension and the variable is positive (red) or negative (blue); color intensity reflects the relative contribution of the variable to the component. c Scatter plots with fitted line and box plots per dimension, depicting the values for the input variables on the Y-axis, for donors ranked according to their score on each dimension on the X-axis. In case of ties, donors were assigned ranks in the order of appearance in the dataset, so that each donor receives a unique rank. The line represents the centered moving average, over a window of 20 observations, with a maximum of 10 missing values. Proportional data is CLR-transformed; unimputed values are shown. Plots are ordered vertically based on the variable’s contribution to the dimension. CLR = centered log ratio; FAMD = factor analysis of mixed data; dim. = dimension; PRS = polygenic risk score; ramif. = ramified; ameb. = ameboid; remyel. = remyelinated (lesion); L. = load; R. = rate; corr. = correlation

First, we explored the relation between input variables and the three dimensions identified in stage 3. Each dimension reflected a distinct neuropathological pattern, in which the variables were represented by and contributing to each dimension to varying extents (Fig. 1b). In addition, we ranked all donors according to their score on each dimension (in ascending order, from low to high) to illustrate how the dimensions were related to the values of the input variables (Fig. 1c). Critically, donors’ scores on one dimension were independent of their scores on the other dimensions—the dimensions are uncorrelated (as shown by the scatter plots in Fig. 1a). A low score on dimension 1 was associated with a predominance of inactive and remyelinated lesions with few HLA+ microglia/macrophages, a high score with a high lesion load, active and mixed lesions populated by foamy microglia, and the presence of nodules and cuffs (perivascular accumulations of leukocytes). Dimension 2 was positively correlated with the proportion of active lesions, the number of reactive sites, and the presence of nodules. Dimension 3 was positively associated with a higher cortical lesion rate and lesion load, a larger proportion of mixed lesions with ramified microglia and inactive lesions, relatively fewer active lesions with foamy microglia, and an absence of cuffs.

Next, we aimed to validate and explore the dimensions. For an overview of the relation between the dimensions and all variables analysed in this study, see Table 1.

Table 1 Overview of the relation between the dimensions and variables analysed in this studyDimensions are associated with year of autopsy and cause of death

The dimensions were not associated with post-mortem delay, pH or brain weight (suppl. Figure 4, Online Resource 1), implying that they were not driven by these technical covariates. Donors who scored high on dimension 2 or 3 were likely to have died in more recent years (dim. 1: p = 0.17, rho = 0.09; 2: p = 3.5 × 10–5, rho = 0.29; 3: p = 0.006, rho = 0.19; suppl. Figure 5, Online Resource 1). There was no association with cause of death, except for dimension 2 (dim. 1: p = 0.12, 2: p = 0.05, 3: p = 0.51; suppl. Figure 6, Online Resource 1); donors who died by legal euthanasia scored higher on this dimension than donors who died of unspecified natural causes (p = 0.02). This finding could reflect an association between year and cause of death, possibly related to Dutch legislation on euthanasia taking effect in 2002. Accordingly, the median year of autopsy was 2013 for donors who died by euthanasia and 2005.5 for those who died from unspecified natural causes.

Dimensions are associated with the clinical manifestations of MS

To investigate the potential clinical correlates of the neuropathological dimensions, we subsequently focused our analysis on data related to demographics and disease experience. There was no significant correlation with sex (dim. 1: p = 0.80; 2: p = 0.80; 3: p = 0.93) or MS clinical phenotype (dim. 1: p = 0.30; 2: p = 0.88; 3: p = 0.30) (Fig. 2a, b), the latter supporting the notion that the historically identified clinical MS phenotypes do not qualitatively differ with regards to the (white matter) pathological features included in our analysis.

Fig. 2

Association between demographic and clinical variables and dimensions. Note that X-axis scale differs among the different plots. a Box plots showing the ranked scores of donors, per dimension, per sex. Sex was known for 214 donors (77 males, 137 females). There are no significant differences between sexes (Mann–Whitney U; dim. 1: p = 0.80; 2: p = 0.80; 3: p = 0.93). b Box plots showing the ranked scores of donors, per dimension, per MS clinical phenotype. MS phenotype was determined for 197 donors (12 relapsing, 119 SP, 66 PP). There are no significant differences (Kruskal–Wallis; dim. 1: p = 0.30; 2: p = 0.88; 3: p = 0.30). In plots in a and b, ties were assigned averaged ranks. c Dot plot (left) showing the correlation between the scores on dimensions 1–3 and age at MS onset for 204 donors, age at death for 214 donors, years till EDSS-6 for 191 donors, and years till death (i.e. duration of disease) for 201 donors. Dot colour indicates the strength of the correlation, dot size the p-value, and a yellow star a significant association. Box plots (right) showing the distribution of the demographic and clinical variables. Significance in c was assessed with Spearman correlation and FDR-adjusted for multiple testing. dim. = dimension; PP = primary progressive; SP = secondary progressive; EDSS = Expanded Disability Status Scale; FDR = False Discovery Rate

Dimension 1 was negatively associated with the years from onset till EDSS-6 (p = 7.7 × 10–8, rho = −0.39), from onset till death (i.e. duration of disease; p = 5.8 × 10–8, rho = −0.39), and age at death (p = 1.6 × 10–12, rho = −0.48); there was no relation with age at MS onset (p = 0.51, rho = −0.06) (Fig. 2c). In addition, the total number of signs and symptoms adjusted for disease duration (the symptom load) was significantly positively correlated with dimension 1 for the general, motor, sensory/autonomic, and cognitive domains (suppl. Figure 7, Online Resource 1). Corresponding with their earlier death, the median age at symptom onset was generally earlier for donors scoring high on dimension 1 (suppl. Figure 8, Online Resource 1). On the whole, a higher score on dimension 1 associates with a more severe disease course of MS.

In contrast, dimension 2 was associated with milder MS, positively correlating with the years till EDSS-6 (p = 0.006, rho = 0.21), till death (p = 4.3 × 10–5, rho = 0.30), and age at death (p = 2.4 × 10–4, rho = 0.26) (Fig. 2c). Age at onset was not associated with the dimension (p = 0.78, rho = 0.02; Fig. 2c). Moreover, for all domains the symptom load was negatively correlated with dimension 2, although not always significantly so (suppl. Figure 7, Online Resource 1). The median age at symptom onset was positively correlated with dimension 2, in line with expectations (suppl. Figure 8, Online Resource 1).

Dimension 3 was not significantly correlated with years till EDSS-6 (p = 0.74, rho = 0.03), age at MS onset (p = 0.31, rho = −0.09), and age at death (p = 0.37, rho = 0.07) (Fig. 2c). However, because donors who scored high on dimension 3 tended to develop MS at a slightly younger age and died at a relatively older age, there was a significant positive correlation with disease duration (p = 0.02, rho = 0.18; Fig. 2c). Interestingly, with regards to the symptom load, the overarching domain might matter: there was a significant negative correlation for the psychiatric and sensory/autonomic domains only (suppl. Figure 7, Online Resource 1). Furthermore, there was a positive association with median age at symptom onset for the cognitive domain (suppl. Figure 8, Online Resource 1).

The identified dimensions correlate with clinical profiles of NBB MS donors, and therefore likely reflect clinically relevant pathological processes. In the following sections, each dimension will be further examined with regards to comorbidities, drug use, cortical neuropathology, potential immune cell involvement, and genetics.

Dimensions correlate with certain comorbidities and use of MS-relevant drugs

To assess associations between our dimensions and comorbidity, we compared the MS donors with and without other diagnoses (i.e. with or without one or more observations of relevant categories and classes); similarly, we compared the donors who did and did not use a particular class of drugs. MS donors with an autoimmune disease scored significantly lower on dimension 1. Donors with a cardiovascular system disease or its subclass hypertension had a higher age at death and scored lower on dimension 1 (suppl. Figure 9, Online Resource 1). Consistently, donors who used cardiovascular system drugs died at an older age and scored lower on dimension 1 (although the latter was not significant after multiple-testing correction; suppl. Figure 10, Online Resource 1), corroborating the quality of our datasets. Type 2 diabetes mellitus was more common in older donors, and donors with this comorbidity scored significantly higher on dimension 2 (suppl. Figure 9, Online Resource 1). Regarding drug therapies relevant for MS, donors using disease-modifying drugs scored higher on dimension 1, died at a younger age, and died more recently. Donors who used MS-relevant immunosuppressive drugs generally scored lower on dimension 2, died at a younger age, and in less recent years (suppl. Figure 11, Online Resource 1).

Dimension 1: demyelination with an active immune system

Microglia morphology was strongly associated with dimension 1, in the form of a positive correlation with the microglial/macrophage activation score (suppl. Figure 12, Online Resource 1); the associated lesion type (active or mixed) seemed to be of less importance. End-stage (remyelinated and inactive) lesions were relatively infrequent in donors with a high score on dimension 1 (Fig. 3a).

Fig. 3

Dimension 1. a Vertical bar graphs showing on the Y-axis the ramified, ameboid and foamy microglia proportions relative to the total number of active and mixed lesions, white matter lesion proportions relative to the total number of white matter lesions, remyelinated and inactive lesion proportions relative to the total number of end-stage (remyelinated plus inactive) lesions, and active and mixed lesion proportions relative to the total number active plus mixed lesions, from top to bottom respectively, with donors ranked according to their score on the X-axis. The lines represent the centered moving average of the ratio remyelinated to end-stage, and the ratio active to active plus mixed, over a window of 20 observations, with a maximum of 10 missing values. b Box plot (top) showing the ranked scores for 175 donors with and 31 donors without (a) cortical lesion(s). Vertical bar graph (bottom) showing the proportions of cortical lesion types on the Y-axis for donors ranked according to their score on dimension 1. In a and b, ties were assigned ranks in the order of appearance in the dataset. c Box plots showing, from top to bottom, the ranked scores of 83 donors with and 49 without lesions in the MOB tissue block used to determine B cell presence; 6 donors with and 126 without CD20+ cells in pch; 21 donors with and 111 without CD20+ cells in pvs; 69 donors with and 63 without CD20+ cells in men; 9 donors with and 123 without CD138+ cells in pch; 14 donors with and 118 without CD138+ cells in pvs; 62 donors with and 70 without CD138+ cells in men; 10 donors with and 16 without CD20+ cells in subc. active lesions; 11 donors with and 25 without CD20+ cells in subc. mixed lesions; 6 donors with and 21 without CD138+ cells in subc. active lesions; and 4 donors with and 32 without CD138+ cells in subc. mixed lesions. d, e Scatter plots with regression lines, depicting on the Y-axis the rank of the CD3+ cell count in subc. NAWM for 53 donors (d) and in subc. active lesions for 22 donors (e), with donors ranked according to their score on dimension 1 on the X-axis. f Scatter plot with regression line, depicting the rank of the MS PRS (polygenic risk score) on the Y-axis for 194 donors, ranked according to their score on dimension 1 on the X-axis. g Two scatter plots with regression lines showing the ranked score of donors on dimension 1 on the Y-axis, and allele dosage on the X-axis. In c–g, ties were assigned averaged ranks. h Overview of the correlation between relevant input and validation variables and dimension 1, created with BioRender.com. Position on the axis is a close approximation of the Spearman correlation between the dimension and the variable(s); variables were grouped when appropriate and positioned based on the average of the correlation coefficients. Significance in b–f was assessed with Mann–Whitney U for binary variables and Spearman correlation for continuous variables and FDR-adjusted for multiple testing; for g see text. ns p > 0.1; + p ≤ 0.1; * p ≤ 0.05; ** p ≤ 0.01; *** p ≤ 0.001; **** p ≤ 0.0001; ramif. = ramified; ameb. = ameboid; remyel. = remyelinated (lesion); MOB = medulla oblongata; pch = parenchyma; pvs = perivascular space; men = meninges; subc. = subcortical; NAWM = normal-appearing white matter; dim. = dimension; EDSS = Expanded Disability Status Scale; CSF = cerebrospinal fluid; NfL = neurofilament light chain; FDR = False Discovery Rate

Donors with (a) cortical lesion(s) scored higher on dimension 1 (p = 3.8 × 10–6). Moreover, donors who scored high on dimension 1 had relatively fewer intracortical lesions (p = 0.006, rho = −0.24). Although dimension 1 was correspondingly positively correlated with the proportion of leukocortical and subpial lesions, this did not reach statistical significance (leukocortical: p = 0.26, rho = 0.09; subpial: p = 0.22, rho = 0.10) (Fig. 3b & suppl. Figure 13, Online Resource 1).

Donors with CD20+ and/or CD138+ cells in the parenchyma and/or perivascular space of the brainstem scored higher on dimension 1 (CD20+ parenchyma (pch): p = 0.01; CD20+ perivascular space (pvs): p = 0.06; CD138+ pch: p = 0.001; CD138+ pvs: p = 0.04). Furthermore, presence of CD20+ cells in active and mixed subcortical lesions was associated with a higher score on dimension 1 (active: p = 0.04; mixed: p = 0.02). Presence of CD20+ or CD138+ cells in the meninges was not associated with dimension 1 (CD20+: p = 0.70; CD138+: p = 0.70), as was presence of CD138+ cells in active and mixed lesions (active: p = 0.86; mixed: p = 0.89). Donors with a lesion in the brainstem tissue block that was used for investigating B cell presence scored higher on dimension 1 (p = 5.3 × 10–7), as expected (Fig. 3c) [13]. There was no significant association with the mean number of CD3+ cells in NAWM (normal-appearing white matter) of the pyramidal tract at the level of the medulla oblongata with dimension 1 (p = 0.54, rho = −0.09; suppl. Figure 14, Online Resource 1). However, CD3+ cell numbers in subcortical perilesional NAWM were positively correlated with dimension 1 (p = 0.04, rho = 0.39), as were CD3+ cell numbers in active but not mixed lesions (active: p = 0.08, rho = 0.50; mixed: p = 0.54, rho = 0.15) (Fig. 3d, e; suppl. Figure 15, Online Resource 1).

Dimension 1 correlated positively with axonal damage (NfL in CSF: p = 0.05, rho = 0.26); consistently, there seemed to be a lower axonal density and more axonal stress in donors with a high score on this dimension, although these associations were not significant (Bielschowsky: p = 0.48, rho = −0.10; APP+: p = 0.43) (suppl. Figure 16, Online Resource 1).

The PRS (polygenic risk score) and two individual variants (rs3135388, highly correlated with the HLA-DRB1*1501 allele [8], and rs10191329, the SNP that was recently associated with disease progression [23]) were not significantly associated with dimension 1 (PRS: p = 0.14, rho = 0.12; rs3135388: p = 0.17; rs10191329: p = 0.31) (Fig. 3f, g). However, the progression-associated SNP has a relatively low minor allele frequency and increased risk is mediated exclusively by the homozygous carriers. The five NBB donors homozygous for the risk allele (A) of rs10191329 (1.9 < dosage < 2.1) all scored high on dimension 1 (their median rank was 163, with the maximum possible rank being 194).

For a comprehensive overview of dimension 1, see Fig. 3h.

Dimension 2: active lesions, ramified microglia, and an association with the HLA region

Dimension 2 is associated with ramified microglia morphology and a corresponding lower microglial/macrophage activation score (suppl. Figure 12, Online Resource 1), as well as a shift from predominantly mixed to more active lesions. The ratio of remyelinated to inactive lesions remains quite stable (Fig. 4a).

Fig. 4

Dimension 2. Legends for a–f as for Figs. 3a–c, f–h, respectively; the only difference being that ranked scores of donors with and without CD20+ and CD138+ cells in subcortical lesions are not shown in c but in suppl. Figure 17 (Online Resource 1). HLA = human leukocyte antigen

Donors with (a) cortical lesion(s) did not score differently than those without (p = 0.38). There was a clear relation with cortical lesion location: scoring high on dimension 2 was associated with relatively fewer leukocortical lesions (p = 3.1 × 10–4, rho = -0.31), and more intracortical and subpial lesions (intracortical: p = 0.03, rho = 0.18; subpial: p = 0.02, rho = 0.21) (Fig. 4b & suppl. Figure 13, Online Resource 1).

There was no association between CD20+ and CD138+ cells in the brainstem and dimension 2 (for CD20+ cells: pch: p = 0.60; pvs: p = 0.12; meninges (men): p = 0.75; for CD138+ cells: pch: p = 0.82; pvs: p = 0.89; men: p = 0.29), and donors with or without a lesion in the brainstem did not score differently on the dimension (p = 0.19) (Fig. 4c). Presence of CD20+ and CD138+ cells in subcortical active and mixed lesions was also not associated with dimension 2 (for CD20+ cells: active: p = 0.94; mixed: p = 0.38; for CD138+ cells: active: p = 0.11; mixed: p = 0.86; suppl. Figure 17, Online Resource 1). Moreover, neither the mean number of CD3+ cells in NAWM nor in subcortical active and mixed lesions was correlated with dimension 2 (brainstem NAWM: p = 0.54, rho = −0.12; subcortical NAWM: p = 0.20, rho = 0.25; active: p = 0.20, rho = 0.35; mixed: p = 0.69, rho = 0.07; suppl. Figure 14 & 15, Online Resource 1).

Regarding axonal damage, density and stress, there was no correlation with dimension 2 (NfL in CSF: p = 0.73, rho = −0.04; Bielschowsky: p = 0.48, rho = −0.10; APP+: p = 0.99; suppl. Figure 16, Online Resource 1).

The PRS and the risk allele (A) of the HLA-DRB1*1501 tag SNP were associated with dimension 2 (PRS: p = 0.002, rho = 0.24; rs3135388: p = 0.002; Fig. 4d, e). After recalculation of the MS PRS without the HLA region, the PRS was no longer significantly correlated with dimension 2 (p = 0.37, rho = 0.09; suppl. Figure 18, Online Resource 1). Dimension 2 did not correlate with the progression SNP rs10191329 (p = 0.53; Fig. 4e), and did not show a particular clustering of the five homozygous carriers of the risk allele.

For an overview of dimension 2, see Fig. 4f.

Dimension 3: loss of (white matter) lesion activity and scar formation

Dimension 3 is negatively correlated with the microglial/macrophage activation score (suppl. Figure 12, Online Resource 1). In contrast to dimension 2, this was related to a relative increase of mixed compared to active lesions. Moreover, scoring higher on dimension 3 was associated with a decrease in the ratio of remyelinated to inactive lesions (Fig. 5a).

Fig. 5

Dimension 3. Legends for a–c and f–h as for Fig. 3a–c, f–h, respectively. d, e Scatter plots with regression lines, depicting on the Y-axis the rank of CSF NfL levels for 64 donors (d) and the rank of the percentage Bielschowsky+ area for 55 donors (e), with donors ranked according to their score on dimension 3 on the X-axis. sens. = sensory/autonomic, psych. = psychiatric

Donors with one or more lesions in the cortex generally scored higher on dimension 3 (p = 0.02). In addition, higher scoring donors had relatively more subpial lesions (p = 0.03, rho = 0.19). There was a non-significant negative correlation with intracortical and leukocortical lesion proportions (leukocortical: p = 0.11, rho = -0.14; intracortical: p = 0.65, rho = -0.03) (Fig. 5b & suppl. Figure 13, Online Resource 1).

Donors with CD138+ cells in the parenchyma but not the perivascular or meningeal regions of the brainstem scored lower on dimension 3 (pch: p = 0.04; pvs: p = 0.95; men: p = 0.18); donors with CD138+ cells in subcortical active and mixed lesions did not score differently (active: p = 0.89; mixed: p = 0.86). There was no significant association with the presence of CD20+ cells at any of the locations in the brainstem (pch: p = 0.17; pvs: p = 0.29; men: p = 0.70), or with the presence of a lesion in the investigated tissue block (p = 0.60). However, the presence of CD20+ cells in subcortical active and mixed lesions was associated with a lower score on dimension 3 (active: p = 0.08; mixed: p = 0.06) (Fig. 5c). Regarding T cells, there was no correlation between CD3+ cells in the NAWM or subcortical lesions and dimension 3 (brainstem NAWM: p = 0.88, rho = −0.02; subcortical NAWM: p = 0.51, rho = 0.13; active: p = 0.56, rho = -0.15; mixed: p = 0.20, rho = 0.29; suppl. Figure 14 & 15, Online Resource 1).

Dimension 3 was negatively correlated with axonal damage (NfL in CSF: p = 0.02, rho = −0.34) and positively with axonal density (Bielschowsky: p = 0.06, rho = 0.31) (Fig. 5d, e). Although the association with axonal stress was not significant, donors with APP+ axonal fragments or bulbs generally scored lower on dimension 3 (p = 0.42; suppl. Figure 16, Online Resource 1).

There was no significant association with the MS PRS (p = 0.60, rho = -0.04), the HLA-DRB1*1501 tag SNP (p = 0.31) or the progression SNP (p = 0.31) (Fig. 5f, g).

For an overview of dimension 3, see Fig. 5h.