The data utilized for our study were acquired from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (https://adni.loni.usc.edu/). The ADNI was established in 2003 as a collaboration between public and private entities, with approval from the institutional review boards at all ADNI sites. The complete list can be found at http://adni.loni.usc.edu. The primary aim of the ADNI was to test whether it is possible to combine different indicators for predicting the progression of MCI to AD. All participants provided explicit consent before taking part in the study. All procedures were carried out in accordance with the applicable rules and regulations. The Publication Committee of the ADNI approved this study.

According to the ADNI protocol, in the study, MCI was identified when a patient or caregiver reported cognitive decline, if the participant showed indications of impairment in the logical memory II subtest of the Wechsler Memory Scale, if the participant achieved a Mini-Mental State Examination (MMSE) score of 24 or higher, and if the participant had a clinical dementia rating of 0.5. Participants with MCI who met the dementia diagnostic criteria were excluded. AD was diagnosed using the National Institute of Neurological and Communicative Disorders and Stroke (NINCDS) and the Alzheimer's Disease and Related Disorders Association (ADRDA) criteria for probable AD (see http://adni.loni.usc.edu/methods/documents for detailed inclusion and exclusion criteria).

We downloaded 744 participants' peripheral blood RNA microarrays (normalized by the robust multi-chip average method) from the ADNI database. The source of the microarray was GPL13667 (Affymetrix GPL platform), Human Genome U219 Array (see https://ida.loni.usc.edu/pages/access/studyData.jsp for details on RNA microarrays and data preprocessing). Among the744 participants, 383 were patients with MCI.

We set a follow-up time of 72 months, with AD conversion as the endpoint event. We defined the occurrence of the endpoint events as the follow-up endpoint and had the follow-up deadline as the endpoint for those who did not experience the endpoint event. A total of 82 patients with MCI had the endpoint events during the follow-up period and were included in the progressive MCI (pMCI) group, and 59 patients with MCI had no endpoint events within 6 years and were included in the stable MCI (sMCI) group. These patients were included in two cohort studies, ADNI 2 and ADNI GO, with the follow-up period ranging from 6 to 72 months.

In our study, the ADNI 2 dataset (51 cases of pMCI and 29 cases of sMCI) was used for the identification of feature genes, the development of a nomogram prognostic model, and internal validation, whereas the ADNI GO dataset (31 cases of pMCI and 30 cases of sMCI) was used for external validation. According to the study protocol, the remaining patients with MCI were excluded for the following reasons: 11 patients with MCI had outlier RNA microarray samples at the time of weighted gene co-expression network analysis (WGCNA), 29 patients with MCI were converted to cognitively normal cohort (CN) at the follow-up endpoint, 169 patients with MCI were followed-up for < 6 years and no endpoint events occurred, and 33 patients with MCI had no follow-up information. Our study flowchart is presented in Figure 1.

WGCNA can be utilized to detect sets of genes with comparable expression patterns (Langfelder and Horvath, 2008). We constructed a WGCNA network for the progression from MCI to AD with all genes from the ADNI 2 dataset using the “WGCNA” R package.

The “limma” R package was used to identify differentially expressed genes (DEGs), with the following screening criteria: a |log2fold change (FC)| of >0.263 and a p-value of < 0.05, where a log2FC of >0.263 and a p-value of < 0.05 was considered Up and a log2FC of < -0.263 and a p-value of < 0.05 was considered Down. The “pheatmap” and “ggplot2” R packages were used to plot the volcano and heatmap of the DEGs.

Subsequently, the obtained DEGs were intersected with the genes in the midnightblue module using WGCNA to obtain DEGs related to the prognosis of MCI.

The least absolute selection and shrinkage operator (LASSO) is a data mining method that achieves an equilibrium between the model variance (the variance of regression coefficients) and the bias (the difference between the predicted value and the true value) by adjusting the parameter lambda (λ). The value of λ with the smallest error was selected as the optimal value using the 10-fold cross-validation method, and the variables included in the model corresponding to this value of λ were significant variables. The “glmnet” R packages were used for the LASSO.

The support vector machine-recursive feature elimination (SVM-RFE) constructs variable ranking coefficients based on the weight vector ω generated by an SVM during training, retains the variables with significant effects, and finally obtains the decreasing ranking of all variable attributes. In the 10-fold cross-validation method, the variables with the minimum root mean square error and the maximum accuracy were considered significant variables. The “e1072” R package was used for the SVM-RFE.

Subsequently, the significant variables obtained by the LASSO were intersected with the significant variables obtained by the SVM-RFE to obtain hub variables related to the prognosis of MCI.

A receiver operating characteristic (ROC) curve was used to evaluate the hub variables and find the optimal cut-off value. The Kaplan–Meier analysis was conducted to compare the predicted individual risk and observe the non-progression proportion. A univariate Cox proportional-hazards regression analysis was conducted to preliminarily assess the impact of the hub variables related to the prognosis of MCI. The multivariable Cox proportional-hazards model was used to screen the independent risk factors. Hazard ratio (HR), 95% confidence interval (CI), and p-values were taken into account.

Based on the results of the above analysis, we found that a low MoCA score and low EBF1 expression are independent risk factors for predicting the progression of MCI to AD; meanwhile, age and gender may influence the reliability of the findings. Consequently, these independent risk factors and potential confounders (age and gender) were applied to construct the nomogram model (Li et al., 2022). The empirical samples of the ADNI 2 dataset, which were obtained by bootstrap resampling, were used as an internal validation set. External validation of nomograms is required to ensure accuracy outside the original patient data (Cote and Grassbaugh, 2024). The ADNI GO dataset was used as an external validation set. The performance of the nomogram was evaluated using the C-index and calibration curve in the training set, the internal validation set, and the external validation set. Decision curve analysis (DCA) was conducted to evaluate the clinical value of the nomogram (Zhang et al., 2018). The “autoReg,” “rms,” “bootstrap,” “pROC,” “survival,” “ggDCA,” and “rmda” R packages were used in our study.

The gene ontology (GO) analysis mainly includes three parts: biological process, molecular function, and cellular component (Ashburner et al., 2000). The Kyoto Encyclopedia of Genes Genomes (KEGG) analysis provides information on gene-related signaling pathways (Kanehisa, 2002). In our study, the “clusterProfiler” R package was used for GO and KEGG functional enrichment analyses. The p-value was set at < 0.05, and the false discovery rate was set at < 0.25.

We used the Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) algorithm to analyze the correlation between the EBF1 (early B cell factor 1) expression and the immune cell infiltration levels. The “CIBERSORT” R package was used in our study. The differences in immune cell abundances between the high EBF1 expression and low EBF1 expression groups were estimated using the Wilcoxon rank-sum test. The correlation between the EBF1 expression and the differential immune cells was analyzed using the Spearman correlation analysis.

A total of 51 pMCI and 29 sMCI cases were included in the training set (ADNI 2 dataset) and 31 pMCI and 30 sMCI cases were included in the external validation set (ADNI GO dataset). The demographic characteristics, results of the MoCA, and 18F-AV-45 PET of the two sets are shown in Table 1. In the two sets, there was a statistically significant difference between the patients with pMCI and those with sMCI in the MoCA score (p < 0.001); therefore, the MoCA score as a variable was included in the subsequent analysis.

After data preprocessing, we obtained the gene expression matrix of 80 samples of the ADNI 2 dataset (14,893 genes). A sample clustering tree, as shown in Figure 2A, was obtained through the clustering of the samples.

Choosing the ideal soft threshold β can enhance the robust connection and diminish the feeble connection among genes, resulting in a constructed network that closely resembles a scale-free network and is more similar to the gene regulatory network in actual biology. Therefore, we selected the soft threshold of β = 7 (using the scale-free topology criterion with R² = 0.85) to achieve a scale-free network (Figures 2B, C).

According to the soft threshold β of 7, the histogram of the distribution of gene connectivity was plotted (Figure 2D). The scale-free network distribution was examined, which showed the number of nodes (k) corresponding to the gene connectivity and was negatively correlated with the probability of node occurrence (p (k)) (correlation coefficient 0.82, slope −1.8), suggesting that the network constructed by the selected soft threshold tended to converge to the scale-free network (Figure 2E).

Subsequently, we created the adjacency matrix and formed a TOM and dissTOM. The correlation heatmap between the genes in the ADNI2 dataset was plotted according to the dissTOM (Figure 2F), and hierarchical clustering was performed to merge the modules with higher similarity (Figure 2G). In the end, a total of 14 modules were discovered and the genes within each module exhibited higher similarity. The smallest module was the midnightblue module, which contained 210 genes. The largest module was the turquoise module, which contained 3,642 genes. The gray module contained all genes that could not be clustered into other modules.

The midnightblue module had a strong association with the prognosis of MCI; therefore, it was chosen as a module of clinical significance for further analysis (Figures 3A, B). Notably, a significant correlation was observed between the MM and GS of the midnightblue module (Figure 3C). By applying MM > 0.8 and GS > 0.3, we identified nine key genes (FCRLA, P2RX5, CD72, MS4A1, EBF1, BLNK, CR2, ABCB4, and FCRL2) in the midnightblue module (Figure 3D).

A total of 178 DEGs were acquired from the ADNI2 dataset, which included 114 downregulated and 37 upregulated genes. The volcano plot and heatmap of the DEGs are shown in Figures 4A, B. A total of 40 DEGs related to the prognosis of MCI were identified by taking the intersection of the DEGs and genes in the midnightblue module using WGCNA (Figure 4C).

The LASSO was performed with the MoCA score and 40 DEGs related to the prognosis of MCI, and eight significant variables were extracted, including the MoCA score, ARHGAP32, P2RX5, EBF1, FCRL2, FCRL5, CELSR1, and SOBP (Figures 4D, E).

In the optimal parameters (minimum root mean square error = 0.2875 and maximum accuracy = 0.7125), the SVM-RFE identified seven significant variables, including the MoCA score, ARHGAP32, P2RX5, EBF1, FCRLA, CELSR1, and CR2 (Figures 4F–H).

The results of the LASSO and SVM-RFE were overlapped by the Venn diagram, and finally, five hub variables were obtained, namely the MoCA score, ARHGAP32, P2RX5, EBF1, and CELSR1 (Figure 4I).

To determine the ability of the MoCA score, ARHGAP32, P2RX5, EBF1, and CELSR1 to distinguish between the sMCI and pMCI and to find the optimal cut-off value, we plotted an ROC curve. The area under the ROC curve (AUC) of the MoCA score, ARHGAP32, P2RX5, EBF1, and CELSR1 were 0.791 (95% CI 0.692–0.891), 0.684 (95% CI 0.551–0.818), 0.688 (95% CI 0.570–0.807), 0.699 (95% CI 0.581–0.817), and 0.691 (95% CI 0.561–0.821) (Figure 5A). The optimal cut-off sensitivity, specificity, Youden index, positive predictive value, and negative predictive values are listed in Table 2.

To assess whether the MoCA score, ARHGAP32, P2RX5, EBF1, and CELSR1 can predict AD dementia in patients with MCI, based on their optimal cut-off value, we divided the patients into the high MoCA score (≥22.5) and low MoCA score (< 22.5) groups, the high ARHGAP32 expression (≥4.565) and low ARHGAP32 expression (< 4.565) groups, the high P2RX5 expression (≥8.577) and low P2RX5 expression (< 8.577) groups, the high EBF1 expression (≥5.845) and low EBF1 expression (< 5.845) groups, and the high CELSR1 expression (≥4.631) and low CELSR1 expression (< 4.631) groups.

The Kaplan–Meier survival curves of the MoCA score, ARHGAP32, P2RX5, EBF1, and CELSR1 showed that the non-progression proportion of the patients in the low groups in the follow-up period was much lower than that of the patients in the high groups (p < 0.001, p = 0.001, p = 0.004, p < 0.001, and p < 0.001, Figures 5B–F).

The results of the univariate Cox proportional-hazards regression analysis showed that a low MoCA score, low ARHGAP32 expression, low P2RX5 expression, low EBF1 expression, and low CELSR1 expression increased the risk of AD dementia in patients with MCI. After the univariable analysis, the five hub variables were entered into the multivariable Cox proportional-hazards model, using stepwise backward selection, and the results demonstrated that the prognosis of MCI was significantly correlated with a low MoCA score (p < 0.001) and low EBF1 expression (p < 0.001) (Table 3). To reduce the impact of potential confounders, age and gender were also included in the final multivariable Cox proportional-hazards model (Figure 5G).

Using the final multivariable model, a nomogram was constructed to predict the progression from MCI to AD (Figure 6A). The nomogram considered two independent risk factors (MoCA score and EBF1 expression), along with two potential confounders (age and gender). A score was assigned to each of the four variables on the point scale axis, and the scores of the variables were used to calculate a cumulative score. By projecting the total score to the total point scale, we were able to estimate the probability of patients with MCI who will progress to AD.

A variety of metrics, including the C-index, calibration, and DCA, were used to evaluate the performance of the nomogram. The C-index of the nomogram was 0.736 in the training set. The empirical samples of the ADNI 2 dataset, which were obtained by bootstrap resampling, were used as an internal validation set, and the C-index was 0.824 in the internal validation set. In the external validation set, the C-index was 0.751. The calibration measured how well the probabilities predicted by our nomogram model compared with the reality. Figures 6B–D present the calibration curves of the nomogram predicting a progression probability of 6 years, which demonstrated a strong agreement between the predicted probabilities by the nomogram and the actual probabilities in the training set, the internal validation set, and the external validation set. Based on the results of the DCA curves, we should avoid using our nomogram when the risk threshold is less than 13% and greater than 96% in the training set and when the risk threshold is less than 11% in the internal and external validation sets. Thus, except for a small range of risk thresholds, the net benefit of predicting progression from MCI to Alzheimer's disease in 6 years using our nomogram is greater than assuming that all patients with MCI will progress to AD or that none will progress to AD (Figure 6E). Therefore, it would be appropriate to use this nomogram in our study to predict progression from MCI to AD in 6 years. The results showed a broad spectrum of alternative threshold probabilities in the nomogram prediction model, demonstrating that our nomogram model worked well as a prediction tool.

The results above showed that, in addition to the MoCA score, EBF1 is valuable in predicting progression from MCI to AD. In the ADNI 2 dataset and ADNI GO dataset, the difference in EBF1 was statistically significant between patients with pMCI and patients with sMCI (p = 0.002 and p = 0.006), and the EBF1 expression level was lower in patients with pMCI than in patients with sMCI (Figures 7A, B).

We divided the patients with pMCI of the ADNI2 dataset into the high EBF1 expression and low EBF1 expression groups. Based on the screening criteria, with a |log2fold change (FC)| of >0.263, a p-value of < 0.05, and a false discovery rate of < 0.05, a total of 276 DEGs related to EBF1 were acquired from the ADNI2 dataset, which included 217 downregulated and 59 upregulated genes (Figures 7C, D).

The KEGG analysis showed that the genes were concentrated in the B cell receptor signaling pathway, primary immunodeficiency, and cellular senescence signaling pathway (Figure 7E).

GO analysis includes three parts: biological process, cellular component, and molecular function. The biological process was enriched in the B cell activation, B cell receptor signaling pathway, B cell proliferation, lymphocyte activation, lymphocyte proliferation, the regulation of the nucleobase-containing compound metabolic process, and intracellular signal transduction (Figure 7F). The cellular component was enriched in the nuclear part, cell surface, biogenesis of lysosome-related organelles (BLOC) complex, chromatin, nuclear chromatin, side of the membrane, cytosol, immunoglobulin complex, chromosomal part, and nuclear chromosome part (Figure 7G). The molecular function was concentrated in the protein homodimerization activity, HMG box domain binding, protein dimerization activity, interleukin-6 receptor binding, DNA binding, lipid binding, and enhancer binding (Figure 7H).

According to the findings of the analysis of the immune cell infiltration, the level of the expression of EBF1 was found to have a positive association with the expression patterns of B cells naïve (r = 0.583, p < 0.001) and mast cells resting (r = 0.234, p = 0.037) and a negative association with the expression patterns of NK cells resting (r = −0.345, p = 0.002) (Figure 8).