Key markers linked to cell cleanup processes identified in sleep deprivation

Three SD datasets were obtained from GEO (http://www.ncbi.nlm.nih.gov/geo/↗), including GSE9442↗, GSE33302↗, GSE9441↗, and GSE37667↗, along with a single-cell RNA sequencing (scRNA-seq) dataset, GSE137665↗ (Jha et al., 2022; Maret et al., 2007; Pellegrino et al., 2012; Vecsey et al., 2012). The GSE9442↗ dataset comprised 9 control samples and 8 SD samples, whereas GSE33302↗ included 12 control and 12 SD samples. The GSE9441↗ dataset, containing nine control and nine SD samples, was excluded from the present analysis because these samples were derived from peripheral tissue (liver). GSE37667↗ contained nine control and nine SD samples, but nine samples collected after recovery sleep were excluded. The GSE137665↗ single-cell dataset included three SD samples and three control samples, with three post-recovery sleep samples excluded (Table 1). The GSE9442↗ and GSE33302↗ datasets were integrated via the removeBatchEffect function from the limma package (ver. 3.60.6). Specifically, the dataset source was defined as the batch effect, while the experimental grouping was included in the design matrix to strictly preserve the biological variance during batch correction. The GSE9441↗ and GSE37667↗ datasets were reserved as independent validation sets for assessing the performance of the diagnostic model. Autophagy-related genes (ARGs) were retrieved from the GeneCards database, and those with a Relevance Score greater than 4 were selected for subsequent analysis as highly associated with autophagy.

Differentially expressed genes (DEGs) linked to SD were identified utilizing “limma” (ver. 3.60.6) via differential expression analysis, with a threshold of adjusted p < 0.05 for gene selection. Subsequently, the DEGs were intersected with 597 ARGs to obtain SD-DE-ARGs. Finally, visualization was conducted using the “ggplot2” (ver. 3.5.1) and “ggvenn” (ver. 0.1.16) packages.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were carried out on SD-DE-ARGs via the “clusterProfiler” package (ver. 4.4.4) to explore their biological functions. A false discovery rate (FDR)-adjusted p < 0.05 denoted statistical significance. Results were visualized via the “circlize” package (ver. 0.4.16) and “ggplot2” package (ver. 3.5.1).

ML algorithms, like Least Absolute Shrinkage and Selection Operator (LASSO) and Support Vector Machine Recursive Feature Elimination (SVM-RFE), were leveraged for SD-DE-ARGs screening. The LASSO algorithm determines the variables by analyzing the minimum classification error λ, which is utilized to select feature variables and develop the optimal classification model. Five-fold cross-validation and regression analysis with tuning/privacy parameters were enabled by the “glmnet” package (ver. 4.1.8) for feature selection, model construction, and handling high-dimensional data (Qiu et al., 2024). SVM-RFE, an approach for feature selection utilizing support vector machines, was conducted using five-fold cross-validation through the “e1071” package (ver. 1.7.16). This method recursively eliminates the least important features to select the most representative subset (Peng et al., 2024). Finally, the intersection of the genes selected by both ML algorithms was defined as the predictor genes for SD. Visualization was performed using the “ggvenn” (ver. 0.1.16) and “ggplot2” (ver. 3.5.1) packages.

Male Sprague-Dawley rats, eight weeks of age and weighing 200–220 g at baseline, were obtained from Guangdong GemPharmatech Co., Ltd. All experimental procedures were conducted in accordance with the guidelines approved by the Institutional Animal Care and Use Committee (IACUC) of the First Affiliated Hospital of Jinan University (Approval No. IACUC-2025080107). Each SD and control group consisted of three randomly assigned rats, determined based on the principle of minimum sample size (n = 3). Randomization was performed to reduce potential selection bias. All cages were maintained on the same rack under identical environmental conditions to control for spatial confounders. However, due to operational constraints, euthanasia and subsequent tissue processing were conducted in a fixed sequence, with the SD group processed before controls. Investigators were not blinded to group allocation during outcome assessment. During acclimation, rats were housed in standard polycarbonate cages (45 cm × 30 cm × 20 cm) at a density of three rats per cage, with bedding replaced twice weekly. The animal facility maintained a 12:12 h light-dark cycle, a temperature of 22 ± 2 °C, and a relative humidity of 50 ± 5%. Environmental enrichment, including PVC tunnels and wooden blocks, was provided and renewed weekly. Standard rodent chow and water were available ad libitum. Rats were acclimated for two weeks before experimentation. The SD model was established according to the method described by Kaushal et al. (2012). The apparatus consisted of a rotatable deprivation bar traversing the cage floor intermittently, compelling the rats to move to avoid contact. When activated, the deprivation bar rotated at ten revolutions per minute, alternating between clockwise 1.5 turns, counterclockwise 1.5 turns, and clockwise 1.5 turns. At the onset of the experiment, rats in both groups were allowed unrestricted movement, feeding, and drinking within the SD chamber. For the SD cohort, the apparatus was activated daily from the fourth hour of the light phase to the twelfth hour of the dark phase for seven consecutive days, restricting sleep to approximately four hours per day. In contrast, the apparatus remained inactive throughout the experimental period in the control cohort. At the end of the deprivation protocol, rats were deeply anesthetized by intraperitoneal injection of 2% pentobarbital sodium (50 mg/kg). Adequate anesthetic depth was confirmed by the absence of pedal reflexes. Euthanasia was performed via cervical dislocation, and death was confirmed by cessation of heartbeat and respiration. The brain was immediately removed and placed on ice-cold PBS-soaked filter paper. For hippocampal dissection, the overlying cortex was carefully peeled away using fine forceps, and the hippocampus was bilaterally excised with removal of the alveus and fimbria. The entire hippocampal formation, including the dentate gyrus and CA1-CA3 subfields, was collected. For prefrontal cortex dissection, the medial prefrontal cortex was targeted based on coordinates from the Paxinos and Watson Rat Brain Atlas, approximately Bregma +4.70 mm to +2.70 mm (Paxinos & Watson, 2006). Coronal slices (∼2 mm) were sectioned using a chilled razor blade. Minor inclusion of the corpus callosum was unavoidable but consistent across all samples. Dissected tissues were immediately snap-frozen in liquid nitrogen and stored at −80 °C. mRNA expression of the predictor genes in these tissues was quantified by RT-qPCR. Animals exhibiting >20% body weight loss, severe lethargy, unhealed wounds, or signs of prolonged distress would have been euthanized before the planned endpoint. Nevertheless, no animals met these humane endpoint criteria. All rats were monitored daily to ensure welfare. Total RNA was extracted from brain tissues using TRIzol reagent (Invitrogen, USA), and RNA concentration was determined spectrophotometrically. Complementary DNA (cDNA) was synthesized using the EntiLink™ 1st Strand cDNA Synthesis Kit (ELK Biotechnology, EQ003) according to the manufacturer’s instructions. Reverse transcription was performed at 25 °C for 5 min, 42 °C for 30 min, and 85 °C for 5 min.

Quantitative analysis of mRNA expression was carried out using the EnTurbo™ SYBR Green PCR SuperMix Kit (ELK Biotechnology, EQ001) on a StepOne™ Real-Time PCR System (Life Technologies, Carlsbad, CA, USA).

The relative gene expression levels were calculated using the 2^−ΔΔCt method, with GAPDH serving as the internal reference. The primer sequences used in the RT-qPCR assays are listed in Table S1.

Expression differences of the predictor genes between the two groups were assessed in the training datasets (GSE9442↗ and GSE33302↗). Two independent datasets, encompassing both target brain and blood samples (GSE9441↗ and GSE37667↗), were employed for validation. Results were visualized via the “ggstatsplot” package (ver. 0.13.0). A diagnostic model for SD was constructed based on multivariate logistic regression using the identified predictor genes. Calibration curves were generated for the diagnostic model. Receiver operating characteristic (ROC) analysis was carried out via “pROC” (ver. 1.18.5), and the diagnostic performance of the model was examined via the area under the curve (AUC) with its 95% confidence interval (95% CI).

The functional relevance of the predictor genes was investigated using gene set enrichment analysis (GSEA). Briefly, Spearman correlation coefficients were computed between each predictor gene and all other genes within the normalized expression matrix. Genes were subsequently ranked in descending order according to these correlation coefficients to generate a pre-ranked gene list. GSEA was then performed using the clusterProfiler R package (ver. 4.12.6). The human-specific GO gene set from the Molecular Signatures Database (MSigDB; “c5.go.v2024.1.Hs.symbols.gmt”) was used as the reference. Pathways with a nominal P value < 0.05 were considered significantly enriched.

To assess SD-induced alterations in the immune microenvironment, single-sample gene set enrichment analysis (ssGSEA) was conducted using the GSVA package in R. The fully normalized, transcriptome-wide expression matrix from the merged SD dataset served as input. A curated gene signature encompassing 28 immune cell subpopulations was applied to estimate the relative infiltration abundance of each immune cell type in individual samples (Charoentong et al., 2017). Pearson correlation analysis was subsequently performed to evaluate associations between immune cell infiltration scores and the expression levels of the predictor genes. The results were visualized using the ggplot2 R package (ver. 3.5.1).

scRNA-seq data obtained from GEO (GSE137665↗) were preprocessed and annotated using the Seurat package (ver. 5.3.0) (Jha et al., 2022). Quality control criteria included detection of 300–2,000 genes per cell and a mitochondrial gene proportion below 10%. Cells passing quality control were normalized and scaled. Highly variable genes were identified via the FindVariableFeatures function. Linear dimensionality reduction was performed via principal component analysis (PCA), followed by nonlinear reduction using UMAP (dims = 1:20). Unsupervised clustering was conducted using the Louvain algorithm with a resolution parameter of 0.5, selected based on clustering stability and biological consistency. Cell type annotation was performed manually based on established marker genes. Differential expression analysis was performed using the FindAllMarkers function (logfc.threshold = 0.25, min.pct = 0.1), and the top 30 marker genes per cell type were selected for functional enrichment analysis. Cell type-specific marker genes and target gene expression were visualized using violin plots and dot plots.

To identify potential candidate drugs targeting key predictor genes, a reverse drug screening was performed using the Drug Signatures Database (DSigDB, Version 1.0; https://dsigdb.tanlab.org/DSigDBv1.0/↗). DSigDB integrates drug-induced gene expression profiles from multiple public databases, including Connectivity Map, LINCS, and DrugBank. Potential drugs and their corresponding predictor gene targets were visualized as interaction networks via Cytoscape (Otasek et al., 2019).

All data analyses were enabled by R (ver. 4.4.1). To evaluate statistical differences across groups, either the Wilcoxon test or the t-test was applied based on data distribution characteristics. Correlation analyses were conducted using Pearson correlation. Each statistical test was two-sided, with p < 0.05 denoting statistical significance.

This animal study gained the approval of the Animal Ethics Committee of the First Affiliated Hospital of Jinan University (Approval No.: IACUC-2025080107) and followed the ARRIVE 2.0 guidelines. All experiments strictly adhered to the 3R principle, minimizing animal use and suffering while ensuring robust statistical outcomes. All experimental procedures and handling protocols were conducted as per the IACUC of the First Affiliated Hospital of Jinan University.

The GSE9442↗ and GSE33302↗ datasets were merged, and batch effect removal was performed, resulting in 20 SD samples and 21 control samples (Figs. 2A–2B). Subsequently, the differential analysis identified 702 upregulated genes and 770 downregulated genes, with heatmaps and volcano plots presented in Figs. 2C–2D. Finally, ARGs were retrieved from GeneCard, and genes with a correlation coefficient greater than 4, strongly linked to autophagy, were selected. The intersection of DEGs and high ARGs yielded 74 SD-DE-ARGs (Fig. 2E).

To elucidate the systemic biological roles of SD-DE-ARGs beyond their canonical functions, GO and KEGG enrichment analyses were performed. In the Biological Process (BP) category, while the genes were primarily clustered in autophagy and macroautophagy regulation, the analysis also revealed their specific involvement in stress adaptation mechanisms, including cellular responses to chemical stress and misfolded proteins. This indicates that SD-DE-ARGs act as critical mediators of cellular quality control under pathological stress. Cellular Component (CC) analysis showed that, in addition to autophagosomes, the gene products were significantly localized to protein kinase complexes and protein-folding chaperone complexes, highlighting their roles in signal transduction and protein stability. Molecular Function (MF) analysis further revealed prominent activities in ubiquitin ligase binding, kinase regulator activity, and SNARE binding, emphasizing their regulatory complexity in protein degradation and vesicle fusion. KEGG pathway analysis uncovered extensive signaling crosstalk and pleiotropic effects. Beyond the expected enrichment in the Autophagy (animal) pathway, SD-DE-ARGs were significantly associated with non-autophagic pathways, including Apoptosis (reflecting regulation of cell death), Lipid metabolism and atherosclerosis (suggesting metabolic implications), and disease-specific pathways such as Prostate cancer and Amyotrophic lateral sclerosis (ALS). These results (Figs. 2F–2I) indicate that SD-DE-ARGs form a complex signaling network connecting autophagy with metabolic regulation and tumorigenesis.

LASSO and SVM-RFE algorithms were employed to further screen the 74 SD-DE-ARGs. Using the LASSO regression algorithm, five diagnostic genes were identified from the SD-DE-ARGs: CDKN1A, HSPA5, NR4A1, OFD1, and PRKAB1 (Figs. 3A–3B). Through the application of the SVM-RFE algorithm, five potential diagnostic biomarkers were pinpointed, including PRKAB1, CDKN1A, NR4A1, HSPA5, and HSPB1 (Fig. 3C). Subsequently, after integrating the feature genes derived from both analytical approaches, four core regulatory genes were identified: CDKN1A, HSPA5, NR4A1, and PRKAB1 (Fig. 3D).

The expression of four predictor genes was investigated across the two groups. In the training sets (GSE33302↗ & GSE9442↗), CDKN1A, HSPA5, and NR4A1 were significantly upregulated in the SD cohort, whereas PRKAB1 was notably downregulated (Fig. 4A). To further validate the differential expression of predictor genes, RT-qPCR was performed on three Sprague-Dawley rats and three normal rats. Compared with the healthy control group, CDKN1A (p < 0.01), HSPA5 (p < 0.05), and NR4A1 (p < 0.01) were significantly upregulated in the SD group, whereas PRKAB11 showed an insignificant difference across the two cohorts (Fig. 4B, Tables S2–S3). Ultimately, based on the significant upregulation confirmed by RT-qPCR and the lack of significance for PRKAB1, we established CDKN1A, HSPA5, and NR4A1 as the final panel of key predictor genes (Fig. 4B).

To evaluate the diagnostic efficacy of the three identified ARGs (NR4A1, HSPA5, and CDKN1A), a diagnostic nomogram was constructed for SD risk assessment (Fig. 4C). The nomogram assigns weighted points according to the regression coefficients of each gene, with NR4A1 contributing the most to the risk score (expression range: 8.9–10.3), followed by HSPA5 (10.5–11.7) and CDKN1A (6.0–6.8). In the training cohort (GSE33302↗ and GSE9442↗), the calibration curve demonstrated strong concordance between predicted probabilities and observed outcomes, with minimal deviation from the ideal diagonal line (Fig. 4D). ROC analysis revealed an AUC of 1.000 (95% CI [1.000–1.000]), indicating excellent discriminatory ability of the three-gene signature within the training dataset (Fig. 4E). External validation using two independent GEO datasets confirmed the generalizability of the diagnostic model. In the first validation cohort (GSE9441↗), the nomogram maintained robust predictive performance, with an AUC of 0.969 (95% CI [0.923–1.000]), and the calibration curve showed good agreement between predicted and observed probabilities (Figs. 4F–4H). In the second validation cohort (GSE37667↗), although performance slightly declined relative to the training and first validation sets, the model still exhibited satisfactory discrimination, with an AUC of 0.778 (95% CI [0.532–1.000]). The calibration curve demonstrated reasonable consistency, despite some deviation from the ideal line in the high-probability region (Figs. 4I–4K).

To unravel the biological roles of the predictor genes, GSEA was performed. The findings indicated that these pivotal genes are likely to be significantly involved in the pathological processes of SD by regulating BP, such as endocytosis, positive regulation of phagocytosis, assisting protein folding, and responding to vascular endothelial growth factor signaling pathways (Figs. 5A–5C).

Correlation analysis of three key predictor genes with immune cells revealed that NR4A1 expression was positively related to effector memory CD8⁺ T cells, myeloid-derived suppressor cells, and macrophages, but negatively linked to effector memory CD4⁺ T and plasmacytoid dendritic cells. HSPA5 showed broad positive correlations with multiple immune cells and a negative relation to effector memory CD4⁺ T and memory B cells. CDKN1A was positively related to effector memory CD8⁺ T cells, but negatively correlated with effector memory CD4⁺ T cells. Therefore, the predictor genes possibly play distinct roles in immune regulation (Fig. 6A).

To characterize the global immune landscape in SD, the relative abundance of the 28 immune cell types was compared between SD and control cohorts. Relative to controls, monocytes, natural killer T cells, CD4⁺ T cells, plasmacytoid dendritic cells, and memory B cells were decreased in the SD group. In contrast, type 2 helper T cells, CD8⁺ T cells, myeloid-derived suppressor cells, central memory CD4⁺ T cells, activated B cells, and macrophages were significantly increased (Fig. 6B).

scRNA-seq data from the SD dataset (GSE137665↗) were processed using a standardized workflow and analyzed with the Seurat package. Quality control metrics confirmed high data fidelity, with a median of ∼500 genes detected per cell and a mitochondrial gene proportion of ∼4% (Fig. S1). Unsupervised clustering via UMAP dimensionality reduction identified 22 transcriptionally distinct clusters, which were annotated into 12 major cell types based on canonical marker gene expression (Figs. 7A–7B). The cellular composition was dominated by endothelial cells and glutamatergic neurons, occupying the lower-left and right regions of the UMAP embedding, respectively.

Cell type identities were validated through established marker genes (Fig. 7C). Endothelial cells specifically expressed Flt1 and Cldn5; pericytes and vascular smooth muscle cells expressed Acta2 and Vtn; astrocytes expressed Aldoc and Slc1a3; oligodendrocytes expressed Cldn11 and Mbp; glutamatergic neurons expressed Snap25 and Slc17a7; GABAergic neurons expressed Dlk1 and Gad1; microglia expressed Cx3cr1 and P2ry12; and macrophages expressed Lyz2 and Cd14.

GO functional enrichment analysis revealed cell type-specific biological signatures (Fig. 7D). Microglia were enriched in chemotaxis, leukocyte migration, and antigen processing and presentation via MHC class II pathways, consistent with their role as CNS immune sentinels. Macrophages similarly displayed enrichment in antigen presentation pathways but exhibited distinct transcriptional profiles reflecting their peripheral origin. Astrocytes were enriched in glial cell differentiation, gliogenesis, and regulation of synapse structure, highlighting neurosupportive functions. Oligodendrocytes showed enrichment in cartilage development and collagen fibril organization, consistent with myelination roles. Glutamatergic neurons were enriched in synaptic structure regulation, cognition, and cilium movement pathways. Endothelial cells exhibited enrichment in vascular endothelial proliferation regulation, angiogenesis, and acetylcholine receptor signaling, reflecting their angiogenic capacity and neurovascular coupling functions.

Notably, the three key predictor genes, CDKN1A, HSPA5, and NR4A1, displayed highly selective, cell type-specific expression patterns (Fig. 7E). CDKN1A demonstrated the most stringent specificity, with robust expression restricted to endothelial cells and negligible expression in all other cell types, including microglia, macrophages, and other immune populations, suggesting a specialized role in vascular homeostasis and endothelial stress responses under SD. HSPA5 exhibited a broader but distinct pattern, with the highest expression in endothelial progenitor cells (EPCs) and moderate expression in oligodendrocytes, neurons, and astrocytes, while remaining absent in microglia and macrophages. NR4A1 was predominantly expressed in neuronal lineages, peaking in glutamatergic neurons and showing moderate expression in GABAergic neurons and other neuronal subtypes, while remaining virtually undetectable in non-neuronal populations, including immune cells.

The three human predictor genes (CDKN1A, HSPA5, and NR4A1) were queried against the Enrichr platform using the DSigDB database to identify potential therapeutic compounds. Drugs with adjusted P < 0.01 were retained, yielding 35 candidate compounds and small molecules targeting these genes, which may represent potential therapeutic strategies for mitigating the adverse effects of SD (Fig. S2, Table S4).

The following information was supplied regarding data availability:

All data generated or analyzed in this study were obtained from GEO: GSE33302↗ and GSE9442↗; GSE9441↗, GSE37667↗; GSE137665↗.

All analysis scripts, including data preprocessing, differential expression analysis, ML pipelines, and visualization routines, are available at GitHub and Zenodo:

– https://github.com/tobejoy613/Autophagy-Genes-as-Sleep-Deprivation-Biomarkers↗.

– tobejoy613. (2026). tobejoy613/Autophagy-Genes-as-Sleep-Deprivation-Biomarkers: Initial release for PeerJ publication (v1.0.0). Zenodo. https://doi.org/10.5281/zenodo.20034930↗.