A total of five datasets were downloaded from the GEO database¹, namely, two AD datasets (GSE33000 and GSE44770), a dataset consisting of altered microRNA expressions in circulating blood after exercise (GSE144627), a skeletal muscle load-related dataset (GSE155933), and a single-cell transcriptome dataset of the skeletal muscle (GSE138826) (Supplementary Table 3). The GSE33000 dataset is based on the GPL4372 platform. This was examined on the Rosetta/Merck Human 44k 1.1 microarray. In this study, 467 human brain tissues were selected; among them, 310 were specimens from AD patients and 157 were specimens from non-demented controls. All brain tissue samples were obtained from the Harvard Brain Tissue Resource Center (HBTRC) (Narayanan et al., 2014). The GSE44770 dataset is also based on the GPL4372 platform. This was examined on the Rosetta/Merck Human 44k 1.1 microarray. A total of 230 human brain tissues were included; among them, 129 specimens were obtained from AD patients and 101 specimens were obtained from non-demented controls (Zhang et al., 2013). The GSE144627 dataset is based on the GPL15520 platform. The miRNA samples of circulating exosomes from 10 older adults (five endurance trainers and five sedentary individuals) were analyzed using the Illumina MiSeq (Homo sapiens) system (Nair et al., 2020). The GSE155933 dataset is based on the GPL24047 platform. This was assayed using the Affymetrix Human Transcriptome Array 2.0. A total of 230 skeletal muscle samples were included and the dataset comprised information on muscle conditions before and after muscle resistance training/unloading (Timmons et al., 2019). The GSE138826 dataset, based on the GPL24247 platform, was screened using Illumina NovaSeq 6000. It is a single-cell dataset from the mouse tibialis anterior muscle tissue (Oprescu et al., 2020). Differential expression analysis was performed using the limma package in R software; the fold change (FC) and false discovery rate (FDR) of all differential genes were noted, and a p-value < 0.05 was the criterion for statistical significance (Costa-Silva et al., 2017).

The necroptosis gene set, M24779.gmt, consists of eight necroptosis genes, all of which were downloaded from the Gene Set Enrichment Analysis (GSEA)² database. Finally, 67 necroptosis-related genes were included in the study (Zhao et al., 2021). After Pearson correlation analysis for all lncRNAs and the 67 necroptosis-related genes in the gene expression matrix, a total of 32 lncRNAs were identified and defined as significant necroptosis-related lncRNAs in AD (Pearson correlation coefficients > 0.4 and p < 0.001).

To ensure the reproducibility of the results, we set the seed at “123456.” The random forest analysis was performed using the “randomForest” package in the R software. The parameters were set using the default function with 500 “trees.” Since the error of the random forest model was minimized when 324 trees were used, accordingly, the importance of each gene was defined. The neural network model was constructed using the “neuralnet” and “NeuralNetTools” packages, with the seed set to “12345678” and the parameter to “hidden = 5.” The GSE33000 dataset was used as the training group, and the GSE44770 dataset was used as the test group. The “pROC” package was used to plot the receiver operating characteristic (ROC) curve for the model.

The visualization function in R software was used to draw the Venn diagram. Prediction of mRNA targets for corresponding miRNAs was based on the database results. miRNA/mRNA interactions were identified using the miRWalk database³; subsequently, the miRNA–mRNA regulatory networks were constructed.

The scRNA-seq analysis was performed using the “Seurat” package. Briefly, the cell clustering and clustering for the GSE138826 dataset followed the methods described previously by Oprescu et al. (2020). SCENIC supports the analysis of positive transcriptional regulation, which implies that it can screen transcription factors that regulate miRNAs during positive transcription. An improved version of the SCENIC method was used to screen the key transcription factors in skeletal muscles, as described previously (Suo et al., 2018; Chen Y. et al., 2021; Lin et al., 2021a,b). Moreover, to quantify the cell-type specificity of a regulon, we adapted an entropy-based strategy that was previously used for the analysis of gene expression data (Cabili et al., 2011).

As an extension of the GSEA, a single sample gene set enrichment analysis (ssGSEA) allows for defining an enrichment score corresponding to the absolute enrichment of a gene set for each sample within a given dataset. This was published in 2009 (Barbie et al., 2009). The “GSVA” and “GSEABase” packages in the R software were used for immuno-infiltration analysis; the “limma” package was used for differential analysis, and the results of the final analyses were visualized using the “ggpubr” package. Immune checkpoint-related genes were extracted to screen the differentially expressed immune checkpoint genes in AD.

The R (version 3.6.3) software was used for GSEA and visualization, as described in previous studies (Wang et al., 2020). The reference gene collection “c2.cp.v7.2.symbols.gmt” and the species Homo sapiens were chosen for this study, and the analysis was carried out according to the clusterProfiler package’s instructions. Pathways meeting a FDR < 0.25 and p.adjust < 0.05 were considered significantly enriched.

The STRING database⁴ was used to calculate the level of interaction likelihood and construct the network of interactions between known and predicted proteins. The more complex was the protein-protein relationship network encoded by a gene, the more was its importance for the gene function. The protein–protein interaction (PPI) network encoded by DEGs was constructed using STRING, and the key genes were screened. Drug sensitivity analysis was based on the results of the “pRRophetic” package; the “limma” package was used for variance analysis. The “pRRophetic” package uses the expression matrix and drug handling information from the Cancer Genome Project (CGP), a database comprising 138 anticancer drugs tested against 727 cell lines. The “cgp2016” dataset was used for the bulk prediction of drug data.

Based on the neural network prediction model, we constructed a new logistic regression-based AD risk prediction model in order to provide a visual basis for the formulation of exercise prescription. Factors with p-values < 0.05 were further screened and used to construct a column line graph prediction model, as described previously (Chen et al., 2020; Ying et al., 2021; Zhou et al., 2021). The GSE44770 dataset was used as an external dataset to validate the findings. These data were used for principal component analysis (PCA) to assess their ability to discriminate between different outcomes as our previous research (Ying et al., 2021). Calibration curves were plotted to evaluate the accuracy of the nomogram. Moreover, the C-index and AUC values were also calculated for further verification of the prediction accuracy. Besides, accuracy, F-value, precision, and recall of each dataset and proposed nomogram were also described as our previous research (Wang et al., 2020). Finally, decision curves were used to evaluate the clinical applicability of the nomogram (Kang et al., 2020).

The research flow schema of this study is shown in Figure 1. Sixty-seven necroptosis-related genes were initially included (Zhao et al., 2021). After Pearson correlation analysis for all lncRNAs in the GSE33000 gene expression matrix and 67 necroptosis-related genes, a total of 32 lncRNAs were identified for subsequent analysis. These lncRNAs were considered as significant necroptosis-associated lncRNAs (Pearson correlation coefficients > 0.4 and p < 0.001) in AD. The PPI network for these lncRNAs is positively associated with mRNAs as shown in Supplementary Figure 1. Our main objective was to screen lncRNAs that may regulate exosomal microRNAs. MicroRNAs negatively regulate mRNAs. Therefore, based on the ceRNA theory, we only screened the lncRNAs–mRNA pairs, wherein the lncRNAs were positively associated with mRNAs. Figure 2A presents a circle diagram of necroptosis-related lncRNAs–mRNAs relationship pairs. The lncRNAs are in the middle of the circle; mRNAs are in the outer periphery of the circle, and the linkage of the outer circle origin with the inner circle origin indicates a significant correlation. To determine whether these lncRNAs–mRNAs relationship pairs are differentially expressed in AD patients and normal patients, we performed the difference analysis. Red connecting lines indicate that these mRNAs are highly expressed in AD, and purple indicates that these mRNAs are lowly expressed in AD. Volcano maps were used to visualize the significantly different lncRNAs and mRNAs. Since data on microRNAs in the GSE33000 dataset are lacking, the dataset of microRNAs is not shown in the volcano plot. The volcano plots for the differentially expressed genes between AD and normal brain tissues are shown in Figures 2B,C. GATA3, ALK, TRAF2, HSPA4, HDAC9, BCL2L11, FASLG, and IDH1 are significantly overexpressed in AD patients. Also, the EGFR expression is downregulated in AD patients. DLEU2, DKFZP434H168, MYCNOS, C22or24, C8orf49, FLJ13224, DGOR5, and HPYR1 are lncRNAs that are also overexpressed in AD.

Genes showing significant differences were used for analysis using the random forest model. In Figure 2D, the green line indicates the training group error; the red line indicates the test group error, and the black line indicates the overall group error. We found that the random forest model exhibited the lowest error at 324 “trees,” and all genes with scores less than 2 were excluded (Figure 2E). Due to the large error, it was necessary to further optimize the random forest model. The ROC curves for the AD neural network prediction model (Figure 2F) are shown in Figures 2G,H; the AUC of the training group was 0.948 (95% CI: 0.928–0.967) and that of the test group was 0.932 (95% CI: 0.896–0.961). This result suggested that mRNAs and lncRNAs based on necroptosis-related genes show good predictive ability.

Since this neural network model did not include basic clinical features such as age and gender, we constructed a forest plot using the age, gender, and neural network risk score (Figure 3A). A forest plot for the multifactor logistical analysis based on age, gender, and neural network risk scores is shown in Figure 3A. The p-values for multifactor regression analysis of age and neural network risk score were significant and were subsequently used to construct the nomogram prediction model (Figure 3B). A PCA suggested that the model could better distinguish AD from normal patients (Figure 3C). The calibration curve was plotted using data from the entire group (Figure 3D). The mean absolute error of this calibration curve was 0.016, while the mean squared error was 0.0049 (Figure 3D). To determine the prediction accuracy of this neural network model, we calculated the AUC for the ROC curve (Figure 3E). The AUC value of the training group was 0.979, while that of the test group was 0.957. In addition, the C-index was calculated and reached 0.979 [95% credibility interval (95% CI; 0.959–0.983)] in the training set and 0.957 (95% CI; 0.968–0.990) in the validation set (Table 1). Besides, we also described the accuracy, F-value, precision, and recall of each dataset and proposed nomogram in Table 2. This indicated that the prediction model was capable of accurately representing the transcriptomic aspects of the brains of AD patients and had high predictive potential. The DCA curve for this nomogram prediction model is shown in Figure 3F, which suggests that the model has a very high predictive power.

Previous studies show that microenvironmental alterations characterized by an imbalance in the immune microenvironment of the brain induced upon damage to the blood–brain barrier are closely related to the AD pathogenesis (Lu et al., 2021). In our investigation, using ssGSEA, we determined the degree of immune cell infiltration. Cytolytic activity, MHC class I, neutrophils, para-inflammation, T cell co-inhibition, Tfh, TIL, Type I, and Type II IFN responses were all found to be significantly associated with the AD risk assessment model (Figure 4A). ICOSLG, PDCD1, and TNFRSF25 were found to be overexpressed in patients in the high-risk group. In the high-risk group, the expressions of TNFRSF18, HAVCR2, and CD276, were found to be downregulated (Figure 4B). Findings from GSEA revealed a substantial association between the gut immune network for IgA production and the AD risk score (Figure 4C). The aforementioned results indicated that this model of AD risk assessment was intimately connected to the local immunological microenvironment. Imbalances in the circulatory system and intracerebral immune microenvironment were found to be strongly related to neuronal necroptosis. Additionally, Supplementary Figure 2 illustrates the differential susceptibility of neuronal cells to medicines in individuals at high and low risk of developing AD. Thus, those medicines may contribute to AD by impairing neuronal necroptosis.

Previous studies suggest that regular exercise leads to the upregulation of three microRNAs, namely, miR-486-5p, miR-215-5p, and miR-941, in the peripheral blood, while that of the exosomal miR-151b is downregulated (Nair et al., 2020). We performed a differential analysis using the GSE144627 dataset and found that miR-1306-3p, miR-215-5p, miR-432-5p, miR-129-5p, miR-370-3p, and miR-197-3p were all upregulated in the exercise group (Supplementary Table 1). The logFC values for miR-215-5p and miR-197-3p expressions were greater than 2. The Venn diagram demonstrated that the microRNA, miR-215-5p, was co-significantly differentially expressed in both previous studies and the GSE144627 dataset (Figure 5A). The miRWalk database was used to construct a miRNA–mRNA regulatory network (Supplementary Table 2). We reasonably hypothesized that the upregulation of exosomal hsa-miR-215 expression in circulating blood after exercise may contribute to the inhibition of AD necroptosis-related genes (e.g., IDH1, SIRT1, and BCL2L11). Among them, SIRT1 and BCL2L11 have been previously validated as the target genes of miR-215-5p.

Based on the TransmiR v2.0 database, a total of 44 TFs, the possible upstream transcriptional regulators regulating hsa-miR-215 transcription, were predicted (Figure 5B). Cell clustering in the GSE138826 dataset was performed based on the findings of Oprescu et al. (2020; Figure 5C). Subsequently, the connection specificity index (CSl) matrix was constructed using the SCENIC method, and a total of ten modules were obtained (Figure 5D). In these modules, 316 key TFs were included. The Venn diagram shows the intersection of key TFs from skeletal muscle scRNA-seq data with those predicted to potentially regulate miR-215 using the TransmiR database (Figure 5E). The regulator rankings for FAPs and skeletal muscle based on the Regulon specificity score (RSS) are shown in Figure 5F. Thus, Cebpe and GATA6 were identified as key TF in FAPs and skeletal muscle cells, respectively. Since Cebpe and GATA6 could potentially regulate miR-215, we speculated that the upregulation of Cebpe and GATA6 expressions in skeletal muscles could elevate miR-215 levels in exosomes.

The tSNE plot was used to demonstrate the expression characteristics of the regulon, CEBPB, and GATA6 from the scRNA-seq data (Figures 6A,B). GATA6 and CEBPB were found to be significantly upregulated in the de-load group (Figure 6C). The tSNE plot shows the average expression characteristics of each gene module (Figure 6D). The distribution of the rankings of each cell type in each module according to the regulon activity score is also shown (Figure 6E). This suggested that CEBPB may exert an important transcriptional regulatory effect in the M2 and M3 cell clusters, while GATA6 may play a critical transcriptional regulatory role in the M4 and M5 cell clusters. CEBPB and GATA6 were found to be significantly upregulated after muscle exercise and de-loading. In a previous study, we found that regular exercise resulted in a significant increase in baseline expression of exosomal microRNA-215. However, acute exercise did not elicit the same effect on exosomal miRNA-215 expression (Figure 7). These suggest that the upregulation of miR-215 expression is not transient, rather it is the result of longterm exercise treatment. Although further validation is needed, we hypothesize that the transcriptional activity of CEBPB and GATA6 in upregulating miR-215 may be a potential negative regulatory feedback mechanism in exosome homeostasis after skeletal muscle unloading.