The study included individuals with sporadic ALS who were autopsied between 1978 and 2023 at the Department of Pathology, Brain Research Institute, Niigata University, Japan, and were diagnosed neuropathologically as having ALS. Eligibility criteria for cases were the presence of TDP-43-positive inclusions in the central nervous system and provision of consent for DNA analysis. Cases in which the onset of dementia preceded the onset of motor symptoms by more than one year were excluded. Among the cases included in the present study, 75 had already been analyzed and included in a previous report [11].

Exome sequencing was performed to identify variants in known ALS causative genes, and APOE genotype was also determined from the exome sequencing data. DNA samples were extracted from autopsied central nervous system tissues (occipital cortex, motor cortex, and cerebellum) using a DNA extraction kit (QIAamp DNA Mini Kit; Bio-Rad Laboratories, Hercules, CA, USA). Exome sequencing was outsourced to Takara Bio (Shiga, Japan) or Macrogen (Seoul, South Korea) and performed using the Illumina NovaSeq 6000 platform. Analysis of GGGGCC repeat expansions in intron 1 of C9orf72 was performed in all cases using previously reported methods [12].

Annotated exome sequencing data were used to detect variants in known ALS-associated genes, including TARDBP, OPTN, FUS, SOD1, TBK1, SQSTM1, MATR3, TUBA4A, NEK1, HNRNPA2B1, VCP, ELP3, SETX, HNRNPA1, CCNF, VAPB, C21orf2, CHCHD10, NEFH, ANG, DCTN1, CHMP2B, UBQLN2, FIG4, PFN1, EWSR1, TAF15, ANXA11, DAO, ERBB4, TIA1, GLE1, PRPH, ALS2, SPG11, SIGMAR1, KIF5A, and DNAJC7. Given that the allele frequency of the well-established pathogenic variant SOD1 p.D91A is 0.001 across all populations in the gnomAD database [13], non-synonymous variants with an allele frequency of < 0.001 were defined as rare variants (RV) and included in the analysis. All detected rare variants were validated by Sanger sequencing regardless of sequencing depth. Allele frequency data were obtained from the Human Genetic Variation Database (HGVD; http://www.hgvd.genome.med.kyoto-u.ac.jp/) and the Exome Aggregation Consortium databases (ExAC All / East Asian populations; https://gnomad.broadinstitute.org/). For autosomal recessive ALS-associated genes (ALS2, SPG11, and SIGMAR1), only homozygous variants were considered.

APOE genotypes were determined using annotated exome sequencing data, with alleles defined by amino acid residues at positions 112 and 158: ε2 (Cys112/Cys158), ε3 (Cys112/Arg158), and ε4 (Arg112/Arg158) [14]. APOE genotypes were assigned based on the combination of these alleles.

Immunohistochemical analyses were performed on 4-μm-thick sections prepared from formalin-fixed, paraffin-embedded blocks obtained from multiple regions of the central nervous system, including the frontal cortex, motor cortex, temporal cortex, hippocampus medulla oblongata, and spinal cord. ALS-associated TDP-43 pathology was evaluated using immunohistochemistry with an antibody against phosphorylated TDP-43 (pTDP-43; pS409/410). We previously demonstrated that the majority of patients with sporadic ALS can be classified into two pathological types (type 1 and type 2, as mentioned above) based on the presence or absence of TDP-43 neuronal cytoplasmic inclusions (NCIs) in the hippocampal dentate gyrus [4]. Accordingly, all cases in the present study were classified using this criterion. Representative pTDP-43 immunohistochemistry images in type 1 and type 2 patients are presented in Fig. 1. In patients with both types, pTDP-43 NCIs and GCIs were observed in the motor cortex (Fig. 1a, d). By contrast, in type 1 patients, they were absent from the hippocampal dentate granule cells (Fig. 1b) and exceedingly rare in the temporal cortex (Fig. 1c). In type 2 patients, however, NCIs and GCIs were present in varying numbers in these regions (Fig. 1e, f). AD-related pathological changes were staged according to the criteria proposed by Thal et al. for senile plaques and by Braak et al. for neurofibrillary tangles (NFT) [15,16,17]. Immunohistochemistry for amyloid-β and phosphorylated tau was performed on temporal lobe sections including the hippocampus, with additional staining of the frontal cortex, motor cortex, and occipital cortex as needed. The primary antibodies used in this study are listed in Supplementary Table 1.

Representative TDP-43 pathology of ALS type 1 and type 2. a—c. Phosphorylated TDP-43 (pTDP-43) immunohistochemistry in a patient with ALS type 1. The motor cortex shows a few neuronal cytoplasmic inclusions (NCIs) and scattered glial cytoplasmic inclusions (a). No pTDP-43-immunoreactive inclusions are evident in the hippocampal granule cells (b) or temporal cortex (c). Arrow, NCI. d—f. pTDP-43 immunohistochemistry in a patient with ALS type 2. Abundant NCIs are evident in the motor cortex (d), hippocampal granule cells (e), and temporal cortex (f). Bar = 50 μm

The primary analytical approaches used in this study were Bayesian SEM and random forest analysis. SEM was employed to simultaneously evaluate the complex relationships among genetic factors, AD-related pathologies, and temporal clinical parameters, and to estimate their direct and indirect effects on the distribution of TDP-43 pathology. In the present study, several endogenous variables were ordinal in nature, and the sample size was relatively small; therefore, Bayesian SEM was adopted, as it is known to provide more stable parameter estimates than conventional maximum likelihood estimation under these conditions. In addition, random forest analysis was used as a complementary approach. This nonparametric machine learning method does not require a priori assumptions regarding causal pathways and is capable of capturing nonlinear relationships and interactions among variables.

Determinants of TDP-43 pathological subtype (type 1 vs. type 2) were evaluated using both SEM and random forest models. The overall analytical workflow is summarized as follows. Two analytical frameworks were constructed by varying the definitions of temporal clinical parameters: model (a) and model (b). In both random forest models, the explanatory variables included Thal Aβ phase, Braak NFT stage, APOE4 status, and the presence of RV. Temporal clinical parameters were defined as follows: in model (a), age at death and entire disease course were included, whereas in model (b), age at onset and survival time were used. “Entire course” was defined as the interval from symptom onset to death, including the period after initiation of mechanical ventilation. “Survival time” was defined as the interval from symptom onset to death or initiation of mechanical ventilation. The hypothesized pathways among variables in the SEM analyses for models (a) and (b) are illustrated in Fig. 3a and b, respectively. In model (a), a final SEM model was constructed by retaining only the paths whose 95% credible intervals did not include zero. In addition, determinants of age at death were evaluated using both SEM and random forest models. In this random forest analysis, the explanatory variables included APOE4 status, Thal Aβ phase, Braak NFT stage, TDP-43 pathology type (type 2), and RV. The corresponding SEM pathways are shown in Fig. 3.

Categorical variables are summarized as counts and percentages, and continuous variables as medians and ranges. Comparisons between the TDP-43 pathology type 1 and type 2 groups were performed using Fisher’s exact test for categorical variables, the Wilcoxon rank-sum test for continuous variables, and the log-rank test for survival time.

For SEM, variables were defined as follows. APOE genotypes ε2/ε4 and ε3/ε4 were coded as APOE4 = 1, ε4/ε4 as APOE4 = 2, and all other genotypes as APOE4 = 0. The presence of at least one RV in known ALS-associated genes was coded as RV = 1, and the absence of such variants as RV = 0. The Thal Aβ phase, Braak NFT stage, and TDP-43 pathological types were classified according to their respective classification systems.

Bayesian SEM was used to accommodate ordinal endogenous variables and a relatively small sample size [18]. Bayesian estimation was performed using the blavaan package in R. All endogenous variables, except for age at death, were treated as ordinals in the blavaan analysis. Standardized path coefficients were estimated, and 95% credible intervals (posterior intervals, indicating a 95% probability that the true value lies within the interval) were calculated. Effects whose 95% credible intervals did not include zero were considered statistically significant. Model fit was evaluated using the posterior predictive p-value (PPP), and PPP values between 0.1 and 0.9 were considered to indicate an acceptable model fit.

To assess the robustness of the results, SEM (a) model was also estimated using the unweighted least squares parameter with mean- and variance-adjusted chi-squared test statistic (ULSMV). Standardized coefficients and p values were calculated. ULSMV estimation was performed using the lavaan package in R. For ULSMV analyses, missing data were imputed using the mice package, with the number of imputations set to m = 20. For missing data imputation using the mice package, the predictive mean matching (pmm) method was applied to continuous variables (age at death and entire course), whereas the proportional odds logistic regression (polr) method was used for ordinal variables.

Random forest models were evaluated using five-fold cross-validation. To prevent data leakage, missing data imputation was performed after splitting the data into training and validation folds. Specifically, the imputation model was fitted using only the training data with the mice package and then applied to both the training and validation folds. The random forest model was trained on the training fold and evaluated on the validation fold. Model performance was assessed using root mean squared error (RMSE), mean absolute error (MAE), and R2 for the analysis of determinants of age at death, and area under the receiver operating characteristic curve (AUC) for the classification of TDP-43 pathological type. The average performance across the five folds was used as the final evaluation metric. For missing data imputation, pmm was applied to continuous variables (age at death, entire course, age at onset, and survival time), whereas polr was used for ordinal variables. The number of imputations was set to m = 20. To visualize the effects of individual variables, partial dependence plots (PDPs) were generated by averaging the effects of other variables. PDPs were constructed using the combined training and validation data, with predictions averaged across the 20 imputed datasets and further averaged over the five folds. Variable importance was assessed using the percent increase in mean squared error (%IncMSE) for the regression task (age at death) and Mean Decrease in Accuracy for the classification task (TDP-43 pathology type).

All statistical analyses were conducted using R software version 4.4.0.