Aggrephagy-related gene signature correlates with survival and tumor-associated macrophages in glioma: Insights from single-cell and bulk RNA sequencing

1.   Introduction

Glioma is the predominant malignant intracranial neoplasm affecting adults, constituting approximately 50% of all brain tumors. Its global incidence ranges from 4.67 to 5.73 per 100,000 individuals ^[1,2]. The quality of life and prognosis for patients with high-grade gliomas such as glioblastoma (GB) remains poor despite advances in clinical research. GB has a median survival time of 14.6 months and a high recurrence rate despite various treatments, including surgery, chemotherapy, and radiotherapy ^[3]. Currently, the treatment choices and prognostic assessment for patients with glioma predominantly rely on the WHO grading system and molecular subtypes characterized by IDH mutation, 1p19q codeletion, MGMT promoter methylation. Nevertheless, the current markers fail to fulfill the demands for guiding personalized therapies and predicting survival outcomes. Consequently, the need to identify more accurate and personalized biomarkers is paramount.

Clearing misfolded proteins and protein aggregates are essential for maintaining cellular homeostasis ^[4]. The formation of protein aggregates caused by misfolded proteins is linked to several human diseases, including neurodegenerative diseases, eye lesions, and type Ⅱ diabetes ^[5,6,7]. These abnormal protein aggregates can disrupt cellular homeostasis and cause cell death. Eukaryotes physically clean protein aggregates through the formation of aggresomes, a process known as aggrephagy ^[8]. The aggresome formation initiates with the ubiquitination of misfolded proteins, followed by their transport to the microtubule-organizing center (MTOC) facilitated by dynein motor proteins. At the MTOC, the cargo is encapsulated by intermediate filament proteins, resulting in the assembly of the aggresome structure. Subsequently, chaperones are recruited to the aggresome, facilitating its autophagic clearance ^[9]. Understanding of aggrephagy has developed over the years. Previous studies suggested that aggrephagy was mediated by the interaction between receptor molecule p62/SQSTM1 and the Atg8/LC3 family of proteins on the autophagosome membrane, forming autophagosome and degraded by fusion with lysosomes ^[8]. However, Ma et al. found that the molecular chaperone protein CCT2 explicitly regulates the clearance of solid aggregates, suggesting the existence of intracellular molecular pathways that target the clearance of different types of aggregates ^[10].

Aggrephagy plays a crucial role in protein quality control (PQC), managing protein synthesis, folding, and degradation. Cancer cells, known for their rapid division and high activity, rely on aggrephagy to manage proteostatic stress ^[11]. Aggresomes, a key component in this process, act as a last defense in PQC, particularly under overwhelming cellular stress ^[12,13]. Aggresome formation is a noted drug resistance mechanism in overcoming proteotoxic stress caused by proteasome inhibition-based therapy in various cancers, including multiple myeloma and pancreatic cancer ^[14,15]. In brain tumors, the presence of aggresomes correlates with tumor malignancy. For instance, Choroid plexus tumors (CPTs) with over 25% aggresome-positive cells exhibit higher aggression ^[16]. Furthermore, Laura et al. indicate that AQP4 disaggregation may potentiate the invasiveness potential of gliomas ^[17]. Significant studies have investigated the association between gliomas and autophagy and associated prognostic models ^[18,19,20]. However, aggrephagy, a specific form of selective autophagy that targets the breakdown of misfolded protein aggregates, remains poorly understood in glioma.

In the current study, we concentrated on the genes involved in regulating aggrephagy, the aggrephagy-related genes (ARGs). We screened out prognostic ARGs to construct an aggrephagy-related gene prognostic signature (ARPS) and evaluated the reliability of ARPS in 2 independent cohorts and characterized the clinical, molecular, and TME profiles. Based on single-cell sequencing data derived from 9 GB tumor samples, we investigated the expression landscape of ARGs across different cell types and the impact of ARGs on macrophages. We identified multiple aggrephagy clusters of tumor-associated macrophages with distinct intercellular communication with malignant cells, immune characteristics, and metabolic activities. Based on our integrated single-cell and bulk RNA sequencing analysis, our study first demonstrates that ARPS is a promising indicator of prognosis, TME characteristics, and immunotherapy response for gliomas. Furthermore, ARGs may direct intercellular communication of macrophages with tumor cells to influence glioma progression.

2.   Materials and methods

2.1.   Data collection

Clinical information and bulk RNA sequencing data were obtained from the Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/) database. After excluding patients with missing survival data, 669 patients from the TCGA database were enrolled for analysis. Furthermore, 313 patients (mRNAseq_325 dataset) and 657 patients (mRNAseq_693 dataset) were extracted from the Chinese Glioma Genome Atlas (CGGA) database (www.cgga.org.cn/), both of which were validation cohorts. Besides, we downloaded the RNA-seq data for 1152 normal brain tissues from the Genotype-Tissue Expression (GTEx; https://gtexportal.org/home/) database. We processed the RNA-seq data by the STAR pipeline and transformed it into transcripts per kilobase million (TPM) value. Each gene's Expression level was log2(TPM+1). Table S1 summarizes the clinicopathological characteristics of three cohorts of gliomas.

2.2.   Acquisition and functional enrichment analysis of ARGs

By searching the Reactome database, we obtained 44 aggrephagy-related genes, namely ARGs (IdentifierR-HSA-9646399, Table S2). Interaction networks between 44 ARGs were analyzed with geneMANIA (http://www.genemania.org). The Gene Ontology (GO) analyses were conducted for the functional annotation of ARGs via the "ClusterProfiler" package ^[21] of R.

2.3.   Development and validation of the aggrephagy-related gene prognostic signature (ARPS)

The ARGs in the training cohort were incorporated into the least absolute shrinkage and selection operator (LASSO) regression using the R package "glmnet". The minimum 10-fold cross-validation was used to select the optimal value of λ, and 12 ARGs were identified. Subsequently, the multivariate proportional hazards (Cox) regression analysis of these 12 ARGs was performed to construct the ARPS. The calculation formula of ARPS is $\sum _{i = 1}^{n}Expi\mathrm{*}coefi$, where xi and Coefi refer to the expression level of selected ARGs and the corresponding coefficient during multivariate Cox regression analysis, respectively. The median ARPS was taken as the cut-off to stratify different ARPS subgroups. Using the R package "survminer", the Kaplan-Meier curve with log-rank test was plotted to compare overall survival (OS) between ARPS subgroups. To assess the predictive capability of the ARPS, we performed receiver operating characteristic (ROC) curve analysis utilizing the 'timeROC' R package. Subsequently, the area under the curve (AUC) value for one year, three years, and five years were calculated. Variables with p-values less than 0.1 from the univariate Cox regression analysis were included in the multivariate Cox regression analysis to evaluate the independent prognostic value of ARPS. The robustness of ARPS was validated in 2 independent cohorts from CGGA cohorts.

2.4.   Construction and validation of a nomogram scoring system

The independent prognostic clinical characteristics and ARPS were developed using the 'rms' package to develop a predictive nomogram. In the nomogram, each variable was matched with a score, and the total score was obtained by adding the scores across all sample variables. Time-dependent ROC curves for 1-, 3-, and 5-year survival were used to assess the nomogram. Calibration plots of the nomogram were used to depict the predictive value between the predicted 1-, 3-, and 5-year survival events and the virtually observed outcomes. The decision curve analysis (DCA) curves were used to compare the predictive accuracy of the nomogram, ARPS, WHO grade, IDH status, and age by the 'ggDCA' package.

2.5.   Comprehensive analysis of molecular and immune characteristics in ARPS subgroups

We carried out Gene set enrichment analysis (GSEA) to create a series of all genes indicated by their association with ARPS. The gene set enrichment analysis (GSEA) method was based on the "c2.cp.kegg.v7.4.symbols.gmt" gene sets. Each patient's immune, stromal, and ESTIMATE scores were calculated using the ESTIMATE algorithm with the "Estimate" package ^[22]. We incorporated the expression data of 669 glioma samples to CIBERSORT (https://cibersort.stanford.edu/) and iterated 1000 times to estimate the relative proportion of 22 types of immune cells ^[23]. To mitigate the limitations of a single algorithm, we employed the IOBR package ^[24] in conjunction with a suite of algorithms, including TIMER, EPIC, MCPcounter, xCell, and quantiseq, to assess the differences in immune infiltration between high and low ARPS groups.

2.6.   Mutation and immunotherapy response of different ARPS subgroups

We used 'Maftools' R package to determine the different mutation profiles, and tumor mutation burden (TMB) of diffent ARPS subgroup patients. compute_gold_standards function of "Estimate Systems Immune Response (EaSIeR)" R package ^[25] was used to compute the different hallmarks of immune response. The hallmarks included in our study were Cytolytic activity (CYT) ^[26], Tertiary lymphoid structures signature (TLS) ^[27], IFNy signature (IFNy) ^[28], Expanded immune signature (Ayers_expIS) ^[28], T-cell inflamed signature (Tcell_inflamed) ^[28], Roh immune score (Roh_IS) ^[29], Davoli immune signature (Davoli_IS) ^[30], Chemokine signature (chemokines) ^[31], Tumor Immune Dysfunction and Exclusion (TIDE) algorithm ^[32]. The specific information and implications of these biomarkers, as well as their association with immunotherapy response, have been included in the Table S3. Then, the correlation between ARPS and above response biomarkers were analyzed.

2.7.   ScRNA-seq data collection and visualization of TME cell types

We obtained Single-cell mRNA sequence (scRNA-seq) data from GSE131928 10X (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc = GSM3828673) ^[33] dataset. Expression values were calculated by RSEM v1.2.3 in paired-end mode using parameters "--estimate-rspd --paired end -sam -p 6", from which TPM values for each gene were extracted, as described in the study. MAESTRO (Model-based AnalysEs of Single-cell Transcriptome and RegulOme) pipeline was adopted to perform quality control, clustering, and cell-type annotation, as described in the TISCH database ^[34]. The Idents and DimPlot functions were applied to annotate and visualize cells of major TME cell types or subtypes.

2.8.   Pseudotime trajectory analysis of aggrephagy regulators for macrophages

We used the "Monocle" R package (Version 2.0) to investigate the relationship between cellular pseudotime trajectories of macrophages and ARGs ^[35]. Highly variable genes were filtered according to mean expression ≥ 0.1 and empirical_dispersion ≥ 1*fit_dispersion. Dimensionality reduction was performed using the DDRTree method. Next, the 'plot_pseudotime_heatmap' function was used to visualize a heatmap showing the dynamic expression of ARGs in pseudotime trajectories of macrophages in gliomas.

2.9.   Non-negative matrix factorization of aggrephagy regulators in macrophages

We first performed a dimension reduction analysis for 41 aggrephagy regulators in macrophages using the "non-negative matrix factorization algorithm" of the NMF R package (Version 0.20.6) to observe the effect of aggrephagy regulator expression on TME cell types. We identified different aggrephagy subtypes depending on the scRNA expression matrix. We performed all these steps like those in published studies ^[36,37].

2.10.   Cell-cell communication analysis for aggrephagy subtypes of macrophages

CellChat contains databases of human and mouse ligand-receptor interactions that can analyze the intercellular communication networks of different cell clusters ^[38]. We first used CellChatDB.human to assess the significant signaling inputs and outputs between all NMF macrophage cell clusters and malignant cells. Then, we visualize the strength or weakness of the cell-cell communication network between different macrophage clusters and malignant tumor cells by the netVisual_circle function. The AddModuleScore function scored genetic features associated with inflammation and M1/M2 polarization. The scMetabolism package was used to quantify the metabolic activity of different cell clusters at a single-cell resolution ^[39].

2.11.   Statistical analysis

We performed comparisons of continuous variables between the two groups using an independent t-test. We tested categorical data with the χ² test. Kaplan-Meier survival analysis with a log-rank test was used for univariate survival analysis. Multivariate survival analysis was conducted with the Cox regression model. A two-sided P < 0.05 was considered significant.

3.   Results

3.1.   The landscape of aggrephagy-related genes in TME cells in glioma

First, the GSE131928 10X dataset was used for the exploration of the ARGs landscape. The dataset contained 13,558 cells, which were categorized into 26 distinct clusters, and annotated with eight cell types, including Astrocyte-like (AC-like), Mesenchymal-like (Mes-like), Neural Progenitor Cell-like (NPC-like), Oligodendrocyte Progenitor Cell-like (OPC-like) Malignant cells, exhausted CD8⁺T cells (CD8Tex), Malignant, Monocyte/Macrophage (Mono/Macro), and oligodendrocyte from 9 GB patients (Figure 1A, B). Besides, the expressions of 41 ARGs detected in the dataset were indeed different among eight cell types in glioblastoma tissue (Figure 1C). GO analysis of the 44 ARGs (Figure 1D) showed an enrichment of GO terms indicative of cytoskeleton and protein binding such as cytoskeleton-dependent intracellular transport (GO:0030705), microtubule (GO:0005874), cytoplasmic dynein complex (GO:0005868), ubiquitin protein ligase binding (GO:0031625). Besides, we confirmed that the corresponding ARGs physically interact using the protein-protein interaction networks from the GeneMANIA tool (Figure 1E).

Figure 1.  Overview of aggrephagy-related genes (ARGs) in the scRNA-seq for glioblastoma. A. Clusters annotations using the Seurat t-distributed stochastic neighbor embedding (t-SNE) plot of 13,558 cells in 9 GB patients. B. t-SNE plot colored by various cell clusters. C. heatmap distribution of 44 ARGs in AC-like malignant cells, CD8Tex cells, malignant cells, Mes-like malignant cells, monocytes/macrophages, NPC-like malignant cells, oligodendrocytes, and OPC-like malignant cells. D. The Gene Ontology (GO) analysis for the functional annotation of ARGs via the "ClusterProfiler" package. E. The protein-protein interaction networks between ARGs using the GeneMANIA tool.

DownLoad: Full-Size Img PowerPoint

3.2.   Establishment and validation of a five-gene prognostic signature based on aggrephagy regulators

Subsequently, we integrated data from the TCGA database and the GTEX database to compare the expression profiles of ARGs in glioma and normal brains. As illustrated in Figure 2A, the majority of ARGs exhibited upregulated expression in glioma tissues, while a minority of ARGs were found to be downregulated. To assess the role of ARGs genes in assessing glioma prognosis, we screened out 12 OS-related genes (TUBA1C, HSF1, TUBB2A, ARL13B, VIM, TUBA4B, DYNLL2, TUBB6, DYNC1I1, TUBB4A, RPS27A, DYNC1H1) using LASSO regression analysis with 10-fold cross-validation under optimal conditions in TCGA-GBMLGG cohort (Figure 2B, C). We then analyzed these 12 genes using multivariate Cox regression to identify independent predictors for OS. Consequently, 5 ARGs (HSF1, DYNC1H1, DYNLL2, TUBA1C, TUBB6) were independent predictors for OS (p < 0.05). Figure S1 reveals the distinct expression patterns of selected ARGs across different cell types, glioma tissues, normal brain tissues, and across different grades of glioma, suggesting their potential role in glioma progression. Then, an aggrephagy-related gene prognostic signature (ARPS) was constructed and calculated as ARPS = Expression_HSF1*(0.29837)-Expression_DYNC1H1*(0.45814) + Expression_DNYLL2*(0.68889) + Expression_TUBA1C*(0.63961) + Expression_TUBB6*(0.18359). With the median ARPS as the cut-off value, we categorized glioma patients into different ARPS-high (n = 335) and ARPS-low (n = 334) subgroups. As shown in Figure 2D, the ARPS-high group had more patients with dead status. The gene expressions in ARPS significantly differ between high and low-risk groups. The time-dependent ROC curves showed the solid ability of ARPS to predict OS (1-year AUC = 0.865, 3-year AUC = 0.894, 5-year AUC = 0.840) (Figure 2E). There was a shorter OS in ARPS-high patients (log-rank p value < 0.05, Figure 2F). We evaluated the prognostic ability of ARPS across different clinical subgroups by stratifying patients into high and low-risk groups based on median ARPS, followed by Kaplan-Meier analysis (Figure S2). The results indicate that patients with high ARPS exhibit poorer prognosis compared to those with low ARPS across different clinical subgroups, except for patients with 1p/19q codeletion (HR = 1.76, p = 0.234).

Figure 2.  Establishment of the aggrephagy-related genes prognostic signature (ARPS) in the TCGA GBMLGG cohort. A. Expression heatmap of ARGs in Glioma Tissues and Normal Brain Tissues from TCGA and GTEX Databases. (B, C) The LASSO regression was performed with the minimum criteria. D. The distribution plots of ARPS, survival status, and expression of 5 selected ARGs. E. Time-dependent ROC curves to predict the 1-, 3-, and 5-year survival according to the ARPS. F. Kaplan-Meier analysis of the OS between high ARPS and low ARPS groups (HR = 5.43, 95%CI: 4.05–7.28, p < 0.001).

DownLoad: Full-Size Img PowerPoint

3.3.   External validation of the robustness of the ARPS

To validate the robustness and accuracy of the ARPS, we calculated ARPS for rach patient in 2 validation cohorts (CGGA-325 and CGGA-693 cohorts) using the same formula above. In each cohort, the patients were divided into a high-ARPS group and a low-ARPS group with the median ARPS as the cut-off value. The distribution plot of the ARPS revealed that survival time decreased with an increase in ARPS (Figure 3A, D). A significantly shorter OS was observed in ARPS-high patients (Figure 3B, E). The ROC curves confirmed the robustness of ARPS to predict OS in both CGGA325(1-year AUC = 0.724, 3-year AUC = 0.820, 5-year AUC = 0.837) and CGGA-693 cohort (1-year AUC = 0.647, 3-year AUC = 0.694, 5-year AUC = 0.708) (Figure 3C, F). Additionally, we conducted a comparative analysis of APRS with commonly utilized prognostic factors such as tumor grade, IDH status, 1p/19q codeletion, age, and MGMT promoter methylation status in both the training and validation cohorts. The concordance (C-index) values of these factors were presented in Table S4. Our findings revealed that APRS exhibited high concordance in predicting patient prognosis, indicating its potential applicability in clinical settings.

Figure 3.  Validation of the ARPS in two independent cohorts. A, D. The distribution plots of ARPS, survival status, and expression of 5 selected ARGs in CGGA325 and CGGA693 cohorts. B, E. Kaplan-Meier curves of ARPS subgroups in CGGA325 and CGGA693 cohorts. C, F. ROC curves for predicting 1-, 3-, and 5-year OS in CGGA325 and CGGA693 cohorts. OS, overall survival; ROC, receiver operating characteristic.

DownLoad: Full-Size Img PowerPoint

3.4.   Clinical correlation analysis and nomogram construction

To investigate the impact of ARPS on clinical characteristics, we explored the correlation between ARPS and clinical features (age, gender, tumor grade, IDH status, 1p/19q codeletion status, and primary therapy outcome) in TCGA-GBMLGG cohort. We observed significantly higher ARPS in grade 4, IDH WT, 1p/19q non-codeletion, age > 60, and patients with poorer therapy outcomes (Figure 4A–F). The Figure S3 illustrates the relationship between ARPS and other clinical features, highlighting its potential to provide distinct classification perspectives for glioma patients. Besides, we performed univariate and multivariate Cox regression analysis, which confirmed ARPS as an independent prognostic factor after adjusting for other clinicopathological characteristics (HR = 1.552, 95%: 1.246–1.932, P < 0.001, Table 1). To provide a more accurate assessment of patient prognosis by incorporating clinical information and the ARPS, we incorporated the ARPS and independent clinicopathological parameters (Age, WHO grade, IDH status) to establish a nomogram predicting the 1-, 3-, 5-year OS (Figure 4G). The ROC curves confirmed the excellent capability of a nomogram to predict OS. In the training cohort, the nomogram's 1-, 3-, and 5-year AUC values were 0.888, 0.930, and 0.883, respectively (Figure 4H). The calibration plot illustrated that the nomogram achieved good consistency between the predicted and observed OS outcomes (Figure 4I). The DCA curves illustrated that the nomogram had the highest net benefits in predicting patient OS (Figure 4J–L). The evaluation of the proportional hazards assumption for the nomogram yielded a p-value greater than 0.05, indicating the nomogram's validity for prognostic predictions.

Figure 4.  Clinical correlation analysis of ARPS and construction of a nomogram. A–F. ARPS in different clinical features: age, sex, primary therapy outcome, IDH status, 1p/19q codeletion, and WHO grade (*p < 0.05, **p < 0.01, ***p < 0.001, ns No significance) G. Nomogram for predicting the 1-, 3-, 5-year OS of glioma patients in TCGA cohort. H. ROC curves of the nomogram for predicting the 1-, 3-, and 5-year ROC curves in the TCGA cohort. I. Calibration curve of the nomogram. J–L. DCA curves showed the net benefits of the nomogram, ARPS, WHO grade, age, and IDH status in predicting 1-, 3-, and 5-year OS.

DownLoad: Full-Size Img PowerPoint

Table 1.  Univariate and Multivariate analysis of ARPS in TCGA GBMLGG cohort.

Characteristics	Univariate analysis		Multivariate analysis
	Hazard ratio (95% CI)	P value	Hazard ratio (95% CI)	P value
ARPS	2.718 (2.428–3.043)	< 0.001	1.552 (1.246–1.932)	< 0.001
WHO grade
G2	Reference
G3	2.999 (2.007–4.480)	< 0.001	1.823 (1.157–2.874)	0.010
G4	18.615 (12.460–27.812)	< 0.001	2.504 (1.359–4.615)	0.003
IDH status
WT	Reference
Mut	0.117 (0.090–0.152)	< 0.001	0.346 (0.218–0.551)	< 0.001
1p/19q codeletion
codel	Reference
non-codel	4.428 (2.885–6.799)	< 0.001	1.246 (0.731–2.125)	0.419
Gender
Female	Reference
Male	1.262 (0.988–1.610)	0.062	1.044 (0.788–1.383)	0.764
Age
< =60	Reference
> 60	4.668 (3.598–6.056)	< 0.001	1.481 (1.081–2.030)	0.015

---

 | Show Table

DownLoad: CSV

3.5.   GSEA analysis and TME characteristics of different ARPS subgroups

To elucidate the potential mechanism of the excellent predictive capability of ARPS, we further investigated biological pathways related to ARPS (Figure 5A, B). The gene set enrichment analysis (GSEA) revealed that the gene sets of patients with high ARPS scores were enriched in DNA replication, cell cycle, mismatch repair, and antigen processing and presentation (p < 0.05, FDR < 0.025). Using the ESTIMATE algorithm, we found that patients with high ARPS scores had higher immune scores, stromal scores, ESTIMATE scores, and lower tumor purity (Figure 5C–F). It indicated that ARPS positively correlated with immune cell infiltration. We then used the CIBERSORT algorithm to assess the infiltration level of various immune cells for glioma patients. Correlation analysis showed that the ARPS was positively related to infiltration of macrophages (M0, M2, M1), neutrophils, regulatory T cells (Tregs), T cells gamma delta, and resting NK cells but was negatively associated with monocytes, CD4 naïve T cells, activated NK cells, eosinophils, and activated mast cells (Figure 5G). Similarly, we found high ARPS patients had a higher fraction of macrophages (M0, M2, M1), neutrophils, Tregs, T cells gamma delta, plasma cells, B cells naïve, T cells CD4 memory resting but a lower fraction of eosinophils, mast cells, monocytes, CD4 naïve T cells, activated NK cells, eosinophils, and activated mast cells (Figure 5H). Furthermore, we utilized various immune infiltration algorithms, including TIMER, EPIC, MCPcounter, xCell, and quantiseq, to validate our findings (Figure S4). We observed that the high ARPS group exhibited an overall increase in immune infiltration, with a consistent augmentation of macrophages across multiple algorithms.

Figure 5.  Molecular characteristic and immune cell infiltration in high and low ARPS subgroups. A. Gene set enriched in the ARPS-high subgroup (p < 0.05, FDR < 0.25). B. Gene set enriched in the ARPS-low subgroup (p < 0.05, FDR < 0.25). C–F. Differences among stromal score, immune score, ESTIMATE score, and tumor purity, respectively. G. Correlation between ARPS and infiltration level of immune cells. H. The fractions of different immune cells between high- and low-ARPS groups. (*p < 0.05, **p < 0.01, ***p < 0.001, ns Not significant).

DownLoad: Full-Size Img PowerPoint

3.6.   Mutation and immunotherapy response of different ARPS subgroups

To get further insight into the immunologic nature of ARPS-defined groups, we compared and found significant different genetic mutation profiles (Figure 6A, B). The top 10 genes that exhibited the most significant differences were IDH1, EGFR, PTEN, CIC, ATRX, RB1, NF1, KEL, TTN, and AHNAK2. Patients with high-ARPS had markedly higher frequencies of EGFR, PTEN, RB1, NF1, KEL, TTN, AHNAK2 mutation compared to ARPS-low patients. However, the exact opposite was observed regarding the mutation levels of IDH1, ATRX, CIC (Figure 6C). Then, we further explored potency of ARPS to be a biomarker of immunotherapy benefits. First, we found significantly higher TMB in ARPS-high patients (Figure 6D, E). Besides, as shown in Figure 6F–G, there was significant correlation between ARPS and known predictive hallmarks of immune response. Among them, ARPS demonstrated a positive correlation with most predictive signatures, except for the chemokine signature. Considering the association between higher ARPS scores and increased tumor mutational burden (TMB), as well as its relevance to predictive signatures, it is postulated that patients with high ARPS scores may have enhanced responsiveness to immunotherapy.

Figure 6.  Mutation profiles and immunotherapy response of different ARPS subgroups. A, B. The mutation profile in ARPS-low and ARPS-high subgroups, respectively. C. The top 10 genes with the most significant mutation differences between ARPS subgroups. D. Spearman correlation analysis of the ARPS and TMB. D. TMB in different ARPS subgroups. E. Correlation between ARPS and known predictive signatures of immune response including Cytolytic activity (CYT), Tertiary lymphoid structures signature (TLS), IFNy signature (IFNy), Expanded immune signature (Ayers_expIS), T-cell inflamed signature (Tcell_inflamed), Roh immune score (Roh_IS), Davoli immune signature (Davoli_IS), Chemokine signature (chemokines), and Tumor Immune Dysfunction and Exclusion (TIDE). F. The Normalized values of different signatures in ARPS-high and ARPS-low groups. (*p < 0.05, **p < 0.01, ***p < 0.001, ns Not significant).

DownLoad: Full-Size Img PowerPoint

3.7.   ARGs-mediated macrophages resemble classical features

As seen in Figure 5G, the ARPS had the highest correlation coefficients with the infiltration fraction of macrophages. As the most abundant immune cells in glioma TME, we further explore the effect of aggrephagy-related genes on tumor-associated macrophages using scRNA-seq data. Pseudotime trajectory analysis demonstrated that the aggrephagy regulators were critical in the trajectory process of macrophages (Figure 7A). We obtained six main aggrephagy clusters for 3679 macrophages by NMF clustering, including five clusters with the expression of aggrephagy regulators (TUBB4B+Mac, n = 1172; PARK7+Mac, n = 549; VIM+Mac, n = 301; TUBA1A+Mac, n = 560; DYNLL1+Mac, n = 332) and one cluster consisting of 765 macrophages exhibited a notable low expression of aggrephagy regulators (defined as Non-Aggre mac), shown in Figure 7B, C. The CellChat analysis revealed different ligand-receptor links between these aggrephagy-related macrophage clusters and malignant cells. VIM+ mac had the highest number of interactions with malignant cells and incoming signals among them. In contrast, the Non-Aggre mac showed the lowest interactions with malignant cells (Figure 7D). Similary, Figure 7E illustrated the key incoming and outgoing signals for different cell types and their coordination with each other, as well as certain signaling pathways to drive communication. We found higher coordination in TUBA1A+ Mac, TUBB4B1+ Mac, VIM+ Mac with malignant cells than the Non-Aggre Mac. Furthermore, using the AddModuleScore function of Seurat for different clusters, we found that VIM+ mac was significantly related to the proinflammatory phenotype of macrophages (M1 upregulated) than Non-Aggre Macrophages (Figure 7F–G). Last, we assessed the metabolic activity of these six clusters using the 'scMetabolism' package. As shown in Figure 7H, VIM+Mac, TUBB4B+Mac, and TUBA1A+Mac had higher metabolic activity than Non-Aggre macrophages (Figure 7H). The study design is summarized in Figure 8.

Figure 7.  ARGs-mediated macrophages resemble classical features. A. Pseudotime trajectory analysis revealed that the aggrephagy regulators played a critical role in the trajectory process of macrophages. B, C. NMF clustering of macrophages by ARGs. D. Cell-Cell communications between major aggrephagy-related macrophage cells and malignant cells by Cell Chat analysis. E. Outgoing and incoming signaling patterns in different macrophage clusters by cell chat analysis. F–G. t-SNE plots and module scores of gene signatures related to inflammation and M1/M2 polarization of different macrophage clusters. (M1up1: M1 upregulated; M1dn2: M1 downregulated) H. Metabolic activity of different macrophage clusters using the 'scMetabolism' package.

DownLoad: Full-Size Img PowerPoint

Figure 8.  The flow diagram of this study.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

Despite accumulating evidence highlighting the potential involvement of aggrephagy in neurodegenerative disorders, further investigation regarding its role in glioma remains limited. The relationship between aggrephagy and prognosis and glioma tumor microenvironment remain largely unknown. In the current study, we identified 44 aggrephagy-related genes (ARGs). These genes have different expression patterns across different cell types of gliomas. Next, we screened out 5 ARGs to construct the prognostic signature, ARPS. The robustness and accuracy of ARPS in predicting prognosis were validated in 2 independent cohorts, with worse survival in the ARPS-high group and better survival in the ARPS-low group. Univariate and multivariate Cox regression analysis proved that ARPS was an independent prognostic factor for OS. Additionally, the ARPS outperforms several published autophagy-related gene prognostic models in assessing the prognosis of glioma patients ^[19,20]. These results confirmed that ARPS has stable power of prognosis and would be widely applicable to glioma patients with different clinicopathological characteristics.

The ARPS comprised 5 ARGs, which were HSF1, DYNC1H1, DYNLL2, TUBB6, and TUBA1C. Upon accumulation of protein aggregates, HSF1 facilitates casein kinase 1-mediated phosphorylation of p62, whereas inhibition of HSF1 impairs protein aggregate-induced autophagosome formation ^[40]. HSF1 is often constitutively active in cancer cells. Its downstream target, inducible heat shock protein 70 (Hsp70), has exhibited tumor-specific localization to lysosomes and promotes autophagy ^[41]. In glioma, HSF1 was found to be upregulated and play a significant role in driving the constitutive overexpression of BAG3 and regulating the expression of FoxM1 ^[42,43]. Targeting HSF1 has shown promise in enhancing apoptotic processes, inhibiting the migratory and invasive capabilities of glioma cells, and suppressing tumor growth in murine models ^[42,43]. The microtubule motor protein encoded by DYNC1H1 is involved in many cellular processes, such as mitosis and intracellular transport ^[44]. DYNC1H1 mutations have been implicated in nervous system diseases and pancreatic cancer ^[45], and these mutations are consistent with a high immune activity of tumor mutation load in various cancer types ^[46]. DYNLL2 gene, also known as Dynein Light Chain LC8-Type 2, is a protein-coding gene that plays a crucial role in various cellular processes. It is primarily associated with short-rib thoracic dysplasia 11 (SRTD11) and Bardet-Biedl syndrome (BBS) ^[47]. Jiang found that DYNLL2 played an oncogenic role and indicated a poor prognosis in osteosarcoma ^[47]. Besides, DYNLL2 exhibits a significant negative correlation with the levels of macrophage M0 infiltration ^[47]. Tubulin beta 6 class V (TUBB6) was recognized as a potential mutation hot spot in human colorectal cancers accompanied by microsatellite instability and served as a biomarker for predicting GC peritoneal metastasis ^[48,49]. In glioma, TUBB6's expression aligns with higher tumor grades and poorer prognoses and is intimately related to M0 macrophage infiltration ^[50]. TUBA1C is a type of α-tubulin associated with microtubules participating in cellular mitosis and division progress. Studies suggested that the expression of TUBA1C predicts poor prognoses of multiple cancers, including hepatocellular carcinoma, lung adenocarcinoma, breast cancer, and osteosarcoma ^{[51,52,53,54]}. Similarly, Zhu et al. found that TUBA1C showed high expression in LGG and is related to a poor prognosis ^[55]. TUBA1C co-expressed with immune-related genes and checkpoints, and its expression positively correlates with the infiltration of B cells, CD8+ T cells, CD4+ T cells, and macrophages ^[55]. In summary, the dysregulation of these proteins can contribute to glioma development by affecting various cellular processes, including stress response, microtubule dynamics, and immune cell interactions within the tumor microenvironment. Further research in this area is ongoing to better understand their roles and potential as therapeutic targets in glioma.

To obtain further biological insight into different ARPS subgroups, we performed GSEA and found that the gene sets of ARPS-high patients mostly enriched in the cell cycle, DNA replication, and immune-related pathways. Thus, we further used different algorithms to understand the landscape of the glioma TME in two ARPS subgroups. By the Estimate algorithm, we found higher immune, stromal, and ESTIMATE scores and lower tumor purity in APRS high subgroup. Besides, the composition of immune cells varies significantly between the two subgroups. Specifically, the proportions of macrophages (including M0, M1, and M2), Treg cells, and neutrophils were higher, while the proportion of activated NK cells and monocytes was lower in the ARPS-high subgroup. In most tumors, including breast, bladder, and ovarian cancers, M2 macrophage cell levels have been proven to be associated with chronic inflammation, mediate an immunosuppressive phenotype, and poor prognosis ^[56,57,58]. In contrast, elevated levels of M1 macrophages are associated with acute inflammation and favor an excellent prognosis ^[56,57,58]. In the present study, the proportion of M1 and M2 macrophages increased in the ARPS-high subgroup, with higher elevation of M2 macrophages. Besides, it is worth noted that parameters such as IFN-y signature and cytotoxic activity were found to be higher in high ARPS subgroup, which may be signs of immune activation. It is important to highlight that these scores are calculated based on gene expression levels, there may be compensatory upregulation of genes associated with immune activation as ARPS increases. However, integrating the results obtained from both bulk RNA-seq and single-cell RNA-seq analyses, we consistently observed a prevalent proportion of tumor-associated macrophages in the immune microenvironment of glioma, with M2 macrophages being the predominant subtype. Therefore, we concluded that the ARPS-high subgroup exhibited a higher degree of immunosuppression.

Furthermore, our study investigated the relationship between ARPS and the immunologic nature of glioma patients and their genetic mutation profiles. The study found that patients with high ARPS had significantly different genetic mutation profiles compared to those with low ARPS. We observed higher frequencies of specific mutations (EGFR, PTEN, RB1, NF1, KEL, TTN, and AHNAK2) and lower mutations (IDH1, ATRX, and CIC) in the ARPS-high group. Besides, we found that patients with high ARPS had significantly higher TMB (tumor mutational burden), a known predictor of immunotherapy response ^[59]. Additionally, the research revealed a noteworthy association between ARPS (the specific factor being studied) and several predictive indicators of the immune response, such as Cytolytic activity (CYT), Tertiary lymphoid structures signature (TLS), IFNy signature (IFNy). These findings suggest that ARPS may be a potential biomarker for predicting the treatment response of ICB in glioma patients. However, when interpreting these results, it is important to acknowledge that the efficacy of immunotherapy in glioblastoma remains uncertain, and there is a lack of a definitive indicator that predicts the response to immunotherapy. Although the indicators used in this study have shown effectiveness in other tumor types, including tumor mutation burden, tumor lymphocyte infiltration, and Tumor Immune Dysfunction and Exclusion (TIDE) score, further validation is required to assess their predictive value in immunotherapy response specifically for glioblastoma.

The rapid development of single-cell sequencing provides new perspectives to investigate the relationship between genes of interest and the tumor microenvironment, as well as identifying new cell clusters ^[60,61,62]. Our study found that the expression of ARGs was highly heterogeneous in different cell types. Through NMF clustering, we found that ARGs-mediated macrophage subtypes exhibited different communication characteristics with tumor cells. Notably, the VIM+ Mac subtype showed M1 characteristics compared to macrophages with low VIM expression. It exhibited significant activation of metabolism-related pathways such as steroid biosynthesis, oxidative phosphorylation, and inositol phosphate metabolism. A previous study revealed that vimentin had abundant expression in activated macrophages and foam cells. The vimentin had a potent suppressive effect on oxidative stress, and Vim-/- mice showed increased vascular inflammation, with increased expression of CD36 on macrophages ^[63]. Similarly, Markovitz et al. found that macrophages secreted vimentin in response to proinflammatory signaling pathways and were probably involved in immune function ^[64]. The above findings revealed the existence of VIM+ macrophages in GB, which have higher interactions with malignant cells and exhibit proinflammatory characteristics.

Our study has several limitations that cannot be ignored. First, the selected aggrephagy-related genes may not represent the complete aggrephagy-related gene landscape. The impact of these genes on the pathological process of glioma through aggrephagy remains elusive and warrants further experimental verification. Second, the prognostic signature ARPS was developed and validated in publicly available databases and required confirmation through prospective clinical research. The low depth of single-cell RNA sequencing and limited sample sizes may result in lower expression levels of certain aggrephagy regulators in glioblastoma, leading to increased zero observations and potential biases in the clustering methodology. Additionally, the conclusions drawn from this study would benefit from validation in a larger cohort of GB patients using single-cell RNA sequencing. Moreover, it should be noted that the ICB response indicators employed in this study have primarily been validated in specific cancer types such as melanoma and non-small cell lung cancer (NSCLC). Therefore, further investigation is necessary to establish the association between ARPS and the potential benefits of ICB therapy in glioma. Nonetheless, integrating bulk RNA sequencing and single-cell RNA sequencing provides a practical approach to understanding the role of aggrephagy regulators in glioma and the tumor microenvironment.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

We would like to acknowledge and thank the following databases for providing the data used in this study: TCGA (https://portal.gdc.cancer.gov/), CGGA (https://www.cgga.org.cn/) websites. We also appreciate the resources provided by the GEO database (https://www.ncbi.nlm.nih.gov/geo) for sharing scRNA-seq data. Without their contribution, this study would not have been possible.

Conflict of interest

The authors declare there is no conflict of interest.