The GEO database (https://www.ncbi.nlm.nih.gov/geo/) provided all of the gene expression datasets used in this study. Four datasets, numbered GSE39445, GSE140829, GSE56931, and GSE63061, were selected for this study. These datasets were profiled using microarrays of peripheral blood samples. Among them, GSE39445 and GSE56931 showed insufficient sleep status in humans. GSE39445 consists of two groups: individuals underwent SD for 39–41 h as the sleep restriction group, and individuals underwent 10 h of sleep opportunities each day for 1 week as the sleep extension group (Moller-Levet et al., 2013). In GSE56931, participants underwent a 24-h normal sleep/wake cycle (baseline) and 38 h of continuous wakefulness (sleep deprivation), followed by subsequent recovery sleep (recovery) (Arnardottir et al., 2014). In this study, the sleep restriction group in GSE39445 and sleep deprivation group in GSE56931 were considered the SD groups, while the sleep extension group in GSE39445, and baseline and recovery groups in GSE56931 were considered the control groups in our analysis. In addition, GSE140829 and GSE63061 were downloaded from the NCBI GEO database, both of which included gene expression data for patients with MCI and healthy controls. The GSE140829 dataset contained 134 MCI samples and 249 normal samples (Nachun et al., 2019), and GSE63061 contained 110 patients with MCI and 134 control individuals (Sood et al., 2015). In this work, hub genes associated with SD and MCI were found using GSE39445 and GSE140829, and validated using GSE56931 and GSE63061, respectively. For background correction and data normalization, we used the R package “limma” to handle the raw data from these datasets.

The coding genes for all proteins found in the mitochondrial membrane, matrix, cristae, and mitochondria-associated endoplasmic reticulum membranes were referred to in this study as mitochondria-related genes (MRGs). A total of 1,136 MRGs was acquired from the human MitoCarta3.0 database (https://www.broadinstitute.org/mitocarta/mitocarta30-inventory-mammalian-mitochondrialproteins-and-pathways) (Rath et al., 2021). Next, using the R package “limma” with p < 0.05, we carried out a t-test difference analysis to find DEGs between SD or MCI and control samples. A heatmap is used to display the resulting differential gene expression data. The intersection of the genes among the MRGs, SD DEGs, and MCI DEGs groups was then shown using the R package “VennDiagram”.

The biological roles of DEGs were further investigated by using Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analyses. GO analysis is a popular bioinformatics technique for gene annotation of cellular component (CC), molecular function (MF), and biological process (BP), while KEGG is extensively utilized to comprehend biological mechanisms and functions. The GOplot program package was used to visualize the GO and KEGG analyses, with a cutoff p < 0.05. Enrichment analysis of target genes in transcription factor targets was performed using the Metascape database (https://metascape.org/gp/index.html#/main/step1).

The key candidate genes for SD and MCI were screened using two machine learning methods: random forest (RF) and least absolute shrinkage and selection operator (LASSO). LASSO is a regularization and variable selection method for statistical models that constructs a penalty function to solve complex collinear data and develop a relatively refined model (Tibshirani, 2011). RF, on the other hand, is developed by Breiman, which uses decision trees to train classification samples and generates predictions based on the classification outcomes (Breiman, 2001). The R packages “glmnet” (Friedman et al., 2010) and “randomForest” (Breiman, 2018) were used to conduct the LASSO regression and RF analysis.

The “limma” R package was utilized to compare the differential expression of the characteristic gene between SD or MCI and control samples in training and validation sets. A boxplot was used to show differences in gene expression. In addition, the diagnostic impact of the characteristic gene in training and validation sets was assessed by computing the area under the curve (AUC) of a receiver-operating characteristic (ROC) curve using the “pROC” R package (Robin et al., 2011).

The functional categories and pathways enriched for the common biomarker were analyzed using GO and KEGG enrichment analysis based on gene set variation analysis (GSVA) and gene set enrichment analysis (GSEA), in order to further investigate the possible role of the common biomarker in SD and MCI. GSEA and GSVA were carried out using the R packages “clusterprofiler” (Yu et al., 2012) and “GSVA” (Hanzelmann et al., 2013). The reference gene set was determined to be h.all.v6.2.symbols.gmt from the Molecular Signatures Database, using a cut-off criterion of < 0.05 as a significant threshold.

The percentage of immune cells in each group was ascertained using the single-sample GSEA (ssGSEA) method. Based on the expression levels of immune cell-related genes generated from the expression profiles, the ssGSEA algorithm was used to examine the infiltration abundance of 28 immune cell types and the correlation between common biomarkers and immune-infiltrating cells. In addition, the packages “limma,” “reshape2,” “tidyverse,” and “ggplot2” were used to evaluate and display the outcomes.

Depending on the features of the data distribution, non-parametric tests or t-tests were used to assess the statistical significance of the differences between the two groups. R 4.2.3 was used for all analyses. A threshold of p < 0.05 was established for statistical significance.

The flow diagram in Figure 1 served as the guide for conducting this study. We used the expression profiling datasets GSE39445 and GSE140829 to screen for DEGs in the SD and MCI groups, respectively. A total of 1926 DEGs was obtained from the SD dataset, including 561 upregulated DEGs and 1,365 downregulated genes. In addition, 3,815 DEGs were identified in the MCI dataset using differential expression analysis, of which 2,151 genes had upregulated expression, and 1,664 genes had downregulated expression. As shown in Figures 2A, B, we employed heatmaps to present the results of DEGs among different samples in the SD and MCI datasets.

To examine the interrelationship between mitochondrial dysfunction and cognitive impairment induced by SD, we employed the MitoCarta 3.0 database to identify MRGs, and then, performed a cross-comparative analysis to obtain mitochondria-related DEGs. The Venn diagram shows that 32 overlapping MRGs were identified in the SD and MCI datasets (Figure 2C). The majority of the mitochondria-related DEGs were found to be correlated with each other and to exhibit a high and substantial degree of correlation in the SD and MCI datasets, according to our correlation analysis conducted to examine the correlation of gene expression between each of the DEGs (Figures 2D, E).

Furthermore, using enrichment analysis of these overlapping genes in transcription factor targets, we identified that the targets of the 32 MitoDEGs were primarily regulated by the zinc finger transcription factor ZNF85 (GO: M30399), transcriptional repressor MXD1 (GO: M40757), and master neural transcription factor BRN2 (GO: M12934), as shown in Figure 2F. Subsequently, we performed GO functional annotation and KEGG pathway enrichment analyses to further understand the biological functions and signaling pathways involved in these DEGs. Figure 2G showed the top 10 GO terms for BP, CC, and MF, respectively. GO analysis showed that mitochondrial transport (GO:0006839) was the most significant biological process, whereas NADH dehydrogenase activity (GO:0003954) and thiolester hydrolase activity (GO:0016790) were the most significant molecular functions. In terms of CC, these genes were mainly distributed in the mitochondrial matrix (GO:0005759).

As shown in Figure 2H, KEGG pathway analysis results revealed the 32 MitoDEGs were mainly enriched in pathways of neurodegenerative diseases, for example, Parkinson’s disease (hsa05012), prion disease (hsa05020), pathways of neurodegeneration-multiple diseases (hsa05022), Huntington’s disease (hsa05016), amyotrophic lateral sclerosis (hsa05014) and Alzheimer disease (hsa05010). The identification of these DEGs could help us obtain a key entry point for studying the cognitive impairment induced by SD, which could lead to an improved understanding of the underlying mechanisms of SD-induced cognitive impairment and facilitate the identification of new therapeutic targets.

Based on the aforementioned analysis, it was found that these mitochondria-related DEGs may affect the disease process of cognitive impairment induced by SD through different functions and pathways. Therefore, based on the 32 MitoDEGs, we combined two different machine learning algorithms, LASSO and RF, to screen key DEGs from the SD and MCI datasets. The LASSO logistic regression algorithm identified 27 candidate genes from the 32 MitoDEGs based on the SD dataset GSE39445, and five candidate genes were identified using the RF algorithm (Figures 3A–C). Similarly, based on the MCI dataset GSE140829, 18 candidate genes were identified using LASSO regression, and five candidate genes were identified using the RF algorithm (Figures 3D–F). A venn diagram was plotted to present the overlapping candidate genes among the aforementioned four groups, and ATP synthase mitochondrial F1 complex assembly factor 1 (ATPAF1), a nuclear gene encoding a 31–32-kDa mitochondrial protein essential for ATP synthase F1 assembly and ATP synthase activity, was identified as the key candidate gene (Figure 3G).

Subsequently, we investigated ATPAF1 expression in the SD and MCI datasets. The analysis of expression differences indicated a significant downregulation of ATPAF1 in SD samples in both the training and validation sets (p < 0.05; Figures 4A, E). In addition, analysis of the MCI datasets showed that the expression of ATPAF1 was significantly downregulated in both the training and validation sets (Figures 4C, G), suggesting that this gene may play an important role in both cognitive impairment and sleep loss. Moreover, ROC curves were displayed and AUC was computed to distinguish patients with SD and MCI from controls. With AUCs of 0.572 and 0.578 in the SD datasets (Figures 4B, F), and 0.635 and 0.587 in the MCI datasets (Figures 4D, H), ATPAF1 showed a favorable diagnostic value.

As ATPAF1 was the key overlapping candidate gene screened by two different machine learning algorithms and its expression was downregulated in both the SD and MCI groups, we speculated that ATPAF1 might be the primary driver among the DEGs and was likely to be an effective target for drug therapy. Therefore, we investigated the biological role of ATPAF1 to explore its potential mechanism in the development of cognitive impairment induced by SD. GSEA was performed to investigate the function and pathways of ATPAF1 in the SD and MCI datasets. According to GSEA analysis results, in the SD dataset, “olfactory receptor activity” and “sensory perception of smell” were the negatively enriched GO terms (Figure 5A), while “olfactory transduction” and “maturity onset diabetes of the young” were the negatively enriched KEGG pathways (Figure 5B). Furthermore, we found that in the MCI dataset, ATPAF1 was positively enriched in “structural constituent of ribosome” and “ribosomal subunit” based on GO terms (Figure 5C), as well as in “ribosome” and “oxidative phosphorylation” based on KEGG pathways (Figure 5D), thus indicating that ATPAF1 is involved in regulating ribosomal protein synthesis and oxidative phosphorylation.

Next, GSVA was performed to investigate the possible biological functions of ATPAF1 (Figure 6). Among the enriched GO terms and KEGG pathways, we focused on “acetylcholine metabolic process” in the SD dataset, and “septin cytoskeleton” and “valine leucine and isoleucine degradation” in the MCI dataset. As one of the first neurotransmitters identified in the central nervous system, acetylcholine (ACh) plays a crucial role in learning and memory (Huang et al., 2022). It promotes the conduction of brain nerves and accelerates information transmission, whereas a marked reduction in Ach levels is one of the major characteristics of people with age-related memory loss and AD (Arendt and Bigl, 1986; De Jaeger et al., 2013). Septins, which are enriched in the mammalian nervous system, have been found to be associated with Tau-based paired helical filament core, and contribute to the formation of neurofibrillary tangle in AD (Kinoshita et al., 1998; Hall and Russell, 2004). Furthermore, the effects of the branched-chain amino acids isoleucine, leucine, and valine on cognitive function have been studied previously. Isoleucine, leucine, and valine are significantly associated with the risk of dementia and AD owing to their regulation of the phosphorylation of Tau protein (Tynkkynen et al., 2018).

Using the ssGSEA method, we examined the percentage of 28 immune cell infiltrations in each sample of the SD dataset in order to examine the differences in immune cell composition between SD and control samples and understand the potential immune mechanisms. The violin plot of immune cell infiltration showed significant differences in T.follicular helper cells, type 17 T helper cells, and effector memory CD8 T cells between the SD and control samples (Figure 7A). Next, the relationship between ATPAF1 expression levels and immune cell abundance was analyzed, which showed that ATPAF1 expression was significantly correlated with 18 types of immune cells, such as type 17 T helper cells and effector memory CD8 T cells (Figure 7B). Similarly, the degree of immune cell infiltration in MCI and control samples was further explored. As shown in Figure 7C, the proportions of activated CD4 T cells, activated CD8 T cells, immature dendritic cells, effector memory CD4 T cells, and effector memory CD8 T cells were significantly lower in the MCI group than in the control, whereas the proportions of activated dendritic cells, MDSC, monocytes, natural killer T cells, and T.follicular helper cells were significantly higher in the MCI group than in the control group. The correlation analysis showed results that ATPAF1 expression was significantly correlated with 11 types of immune cells (Figure 7D). Notably, ATPAF1 was positively correlated with cells such as activated CD4 T cells, activated CD8 T cells, and effector memory CD4 T cells that were significantly lower in number in MCI samples, and negatively correlated with cells such as activated dendritic cells that were significantly higher in number in the MCI group.