The ALDH Family Contributes to Immunocyte Infiltration, Proliferation and Epithelial-Mesenchymal Transformation in Glioma.

Gliomas are malignant tumors that originate from the central nervous system. The aldehyde dehydrogenase family has been documented to affect cancer progression; however, its role in gliomas remains largely unexplored. Bulk RNA-seq analysis and single-cell RNA-Seq analysis were performed to explore the role of the aldehyde dehydrogenases family in gliomas. Training cohort contained The Cancer Genome Atlas data, while data from Chinese Glioma Genome Atlas and Gene Expression Omnibus were set as validation cohorts. Our scoring system based on the aldehyde dehydrogenases family suggested that high-scoring samples were associated with worse survival outcomes. The enrichment score of pathways were calculated by AUCell to substantiate the biofunction prediction results that the aldehyde dehydrogenases family affected glioma progression by modulating tumor cell proliferation, migration, and immune landscape. Tumor immune landscape was mapped from high-scoring samples. Moreover, ALDH3B1 and ALDH16A1, two main contributors of the scoring system, could affect glioblastoma cell proliferation and migration by inducing cell-cycle arrest and the epithelial-mesenchymal transition. Taken together, the aldehyde dehydrogenases family could play a significant role in the tumor immune landscape and could be used to predict patient prognosis. ALDH3B1 and ALDH16A1 could influence tumor cell proliferation and migration.

Introduction

Gliomas are the most common cancer type of the central nervous system, representing approximately 80% of malignant brain tumors. Glioblastoma (GBM) is the most aggressive type of gliomas (1, 2), with a reported overall survival time of only 15-18 months (3, 4). However, the current classification based on histology has major limitations and cannot precisely predict patient outcomes (5). Over the past decades, the identification of molecular signatures has provided new insights into glioma prognosis prediction (6). For instance, isocitrate dehydrogenase (IDH) 1/2 mutations and 1p19q co-deletion have been considered glioma prognosis-related biomarkers (7). Interestingly, classification based on these signatures can be used to subdivide GBM into three subtypes (classical, proneural and mesenchymal) and assist during clinical practice for treatment selection and prognosis prediction (8). Since gliomas are highly heterogeneous tumors, it is vital to identify novel biomarkers to help clinicians during the differential diagnosis and make treatment decisions.

Aldehyde dehydrogenases (ALDHs) are a group of enzymes that catalyze the oxidation of aldehydes to less toxic carboxylic acids (9). The human ALDH superfamily consists of 19 putative members distributed in different organs. Multiple studies reported that ALDHs were associated with tumorigenesis and tumor resistance to therapy. In this regard, some ALDHs have been discovered to be biomarkers of cancer stem cells (10, 11). For instance, ALDH1A3 has been reported to participate in extracellular matrix reorganization and cell adhesion in gliomas by inducing mesenchymal transformation (12–14) and is responsible for the aggressiveness of mesenchymal-like glioma stem cells (15, 16). Furthermore, ALDH1A1 has been associated with glioblastoma proliferation and migration (17, 18), while decreased ALDH1A1 expression could restore GBM cell sensitivity to temozolomide, a common and efficient chemotherapeutic drug for GBM (19). Accordingly, the association between ALDHs and cancer is complicated, warranting further studies.

In the present study, the expression profile of ALDHs in gliomas and their different clinical features were explored. A scoring system based on the least absolute shrinkage and selection operator (LASSO) regression analysis calculated the risk score for each sample. The high-risk group showed worse overall survival outcomes relative to the low-risk group. Besides, GO/KEGG enrichment analysis and AUCell analysis from single-cell RNA-Seq analysis corroborated that ALDHs promoted glioma progression by affecting cell proliferation and migration. Notably, we found that two contributors of that system, ALDH3B1 and ALDH16A1, could induce cell cycle arrest at the G2/M phase and epithelial-mesenchymal transition by performing a CCK8 assay, colony-forming assay, transwell assay and flow cytometry analysis. Finally, a prognostic risk score model was established to assist in clinical decision-making.

Materials And Methods

mRNA Data Acquisition and Processing

mRNA data were downloaded from public databases, including The Cancer Genome Atlas (TCGA), Chinese Glioma Genome Atlas (CGGA) and Gene Expression Omnibus (GEO). Samples were separated as training (TCGA dataset) and validation cohorts (CGGA dataset and GSE108474 dataset); CGGA sequencing data and microarray data were termed CGGA1 and CGGA2 datasets, respectively. Detailed information is listed in .

mRNA data for Single-cell RNA-Seq analysis was downloaded from the GSE139448 dataset (20). GSE139448 dataset consisted of three GBM samples, GBM27, GBM28 and GBM29. R packages “Seurat”, “NormalizeData” and “FindVariableGenes” were applied to normalize the expression matrix and identify markers for each cell type. Cluster analysis was performed using PCA and visualized with the UMAP algorithm. Six cell types were clustered, including neoplastic (tumor cells), oligodendrocyte progenitor cells (OPC), macrophage, vascular, oligodendrocyte and NK cells. Among these cell types, OPC, vascular and oligodendrocyte were defined as tumor stromal cells.

Single Nucleotide Polymorphisms and Copy Number Variations

The single nucleotide polymorphisms (SNPs) and copy number variations (CNVs) were generated by performing GISTIC 2.0 (21). The nomogram and corresponding ROC curve were constructed.

The Risk System and Biofunction Prediction

A risk score of each sample was calculated by LASSO regression analysis as described in our previous work (22). High and low-risk groups were separated based on the median value of the risk score. Genes enriched in three GO pathways: GO_REGULATION_OF_EXTRACELLULAR_MATRIX_DISASSEMBLY, GO_REGULATION_OF_CELL_CYCLE_CHECKPOINT and GO_T_CELL_MEDIATED_IMMUNITY were downloaded from http://www.gsea-msigdb.org/gsea/msigdb. The correlation between contributors of the scoring system and genes involved in those three pathways was calculated using the Pearson correlation coefficient.

Then, we examined the activation of those pathways in GBM samples with single-cell RNA-seq analysis. GO enrichment analysis and the AUC score of three GO pathways, including T cell-mediated immunity, regulation of cell cycle arrest and extracellular matrix disassembly, were performed by using the “clusterProfiler” and “AUCell” package, respectively (23, 24). Only tumor cells were analyzed for this part.

Tumor Immune Landscape

The immune landscape of glioma between two groups was evaluated by the ESTIMATE algorithm (25), the CIBERSORT algorithm (26) and 28 subpopulations of immunocytes (27) with R packages “ESTIMATE”, “CIBERSORT” and “ssGSEA”. The expression profile of tumor immune escape-related genes was also illustrated (28, 29).

Cell Culture and Transfection

Human glioma cells, U251 and U87MG, were purchased from the Chinese Academy of Sciences. Cells were cultured in a DMEM medium with 10% fetal bovine serum and 1% penicillin-streptomycin. Cells were subdivided into the control group, the siRNA negative control group (the si-NC group), the transfection group (ALDH3B1: siRNA-462, siRNA-1155; ALDH16A1: siRNA-1192, siRNA-536). The siRNA was diluted to 100nM for transfection, using Lipofectamine® 2000 Reagent according to the manufacturer’s protocol.

Real-Time Quantitative PCR

RNA was extracted using Trizol reagent. Agarose gel electrophoresis was used to demonstrate the high integrity of extracted RNA. The cDNA was then synthesized and added to a PCR mix along with forward and reverse primers of each gene. The PCR mix consisted of 2 ul cDNA, 10 ul SYBR Green, 0.8 ul of the former primer, 0.8 ul latter primer and 6.4 ul ddH2O, as the description of TSINGKE Master qPCR Mix (Tsingke Biotechnology Co., Ltd., Beijing). The primers were designed by primer5 and synthesized by Sangon Biotech (Shanghai, China). The cycling conditions were set as follows: 95°C for the first 10 minutes, followed by 40 cycles at 95°C for 15 seconds and, finally, 60°C for 30 seconds. The primer sequences used are listed below:

β-actin: Forward: ACCCTGAAGTACCCCATCGAG; Reverse: AGCACAGCCTGGATAGCAAC

ALDH3B1: Forward: GTTTCCTGCCTGACAACCCTT; Reverse: AGCAGCCTCTCATTGTCCA

ALDH16A1: Forward: CCTGAGCGCCACCTACTGCAT; Reverse: CTGCCAGAAGTTGCTACACCCTT

Western-Blot Assay

DTT and PMSF were added to RIPA before lysing cells and the supernatant was collected after cell lysis with RIPA (Shanghai, China). The protein concentration was determined by the bicinchoninic acid method. Samples were run on 10% SDS-PAGE gel at 75mV. The bolts were transformed to membrane at 300mA. The membrane was incubated with 5% non-fat milk for 2 hours, and then primary antibodies were added. Primary antibodies including ALDH3B1 (1:1000, Proteintech, Cat# 19446-1-AP, RRID: AB_2861362), ALDH16A1 (1:2000, Abcam, Cat# ab137082, RRID: AB_2861363), E-cadherin (1:5000, Proteintech Cat# 20874-1-AP, RRID : AB_10697811), N-cadherin (1:1000, Abcam Cat# ab76011, RRID : AB_1310479), snai1 (1:1000, Proteintech Cat# 13099-1-AP, RRID : AB_2191756), β-catenin (1:5000, Proteintech Cat# 51067-2-AP, RRID : AB_2086128), cyclin A (1:2000, Proteintech Cat# 18202-1-AP, RRID : AB_10597084), cyclin B1 (1:2000, Proteintech Cat# 55004-1-AP, RRID : AB_10859790), aurora A (1:500, Proteintech Cat# 10297-1-AP, RRID : AB_2061337) and β-actin (1:5000, Proteintech Cat# 66009-1-Ig, RRID: AB_2687938) were diluted by antibody dilution buffer.

On the second day, the membrane was incubated with horseradish peroxidase-conjugated secondary antibodies (1:2000, Proteintech Cat# SA00001-1, RRID: AB_2722565; 1:2000, Proteintech Cat# SA00001-2, RRID: AB_2722564) for 1 hour at room temperature. Bands were visualized using enhanced chemiluminescence (ECL) solution, and pictures were captured by AI600 Chemiluminescent Imager (Cytiva, USA).

CCK8 Assay

Cells were plated into 96-well plates at a concentration of 1×103 cells per well. Then, cells were cultured at 37°C with 5% CO2. The OD value of 450nm was measured at 0h, 12h, 24h, 48h and 72h after adding culture medium containing 10% cck8 (Dojindo, Japan).

Transwell Assay

100 μl single-cell suspensions (1x106 cells per milliliter) with empty DMEM were added to the upper chamber (Corning), while DMEM with 20% bovine serum was added to the lower chamber. After incubating cells at 37°C with 5% CO2 for 48 hours, the chamber was washed with PBS and cells on the upper chamber side surface of the membrane were removed with a cotton swab. Cells attaching to the lower surface were fixed and stained. The OD value at 550nm was measured after cell discoloration with acetic acid.

Flow Cytometry Assay

1×106 cells per milliliter single-cell suspensions were collected and fixed with pre-cooled 75% ethanol. Each sample was stained with Propidium Iodide Solution consisting of 500ul PI stain buffer, 25ul PI stain and 10ul RNase A Solution. After incubating cells at 37°C in the dark for 30 minutes, samples were uploaded on a flow cytometer. The percentage of cells in the G0, G1 and S phases were analyzed using ModFit LT (version 5.0).

Colony-Forming Assay

Cells were digested with 0.25% trypsin, then were plated into a 6-well plate at a concentration of 200 cells per well. Afterward, cells were incubated in cell incubators for 14 days. Then, 4% methanol was applied to fix the cells colony, and crystal violet was used to stain those colonies. The OD value at 550nm was measured after cell discoloration.

Statistical Analysis

Differences between the two groups were analyzed by the Mann-Whitney U test, while one-way ANOVA analysis was conducted for multiple groups. A two-way ANOVA analysis was conducted to ananlyze the CCK8 assay. Kaplan-Meier survival curves were generated and compared by using the log-rank test. Statistical analysis for the colony-forming assay, the CCK8 assay and the transwell assay were carried out by GraphPad Prism (version 8.0). A P-value < 0.05 was considered to be statistically significant. All bioinformatics analyses were conducted by R (version 3.6.2).

Results

Abnormal ALDH Expression Was Associated With Gliomas Pathological Grades

TCGA data (522 lower-grade glioma (LGG) samples and 150 GBM samples) was used as the training cohort, while the CGGA1 (180 LGG samples and 136 GBM samples), CGGA2 (173 LGG samples and 124 GBM samples) and GSE108474 (170 LGG samples and 124 GBM samples) datasets were set as validation cohorts.

The correlation between ALDH expression and clinicopathological parameters of gliomas, including pathological glioma grade and IDH status, was first investigated. Analysis of the TCGA dataset showed higher expression levels of ALDH16A1, ALDH3B1, ALDH7A1, ALDH1A2, ALDH3A1 and ALDH1A3 in GBM compared to LGG, while ALDH4A1, ALDH8A1, ALDH1A1, ALDH6A1, ALDH2, ALDH5A1, ALDH1L2 expressions were lower in GBM ( ). In the validation cohorts, lower expressions of ALDH1A1, ALDH6A1, ALDH2, ALDH5A1 were observed in GBM compared to LGG, while ALDH16A1, ALDH3B1 and ALDH3A1 were preferentially enriched in GBM ( ).

Figure 1

Expression profile of the ALDH family. The ALDH family expression of normal brain tissue and glioma. The expression profile of ALDHs between WHO grades of glioma from the TCGA (A), CCGA1 (B), CGGA2 (C) and GEO (D) database. (E) Genomic alteration of ALDHs in glioma. OPC: oligodendrocyte precursor cell. (F–H) The construction of the scoring model based on the LASSO algorithm. NS, no significantly statistical; *P < 0.05; **P < 0.01; ***P < 0.001.

In the LGG cohort from the TCGA dataset, ALDHs including ALDH2, ALDH6A1, ALDH5A1, ALDH9A1, ALDH1A1 were upregulated in WHO grade II gliomas while the increased expression of ALDH16A1, ALDH3B1, ALDH3A2, ALDH1A2, ALDH3A1 was found in WHO III grade gliomas ( ). Similar alterations of ALDH2, ALDH6A1, ALDH5A1, ALDH1A1, ALDH16A1, ALDH3B1, ALDH3A2, ALDH3A1 expression were observed in CGGA and GEO108474 validation datasets ( ).

The IDH status was closely associated with glioma progression, and IDH mutant gliomas had better clinical outcomes. Analysis of the TCGA dataset showed that ALDH16A1, ALDH3B1, ALDH1A3, ALDH3A1, ALDH7A1 were significantly increased in IDH wildtype gliomas, while ALDH2, ALDH6A1, ALDH5A1, ALDH1A1, ALDH1L2, ALDH18A1, ALDH1B1 expression were higher in IDH mutant gliomas than IDH wildtype gliomas ( ). Similar gene expression dysregulation between IDH wildtype and IDH mutant gliomas were verified in validation datasets ( ).

Moreover, the incidence of somatic ALDH mutations in gliomas was 5.33% (48/901), and most of them were missense mutations. This result further substantiated the presence of ALDH dysregulations in gliomas ( ). Taken together, the expression of most ALDHs was altered, and this dysregulation was associated with glioma progression.

The Risk Scoring System Can Predict Glioma Patient Prognosis

To explore the prognostic value of different ALDHs in glioma patients, univariate Cox regression was applied to assess the relationship between individual ALDHs and the overall survival outcome of glioma patients from the TCGA dataset. A total of 16 ALDHs was significantly associated with overall survival (OS) with p-values less than 0.05 ( ). ALDH2, ALDH5A1, ALDH6A1, ALDH1L2, ALDH1B1, ALDH18A1, ALDH1A1, ALDH8A1, ALDH1B1 were found to be favorable factor for glioma patients with a hazard ratio (HR) <1. ALDH16A1, ALDH3B1, ALDH3A1, ALDH7A1, ALDH1A3, ALDH1L1, ALDH1A2 were risk genes with HR > 1. Moreover,overall survival analysis based on ALDHs further confirmed their prognostic roles in gliomas ( ).

Subsequently, a risk model based on the LASSO regression analysis was constructed, which consisted of five genes ( ). The risk score for each sample was calculated based on these genes. The performance of our risk model was also assessed using data from the validation cohort. Heatmaps were built to illustrate the correlation between the risk score and ALDHs expression ( ). We found that the expression of ALDH2 and ALDH5A1 were inversely correlated with the risk score, while ALDH3A1, ALDH 3B1 and ALDH16A1 expression positively correlated with the risk score.

Figure 2

The risk model based on the ALDH family. Association between ALDHs expression and risk in TCGA (A), CCGA1 (B) and CGGA2 (C) database. The distribution of risk in gliomas grades (D, E), the IDH status gliomas (F), tumor subtype (G) and tumor response to treatment (H). Overall survival analysis based on risk in TCGA (I), CCGA1 (J) CGGA2 (K) and GEO (L) database. G2, grade II; G3, grade III; CL, classical; ME, mesenchymal; NE, neural; PN, proneural; CR, complete remission/response; PR, partial remission/response; PD, progression diseases; SD, stable diseases. ***P < 0.001.

The risk score distribution for different clinical features was also analyzed. High-risk samples were associated with malignant glioma features, including high-grade ( ), IDH wildtype ( ), and classical or mesenchymal ( ). Besides, high-risk samples also exhibited worse treatment responses than low-risk samples ( ). A similar tendency could also be noticed in validation datasets ( ). Taken together, high-risk samples were associated with worse overall survival outcomes than low-risk samples.

Subsequently, the prognostic value of different risk scores obtained from glioma patients was evaluated by performing survival analysis. High risk patients were strongly associated with poor clinical outcome in TCGA (P < 0.0001, ), CGGA1 (P < 0.0001, ), CGGA2 (P < 0.0001, ) and GSE108474 datasets (P < 0.0001, ). Moreover, our risk model exhibited a good performance in predicting survival outcomes of LCG patients but performed poorly for GBM patients ( ). Nevertheless, a significant difference in the survival outcome of GBM patients in the CGGA1 (P = 0.025) and CGGA2 (P = 0.027) cohorts was found; a similar tendency was noticed in the GSE108474 dataset; however, the difference was not statistically significant. Therefore, the poor predictive performance of our risk model could result from the limited number of GBM samples.

Somatic Mutations and Copy Number Variations Between High and Low-Risk Groups

Next, we explored the difference in somatic mutations and copy number variations between high and low-risk groups. The overall mutation ratio was 89.9% and 97.85% in groups with high and low-risk scores. A higher frequency of IDH1, TP53, ATRX mutation was found in the low-risk group than in the high-risk one (IDH1 91% vs. 30%, TP53 47% vs. 39%, ATRX 33% vs. 21%). In contrast, a lower mutation frequency of TTN, MUC16, and PIK3CA was found in the high-risk group than in the low-risk group (TTN 8% vs. 21%, MUC16 6% vs. 11%, PIK3CA 7% vs. 8%). In addition, mutations of CIC (27%),FUBP1 (12%), and NOTCH1 (10%) were observed in the lowriskgroup, while EGFR (21%), PTEN (18%), and NF1 (12%) mutations were identified in the high-risk score group ( ).

Figure 3

The SNPs and CNVs between high and low-risk groups. SNPs in high (A) and low-risk (B) groups. (C) The overview of the incidence of CNVs ranked by risk. Amplification (D) or deletion (E) regions in samples from high or low-risk group.

CNVs analysis revealed that glioma patients with high-risk scores featured frequent amplification of Chr 7 and deletion of Chr 9 and Chr 10 ( ). GISTIC 2.0 analysis revealed that in gliomas with high-risk scores, the amplification mainly occurred in chromosomal regions 1q32.1, 3q26.33, 4q12, 7p11.2, and 12q14.1, where several key oncogenes (PDGFRA, EGFR, MDM2, SOX2, and CDK4) were located. Focal deletions in chromosomal regions 9p21.3 and 10q23.31, which encode tumor suppressors including PTEN, FAS, and CDKN2A, were observed ( ).

The Clinical Application of the Risk Scoring System

Independent prognostic indicators were identified by employing univariate and multivariate Cox analysis ( ). Schoenfeld individual tests were performed to identify biomarkers for the construction of nomogram ( ). A nomogram was built to predict the survival of patients using those variables, and each variable was assigned a score ( ). Calibration curves proved that this model could precisely predict a patient’s survival outcome ( ). The overall survival analysis suggested shorter survival of patients with high nomogram scores than patients with low nomogram scores. Similar results were obtained in the CGGA datasets ( ). The area under the curve (AUC) of receiver operating characteristic (ROC) curves predicting the overall survival probability at 3 and 5 years was 0.882 and 0.816, respectively, in the TCGA dataset. Validation results from the CGGA1 (AUC at 3 years OS: 0.876; AUC at 5 years OS: 0.890) and CGGA2 (AUC at 3 years OS: 0.881; AUC at 5 years OS: 0.812) also substantiated this conclusion ( ).

ALDHs May Affect Tumor Immune Landscape, Tumor Proliferation and Migration

GSVA analysis was applied to perform GO and KEGG enrichment analysis ( ). Pathways related to cell cycle, cell adhesion and tumor immune landscape were enriched in the high-risk group, including regulation of transcription involved in G1/S transition of the mitotic cell cycle, cell cycle, extracellular matrix disassembly, regulation of extracellular matrix organization, T-cell mediated immunity, regulation of T cell chemotaxis, regulatory T cell differentiation. The GSEA findings were consistent with the GO/KEGG enrichment analysis results ( ).

Figure 4

Biofunction prediction based on bulk RNA-Seq analysis and single-cell analysis. (A) The GO and KEGG enrichment analysis based on the GSVA algorithm. (B) The GO enrichment analysis based on the GSEA algorithm. (C) The cell components of single cell RNA-seq analysis. (D) The distribution of ALDH2, ALDH16A1, ALDH3B1, ALDH5A1, ALDH3A1 and risk in different cell components. The AUC score of three pathways, including T cell mediated immunity (E), positive regulation of cell cycle arrest (F) and extracellular matrix disassembly (G), from the GO enrichment analysis were calculated. (H) The GO enrichment analysis based on single-cell analysis.

Moreover, we examined the relevance of these pathways with pure tumor cells using single-cell RNA-Seq analysis. Six cell components were identified, including tumor cells, OPC, macrophage, vascular, oligodendrocyte and NK cells ( ). The distribution of the components of our risk score model and risk score associated with these cell components was mapped; and ALDHs were mostly enriched in neoplastic cells and macrophage ( ). The activation of three GO pathways, T cell-mediated immunity ( ), positive regulation of cell cycle arrest ( ) and extracellular matrix disassembly ( ), were evaluated by calculating the AUC score. Importantly, results from the GO enrichment analysis of single-cell RNA-seq based on R package “clusterProfiler” ( ) and “GSVA” both proved those results ( ).

Figure 5

Biofunction prediction based on bulk RNA-Seq and single-cell RNA-Seq analysis. (A) The enrichment score distribution of GO pathways based on single-cell analysis. Association between contributors of the scoring system and genes enriched in T cell mediated immunity (B), cell cycle checkpoint (C) and extracellular matrix disassembly (D) in TCGA, CGGA1, CGGA2 and GSE108474 datasets. NS, no significantly statistical; ***P < 0.001.

To improve our understanding of the relationship between components of our risk model and the tumor immune landscape, cell cycle and migration, the correlation between different ALDHs and genes from GO pathways, including T cell-mediated immunity ( ), cell cycle checkpoint ( ) and extracellular matrix disassembly pathways ( ) was analyzed. ALDH16A1, ALDH3A1, ALDH3B1 positively correlated with genes in these pathways, while the opposite finding was observed for ALDH2 and ALDH5A1. Thus, these results suggested that ALDHs could affect the tumor immune landscape, tumor cell cycle and migration.

ALDHs Can Affect Immunocyte Infiltration and the Expression of Immune Checkpoints

First, tumor purity was calculated by applying the ESTIMATE algorithm to evaluate the immunocytes infiltration ratio in high and low-risk groups. The increased estimate score (a combination of stromal score and immune score) and decreased tumor purity in high-risk score samples from TCGA suggested that high-risk score samples possessed a complicated tumor immune landscape with more immunocytes and stromal cells ( ). Results from the validation datasets also drew similar conclusions ( ). The composition of immunocytes in the high and low-risk groups was generated by applying the CIBERSORT algorithm and 28-immune cell lineage analysis in the training ( ) and validation datasets ( ). Besides, the correlation between the enrichment score of immunocytes and risk score was also calculated ( and ). Results showed that the proportion of immunocytes varied significantly between high and low-risk groups. The high-risk group contained more M2 macrophages, regulatory T cells and CD4 memory resting T cells, whereas activated mast cells were relatively lower.

Figure 6

Tumor immune landscape and the expression of immune checkpoint genes based on the scoring system. (A) Tumor purity and estimate score were calculated by the ESTIMATE algorithm. (B) The proportion of immunocytes based on the CIBERSORT algorithm. (C) The enrichment score profile of 28 immunocytes. (D) Correlation between 28 immunocytes enrichment score and risk. (E–K) The expression profile of immune checkpoint genes in the risk model. NS, no significantly statistical; *P < 0.05; **P < 0.01; ***P < 0.001.

Interestingly, significantly infiltration of immunotherapy-related immunocytes was found in the high-risk group, including activated CD8 T cells, regulatory T cells, helper T cells et al. Accordingly, we explored the expression of tumor immune checkpoint genes (ICGs) for a better understanding of the role of ALDHs in immunotherapy ( ). As illustrated, higher expression of ICGs, such as CD40LG, CD40, CTLA4, PDCD1, PDCD1LG2, ICAM1, CD80 and GZMA, was found in the high-risk group, compared to the low-risk group. In summary, ALDHs could participate in reshaping the tumor immune landscape and weakening immunocyte function.

ALDH3B1 and ALDH16A1 Promote Proliferation and Migration of Glioma Cells

Three out of the five components of our risk score, ALDH3A1, ALDH3B1 and ALDH16A1, were identified as tumor promoters. A previous study reported that ALDH3A1 could modulate glioma proliferation and migration (30), whereas the other two genes have been poorly investigated. Therefore, we silenced the expression of ALDH3B1 or ALDH16A1 in two GBM cell lines, U251 and U87MG, to verify their role. RT-qPCR and Western-blotting assay were performed to detect alterations of ALDH3B1 and ALDH16A1 expression. We observed that both mRNA and protein levels of ALDH3B1 and ALDH16A1 were significantly decreased in the silenced group ( ). Furthermore, the CCK8 assays showed that the proliferation ratio of cells with silenced ALDH3B1 or ALDH16A1 was significantly reduced ( ). The colony-forming ability of tumor cells was inhibited by decreasing ALDH3B1 and ALDH16A1 expression ( , statistical results in ).

Figure 7

Silenced ALDH3B1 or ALDH16A1 expression inhibits glioma cells progression. The expression of ALDH3B1 was silenced by siRNA-462, siRNA-1155, while ALDH16A1 was silenced by siRNA-536, siRNA-1192. The effect of ALDH depletion by siRNA transfection was detected by qRT-PCR (A) and Western blot (B). (C) Cell proliferation was measured at the indicated time points by CCK8 assay. (D) Representative images of colony formation. All experiments have been independently replicated three times. NS, no significantly statistical; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001. n=3. Data are represented as mean ± SD.

In addition, flow cytometry results suggested that cells were accumulated at the S phase and G2 phase, implying cell cycle arrest at the G2/M phase (U251 in , U87MG in ). A previous study reported that cyclin A, cyclin B1 and aurora A were G2/M cell cycle checkpoints (31). In this work, cyclin B1 and aurora A expressions were decreased after silenced ALDH3B1 and ALDH16A1 expression; however, no difference in cyclin A expression was observed (U251 in , U87MG in ). Based on these findings, it can be concluded that ALDH3B1 and ALDH16A1 may affect the cell cycle via modulation of cyclin B1 and aurora A.

Figure 8

ALDH3B1 and ALDH16A1 affect glioma cells proliferation and migration in U251. The expression of ALDH3B1 was silenced by siRNA-462, siRNA-1155, while ALDH16A1 was silenced by siRNA-536, siRNA-1192. (A) Cell cycle was arrested at the G2/M phase by decreasing ALDH3B1 and ALDH16A1 expression. (B) Representative images of migration assays for U251 transfected with indicated siRNA. (C) The expression of tumor cell cycle and migration associated proteins in U251, including E-cadherin, N-cadherin, β-catenin, snai1, aurora A, cyclin A, cyclin B1 and β-actin. Scale bar = 560μm. All experiments have been independently replicated three times. NS, no significantly statistical; *P < 0.05; **P < 0.01; ***P < 0.001. n=3. Data are represented as mean ± SD.

In contrast, the transwell assay showed that tumor cell migration was inhibited by suppressing ALDH3B1 or ALDH16A1 expression (U251 in , U87MG in , statistical results in ). E-cadherin expression was elevated while N-cadherin, β-catenin and snai1 expression were decreased, implying that epithelial-mesenchymal transition was inhibited after silencing ALDH3B1 and ALDH16A1 expression (U251 in , U87MG in ). Therefore, ALDH3B1 and ALDH16A1 could promote glioma progression by regulating cell proliferation and migration.

Discussion

Gliomas are the most common primary malignant brain tumors. GBM has been documented with the worst outcomes, with an average survival of less than two years (32). Therefore, it is necessary to explore the underlying biological mechanisms of gliomas to develop new classification systems and design novel treatment approaches (33). The role of ALDHs in gliomas remains elusive; exploring its role may be beneficial in constraining glioma progression.

In this study, we used LASSO regression analysis to construct a risk model of gliomas based on ALDHs. ALDH5A1, ALDH3B1, ALDH3A1, ALDH2 and ALDH16A1 were screened as main contributors of that scoring system to calculate a risk score. Samples with high risk were associated with a worse prognosis. Bioinformatic analysis and in vitro experiments suggested that ALDHs influenced the tumor cell cycle, epithelial-mesenchymal transition, and infiltration of immunocytes. Moreover, our constructed prognostic model, nomogram, based on risk and gliomas type, could precisely predict overall survival outcome, indicating its potential clinical application value.

The dysregulated expression of ALDHs has been widely noticed in gliomas and is reported tightly associated with glioma progression. Recently, several studies reported the function of ALDHs in gliomas. An increasing body of evidence suggests that ALDHs target inhibitors could inhibit stem-like cell proliferation of gliomas (18, 34, 35). Furthermore, ALDH1A3 could reportedly affect glioma glucose metabolism and is recognized as a stem-like cells biomarker (36, 37). Besides, ALDH1A1 and ALDH1A2 have been documented as glioma migration-associated genes (38–40), while ALDH1A1 also participated in GBM resistance to temozolomide (17, 19). To the best of our knowledge, this is the first study that identified the role of ALDH16A1 and ALDH3B1 in gliomas.

Based on our scoring system, three out of five ALDHs (ALDH3A1, ALDH3B1 and ALDH16A1) were identified as unfavorable biomarkers for gliomas. Previous studies reported that ALDH3A1 was involved in the Wnt/β-catenin signaling pathway mediated glioblastoma resistance to temozolomide (30). Accordingly, we focused on the role of ALDH3B1 and ALDH16A1 in gliomas.

In the present study, ALDH3B1 and ALDH16A1 were positively correlated with genes involved in cell cycle checkpoint-related pathways and extracellular matrix disassembly-related pathways. Previous studies reported the epithelial-mesenchymal transition as one of the glioma’s malignant factors, suggesting its critical role in glioma progression (41, 42). Importantly, silencing ALDH3B1 and ALDH16A1 expression induced cycle arrest at the G2/M phase and inhibited the epithelial-mesenchymal transition in the glioma cells. Critically, ALDH3B1 and ALDH16A1 could also affect the expression of cyclin B1 and aurora A, G2/M cell cycle checkpoints that modulate glioma cells proliferation. Moreover, the expression of epithelial-mesenchymal transition-associated proteins was affected after interfering with ALDH3B1 and ALDH16A1 expression, leading to increased E-cadherin expression and decreased expression of N-cadherin, β-catenin and snai1.

Interestingly, based on single-cell RNA-seq analysis, we found that the system’s main contributors were mostly expressed in tumor cells, OPC and macrophages. As previously mentioned, the influence of ALDHs on tumor cells and their role in glioma stem-like cells have been established by various studies. Apart from that, distribution map from single cell RNA-seq analysis also suggested that ALDHs may affect the tumor immune landscape by modulating immunocytes. Previous studies reported that elevated ALDHs activity could modulate macrophage polarization (43), and tumor-infiltrating macrophages were positively correlated with glioma cell proliferation and invasion (44). In addition, ALDHs also played a vital role in the induction and function of regulatory T cells (45). Taken together, those findings suggested that ALDHs could affect the tumor immune landscape by regulating the communication between tumor cells and immunocytes and tumor response to the immune system.

Several ALDHs have also been associated with glioma sensitivity to temozolomide (30, 46). In our study, higher expression of genes associated with antigen representation, immunocyte inhibitor or stimulator, ligand-receptor were found in the high-risk group compared to the low-risk group. This finding suggested that immunocytes may also be affected the ALDHs. For instance, GZMA has been reported to trigger adaptive immunity by stimulating dendritic cells (47) and modulates T cell expansion in immunotherapy (48). In the meantime, the elevated expression of ICGs such as PDCD1, CD274, CTLA4 may also modulate tumor response to immunotherapy (49). In a nutshell, ALDHs are vastly involved in determining tumor response to treatment, and targeting them could have huge prospects in improving current treatment efficiency.

In conclusion, this study systematically illustrated the role of ALDHs in the tumorigenesis of gliomas. The expression pattern, prognostic value and potential mechanism of ALDHs were analyzed. A clinical model was constructed to predict clinical outcomes. Moreover, our in vitro experiments substantiated that ALDH3B1 and ALDH16A1 could affect tumor cells proliferation and epithelial-mesenchymal transition.

Data Availability Statement

The original contributions presented in the study are included in the article/ . Further inquiries can be directed to the corresponding authors.

Author Contributions

Writing – original draft: ZYW, YM, and YT. Writing-review & editing: ZYW, YM, ZHW, and ZD. Data curation: ZYW. Formal analysis: ZYW. Validation: YM, YT, and HZ. Visualization: XZ and SF. Methodology: XL, and QC. Project administration: TS and QC. Funding acquisition: QC. Supervision: TS and QC. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Nature Science Foundation of China (NO.82073893, 81703622); the China Postdoctoral Science Foundation (NO.2018M633002); Hunan Provincial Natural Science Foundation of China (NO.2018JJ3838, 2020JJ8111); and grants from the research projects of Hunan Provincial Health Commission of China (NO.C2019186); Xiangya Hospital Central South University postdoctoral foundation; and the Fundamental Research Funds for the Central Universities of Central South University (No. 2021zzts1027).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.