- 1Department of Neurosurgery, Affiliated Wuxi People’s Hospital of Nanjing Medical University, Wuxi, China
- 2Department of Gastroenterology, Affiliated Hospital of Jiangsu University, Zhenjiang, China
- 3Department of Lung Cancer Surgery, Tianjin Medical University Cancer Institute and Hospital, Tianjin, China
- 4Department of General Surgery, Affiliated Wuxi People’s Hospital of Nanjing Medical University, Wuxi, China
- 5Department of Plastic Surgery, Xiangya Hospital, Central South University, Changsha, China
- 6Suzhou Kowloon Hospital, Shanghai Jiao Tong University School of Medicine, Suzhou, China
Background: Glioblastoma (GBM) is a malignant primary brain tumor. This study focused on exploring the exosome-related features of glioblastoma to better understand its cellular composition and molecular characteristics.
Methods: Single-cell RNA sequencing (scRNA-seq) and spatial transcriptome RNA sequencing (stRNA-seq) were used to analyze the heterogeneity of glioblastomas. After data integration, cell clustering, and annotation, five algorithms were used to calculate scores for exosome-related genes(ERGs). Cell trajectory analysis and intercellular communication analysis were performed to explore exosome-related communication patterns. Spatial transcriptome sequencing data were analyzed to validate the findings. To further utilize exosome-related features to aid in clinical decision-making, a prognostic model was constructed using GBM’s bulk RNA-seq.
Results: Different cell subpopulations were observed in GBM, with Monocytes/macrophages and malignant cells in tumor samples showing higher exosome-related scores. After identifying differentially expressed ERGs in malignant cells, pseudotime analysis revealed the cellular status of malignant cells during development. Intercellular communication analysis highlighted signaling pathways and ligand-receptor interactions. Spatial transcriptome sequencing confirmed the high expression of exosome-related gene features in the tumor core region. A prognostic model based on six ERGs was shown to be predictive of overall survival and immunotherapy outcome in GBM patients. Finally, based on the results of scRNA-seq and prognostic modeling as well as a series of cell function experiments, BARD1 was identified as a novel target for the treatment of GBM.
Conclusion: This study provides a comprehensive understanding of the exosome-related features of GBM in both scRNA-seq and stRNA-seq, with malignant cells with higher exosome-related scores exhibiting stronger communication with Monocytes/macrophages. In terms of spatial data, highly scored malignant cells were also concentrated in the tumor core region. In bulk RNA-seq, patients with a high exosome-related index exhibited an immunosuppressive microenvironment, which was accompanied by a worse prognosis as well as immunotherapy outcomes. Prognostic models constructed using ERGs are expected to be independent prognostic indicators for GBM patients, with potential implications for personalized treatment strategies for GBM. Knockdown of BARD1 in GBM cell lines reduces the invasive and value-added capacity of tumor cells, and thus BARD1-positively expressing malignant cells are a risk factor for GBM patients.
1 Introduction
Glioblastoma (GBM) is a highly invasive primary brain tumor that originates from glial cells and is considered one of the most common and aggressive brain tumors in adults. It typically occurs in the central regions of the brain, including the cerebral cortex, basal ganglia, and white matter areas (1). The pathogenesis of GBM is complex and involves multiple factors such as genetic variations, epigenetic changes, and environmental factors. Some common biomarkers include the expression levels of glioblastoma surface molecules (such as EGFR, IDH1, MGMT), 1p/19q chromosomal deletions, and overexpression of the MDR1 gene (2, 3). The detection of these biomarkers provides valuable information about the molecular subtypes, prognosis, and treatment response of the tumor. The treatment of GBM is a comprehensive approach that includes surgical resection, radiation therapy, and chemotherapy. However, due to the invasive and heterogeneous nature of GBM, the treatment outcomes are often limited (4, 5). In recent years, the introduction of personalized therapy and immunotherapy has brought hope for the treatment of GBM, but many challenges remain. With the increasing importance of molecular alterations in the classification and grading of gliomas, the search for new biomarkers and the establishment of effective molecular subtyping systems are crucial in helping clinicians select the most suitable treatment strategies for GBM patients (6, 7).
Exosomes are a class of small vesicles secreted by cells that are approximately 30 to 150 nm in diameter (8). Exosomes contain biomolecules such as proteins, nucleic acids (e.g. DNA and RNA), lipids, and cell membrane proteins whose composition reflects the origin and state of their parent cells (9). Exosomes have a dual role in tumor immunotherapy. They can either activate the immune response and enhance anti-tumor immunity or they may inhibit the activity of immune cells and provide an escape mechanism for the tumor (10, 11). The complexity and diversity of exosomes allow them to possess different immunomodulatory capacities, and their role depends on the source, environment, and immune status. Furthermore, the presence of exosomes in the tumor microenvironment may inhibit anti-tumor immune responses. Therefore, an in-depth study of the properties and regulatory mechanisms of exosomes is essential to optimize their use in tumor immunotherapy (12, 13).
PD-L1, known for its expression on tumor cells’ surfaces, has been identified as a promoter of tumor immune evasion. Surprisingly, apart from its cell surface expression, PD-L1 is released from tumor cells into the extracellular space in the form of free PD-L1, including exosomal PD-L1 (14). Emerging research has highlighted the crucial involvement of exosomal PD-L1 in tumor-induced immunosuppression. Notably, exosomal PD-L1 exhibits higher resistance to degradation by extracellular protein hydrolases compared to soluble PD-L1. It also contributes to T cell dysfunction and enhances stability, emphasizing its significant role in tumor immune modulation (15). In a recent glioma study, GBM cells were found to excrete oncogenic miRNAs extracellularly via exosomes, promoting the conversion of immune cells’ immunosuppressive phenotypes in the tumor microenvironment, thus achieving the dual role of promoting malignant tumor progression (16).
However, despite significant advances in exosome research using high-resolution technologies such as single-cell sequencing, our detailed understanding of exosomes in specific tumor types like glioblastoma remains limited at present (17). These specific tumors often have complex cellular compositions and microenvironmental features, where the role of exosomes may involve more diverse and complex mechanisms. By studying the characteristics of exosome-related genes, we can unveil their crucial role in regulating immune cell functions, modulating the tumor microenvironment, and influencing tumor immune escape. In this study, we comprehensively characterized the microenvironment of exosome-related genes (ERGs) in GBM through single-cell sequencing and spatial transcriptome sequencing. We combined this data with bulk sequencing to construct a validated prognostic index to aid clinicians in making better treatment decisions (18). Finally, we identified a new immunotherapeutic target, BARD1, and validated our findings through various bioinformatics and cellular assays. These additional studies will further unravel the regulatory mechanisms of exosomes and provide a theoretical basis for developing targeted therapeutic and immunotherapeutic strategies against exosomes.
2 Materials and methods
2.1 Source of raw data
Bulk RNA-seq data, mutation data, and clinicopathological features for TCGA-GBM were downloaded from the UCSC Xena website (https://xena.ucsc.edu/). There are 168 glioma samples in the TCGA cohort. Gene expression profiling data for 374 GBM patients in the validation model were obtained from the China Glioma Genome Atlas (CGGA) data portal (http://www.cgga.org.cn/). All expression profile data were in TPM format. Batch correction and integration of the two sets of gene expression data were performed using the “limma” and “sva” (19) software packages. Single-cell RNA sequencing(scRNA-seq) data for GBM were downloaded from GSE84465 (20) and contained a total of 3589 cells within and near the tumor. Spatial transcriptome RNA sequencing(stRNA-seq) for primary GBM from the 10x Visium platform was downloaded from GSE194329 (21). 121 exosome-related genes were downloaded from the ExoBCD database (https://exobcd.liumwei.org) (22).
2.2 Processing of single-cell sequencing data
We conducted an analysis of single-cell RNA sequencing data using the R packages ‘Seurat’ and ‘SingleR’ (23). To ensure the inclusion of high-quality cellular data, we considered genes expressed in a minimum of three single cells. In addition, we excluded cells with gene counts fewer than 200 or exceeding 10,000, fewer than 1000 counts, and over 20% of mitochondrial and ribosomal genes. To address batch effects between cancer and paracancer samples, we employed the ‘harmony’ R package (24). The scRNA-seq data were normalized using the ‘Seurat’ R package with the ‘NormalizeData’ function. Subsequently, the normalized scRNA-seq data were transformed into Seurat objects, and the top 2000 highly variable genes were identified using the ‘FindVariableFeatures’ function. To reduce the dimensionality of the scRNA-seq data, we performed principal component analysis (PCA) using the ‘RunPCA’ function of the ‘Seurat’ R package. Significant principal components (PCs) were identified using JackStraw analysis, and appropriate PCs for cell clustering analysis were selected based on the proportion of variance. For the clustering of integrated data, we employed the ‘FindNeighbors’ and ‘FindClusters’ functions, and the resultant cells were visualized using UMAP or t-SNE methods. To identify genes specifically expressed in each cluster, we conducted Wilcoxon tests between pairs of cell clusters using the ‘FindAllMarkers’ and ‘FindMarkers’ functions of the ‘scran’ R package. The expression of specific genes was depicted using the ‘featureplot’ function (25). Cell type annotations were based on information from the original text and the tumor single-cell transcriptome database TISCH (http://tisch.comp-genomics.org/).
Cell developmental trajectories of malignant cells were analyzed by pseudotime using the “Monocle” R package (26). After transforming Seurat objects into cellular dataset objects, developmental difference genes were selected using unsupervised analysis. Utilizing the expression of 121 exosome-related genes, we used five commonly used algorithms to score gene sets from single-cell data (AddModuleScore, ssGSEA, AUCell, UCell, and singscore). The “AddModuleScore” algorithm, found within the “SingleR” R package, is used for scoring gene sets (27). Its essence lies in first calculating the mean of all genes in the gene set. Then, the expression matrix is divided into several segments based on the mean value, and a set of control genes is randomly sampled from each segment as background values. “ssGSEA” is a single-sample gene set enrichment analysis method used to evaluate the enrichment level of gene sets in a single sample or cell. It relies on the rank-based ordering of gene expression within a sample and calculates enrichment scores for each gene set. “AUCell” evaluates whether an input gene set is enriched among the top 5% expressed genes in a single sample, based on the ranking of gene expression in individual samples. The distribution of AUC scores across all cells enables the exploration of relative expression features. Given its rank-based scoring method, AUCell remains uninfluenced by gene expression units and normalization procedures. “UCell” (https://github.com/carmonalab/UCell) (28) is an unsupervised cell type identification method used to identify and classify the cell types of individual cells. UCell’s signature score is based on the Mann-Whitney U statistic, which is robust to dataset size and heterogeneity. Compared to other available methods, UCell requires less computational time and memory. “singscore” is a cell state assessment method used to quantify the activity level of specific functions or biological processes in a single sample or cell. It relies on gene sets from the gene expression profile and calculates cell state scores for samples or cells by considering gene weights and directions. These tools and methods have significant applications in single-cell transcriptomics and gene set enrichment analysis, helping researchers uncover changes in cell function, biological processes, and disease relevance, thus providing a deeper understanding of the complexity of biological systems (29).
The “CellCall” R package is a toolkit for inferring intercellular communication networks and internal regulatory signals by integrating intracellular and intercellular signals (30). The most notable feature lies in its ability to combine intercellular ligand-receptor communication with intracellular transcription factor expression, forming a ligand-receptor-transcription factor axis (L-R-TF axis). Simultaneously, it also encompasses pathway activity analysis, allowing for the analysis of receptor-cell pathway changes resulting from communication between two specific cell types.
2.3 Processing of spatial transcriptome sequencing data
The data is then analyzed using Seurat in R. The UMI counts are normalized, scaled, and the most variable features are determined by the “SCTransform” function. Dimensionality reduction is then performed using “RunPCA” for unsupervised clustering analysis. The “FindNeighbors” and “FindClusters” functions were executed with default parameters and the 30 most significant principal components. The “SpatialFeaturePlot” function was utilized for visualizing subgroups and genes.
The “scMetabolism” R package, developed by Fudan University, quantifies metabolic activity at the single-cell level. It operates based on a conventional single-cell matrix file and employs the VISION algorithm to score each cell, deriving a final activity score for each metabolic pathway (31).
For Python, the Scanpy and stlearn packages are employed. Scanpy is a Python-based package designed for analyzing single-cell data, encompassing pre-processing, visualization, clustering, proposed time series analysis, and differential expression analysis. The Institute of Molecular Biosciences at the University of Queensland has introduced an integrated analysis method, the stlearn package (https://github.com/BiomedicalMachineLearning/stLearn). This tool employs gene expression data, tissue morphology data, and spatial location information to initially identify cell types and subsequently reconstruct tissue cell types within tissues. It also infers evolutionary pathways and identifies tissue regions with high cell-to-cell interactions. stLearn integrates analyses to deduce interactions influenced by information on ligand pairs, gene expression, spatial location, and spatial cell type distribution.
Reverse Compositional Transcriptomics Deconvolution (RCTD) is a method utilized to deduce cellular composition in spatial transcriptomic data. This methodology involves analyzing the comprehensive gene expression profile of an entire tissue or sample to reverse-calculate the spatial distribution and relative abundance of each cell type. RCTD’s key advantage lies in its capacity to infer cellular composition from the overall gene expression profile, eliminating the necessity for the isolation and sequencing of individual cells. This attribute is particularly advantageous for spatial transcriptomics research, as it uncovers the relative distribution and interactions of distinct cell types within tissues.
2.4 Construction and validation of the exosome-related index
The TCGA-GBM cohort was used as the training set, while the CGGA dataset was used as the validation set. First, 121 exosome-related genes were used as candidate predictors for univariate cox analysis to identify genes that were statistically associated with patient OS (p< 0.05). Next, we performed LASSO and multivariate regression analyses to further screen for genes and risk factors that were strongly associated with prognosis (32). The exosome-related index (ERI) was calculated for each GBM patient based on the coefficients determined by multivariate cox analysis. The risk score/ERI is calculated as follows:
Here, h0(t) represents the baseline hazard function, indicating the risk level when all predictor variables are at 0, serving as the increment over the risk at time t based on a zero value of predictor variables. b1, b2, …, bp are the regression coefficients of the Cox model, and x1, x2, …, xp are the corresponding predictor variables (gene expression values) at time t, the observed values. Based on the median value of ERI, the patients in the training and validation sets were divided into high and low-scoring groups. Survival curves were also plotted using the Kaplan-Meier method and log-rank tests were used to determine their statistical significance.
2.5 The creation of Nomograms and the analysis of mutations
To calculate the probability of patient survival, we created a nomogram combining factors such as ERI, age, and IDH mutation status as prognostic factors. Using consistency index analysis and decision curve analysis (DCA), we further assessed the net benefit of the line graph and clinical characteristics alone. The ‘oncoplot’ functional waterfall plots from the R software ‘maftools’ package were used to explore detailed mutation characteristics. In addition, we investigated the correlation between ERI and tumor mutational load (TMB) and visualized this using the R software ‘ggplot2’ package (33).
2.6 Inference of the immune microenvironment and prediction of response to immunotherapy
Using the expression profile data, we employed the R package “estimate” to estimate the abundance of interstitial and immune cells, as well as tumor purity, in malignant tumor tissue (34). Subsequently, we assessed the extent of immune infiltration in GBM patients using the TIMER 2.0 database, which encompasses results from seven evaluation methods. These data were utilized to generate a heat map visualizing the relative fraction of immune cell infiltration within the tumor microenvironment (TME). To quantify the relative fraction of infiltrating immune cells and immune-related functions, we utilized the “ssGSEA” R package (35). Xu et al. created a website that provides us with a collection of genes related to cancer and immunology (36), as well as a collection of genes related to Mariathasan’s findings and a list of genes with favorable responses to anti-PD-L1 drugs (37). We utilized the GSVA approach to quantify both gene sets and evaluate their correlation with ERI. Furthermore, we explored the association between the model genes and the 51 immune genes, presenting the results in a circular heat map (38). To investigate variations in biological function across populations, we conducted Gene Set Enrichment Analysis (GSEA) using the MsigDB database “c2.cp.kegg.v6.2.” Additionally, we employed TIDE (http://tide.dfci.harvard.edu/), which stands for Tumor Immune Dysfunction and Rejection, as a computational framework to assess the potential for tumor immune escape based on the gene expression profile of tumor samples. The Cancer Immunome Atlas (TCIA) web tool provided comprehensive immunogenomic analysis results. The Immunophenotype Score (IPS), a quantitative tumor immunogenicity score ranging from 0 to 10, was used to denote the immunophenotype score (39). IPS can be used to predict response to immune checkpoint inhibitors.
2.7 Transfection of cells and real-time PCR
U251MG, LN229, and SW1783 human glioma cells and human astrocytes (NHA) were cultured in Dulbecco’s Modified Eagle’s Medium (DMEM, Gibco, C11995500BT, Canada) supplemented with 10% fetal bovine serum (FBS, Gibco, 10091148, Canada) and 1x penicillin/streptomycin (Gibco, 15140-122, Canada). All cultures were maintained in a CO2 incubator (TFS3111, USA) at 37°C with 5% CO2. BARD1 gene knockdown was achieved using small interfering RNA (siRNA). The specific BARD1 siRNA sequences can be found in Supplementary Table 1. In brief, cells were seeded at 50% confluency in 6-well plates and transfected with negative control (NC) and siBARD1 using Lipofectamine 3000 (Invitrogen, USA).
Total RNA was extracted from cell lines and tissues using TRIzol (Sigma-Aldrich, T9424, America) according to the manufacturer’s instructions. cDNA was synthesized using the PrimeScriptTM RT Reagent Kit (Takara, RR047, Japan). Real-time polymerase chain reaction (RT-PCR) was performed using SYBR Green Master Mix (Q111-02, Vazyme) to quantify mRNA expression levels normalized to GAPDH mRNA levels. The 2−ΔΔCt method was used to calculate the expression levels. All primers were provided by Qingdao BioScience (Beijing, China), and the primer sequences can be found in Supplementary Table 1.
2.8 Cell counting Kit8 assay, Wound-healing assay, and Transwell assay
First, cells (1000 cells per well) were seeded into a 96-well plate and incubated at 37°C for 4 hours with CCK-8 reagent (10 μL) (Dojindo, CK18, Japan). The absorbance was measured at a wavelength of 450 nanometers using an ELx800 plate reader (Thermo, Multiskan Spectrum, USA) to count the cells. Cell growth was represented as fold change from day 0 to day 4 and presented in a graph.
In a 6-well plate, transfected cells were electroplated and cultured in a cell culture incubator until 95% fusion. In each culture well, a straight line was scraped using a sterile 20 μL plastic pipette tip, and the unattached cells and debris were gently rinsed twice with PBS. Finally, the scratch wounds were photographed at 0 hours and 48 hours using Image J software, and the wound width was measured.
Cell invasion and migration studies were performed using a transwell assay. The upper chambers of a 24-well plate were filled with treated SW1783 cells (2×10^5 cells) and incubated for 48 hours. To evaluate the invasive and migratory abilities of the cells, the top surface of the plate was pre-coated with a matrix gel solution (BD Biosciences, USA) or left uncoated. The remaining cells at the bottom layer were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet (Solarbio, China) after removing the surface cells.
2.9 Statistical analysis
The statistical analyses were conducted using R version 4.1.3, 64-bit, along with its support packages. The Pycharm (3.9) integrated development environment for Python was also utilized. Prognostic values and comparisons of patient survival across different subgroups in each dataset were calculated using Kaplan-Meier survival analysis and the log-rank test. The non-parametric Wilcoxon rank sum test was employed to assess the relationship between two groups for continuous variables. Prognostic variables within the clinical characteristics of the different subgroups were identified through univariate and multivariate Cox regression analysis using the R package “survival.” Spearman correlation analysis was conducted to examine correlation coefficients. A significance level of P<0.05 was considered statistically significant for all statistical investigations.
3 Results
3.1 Scoring of exosome-related features in scRNA-seq
A brief flowchart depicting this study is displayed in Figure 1. Single-cell data from four samples were acquired based on the scRNA-seq data of GSE84465, aiming to explore the heterogeneity of glioblastoma and appraise disparities among cells situated within and adjacent to the tumor. Following stringent quality control filtration for the removal of low-quality cells, a total of 3533 cells were encompassed in subsequent analysis (Supplementary Figure 1A). After the elimination of batch effects and normalization of data, the top 2000 highly variable genes were chosen (Supplementary Figure 1B). For dimensionality reduction, the PCA technique was employed, with the top 20 principal components selected for further scrutiny based on P values (Supplementary Figure 1C). Data integration and batch effect elimination for the four single-cell samples were carried out using the Harmony algorithm, and t-SNE was employed to exhibit the distribution of cells before and after integration (Supplementary Figures 1D, E). Following the identification of 15 clusters, established marker genes were employed to annotate cell subpopulations. The distribution of different samples, tissue types, clusters, and cell subpopulations after annotation was demonstrated using t-SNE visualization (Figures 2A–D). Figure 2F delineates the proportional representation of distinct cell types within the four samples. The expression patterns of these genes across each cell type were showcased, predicated upon commonly used cellular marker genes (Figure 2E). The heatmap portrays the relative expression of marker genes within each cell subpopulation, with the top five marker genes calculated (Figure 2G). Grounded in the expression data of 121 exosome-related genes (ERGs), five prevalent algorithms (AddModuleScore, GSVA, AUCell, UCell, and singscore) were employed to assign scores to the gene sets within the single-cell data. The outcomes, depicted in Figures 2H, I, demonstrated relatively heightened exosome-related scores (ERS) in monocytes/macrophages and malignant cells. Moreover, a comparison between exosome-related scores in tumor and adjacent tissues unveiled intriguing observations: monocytes/macrophages, malignant cells, and oligodendroglial within tumors exhibited notably high ERS (Figure 2J). Following the integration of all malignant cells, a differential analysis of malignant cells within the tumor and adjacent tissues was conducted, revealing 33 ERGs with significant differential expression in both contexts (Supplementary Table 2).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1263329/fimmu-14-1263329-HTML/image_m/fimmu-14-1263329-g002.jpg)
Figure 2 Categorization of cellular subpopulations in GBM and enrichment scores of exosome-related genes. (A–D) t-SNE plots of different samples, tissue sources, clusters of cells, and cell subpopulations associated by color. (E) Expression of common marker genes used for cellular annotation in these cell subpopulations. (F) Histogram of cell type content for different samples. (G) Heatmap showing the relative expression of marker genes in eight cell subpopulations. Red color represents highly expressed genes and blue color represents lowly expressed genes. (H) Enrichment scores of exosome-related genes for each cell type in GBM are shown by bubble plots. (I) Enrichment scores of exosome-associated genes in each cell type are shown by t-SNE plots, with darker purple color having higher scores. (J) The difference in enrichment scores of exosome-related genes in each cell type in tumor and peripheral tissues. ns, Not significant; ** p< 0.01; **** P< 0.0001.
3.2 Pseudotime analysis and intercellular communication analysis
The evolution and differentiation of cells at the single-cell level can be inferred through cell trajectory analysis, accomplished by constructing intercellular trajectories to reshape cellular processes over time. The determination of cell trajectories and pseudotime distributions of malignant cells was carried out using the “monocle” R package, revealing the existence of seven cellular states during malignant cell development. For instance, the early state of cell development corresponds to cluster 5 (Figure 3A). The expression of 33 differentially expressed exosome-related genes at various developmental stages is depicted in the heatmap of Figure 3B, where BRAD1 is notably highly expressed in the early stages of malignant cell development. The “CellCall” tool was utilized to infer intercellular communication networks and internal regulatory signals by integrating intracellular and intercellular cues. By leveraging the median exosome-related scores across all malignant cells, we classified these cells into ERS high and ERS low groups. Bubble plots were employed to present the results of signaling pathway activity analysis, revealing that malignant cells with high scores exhibited intensified JAK-STAT signaling pathway activity in relation to monocytes/macrophages (Figure 3C). The JAK-STAT signaling pathway is implicated in pivotal biological processes like cell proliferation, differentiation, apoptosis, and immune regulation. The circle plot showcased in Figure 3D illustrates the intensity of direct ligand-receptor signaling across distinct cell types. Notably, a more pronounced CTF1-IL6ST ligand-receptor pair relationship was observed between malignant cells with high scores and monocytes/macrophages (Figure 3F). Subsequently, we shifted our focus to inferring the presence of ligand-receptor pairs and corresponding transcription factors (TFs) between high-scoring malignant cells and monocytes/macrophages (Figure 3E).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1263329/fimmu-14-1263329-HTML/image_m/fimmu-14-1263329-g003.jpg)
Figure 3 pseudotime analysis and cell communication analysis. (A) Cell trajectory and pseudo-time analysis for the malignant cells. (B) Heat map showing the expression of 33 differentially expressed exosome-related genes during cell development. Red represents high expression and blue represents low expression. (C) Bubble plots present the activity analysis of signaling pathways in different cell types. (D) The intensity of ligand receptors between different cell types is shown by circle plots. (E) Ligand receptor pairs and associated transcription factors between malignant cells with high exosome scores and monocytes/macrophages. (F) Ligand receptor intensities between different cell types.
3.3 Exosome-related features in spatial transcriptome sequencing
Spatial transcriptome sequencing data from a GBM patient was sourced from GSE194329, and following the exclusion of ribosomal and mitochondrial genes, the SCTransform method was applied to rectify sequencing depth and implement a series of normalizations. Through downscaled clustering, a total of 11 cellular subgroups were delineated within the spatial context (Figures 4A, B). The 33 exosome-related genes investigated in the prior study exhibited generally elevated expression across these 11 subgroups (Figure 4C). Guided by the original literature’s annotations, subpopulations 0 and 1 primarily occupied the GBM tumor core. Consequently, the “scMetabolism” R package was utilized for the assessment of metabolic activity within distinct cell subpopulations. Notably, subpopulations 0 and 1 exhibited close ties to Folate biosynthesis and tryptophan metabolic activity (Figure 4D). In this context, tryptophan, an essential amino acid, stands as a pivotal microenvironmental factor influencing the immunobiology of various tumor types. Its metabolism is integral to fostering immunosuppression, invasion promotion, and growth facilitation in malignant gliomas (40). As depicted in Figure 4E, heightened tryptophan metabolic activity was prominently enriched within the tumor’s core region. Leveraging Python’s Scanpy and stlearn packages, we conducted a trajectory analysis of cell subpopulations within the spatial domain. Employing pycharm, we normalized and clustered the spatial transcriptome data, yielding an additional set of 11 cell subpopulations (Figure 4F). Of interest, cluster 0, located within the tumor’s core, exhibited a differentiation trend toward cluster 5 situated in the periphery (Figure 4G). Ultimately, we employed the RCTD method to extrapolate the annotated cell types from single-cell data into the spatial dataset, inferring predominant cell types at each spatial location. Notably malignant cells with elevated exosome-related scores are predominantly located in the core of the tumor (Figure 4H).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1263329/fimmu-14-1263329-HTML/image_m/fimmu-14-1263329-g004.jpg)
Figure 4 correlation analysis of the spatial transcriptome. (A) UMAP plot showing the 11 clusters identified by stRNA-seq. (B) Spatial plot showing the 11 clusters identified by stRNA-seq. (C) Bubble plot presenting the expression intensity of exosome-related genes for different clusters. Red color represents high expression and blue color represents low expression. (D) Bubble plots presenting the metabolic intensities of different clusters. Red represents high expression, blue represents low expression. (E) The spatial plot of tryptophan metabolic intensity. (F) The spatial plot of 11 cellular clusters identified in python. (G) Spatial map showing the developmental trajectory of cluster 0 to cluster 5. (H) The distribution of different cell types in the spatial map was identified by the algorithm of RCTD.
3.4 Calculation of exosome-related index and construction of the prognostic model
To further utilize exosome-related features to aid clinical decision-making, 121 exosome-related genes (ERGs) were used to construct prognostic models in the bulk-seq of GBM. 168 GBM samples with survival information from the TCGA dataset were used as the training set for constructing prognostic risk models, and 388 GBM samples with survival information from the CGGA dataset were samples in the CGGA dataset used for external validation. First, the univariate cox analysis method was used to screen 17 ERGs affecting the overall survival (OS) of GBM patients (P<0.05, Figure 5A). To avoid prognostic signature overfitting and narrowing down the genes predicting OS, a 6-gene signature was constructed by Lasso-Cox regression analysis (Figures 5B, C). Risk scores/ERI were calculated for each patient based on the expression profile of each gene and the corresponding regression coefficients (Figure 5D). By weighting the estimated cox regression coefficients, our model yielded an exosome-related index (ERI) for each patient in the TCGA cohort. Based on the scoring formula, patients were divided into a low-ERI group and a high-ERI group using the median ERI as the cut-off point. The risk graph shows the detailed survival outcomes for each patient in the TCGA cohort and the CGGA external validation cohort (Figures 5E, F). Survival curves for both the training and validation groups showed that patients in the high-ERI group had worse OS compared to those in the low-ERI group (Figures 5G, H). In addition, the ROC curves showed that ERI performed well in predicting OS for these individuals in the TCGA cohort (Figure 5I).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1263329/fimmu-14-1263329-HTML/image_m/fimmu-14-1263329-g005.jpg)
Figure 5 Calculation of exosome-related index and construction of the prognostic model. (A) Forest plot showing the univariate cox analysis obtained for 17 prognosis-related genes. (B) LASSO coefficient profiles. (C) 10-time cross-validation for tuning parameter selection in the Lasso model. (D) Model genes and coefficients were determined by lasso regression and multivariate cox analysis. (E) Distribution of scores between the low ERI and high ERI groups in the TCGA cohort and patient survival. (F) Distribution of scores between the low-ERI and high-ERI groups in the CGGA cohort and the survival status of patients. (G) Kaplan-Meier survival curves for OS of patients in the low ERI and high ERI groups in the TCGA cohort. (H) Kaplan-Meier survival curves for OS of patients in the low-ERI and high-ERI groups in the CGGA cohort. (I) AUC values of ERI at 1, 2, and 3 years in the TCGA cohort.
3.5 Validation of the clinical value of the exosome-related index and construction of nomograms
In the TCGA cohort, ERI could be an independent prognostic indicator for patients compared to other common clinical features (age, grade, IDH mutation status, etc.) according to univariate and multivariate Cox analysis (Figures 6A, B). In addition to that, the area under the curve (AUC) of ERI at 1 year was much higher than other clinicopathological features (Figure 6E). To help clinicians make better clinical decisions, based on the correlation between the above clinicopathological features and ERI, we created a nomogram for predicting 1-year, 1.5-year, and 2-year survival rates in GBM patients (Figure 6C). The calibration curves were also able to show that the nomogram was able to make accurate predictions (Figure 6D). 1-year DCA curves (Figure 6F) and C-index values (Figure 6G) both showed that our constructed nomogram and ERI had the highest net benefit and that the risk model constructed based on the six ERGs was more influential in clinical decision-making than the traditional model in clinical decision making. The results of the chi-square test showed that grouping was associated with the IDH mutation status of patients (Figure 6H). The proportion of IDH mutations was higher in patients in the low-ERI group (Figure 6I). Based on the analysis of the results, we are more confident that ERI and nomograms are reliable clinical predictive scoring systems.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1263329/fimmu-14-1263329-HTML/image_m/fimmu-14-1263329-g006.jpg)
Figure 6 Independent prognostic analysis of exosome-related index (ERI) and clinicopathological factors in the TCGA cohort. (A, B) Univariate and multivariate Cox regression analyses of clinicopathological variables and OS risk scores in the TCGA training cohort. (C) Combined nomograms of age, sex, and IDH mutation status for predicting 1-, 1.5-, and 2-year OS in GBM patients. (D) Correction curves of nomograms. (E) AUC values for ERI and clinical characteristics at 1 year in the TCGA cohort. (F) DCA curves of ERI, nomogram scores, and other clinical features. (G) The predictive effect of different clinical features, nomogram scores, and ERI were evaluated using c - index curves. (H) Heat map of clinical characteristics associated with subgroups as determined by chi-square test. (I) The proportion of IDH mutation status in different score subgroups. *P< 0.05; ***P< 0.001.
3.6 Mutational landscape, analysis of immune infiltration and immune function
Given the pivotal role of genetic mutations in tailoring cancer patient treatments, the somatic mutation profiles of 121 ERGs were examined. As depicted in Supplementary Figures 2A, B, HRNR emerged as the most frequently mutated gene, predominantly featuring missense mutations. Additionally, the distribution of the most frequently mutated genes in gliomas within ERI subgroups was scrutinized, revealing a higher prevalence of TP53 and EGFR mutations in high-ERI subgroups (Supplementary Figure 2C). Moreover, an exploration of the co-mutation patterns among model genes unveiled the co-mutation of HOXC6 and DCT (P<0.05, Supplementary Figure 2D). However, no substantial discrepancy in TMB was observed between patients in the high-ERI and low-ERI groups (Supplementary Figures 2E, F). Subsequently, patients were categorized into four groups based on median TMB values and median ERI (H-TMB+high ERI, H-TMB+low ERI, L-TMB+high ERI, and L-TMB+low ERI). The outcomes indicated that patients with low ERI and high mutation burdens exhibited relatively improved prognoses (Supplementary Figure 2G).
The clinical outcome of patients and their response to therapy is influenced by the tumor microenvironment (TME). Among the factors within the TME, tumor-infiltrating immune cells play a significant role in impacting tumor progression and the effectiveness of antitumor therapies. Tumor-infiltrating immune cells (TIICs) constitute a crucial element of the TME, and their composition and distribution are intimately linked with tumorigenesis and progression (41). Consequently, an investigation was conducted into the immune landscape of high and low-ERI groups using algorithms from XCELL, TIMER, QUANTISEQ, MCPCOUNTER, CIBERSORT, CIBERSORT-ABS, and EPIC platforms (Figure 7A). To further delve into the relationship between the exosome-related index and immune cells and their functions, the enrichment scores of diverse immune cell subpopulations and immune functions were quantified through the “ssGSEA” method. The findings demonstrated that the high ERI group exhibited elevated scores for immune cell infiltration and immune pathway activation (Figure 7C). This encompassed certain immunosuppressive cells, including regulatory T cells (Tregs) and tolerogenic dendritic cells. Given the considerable influence of abnormal expression and function of immune checkpoint molecules on tumor immunotherapy, the expression of immune checkpoints (ICs) was evaluated across different ERI subgroups. Nearly all ICs showed heightened expression in the high ERI group (Figure 7E). Additionally, ESTIMATE was employed to compute the ratio of stromal and immune cells within various ERI subgroups, providing an estimate of tumor purity (Figure 7F). Heat maps were employed to visualize the expression of ICs, TME scores, and immune cell infiltration patterns in distinct subgroups (Figure 7B). Due to the contrasting prognostic outcomes and immune infiltration patterns in patients from the high-ERI and low-ERI groups, a Gene Set Enrichment Analysis (GSEA) was carried out to uncover potential disparities in biological functions between these two groups. For each group, the four most prominent signaling pathways were selected (Figure 7D). The high ERI group was notably associated with various cytokine-related pathways, whereas the low-ERI group exhibited enrichment in signaling pathways related to cell cycle and division. Based on the aforementioned findings, it is speculated that patients within the high ERI group may face a less favorable prognosis, although they exhibit more robust immune function. Patients with higher ERI might correspond to a tumor microenvironment in GBM characterized by immunosuppression, which could contribute to a reduced response rate to immunotherapy.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1263329/fimmu-14-1263329-HTML/image_m/fimmu-14-1263329-g007.jpg)
Figure 7 Immune microenvironment analysis of different ERI subgroups. (A) Differences in immune infiltration in different ERI subgroups were assessed using seven algorithms. (B) Heat map illustrating the differences in TME scores, immune checkpoint expression, and immune cell infiltration in different ERI subgroups (based on ssGSEA). (C) Radar plot illustrating the differences in immune cell infiltration and immune-related pathways calculated by ssGSEA between patients in different ERI subgroups. (D) GSEA analysis of the high-ERI and low-ERI groups focusing on the different enrichment of the KEGG pathway. (E) Differences in immune checkpoint expression between different ERI subgroups. (F) Differences in TME scores between patients in different ERI subgroups. * P<0.05, ** P<0.01, *** P<0.001.
3.7 Prediction of the effectiveness of immunotherapy
A recent study (42) has validated the remarkable advancements made in cancer treatment through immune checkpoint blockade (ICB) therapy. Gaining a deeper comprehension of the mechanisms underpinning cancer immunotherapy and translating this understanding into therapeutic strategies could potentially lead to prolonged survival for patients with limited treatment options. However, the efficacy of ICB therapy remains limited due to primary resistance, resistance development, and associated toxicities. To gain further insights into the role of risk scores in immunotherapy, the TIDE score was utilized to assess patients with potential abnormalities in immune function within tumors and regional lymph nodes. This approach aimed to judiciously identify candidates for immunotherapy. It is recognized that the impact of immunotherapy may fluctuate with tumor progression owing to variations in the degree of immune infiltration. Consequently, we probed whether prognostic models could predict responses to ICB in GBM patients. The utilization of the Immune Profile Score (IPS) can identify individuals who stand to benefit from immunotherapy. The violin plots illustrate the correlation between IPS values and risk groups, with higher IPS values indicating enhanced responses to PD-1 and CTLA-4 inhibitors. We anticipated a favorable immune response to CTLA-4 inhibitors in individuals from the low-ERI group (Figure 8A). Given the pivotal role of the immune microenvironment in mediating ICB responses, we further delved into the correlation between risk scores and ICB response characteristics. Notably, the ERI exhibited a significant negative correlation with these ICB response attributes (Figure 8B). In order to comprehensively explore the variations in immune responses across different subgroups, correlation analyses were conducted using six model genes and classical immune-related genes (Figure 8C). The TIDE score reflects tumor immune dysfunction and rejection, offering a computational framework to assess the likelihood of tumor immune evasion based on gene expression profiles from tumor samples. Elevated tumor TIDE prediction scores correlated with diminished ICB responses and inferior patient survival. Our findings revealed that patients in the high-ERI group exhibited heightened dysfunction scores, while those in the low-ERI group demonstrated elevated exclusion scores. Notably, there was no substantial discrepancy in TIDE scores between these two groups (Figure 8D).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1263329/fimmu-14-1263329-HTML/image_m/fimmu-14-1263329-g008.jpg)
Figure 8 Prediction of the effect of immunotherapy. (A) Comparison of the relative distribution of immune fraction (IPS) in the high-ERI and low-ERI groups. (B) Correlation of ERI with ICB response characteristics and each step of the tumor-immune cycle. (C) Heat map of model gene and immune gene correlations. (D) TIDE between GBM patients in the high-ERI and low-ERI groups. ns, not significant, * P<0.05, ** P<0.01.
3.8 BARD1+malignant cells as a target for prognosis and immunotherapy
Based on the intersection of differential and model genes in Supplementary Table 1, a total of three ERGs are thought to influence GBM progression and treatment (Figure 9A). Interestingly BARD1, CTSB, and GSTP1 were all highly expressed in tumor samples from the TCGA dataset (Figure 9B). Based on whether these three genes were expressed in malignant cells, we divided the malignant cells in single-cell sequencing into expression positive and negative groups. Based on the marker genes of BARD1, CTSB, and GSTP1 expressing positive cells, the algorithm of GSVA was used to impute the content of these cells in the bulk sequencing data. We categorized them into high and low groups based on the optimal cutoff value, and the group with high expression of BARD1 and CTSB+ malignant cells had a poorer prognosis (P<0.05, Figure 9C and Supplementary Figure 3A), whereas the difference in the overall survival curves between the two groups of GSTP1+ malignant cells was not statistically significant (P>0.05, Supplementary Figure 3D). Patients with a lower percentage of BARD1+ malignant cells may respond better to immunotherapy (P<0.05, Figure 9D). Meanwhile, the proportion of BARD1+ malignant cells was significantly positively correlated with cancer-associated fibroblasts (CAF) (Figure 9E). In contrast, patients with lower proportions of CTSB+ and GSTP1+ malignant cells may respond poorly to immunotherapy, and the comparison was not statistically significant (P>0.05, Supplementary Figures 3B, E). In the spatial transcriptome, it is clear that BARD1 expression is overall low, while CTSB and GSTP1 are mainly expressed in the core region of the tumor (Figure 9F, Supplementary Figures 3C, F). Thus, malignant cells with positive BARD1 expression may be a risk factor for GBM. To understand the interactions of these cells in space, we again performed deconvolution of BARD1+ malignant cells, BARD1- malignant cells, and other cells according to the RCTD method. The recipient-ligand interactions were inferred by Python’s stlearn, and Figure 9G shows the top 50 reciprocal ligand-receptor pairs. Among them, there is strong ligand-receptor communication between BARD1+ malignant cells and vascular cells (Figure 9H). Figure 9I demonstrates the spatial enrichment score of COL4A2-CD93 ligand-receptor pairs.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1263329/fimmu-14-1263329-HTML/image_m/fimmu-14-1263329-g009.jpg)
Figure 9 BARD1 is a marker of poor prognosis in GBM patients. (A) Three genes common to exosome-related and model genes were obtained by differential analysis in single-cell sequencing. (B) Expression of three genes in the TCGA cohort. (C) Kaplan-Meier survival curves for OS of patients in the high and low BARD1+ malignant cell expression groups. (D) Prediction of response to immunotherapy by TIDE based on the proportion of BARD1+ malignant cells. (E) Correlation between the proportion of BARD1+ malignant cells and CAF cell content. (F) Spatial map of BARD1 expression. (G) Interacting ligand-receptor pairs of top50. (H) Ligand-receptor communication of different cell types in space. (I) The spatial plot of COL4A2-CD93 ligand-receptor pair scores. *** p< 0.001.
3.9 Experimental validation of BARD1
The expression levels of BARD1 in HA cells and three GBM cell lines were initially compared through cell line experiments, revealing a notable upregulation of BARD1 in tumor cells (Figure 10A). Subsequently, the expression level of BARD1 was assessed 5 days post-transfection using qRT-PCR, to validate the effectiveness of siRNA-mediated BARD1 knockdown in the SW1783 cell line (Figure 10B). Following this, a CCK-8 cell assay demonstrated that the reduction of BARD1 resulting from knockdown significantly hindered the proliferative capacity of the SW1783 cell line (Figures 10C, F). Furthermore, the wound healing assay unveiled that the knockdown of the BARD1 gene notably impeded the migratory and invasive potential of SW1783 cells (Figures 10D, G). In alignment with the wound healing assay outcomes, GBM cells transfected with si-BARD1 exhibited diminished migratory and invasive capabilities in the transwell assay (Figures 10E, H). Collectively, these findings collectively indicate that BARD1 serves as a pro-carcinogenic factor in GBM.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1263329/fimmu-14-1263329-HTML/image_m/fimmu-14-1263329-g010.jpg)
Figure 10 Role of BARD1 in GBM. (A) BARD1 was highly expressed in GBM cell lines compared to normal human astrocyte NHA cell lines. (B) RT-qPCR was performed to detect the relative expression of BARD1 in GBM cells transfected with si-RNAs or negative control (NC). (C, F) CCK8 assay showed that SW1783 cells with reduced BARD1 expression had significantly reduced proliferative capacity compared to the NC group. (D, G) Scratch healing assay. the migration rate of SW1783 cells with reduced BARD1 gene expression was significantly reduced. (E, H) Transwell assay showed that down-regulation of BARD1 expression inhibited the migration and invasion ability of SW1783 cells. All data are expressed as mean ± SD of three independent experiments. * p< 0.05, ** p< 0.01, *** p< 0.001.
4 Discussion
Exosomes are extracellular vesicles that were first discovered 30 years ago. Since then, they have been used for intercellular communication, disease transmission, and drug development (43). Studies have reported that a single glioma cell can secrete approximately 10,000 exosomes within 48 hours (44). Numerous studies have shown that exosomes can influence the process of glioma development, growth, and metastasis by interacting with tumor cells and their surrounding microenvironment. Firstly, bioactive molecules within exosomes can impact tumor cell signaling pathways and gene expression through their transport and delivery mechanisms, consequently enhancing the proliferation, invasion, and metastatic potential of tumor cells (45). For instance, a study investigating the impact of glioma cells on the angiogenic process revealed that these cells could induce angiogenesis by transferring LncRNA-CCAT2 to endothelial cells through exosomes (46). Secondly, miRNAs and other nucleic acid molecules in exosomes can alter the gene expression pattern of tumor cells by horizontal transfer, affecting the phenotype and characteristics of tumors. For example, epithelial-mesenchymal transition (EMT) is an important factor that plays an important role in glioma progression, and inhibition of microRNA-708 in exosomes has been reported to increase cell proliferation and EMT in gliomas by promoting the SPHK2/AKT/catenin pathway (47). In addition, exosomes modulate the immune system response and inhibit immune cell activity, thus helping tumor cells to evade immune surveillance (48). Nonetheless, there remain unanswered questions and challenges concerning the role of exosomes in gliomagenesis. The precise mechanisms of exosome action and their interaction networks within glioma remain incompletely understood, underscoring the need for further investigations to unveil these aspects.
In recent years, the development of single-cell sequencing and spatial transcriptome technologies has brought new breakthroughs in tumor research and provided a deeper understanding of tumor types such as glioblastoma (GBM). However previous bioinformatics analyses rarely integrate them can compensate for each other and provide more comprehensive information (49). Single-cell sequencing technology allows us to perform high-throughput analysis of individual cells, revealing the heterogeneity and subpopulation structure within tumors. Single-cell sequencing allows us to identify different types of tumor cells, immune cells, and other cellular components and to study their functions and interactions in GBM. Nevertheless, single-cell sequencing does not furnish information about the spatial distribution of cells, a crucial aspect in the intricate organization of GBM. In contrast, spatial transcriptome technology provides the spatial distribution pattern of cells in tissues by analyzing gene expression on tissue sections. The spatial transcriptome allows us to understand the positional relationships of different cell types in GBM tissues and to explore the network of cellular interactions. Thus, the integration of single-cell sequencing and spatial transcriptome technologies assumes paramount significance in GBM research. Harnessing their combined potential enables us to acquire a more comprehensive and precise understanding of tumor cells and tissues, thus elucidating the pathogenesis of GBM, the tumor microenvironment, and potential therapeutic targets (50). This integrated analysis will drive tumor research in a deeper and more detailed direction and provide stronger support for future individualized therapy.
In this study, we employed single-cell sequencing and spatial transcriptomic techniques to characterize exosome-related genes in glioblastoma (GBM), thereby providing insight into the potential role of exosomes in GBM development and immunotherapy. Our results demonstrate that at the single-cell level, exosome-related genes exhibit significant heterogeneous expression patterns in GBM and show expression differences in specific cellular subpopulations and microenvironmental locations in spatial transcriptomic data. Single-cell sequencing data revealed heterogeneous expression of these genes in GBM cell populations, suggesting that exosome synthesis and function may be influenced by cellular heterogeneity. Our study found that gene enrichment scores associated with exosomes were significantly higher in GBM samples than in cells from peripheral samples. Thus due to the possible presence of higher levels of exosome release in tumor tissues (51). Cancer-associated cells secrete more exosomes than healthy cells because of the need to exchange information or nutrients between cells. It is estimated that the blood of cancer patients contains twice the amount of exosomes than that of healthy humans. It is reasonable to assume that in tumors, more exosomes are required to meet the intercellular communication needs due to the complex hypoxic environment in which tumors form (52, 53). Interestingly, when we brought this feature into the spatial transcriptome by means of deconvolution, we found that exosome features at the core locations of GBM tumors scored relatively higher. Thus the function of tumor cells in the core region of GBM is affected, leading to a decrease in their exosome secretion. We speculate that hypoxia in the tumor core region affects exosome secretion. Hypoxia is usually present in the tumor core region due to rapid growth and irregular blood supply. Hypoxic conditions alter the metabolic pathways and cellular signaling pathways of tumor cells, thus affecting exosome formation and release. Studies have shown that a hypoxic environment can lead to more exosome production by tumor cells and that the composition and function of exosomes may be altered (54). Under hypoxic conditions, tumor cells may increase exosome release in response to environmental stress and facilitate tumor cell interactions through exosome messaging. In addition, the hypoxic environment may lead to changes in the number and type of specific proteins, RNAs, or other biomolecules in exosomes, which may be associated with characteristics such as tumor aggressiveness, metastatic ability, and treatment resistance (55). Subsequently, we used spatial transcriptomic techniques to combine the expression of these exosome-related genes with cell type and tissue structure to reveal the spatial heterogeneity of exosomes in the GBM microenvironment. Of particular interest, we discovered a remarkably robust intercellular communication between tumor cells exhibiting higher exosome scores and vascular endothelial cells. Notably, angiogenesis, a critical process in glioma development, is known to be influenced by glioma-derived exosomes. These exosomes have been reported to play a significant role in driving angiogenesis during glioma progression. For instance, a study has suggested that glioma cells can induce angiogenesis by transferring Linc-CCAT2 to endothelial cells through exosomal transport. This finding highlights the involvement of exosomes in mediating angiogenesis and its potential impact on glioma pathogenesis (56). In addition, we observed that CD93 has a key role in tumor vascular maturation and extracellular matrix organization and is a potential therapeutic target, and a previous study showed that CD93 regulates Integrin β signaling activation and fibronectin fibril organization during tumor angiogenesis (57).
Considering the role of exosomes in GBM progression as immunosuppressive and aiding tumor metastasis, we developed a new exosome-related scoring system and associated index for risk stratification and prediction of personalized therapy (58). After analysis by lasso-cox regression, we classified patients into two subgroups based on ERI expression. Similar to many previous studies, patients in the high ERI group had a worse prognosis and were accompanied by higher TIICs and TME scores. The exosome-related signature we constructed proved to be effective for risk stratification of GBM patients and as an independent predictive variable to assess the survival of GBM patients. In addition, higher tumor microenvironment TME scores were associated with higher ERI. In several previous bioinformatics studies, mesenchymal and immune scores were found to be significantly higher with malignancy progression and to suggest an extremely poor prognosis (59), and since immunosuppressive cells within the TME render immunotherapy ineffective, TME is considered a red flag for GBM patients (60). Of course, in addition to the irreversible suppressive role played by TME in GBM, hypoxic conditions may also protect tumors from immune responses through various mechanisms, inhibit the activity of natural killer cells and connective tissue cells, and promote the release of various immunosuppressive cytokines and the enhancement of immunosuppressive cell function. Our study suggests that higher ERIs are associated with an immunosuppressive microenvironment, which explains why patients in the low-scoring subgroup showed an advantage in terms of OS and effective response to immunotherapy.
The therapeutic effect of immune checkpoint inhibitors (ICIs) alone for GBM is relatively limited due to the immune escape mechanism and immunosuppressive microenvironment of GBM. However, in recent years, researchers have begun to explore the use of ICIs in combination with other therapeutic approaches (e.g., radiotherapy, chemotherapy, vaccines, etc.) to improve the therapeutic efficacy of GBM (61). Currently, there are several clinical trials exploring the potential of ICIs in the treatment of GBM. Although partial responses and prolonged survival have been observed in some patients, overall, the efficacy of ICIs in the treatment of GBM remains relatively limited. Researchers are working to find more effective treatment strategies and combination regimens to overcome the resistance of GBM to immunotherapy (62). We investigated the distribution of common ICs between the high-ERI and low-ERI groups. Most ICs, including PD-1, CTLA-4, IDO, LAG-3, and TIM-3 were expressed higher in the low subgroup. The interaction of programmed cell death factor 1 (PD-1) protein and programmed cell death ligand 1 (PD-L1) protein generates an immunoregulatory axis that promotes GBM cell invasion in brain tissue (63). PD-L1 elevated in glioma cells binds to PD-1 on Tumor-Associated Macrophages (TAMs) and Tumor-Infiltrating Lymphocytes (TILs), induces a suppressive immune microenvironment, and is associated with poor prognosis in GBM patients (64, 65). The study of the GBM tumor immune cycle and ICB response revealed that ERI showed a significant negative correlation with ICB-related positive signals, while it showed a positive correlation with the suppressive tumor immune cycle. The above findings further support that immunosuppression is characteristic in the high ERI group. Finally, the calculation of the IPS score showed that a higher response rate to immunotherapy was associated with a lower ERI score and also corresponded to a better prognosis.
Interestingly, we finally identified a risk gene for GBM, BARD1, through differential screening at the single-cell sequencing level and selection of a model gene. In recent decades, researchers have investigated the role of the BARD1 gene in cancer progression and its use as a prognostic biomarker and potential candidate for targeted cancer therapy (66). BARD1 is a gene that encodes a protein that plays a critical role in the development and progression of different types of cancer. This protein plays a dual role in cancer as both a tumor suppressor and an oncogene (67). Notably, BARD1 shares a structurally homologous domain with BRCA1, and these two proteins interact to impede the progression of different cancers, including breast and ovarian cancers, through the BRCA1-dependent pathway (68). Furthermore, BARD1 has been implicated in other tumor suppression pathways, such as the tp53-dependent apoptosis signaling pathway (69). Mutations in BARD1 have been associated with susceptibility to various cancers, including lung, breast, and cervical cancers (70). These mutations may result in the generation of distinct BARD1 isoforms that differ from the full-length BARD1 protein. These isoforms may have a dominant negative effect, meaning that they interfere with the function of the full-length protein and promote tumor growth (71). BARD1 is now included in the clinical genome for cancer susceptibility testing. This means that mutations in BARD1 can be used as a diagnostic tool to identify individuals at higher risk of developing certain types of cancer (72, 73). In addition, BARD1 is being investigated as a potential therapeutic target for cancer treatment. However, there is still a gap in the study of BARD1 in GBM, and according to our findings, BARD1-positive GBM cells represent a poor prognosis and poorer immunotherapeutic efficacy, and subsequently, based on a series of knockout experiments, it is more confirmed that aberrant expression of BARD1 in gliomas can contribute to malignant transformation of cells and lead to proliferation of tumor cells, avoidance of apoptosis, and promote abnormal biological behaviors such as invasion and metastasis.
In this study, we performed a comprehensive analysis of exosome-related gene features expressed at the single-cell sequencing level and the spatial transcriptome level. Our study revealed differences in the expression of exosome-related genes in different cell subpopulations by single-cell sequencing, which provides clues for a deeper understanding of the role of exosomes in tumor immune escape and drug resistance mechanisms. On the other hand, through spatial transcriptome technology, we were able to explore the expression patterns of exosome-related genes at the tissue structure and to study their distribution characteristics within the tumor and surrounding microenvironment. In addition to this, we highlighted significant differences in anti-tumor immune response and immune status among different ERI groupings in GBM. Furthermore, we created a novel exosome-related index that provides some new data and novel findings for GBM biomarkers and their clinical applications. This signature not only correctly predicts the prognosis of GBM patients, but also provides additional benefits in terms of prognosis and personalized treatment for these high-risk patients.
However, despite the important application potential of single-cell sequencing and spatial transcriptome technologies in revealing exosome-related gene features, there are still some limitations (74). First, single-cell sequencing technologies may have certain noise and errors during the detection process, and have limitations on the capture rate and coverage of cells, which may result in the expression characteristics of certain cell subpopulations being overlooked or underestimated. Second, spatial transcriptome technologies are still in the developmental stage, and challenges remain for the analysis and interpretation of complex tissue structures. In addition, these technologies require more refined and efficient methods of data processing and analysis. And all our analyses were performed only on data from public databases, large prospective studies as well as additional in vivo and in vitro experimental studies are needed to confirm our findings, and most importantly this exosome-related signature needs to be validated in more real cohorts.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Author contributions
SZ: Conceptualization, Data curation, Methodology, Software, Visualization, Writing – original draft. QW: Formal Analysis, Methodology, Software, Visualization, Writing – original draft. KN: Investigation, Writing – original draft, Resources. PZ: Investigation, Writing – original draft, Methodology. YL: Investigation, Writing – original draft. JX: Investigation, Writing – original draft. WJ: Project administration, Supervision, Validation, Writing – review and editing. CC: Funding acquisition, Project administration, Supervision, Validation, Writing – review and editing. QZ: Funding acquisition, Supervision, Validation, Writing – review and editing.
Funding
This work was supported by the Wuxi Taihu Lake Talent Plan, Supports for Leading Talents in Medical and Health Profession(2020THRC-DJ-SNW) and Basic Research Project of Wuxi Science and Technology Bureau(K20221024).
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.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1263329/full#supplementary-material
Supplementary Figure 1 | quality control and principal component analysis of scRNA-seq. (A) Quality control of scRNA-seq data from GBM samples. (B) Variance plot showing 14,263 genes in all cells, with red dots representing the top 2,000 highly variable genes. (C) Principal component analysis was used to downscale and select the top 20 PCs.(D) t-SNE plots of four samples before over harmony to remove batch effects. (E) t-SNE plots of four samples after performing harmony to remove batch effects.
Supplementary Figure 2 | Mutation landscape of GBM samples. (A) Description of the details of mutations in all 121 exosome-related genes, where the most common mutation type is a missense mutation. SNPs account for an absolute proportion of mutations compared to insertions or deletions, and C>A occurs more frequently than other types. The horizontal histogram lists the top 10 mutated genes in GBM. (B) Mutation landscape of all 121 exosome-related genes. (C) Mutation landscape of the top 10 genes in mutation frequency in the two ERI subgroups. (D) Co-mutation or co-exclusion relationships between model genes. (E) Comparison of tumor mutational load (TMB) between different subgroups. (F) Correlation analysis of ERI and TMB. (G) Survival differences among the four subgroups of H-TMB+ high ERI, H-TMB+ low ERI, L-TMB+ high ERI, and L-TMB+ low ERI.
Supplementary Figure 3 | Analysis of CTSB and GSTP1. (A) Kaplan-Meier survival curves for OS in patients in the high and low CTSB+ malignant cell expression group. (B) Prediction of response to immunotherapy based on the proportion of CTSB+ malignant cells. (C) Spatial map of CTSB expression. (D) Kaplan-Meier survival curves for OS of patients in the high and low GSTP1+ malignant cell expression group. (E) Prediction of response to immunotherapy based on the proportion of GSTP1+ malignant cells. (F) Spatial plot of GSTP1 expression.
Supplementary Table 1 | siRNA and primer sequence information.
Supplementary Table 2 | 33 ERGs differentially expressed in malignant cells from the tumor and peripheral samples.
References
1. Wirsching HG, Galanis E, Weller M. Glioblastoma. Handb Clin Neurol (2016) 134:381–97. doi: 10.1016/B978-0-12-802997-8.00023-2
2. Sasmita AO, Wong YP, Ling APK. Biomarkers and therapeutic advances in glioblastoma multiforme. Asia Pac J Clin Oncol (2018) 14(1):40–51. doi: 10.1111/ajco.12756
3. Mathios D, Phallen J. Circulating biomarkers in glioblastoma: ready for prime time? Cancer J (2021) 27(5):404–9. doi: 10.1097/PPO.0000000000000541
4. Lynes JP, Nwankwo AK, Sur HP, Sanchez VE, Sarpong KA, Ariyo OI, et al. Biomarkers for immunotherapy for treatment of glioblastoma. J Immunother Cancer (2020) 8(1):e000348. doi: 10.1136/jitc-2019-000348
5. Pan S, Sun S, Liu B, Hou Y. Pan-cancer landscape of the runx protein family reveals their potential as carcinogenic biomarkers and the mechanisms underlying their action. J Transl Int Med (2022) 10(2):156–74. doi: 10.2478/jtim-2022-0013
6. Zhao S, Ji W, Shen Y, Fan Y, Huang H, Huang J, et al. Expression of hub genes of endothelial cells in glioblastoma-a prognostic model for gbm patients integrating single-cell rna sequencing and bulk rna sequencing. BMC Cancer (2022) 22(1):1274. doi: 10.1186/s12885-022-10305-z
7. Wang Z, Gao L, Guo X, Feng C, Lian W, Deng K, et al. Development of a nomogram with alternative splicing signatures for predicting the prognosis of glioblastoma: A study based on large-scale sequencing data. Front Oncol (2020) 10:1257. doi: 10.3389/fonc.2020.01257
8. Gong X, Chi H, Strohmer DF, Teichmann AT, Xia Z, Wang Q. Exosomes: A potential tool for immunotherapy of ovarian cancer. Front Immunol (2022) 13:1089410. doi: 10.3389/fimmu.2022.1089410
9. Dai J, Su Y, Zhong S, Cong L, Liu B, Yang J, et al. Exosomes: key players in cancer and potential therapeutic strategy. Signal Transduct Target Ther (2020) 5(1):145. doi: 10.1038/s41392-020-00261-0
10. Mashouri L, Yousefi H, Aref AR, Ahadi AM, Molaei F, Alahari SK. Exosomes: composition, biogenesis, and mechanisms in cancer metastasis and drug resistance. Mol Cancer (2019) 18(1):75. doi: 10.1186/s12943-019-0991-5
11. Chen X, Chi H, Zhao X, Pan R, Wei Y, Han Y. Role of exosomes in immune microenvironment of hepatocellular carcinoma. J Oncol (2022) 2022:2521025. doi: 10.1155/2022/2521025
12. Ge Y, Mu W, Ba Q, Li J, Jiang Y, Xia Q, et al. Hepatocellular carcinoma-derived exosomes in organotropic metastasis, recurrence and early diagnosis application. Cancer Lett (2020) 477:41–8. doi: 10.1016/j.canlet.2020.02.003
13. Xiong J, Chi H, Yang G, Zhao S, Zhang J, Tran LJ, et al. Revolutionizing anti-tumor therapy: unleashing the potential of B cell-derived exosomes. Front Immunol (2023) 14:1188760. doi: 10.3389/fimmu.2023.1188760
14. Chen G, Huang AC, Zhang W, Zhang G, Wu M, Xu W, et al. Exosomal pd-L1 contributes to immunosuppression and is associated with anti-pd-1 response. Nature (2018) 560(7718):382–6. doi: 10.1038/s41586-018-0392-8
15. Niu M, Liu Y, Yi M, Jiao D, Wu K. Biological characteristics and clinical significance of soluble pd-1/pd-L1 and exosomal pd-L1 in cancer. Front Immunol (2022) 13:827921. doi: 10.3389/fimmu.2022.827921
16. Qi Y, Jin C, Qiu W, Zhao R, Wang S, Li B, et al. The dual role of glioma exosomal micrornas: glioma eliminates tumor suppressor mir-1298-5p via exosomes to promote immunosuppressive effects of mdscs. Cell Death Dis (2022) 13(5):426. doi: 10.1038/s41419-022-04872-z
17. Wu J, Zeng D, Zhi S, Ye Z, Qiu W, Huang N, et al. Single-cell analysis of a tumor-derived exosome signature correlates with prognosis and immunotherapy response. J Transl Med (2021) 19(1):381. doi: 10.1186/s12967-021-03053-4
18. Zheng S, Liang JY, Tang Y, Xie J, Zou Y, Yang A, et al. Dissecting the role of cancer-associated fibroblast-derived biglycan as a potential therapeutic target in immunotherapy resistance: A tumor bulk and single-cell transcriptomic study. Clin Transl Med (2023) 13(2):e1189. doi: 10.1002/ctm2.1189
19. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics (2012) 28(6):882–3. doi: 10.1093/bioinformatics/bts034
20. Darmanis S, Sloan SA, Croote D, Mignardi M, Chernikova S, Samghababi P, et al. Single-cell rna-seq analysis of infiltrating neoplastic cells at the migrating front of human glioblastoma. Cell Rep (2017) 21(5):1399–410. doi: 10.1016/j.celrep.2017.10.030
21. Ren Y, Huang Z, Zhou L, Xiao P, Song J, He P, et al. Spatial transcriptomics reveals niche-specific enrichment and vulnerabilities of radial glial stem-like cells in Malignant gliomas. Nat Commun (2023) 14(1):1028. doi: 10.1038/s41467-023-36707-6
22. Wang X, Chai Z, Pan G, Hao Y, Li B, Ye T, et al. Exobcd: A comprehensive database for exosomal biomarker discovery in breast cancer. Brief Bioinform (2021) 22(3):bbaa088. doi: 10.1093/bib/bbaa088
23. Aran D, Looney AP, Liu L, Wu E, Fong V, Hsu A, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol (2019) 20(2):163–72. doi: 10.1038/s41590-018-0276-y
24. Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, et al. Fast, sensitive and accurate integration of single-cell data with harmony. Nat Methods (2019) 16(12):1289–96. doi: 10.1038/s41592-019-0619-0
25. Zhang P, Pei S, Gong Z, Feng Y, Zhang X, Yang F, et al. By integrating single-cell rna-seq and bulk rna-seq in sphingolipid metabolism, cacybp was identified as a potential therapeutic target in lung adenocarcinoma. Front Immunol (2023) 14:1115272. doi: 10.3389/fimmu.2023.1115272
26. Borcherding N, Voigt AP, Liu V, Link BK, Zhang W, Jabbari A. Single-cell profiling of cutaneous T-cell lymphoma reveals underlying heterogeneity associated with disease progression. Clin Cancer Res (2019) 25(10):2996–3005. doi: 10.1158/1078-0432.CCR-18-3309
27. Tirosh I, Izar B, Prakadan SM, Wadsworth MH 2nd, Treacy D, Trombetta JJ, et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell rna-seq. Science (2016) 352(6282):189–96. doi: 10.1126/science.aad0501
28. Andreatta M, Carmona SJ. Ucell: robust and scalable single-cell gene signature scoring. Comput Struct Biotechnol J (2021) 19:3796–8. doi: 10.1016/j.csbj.2021.06.043
29. Zhao Z, Ding Y, Tran LJ, Chai G, Lin L. Innovative breakthroughs facilitated by single-cell multi-omics: manipulating natural killer cell functionality correlates with a novel subcategory of melanoma cells. Front Immunol (2023) 14:1196892. doi: 10.3389/fimmu.2023.1196892
30. Zhang Y, Liu T, Hu X, Wang M, Wang J, Zou B, et al. Cellcall: integrating paired ligand-receptor and transcription factor activities for cell-cell communication. Nucleic Acids Res (2021) 49(15):8520–34. doi: 10.1093/nar/gkab638
31. Wu Y, Yang S, Ma J, Chen Z, Song G, Rao D, et al. Spatiotemporal immune landscape of colorectal cancer liver metastasis at single-cell level. Cancer Discovery (2022) 12(1):134–53. doi: 10.1158/2159-8290.CD-21-0316
32. Zou Y, Xie J, Zheng S, Liu W, Tang Y, Tian W, et al. Leveraging diverse cell-death patterns to predict the prognosis and drug sensitivity of triple-negative breast cancer patients after surgery. Int J Surg (2022) 107:106936. doi: 10.1016/j.ijsu.2022.106936
33. Yuan K, Zhao S, Ye B, Wang Q, Liu Y, Zhang P, et al. A novel T-cell exhaustion-related feature can accurately predict the prognosis of oc patients. Front Pharmacol (2023) 14:1192777. doi: 10.3389/fphar.2023.1192777
34. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612
35. Li XY, Zhao ZJ, Wang JB, Shao YH, Hui L, You JX, et al. M7g methylation-related genes as biomarkers for predicting overall survival outcomes for hepatocellular carcinoma. Front Bioeng Biotechnol (2022) 10:849756. doi: 10.3389/fbioe.2022.849756
36. Xu L, Deng C, Pang B, Zhang X, Liu W, Liao G, et al. Tip: A web server for resolving tumor immunophenotype profiling. Cancer Res (2018) 78(23):6575–80. doi: 10.1158/0008-5472.CAN-18-0689
37. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. Tgfbeta attenuates tumour response to pd-L1 blockade by contributing to exclusion of T cells. Nature (2018) 554(7693):544–8. doi: 10.1038/nature25501
38. Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang TH, et al. The immune landscape of cancer. Immunity (2019) 51(2):411–2. doi: 10.1016/j.immuni.2019.08.004
39. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep (2017) 18(1):248–62. doi: 10.1016/j.celrep.2016.12.019
40. Platten M, Weller M, Wick W. Shaping the glioma immune microenvironment through tryptophan metabolism. CNS Oncol (2012) 1(1):99–106. doi: 10.2217/cns.12.6
41. Low V, Li Z, Blenis J. Metabolite activation of tumorigenic signaling pathways in the tumor microenvironment. Sci Signal (2022) 15(759):eabj4220. doi: 10.1126/scisignal.abj4220
42. Wang Y, Zhang H, Liu C, Wang Z, Wu W, Zhang N, et al. Immune checkpoint modulators in cancer immunotherapy: recent advances and emerging concepts. J Hematol Oncol (2022) 15(1):111. doi: 10.1186/s13045-022-01325-0
43. Edgar JR. Q&A: what are exosomes, exactly? BMC Biol (2016) 14:46. doi: 10.1186/s12915-016-0268-z
44. Simundic AM, Baird G, Cadamuro J, Costelloe SJ, Lippi G. Managing hemolyzed samples in clinical laboratories. Crit Rev Clin Lab Sci (2020) 57(1):1–21. doi: 10.1080/10408363.2019.1664391
45. Gourlay J, Morokoff AP, Luwor RB, Zhu HJ, Kaye AH, Stylli SS. The emergent role of exosomes in glioma. J Clin Neurosci (2017) 35:13–23. doi: 10.1016/j.jocn.2016.09.021
46. Lang HL, Hu GW, Zhang B, Kuang W, Chen Y, Wu L, et al. Glioma cells enhance angiogenesis and inhibit endothelial cell apoptosis through the release of exosomes that contain long non-coding rna ccat2. Oncol Rep (2017) 38(2):785–98. doi: 10.3892/or.2017.5742
47. Chen Y, Deng X, Chen W, Shi P, Lian M, Wang H, et al. Silencing of microrna-708 promotes cell growth and epithelial-to-mesenchymal transition by activating the sphk2/akt/beta-catenin pathway in glioma. Cell Death Dis (2019) 10(6):448. doi: 10.1038/s41419-019-1671-5
48. Hellwinkel JE, Redzic JS, Harland TA, Gunaydin D, Anchordoquy TJ, Graner MW. Glioma-derived extracellular vesicles selectively suppress immune responses. Neuro Oncol (2016) 18(4):497–506. doi: 10.1093/neuonc/nov170
49. Mirzaei R, D'Mello C, Liu M, Nikolic A, Kumar M, Visser F, et al. Single-cell spatial analysis identifies regulators of brain tumor-initiating cells. Cancer Res (2023) 83(10):1725–41. doi: 10.1158/0008-5472.CAN-22-3004
50. Xiong A, Zhang J, Chen Y, Zhang Y, Yang F. Integrated single-cell transcriptomic analyses reveal that gpnmb-high macrophages promote pn-mes transition and impede T cell activation in gbm. EBioMedicine (2022) 83:104239. doi: 10.1016/j.ebiom.2022.104239
51. Wei Y, Wang D, Jin F, Bian Z, Li L, Liang H, et al. Pyruvate kinase type M2 promotes tumour cell exosome release via phosphorylating synaptosome-associated protein 23. Nat Commun (2017) 8:14041. doi: 10.1038/ncomms14041
52. Yu D, Li Y, Wang M, Gu J, Xu W, Cai H, et al. Exosomes as a new frontier of cancer liquid biopsy. Mol Cancer (2022) 21(1):56. doi: 10.1186/s12943-022-01509-9
53. Olejarz W, Kubiak-Tomaszewska G, Chrzanowska A, Lorenc T. Exosomes in angiogenesis and anti-angiogenic therapy in cancers. Int J Mol Sci (2020) 21(16):5840. doi: 10.3390/ijms21165840
54. Yu F, Liang M, Huang Y, Wu W, Zheng B, Chen C. Hypoxic tumor-derived exosomal mir-31-5p promotes lung adenocarcinoma metastasis by negatively regulating satb2-reversed emt and activating mek/erk signaling. J Exp Clin Cancer Res (2021) 40(1):179. doi: 10.1186/s13046-021-01979-7
55. Sun H, Meng Q, Shi C, Yang H, Li X, Wu S, et al. Hypoxia-inducible exosomes facilitate liver-tropic premetastatic niche in colorectal cancer. Hepatology (2021) 74(5):2633–51. doi: 10.1002/hep.32009
56. Cheng J, Meng J, Zhu L, Peng Y. Exosomal noncoding rnas in glioma: biological functions and potential clinical applications. Mol Cancer (2020) 19(1):66. doi: 10.1186/s12943-020-01189-3
57. Sun Y, Chen W, Torphy RJ, Yao S, Zhu G, Lin R, et al. Blockade of the cd93 pathway norMalizes tumor vasculature to facilitate drug delivery and immunotherapy. Sci Transl Med (2021) 13(604):eabc8922. doi: 10.1126/scitranslmed.abc8922
58. Che Mat MF, Abdul Murad NA, Ibrahim K, Mohd Mokhtar N, Wan Ngah WZ, Harun R, et al. Silencing of pros1 induces apoptosis and inhibits migration and invasion of glioblastoma multiforme cells. Int J Oncol (2016) 49(6):2359–66. doi: 10.3892/ijo.2016.3755
59. Zhang C, Cheng W, Ren X, Wang Z, Liu X, Li G, et al. Tumor purity as an underlying key factor in glioma. Clin Cancer Res (2017) 23(20):6279–91. doi: 10.1158/1078-0432.CCR-16-2598
60. Alghamri MS, Banerjee K, Mujeeb AA, Mauser A, Taher A, Thalla R, et al. Systemic delivery of an adjuvant cxcr4-cxcl12 signaling inhibitor encapsulated in synthetic protein nanoparticles for glioma immunotherapy. ACS Nano (2022) 16(6):8729–50. doi: 10.1021/acsnano.1c07492
61. Gatto L, Di Nunno V, Franceschi E, Tosoni A, Bartolini S, Brandes AA. Pharmacotherapeutic treatment of glioblastoma: where are we to date? Drugs (2022) 82(5):491–510. doi: 10.1007/s40265-022-01702-6
62. Qi Y, Liu B, Sun Q, Xiong X, Chen Q. Immune checkpoint targeted therapy in glioma: status and hopes. Front Immunol (2020) 11:578877. doi: 10.3389/fimmu.2020.578877
63. Grujic SD, O'Sullivan DD, Wehrmacher WH. Organizational control of hospital infrastructure determines the quality of care. Qual Assur Util Rev (1989) 4(1):19–24. doi: 10.1177/0885713x8900400105
64. Zhou J, Pei X, Yang Y, Wang Z, Gao W, Ye R, et al. Orphan nuclear receptor tlx promotes immunosuppression via its transcriptional activation of pd-L1 in glioma. J Immunother Cancer (2021) 9(4):e001937. doi: 10.1136/jitc-2020-001937
65. Gordon SR, Maute RL, Dulken BW, Hutter G, George BM, McCracken MN, et al. Pd-1 expression by tumour-associated macrophages inhibits phagocytosis and tumour immunity. Nature (2017) 545(7655):495–9. doi: 10.1038/nature22396
66. Tarsounas M, Sung P. The antitumorigenic roles of brca1-bard1 in DNA repair and replication. Nat Rev Mol Cell Biol (2020) 21(5):284–99. doi: 10.1038/s41580-020-0218-z
67. Hawsawi YM, Shams A, Theyab A, Abdali WA, Hussien NA, Alatwi HE, et al. Bard1 mystery: tumor suppressors are cancer susceptibility genes. BMC Cancer (2022) 22(1):599. doi: 10.1186/s12885-022-09567-4
68. Wu LC, Wang ZW, Tsan JT, Spillman MA, Phung A, Xu XL, et al. Identification of a ring protein that can interact in vivo with the brca1 gene product. Nat Genet (1996) 14(4):430–40. doi: 10.1038/ng1296-430
69. Westermark UK, Reyngold M, Olshen AB, Baer R, Jasin M, Moynahan ME. Bard1 participates with brca1 in homology-directed repair of chromosome breaks. Mol Cell Biol (2003) 23(21):7926–36. doi: 10.1128/MCB.23.21.7926-7936.2003
70. Hawsawi YM, Al-Numair NS, Sobahy TM, Al-Ajmi AM, Al-Harbi RM, Baghdadi MA, et al. The role of brca1/2 in hereditary and familial breast and ovarian cancers. Mol Genet Genomic Med (2019) 7(9):e879. doi: 10.1002/mgg3.879
71. Zhang YQ, Bianco A, Malkinson AM, Leoni VP, Frau G, De Rosa N, et al. Bard1: an independent predictor of survival in non-small cell lung cancer. Int J Cancer (2012) 131(1):83–94. doi: 10.1002/ijc.26346
72. Parvin JD. Brca1 at a branch point. Proc Natl Acad Sci U.S.A. (2001) 98(11):5952–4. doi: 10.1073/pnas.121184998
73. Fabbro M, Rodriguez JA, Baer R, Henderson BR. Bard1 induces brca1 intranuclear foci formation by increasing ring-dependent brca1 nuclear import and inhibiting brca1 nuclear export. J Biol Chem (2002) 277(24):21315–24. doi: 10.1074/jbc.M200769200
Keywords: exosome, glioblastoma, scRNA-seq, stRNA-seq, prognostic model, tumor microenvironment, immunotherapy, BARD1
Citation: Zhao S, Wang Q, Ni K, Zhang P, Liu Y, Xie J, Ji W, Cheng C and Zhou Q (2023) Combining single-cell sequencing and spatial transcriptome sequencing to identify exosome-related features of glioblastoma and constructing a prognostic model to identify BARD1 as a potential therapeutic target for GBM patients. Front. Immunol. 14:1263329. doi: 10.3389/fimmu.2023.1263329
Received: 19 July 2023; Accepted: 17 August 2023;
Published: 31 August 2023.
Edited by:
Zhirui Zeng, Guizhou Medical University, ChinaReviewed by:
Guoliang Ye, The Affiliated Hospital of Medical School of Ningbo University, ChinaJindong Xie, Sun Yat-sen University Cancer Center (SYSUCC) Guangzhou, China
Copyright © 2023 Zhao, Wang, Ni, Zhang, Liu, Xie, Ji, Cheng and Zhou. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Qiang Zhou, zhouqiang@kowloonhospital.com; Chao Cheng, Mr_chengchao@126.com; Wei Ji, wxrmyyjw@163.com
†These authors have contributed equally to this work