Three Liquid-Liquid Phase Separation-Related Genes Associated with Prognosis in Glioma

Introduction

Liquid-liquid phase separation (LLPS) is an intrinsic biophysical process within the cellular environment, leading to the formation of membrane-less compartments characterized by concentrated proteins and RNA molecules.¹ LLPS underlies the genesis of numerous membraneless organelles, including stress granules, processing bodies (P-bodies), and nuclear speckles. By orchestrating specific biological macromolecules into discrete membraneless compartments, cells can adeptly execute vital functions such as chromatin restructuring, and the regulation of gene transcription and translation.² LLPS is instrumental for many physiological functions, and its dysregulation can lead to pathological conditions. For example, pharmacological modulation of LLPS within the core regulatory circuitry attenuates metastatic potential and reverses chemoresistance in osteosarcoma.³ The induction of ILF3 liquid–liquid phase separation promotes esophageal cancer progression.⁴ Aberrant phase separation has been implicated in several diseases, including neurodegenerative disorders and cancers.^5–7 Thus, we believe that investigating the role of LLPS in oncology could enhance our comprehension of tumor pathogenesis, improve prognosis prediction, and refine personalized treatment choices.

Glioma represents the most prevalent malignant primary brain tumor, constituting over 81% of central nervous system malignancies.^8,9 Originating from glial or precursor cells, it typically evolves into astrocytoma, oligodendroglioma, ependymoma, or oligoastrocytoma.¹⁰ Significant heterogeneities in glioma patient prognosis and therapeutic responses consistently pose major clinical challenges.¹¹ Considering these challenges, our study sought to pinpoint these heterogeneities by examining the LLPS patterns in glioma patients. The variety of molecular subtypes is often related to the treatment and prognosis of patients. In this study, three molecular subtypes of glioma patients based on LLPS-related genes were identified with distinct prognosis and tumor immune microenvironment patterns. The prognostic model we developed demonstrated strong predictive capabilities for both prognosis and immune response. The present study is to reveal the correlation between LLPS-related genes and glioma progression. Our discoveries could pave the way for personalized prognosis assessment and optimized immunotherapy choices for glioma patients.

Materials and Methods

Data Collection

Gene expression profiling data of glioma was obtained from the dataset GSE50161 in the Gene Expression Omnibus (GEO) database.¹² This dataset comprises expression data from 34 glioma samples and 13 normal samples. The data was generated on the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 array) that completely covers the Human Genome U133 Set plus 6500 additional genes for analysis of over 47,000 transcripts. The annotations for the probes were based on the “SOFT formatted family file” associated with the GPL570 platform. Transcriptomic expression profiles and the corresponding survival time data for glioma patients were sourced from the UCSC Xena database (http://xenabroswer.net/hub).¹³ This dataset encompassed a total of 529 patient samples. The obtained data was subsequently processed and integrated with the GSE50161 dataset for subsequent analyses referring to reported bioinformatics workflows.^14,15

Screening for Differentially Expressed Genes (DEGs)

Gene expression data from the GSE50161 dataset was subjected to differential expression analysis using the limma package¹⁶ in R software version 4.2.2. Genes were deemed differentially expressed based on the criteria of a p-value < 0.05 and an absolute log2 fold change (|logFC|) > 1, with a Benjamini–Hochberg correction. LLPS-associated genes were curated from two databases: DrLLPS (http://llps.biocuckoo.cn/) and PhaSepDB (http://db.phasep.pro/). Duplicates were removed after the extraction of gene lists from both databases. To discern the genes that were both differentially expressed in the context of glioma and associated with LLPS, an intersection of the DEGs from GSE50161 and the LLPS-related genes was performed and visualized by Venn diagram.

Molecular Subtyping of Glioma

Based on the identified LLPS-related DEGs, the TCGA glioma dataset was subjected to molecular subtype identification using the ConsensusClusterPlus package¹⁷ in R version 4.2.2. ConsensusClusterPlus package has specialized application in consensus clustering analysis, particularly suited for high-dimensional data. The package offers simplicity in usage and provides various clustering algorithm options, including k-means, hierarchical clustering, and PAM. Notably, it supports running consensus clustering across different cluster numbers (k values), automatically selecting the optimal number, and offers visual tools such as trace plots and consensus matrices for result analysis. ConsensusClusterPlus package has been used for clustering analysis of glioma in previous studies.^18,19 In our analysis, we set parameters as follows: maxK = 10, reps = 5, pItem = 1, pFeature = 1, clusterAlg = “km” (k-means), distance = “euclidean”, and seed = 2. As a result, three distinct molecular subtypes, C1, C2, and C3, were delineated. The association between these subtypes and the expression levels of PD1 and PDL1 was further explored.

Immune Cell Infiltration Analysis

CIBERSORT, grounded on the linear support vector regression algorithm, is adept at inferring the proportions of distinct immune cell types within a sample based on gene expression data.²⁰ Utilizing the CIBERSORT algorithm in R version 4.2.2, we estimated the immune cell infiltration ratios for each sample in the TCGA glioma dataset. The immune cell signature genes used in our analysis were sourced from the specified benchmark database file (LM22.txt), which contains marker genes for 22 immune cell types.²¹ Additionally, we examined the disparities in immune cell infiltration levels across the three identified subtypes: C1, C2, and C3. Differences in immune cell infiltration levels across molecular subtypes were evaluated using the t-test statistical method. A p-value less than 0.05 was considered as indicative of a significant difference in immune cell infiltration between the two subtypes.

Selection of Survival-Associated Genes

The overall survival data utilized in our analysis were obtained from TCGA-LGG. Additional survival data, including progression-free interval, disease-specific survival, and more, can be accessed and downloaded from the UCSC Xena Hub (https://tcga-xena-hub.s3.us-east-1.amazonaws.com/download/survival%2FLGG_survival.txt). Incorporating patient survival data, a univariate Cox regression analysis was carried out on the LLPS-related DEGs genes using the survival package in R version 4.2.2. Applying a Log rank test with a threshold of p < 0.001, we pinpointed 312 genes that exhibited a significant association with survival duration. The stringent threshold of p < 0.001 was chosen to focus on the most compelling candidates and to mitigate potential false positives, contributing to the reliability and significance of our gene selection process.

Prognostic Model Construction and Validation

To avert the risk of overfitting in the constructed model, we randomly partitioned the TCGA glioma dataset, selecting 70% of the samples to serve as the training set and designating the remaining 30% as the validation set. The glmnet package in R version 4.2.2 was employed for lasso Cox multivariate regression, with the Log rank test for correction. We opted for the glmnet package due to its efficiency in adapting to various regression models, particularly in lasso Cox multivariate regression. The algorithm, driven by the regularization parameter lambda, swiftly computes parameters on the logarithmic scale. Glmnet’s speed is attributed to cyclic coordinate descent and Fortran subroutines. Parameter sets (n.seed = 4, nfolds = 7, split_train = 0.7, lasso.n = 11) were tailored for optimal performance in Cox regression analysis, considering the package’s versatility and computational efficiency. Leveraging the results of 7-fold cross-validation, the regularization parameter λ was set to lambda.1se. Subsequent Kaplan-Meier (KM) survival analysis and time-dependent ROC analysis were conducted on both the training and validation sets to validate the prognostic model.

Cell Culture

Human normal astrocytes HA1800 (Cat# YS1665C) and glioblastoma cell lines U87 (Cat# YS190C), U251 (Cat# YS301C), and SHG44 (Cat# YS266C) were purchased from the Y-J Biological, Shanghai, China. Cells were maintained in Dulbecco’s Modified Eagle Medium (DMEM) supplemented with 10% fetal bovine serum, 100 U/mL penicillin, and 100 μg/mL streptomycin. Cultures were kept in a humidified incubator at 37°C with 5% CO₂. The medium was replaced every 2 days. Upon reaching an approximate confluence of 80–90%, cells were routinely passaged.

RT-qPCR

Total RNA was isolated from cells using the Trizol reagent (Invitrogen, USA), and then converted into complementary DNA (cDNA) using FastKing One-step Genomic cDNA First Strand Synthesis Premix Kit (TIANGEN, China). RT-qPCR was conducted in a real-time PCR machine using the SYBR Green PCR Master Mix (Lifeint, China) and using GAPDH as an internal reference. The sequence information of amplification primers is presented in Table 1. The relative transcription level of the target gene was calculated using the 2^−ΔΔCT method. Experiments were performed in triplicate.

Table 1 Primer Sequences

Statistical Analysis

Experimental data were presented as mean ± standard deviation. Differences between groups were analyzed using one-way analysis of variance, followed by Tukey’s test. Statistical analyses were conducted using GraphPad 7.0 software, with p-values < 0.05 considered statistically significant.

Results

LLPS-Related Genes Differentially Expressed in Glioma

We conducted differential gene expression analysis on dataset GSE50161 of glioma. Based on a threshold of p < 0.05 and |logFC| > 1, we identified 3227 differentially expressed genes, among them, 1477 genes were down-regulated and 1750 were up-regulated in glioma compared to the normal samples (Figure S1A). These DEGs showed clearly different expression patterns in glioma and normal samples (Figure S1B). From the DrLLPS and PhaSepDB databases, a total of 3773 LLPS-related genes were recognized after removal of duplicates, including 356 upregulated and 317 downregulated genes (Figure S1C-D). By intersecting the set of LLPS-related genes with DEGs, we identified a common set of 673 genes (Figure 1).

Figure 1 Intersection of liquid-liquid phase separation (LLPS)-related genes with glioma-related differentially expressed genes (DEGs), identifying a common set of 673 genes.

Molecular Subtyping of Glioma

According to the identified differentially expressed LLPS-related genes, a refined disease subtyping was performed on the TCGA-glioma dataset. Setting the cluster count (k) to three, three distinct molecular subtypes (C1, C2, and C3) of glioma were obtained (Figure 2A). C1, C2, and C3 contains 263, 109, and 157 glioma samples, respectively. The optimal number of clusters was determined by the cumulative distribution function (CDF), and we observed relatively stable clustering when k = 3 (Figure 2B). Principal component analysis (PCA) indicated a significant difference among three clusters (Figure 2C). Further analysis of the prognostic characteristics of these three molecular subtypes showed that the prognosis of C1 and C3 was poor, and the prognosis of C2 was good, with significant differences (p < 0.0001, Figure 2D). We further delved into the interplay between these subtypes and the expression levels of the immune checkpoint genes PD1 and PDL1. A marked variance in the expression of PD1 and PDL1 was evident across the three molecular subtypes (t-test, p < 0.001). The expression of PD1 and PDL1 was significantly higher in C2 than those in C1 and C3 (p < 0.001, Figure 2E and F).

Figure 2 Molecular subtyping and analysis of glioma samples. (A) Three distinct molecular subtypes (C1, C2, and C3) derived from the TCGA-glioma dataset. (B) Cumulative distribution function (CDF) displaying the optimal clustering stability at k = 3. (C) Principal component analysis (PCA) demonstrating significant differentiation among the C1, C2, and C3 glioma molecular subtypes. (D) Prognostic characteristics of C1, C2, and C3 molecular subtypes. (E-F) Expression levels of immune checkpoint genes PD1 and PDL1 across the three molecular subtypes. ***p < 0.001, ****p < 0.0001.

Immune Infiltration Among Three Molecular Subtypes

We further analyzed the immune cell infiltration level and immune function level of the three subclusters. As expected, we observed that C2 had more immune cell infiltration and immune response (Figure 3). CIBERSORT analysis results showed that the infiltration level of resting dendritic cells, eosinophils, M0 macrophages, M1 macrophages, monocytes, neutrophils, plasma cells, naive CD4 T cells, and follicular helper T cells in C2 was higher than those in C1 and C3 (p < 0.05, Figure S2). The infiltration levels of naive B cells, activated mast cells, resting mast cells, activated NK cells, resting NK cells, activated memory CD4 T cells, resting memory CD4 T cells, CD8 T cells, delta gamma T cells, and regulatory T cells (Tregs) in C2 was lower than those in C1 and C3 (p < 0.0001, Figure S3).

Figure 3 Immune cell infiltration and immune response levels among the three glioma molecular subtypes.

Construction and Performance Evaluation of the Prognostic Model

By integrating patient survival time data, a univariate Cox regression analysis was performed on the 673 intersecting genes. Based on a Log rank test with a threshold of p < 0.001, 312 genes were identified as significantly correlated with survival duration. Through the lasso Cox regression, we discerned a composite prognostic signature comprising three genes (TAGLN2, NTNG2, and IGF2BP2) from the 312 survival-correlated genes (Figure 4A). The univariate Cox regression analysis for 312 survival-correlated genes was displayed in Figure S4. The risk score equation derived from the three genes is:

Figure 4 Prognostic evaluation based on selected genes and risk assessment. (A) Prognostic signature derived from three genes (TAGLN2, NTNG2, and IGF2BP2) identified from the intersecting gene set. (B-C) Kaplan-Meier survival curves representing high and low-risk groups from both training and validation sets. The accompanying Receiver Operating Characteristic (ROC) analysis highlights the prognostic accuracy of the model.

risk_score = 0.025*TAGLN2 - 0.206*NTNG2 + 0.069*IGF2BP2.

Based on this risk score equation, samples with scores above the median were categorized as high-risk, while the remainder were designated as low-risk. To mitigate the risk of overfitting in the predictive model, we randomly selected 70% of the samples from the glioma TCGA dataset to serve as the training set for model construction, with the remaining 30% designated as the validation set. The KM showed significant survival differences between high and low-risk groups in both training and validation sets (Log rank test, p < 0.0001). High-risk group presented significantly lower survival probability than low-risk group (p < 0.05, Figure 4B and C). ROC analysis showed that the prognostic model had good prognostic properties for 1-year (AUC: 0.842), 3-year (AUC: 0.859), and 5-year (AUC: 0.757) survival in the training set (Figure 4B). The good prognostic properties also presented for 1-year (AUC: 0.932), 3-year (AUC: 0.843), and 5-year (AUC: 0.684) survival in the validation set (Figure 4C).

Immune Infiltration and Expression Analysis of Three Prognostic Genes in Glioma

We further analyzed the correlation between three prognostic genes (TAGLN2, NTNG2, and IGF2BP2) and immune cells in glioma samples (Figure 5A). The obtained results showed that some immune cells strongly correlated with three prognostic genes (TAGLN2, NTNG2, and IGF2BP2). TAGLN2 was significantly negatively correlated with regulatory T cells (Tregs), resting mast cells, and activated mast cells and significantly positively correlated with follicular helper T cells, naive CD4 T cells, plasma cells, neutrophils, and eosinophils. NTNG2 was significantly negatively correlated with follicular helper T cells, naive CD4 T cells, neutrophils, and eosinophils and significantly positively correlated with resting mast cells and activated mast cells. IGF2BP2 was significantly negatively correlated with regulatory T cells (Tregs), resting mast cells, and activated mast cells and significantly positively correlated with follicular helper T cells and eosinophils. In addition, of these genes, the expression of TAGLN2 and IGF2BP2 in glioma tumor samples was significantly higher than those in normal samples, while the expression of NTNG2 was significantly lower based on the GSE50161 dataset (Figure 5B). Consistently, RT-qPCR confirmed the upregulation of TAGLN2 and IGF2BP2 and the downregulation of NTNG2 in glioma cells (p < 0.01, Figure 5C).

Figure 5 Correlation and expression analysis of prognostic genes in glioma. (A) Analysis illustrating the correlation between the three prognostic genes (TAGLN2, NTNG2, and IGF2BP2) and immune cells in glioma samples. (B) Differential expression levels of TAGLN2, NTNG2, and IGF2BP2 in glioma tumor samples in comparison to normal samples using the GSE50161 dataset. *p < 0.05, ****p < 0.0001. (C) RT-qPCR results confirming the altered expression patterns of TAGLN2, NTNG2, and IGF2BP2 in glioma cells. **p < 0.01 vs HA1800.

Discussion

The LLPS mechanism is crucial in both the initiation and advancement of tumors.²² LLPS leads to the dysfunction of biomolecular condensates, thereby contributing to the tumorigenesis and progression. For instance, MALR/ILF3-mediated LLPS indicates poor prognoses for patients with esophageal cancer and promote cancer progression.⁴ The emergence of distinct tumor immune microenvironment (TIME) profiles is also linked to LLPS, given its pivotal role in modulating immune signaling pathways.²³ The induction of LLPS of mutant NF2 imprisons the cGAS-STING machinery to abrogate antitumor immunity.²⁴ We postulated that an in-depth investigation of LLPS-related biomarkers could offer valuable insights for distinguishing novel tumor subtypes and predicting prognosis as well as immunotherapy responsiveness. Predicting the prognosis of gliomas has been considered more challenging than for other malignancies.²⁵ A previous study provided a classification of lower-grade glioma (LGG) patients based on LLPS. The LLPS-related prognostic risk score was correlated with clinicopathological features, cancer hallmarks, genomic alterations, and TIME patterns of LGG patients.²⁶ In this study, we focused on identifying LLPS-related DEG in glioma. A total of 673 LLPS-related DEGs in glioma were obtained by intersecting 3227 glioma-related DEGs from the GSE50161 dataset and 3773 LLPS-related genes from the DrLLPS and PhaSepDB databases. Based on the expression profiles of 673 LLPS-related DEGs, we identified three molecular subtypes of glioma patients in the TCGA database. Then, significant differences among three molecular subtypes were observed regarding prognosis and immune microenvironment. For personalized integrative evaluations, a prognostic model containing three LLPS-related DEGs (TAGLN2, NTNG2, and IGF2BP2) was constructed via the LASSO Cox regression analysis. Results revealed that the prognostic model had good predictive properties for prognosis and immune microenvironment of glioma patients.

Our analysis developed a novel prognostic model based on LLPS-related DEGs (TAGLN2, NTNG2, and IGF2BP2) for predicting the prognosis of glioma. A study revealed a multi-stage LLPS process during osmotic stress, in which rapid clustering of IGF2BP proteins initiates stress granules assembly.²⁷ MYC translation mediating PABP1/IGF2BPs liquid–liquid phase separation contributes to tumor progression.²⁸ Mechanism studies of TAGLN2 and NTNG2 in LLPS are few; however, their functions in brain tumors, including glioma, have been reported. This is helpful for our discovery of candidate key molecules of glioma from the perspective of LLPS. It has been reported that TAGLN2 is a candidate prognostic biomarker promoting tumorigenesis in human gliomas.²⁹ Succinylation of lysine-40 on TAGLN2 promotes angiogenesis and tumor growth in glioma.³⁰ Inhibiting TAGLN2 suppresses the malignant progression of glioma and augments sensitivity to temozolomide.³¹ Homozygous missense variants in NTNG2, encoding netrin-G2 involved in neuronal circuitry and synaptic organization in vertebrates, result in a unique neurodevelopmental disorder.³² The pan-cancer landscape of netrin family reveals that the participation of NTNG2 in various cancers, including glioma, shows potential as tumor markers and therapeutic targets.³³ A pan-cancer landscape demonstrates the IGF2BP3 is a potential negative regulator of glioma tumorigenesis by modulating stemness.³⁴ IGF2BP2 SUMOylation stimulates vasculogenic mimicry formation in glioma cells and boosts xenograft growth in nude mice, suggesting a novel approach for glioma molecular targeted therapy.³⁵ IGF2BP2 can induce drug resistance in glioblastoma cells and promote glioma progression.³⁶ The prognosis model consisting TAGLN2, NTNG2, and IGF2BP2 showed good prognostic properties for 1-year (AUC: 0.842), 3-year (AUC: 0.859), and 5-year (AUC: 0.757) survival for patients with glioma. By bioinformatics analysis and in vitro experimental validation, we found that TAGLN2 and IGF2BP2 were upregulated in glioma tumor tissues and cells, while NTNG2 was downregulated. These studies suggest the prognostic roles of TAGLN2, NTNG2, and IGF2BP2 in glioma. However, their involvement in the LLPS mechanism remains to be further investigated.

We further investigated the association between LLPS-related DEGs and TIME. Within the tumor stroma, immune cells serve as fundamental components and play a vital role in its intricate dynamics. An accumulating body of research indicates that both innate immune cells (macrophages, neutrophils, dendritic cells, innate lymphoid cells, myeloid-derived suppressor cells, and NK cells) and adaptive immune cells (T cells and B cells) can foster tumor progression when situated within the tumor microenvironment.³⁷ Based on the LLPS-related DEGs, we identified three molecular subtypes of glioma (C1, C2, and C3). Different LLPS subtypes of glioma patients exhibited distinct survival probability. The C2 subtype had significantly lower survival probability than C1 and C3 subtypes. There are also significant differences in the expression levels of the PD1 and PDL1 immune checkpoint genes across the three molecular subtypes. Furthermore, the degree of immune cell infiltration varies significantly among the three molecular subtypes, potentially leading to disparities in tumor progression and the efficacy of immunotherapy. In addition, we found that these three prognostic biomarkers (TAGLN2, NTNG2, and IGF2BP2) were closely related to immune microenvironment in glioma. TAGLN2 and IGF2BP2 was significantly negatively correlated with resting and activated mast cells, while NTNG2 was significantly positively correlated with resting and activated mast cells. A large-scale bulk RNA-seq analysis defines immune evasion mechanism related to mast cell in gliomas.³⁸ Mast cells modulate proliferation, migration, and stemness of glioma cells through downregulation of GSK3β expression and inhibition of STAT3 activation.³⁹ Therefore, the LLPS-related DEGs (TAGLN2, NTNG2, and IGF2BP2) may regulate glioma progression via affecting mast cells.

Conclusion

Our research underscores the significance of LLPS-related genes in glioma prognosis. We identified three distinct glioma molecular subtypes, with subtype C2 demonstrating a favorable prognosis and increased immune response. A prognostic model was formulated based on three pivotal LLPS-related genes (TAGLN2, NTNG2, and IGF2BP2), showing strong predictive abilities for patient outcomes. However, the LLPS-related prognostic model was proposed based on a public database, which may exist potential biases during model construction and validation. This model needs to be validated by prospective clinical experiments. It also requires molecular function experiments to investigate prognostic biomarkers (TAGLN2, NTNG2, and IGF2BP2) in LLPS mechanisms and glioma.