Molecular landscapes of glioblastoma cell lines revealed a group of patients that do not benefit from WWOX tumor suppressor expression

Introduction Glioblastoma (GBM) is notorious for its clinical and molecular heterogeneity, contributing to therapeutic failure and a grim prognosis. WWOX is one of the tumor suppressor genes important in nervous tissue or related pathologies, which was scarcely investigated in GBM for reliable associations with prognosis or disease progression despite known alterations. Recently, we observed a phenotypic heterogeneity between GBM cell lines (U87MG, T98G, U251MG, DBTRG-05MG), among which the anti-GBM activity of WWOX was generally corresponding, but colony growth and formation were inconsistent in DBTRG-05MG. This prompted us to investigate the molecular landscapes of these cell lines, intending to translate them into the clinical context. Methods U87MG/T98G/U251MG/DBTRG-05MG were subjected to high-throughput sequencing, and obtained data were explored via weighted gene co-expression network analysis, differential expression analysis, functional annotation, and network building. Following the identification of the most relevant DBTRG-distinguishing driver genes, data from GBM patients were employed for, e.g., differential expression analysis, survival analysis, and principal component analysis. Results Although most driver genes were unique for each cell line, some were inversely regulated in DBTRG-05MG. Alongside driver genes, the differentially-expressed genes were used to build a WWOX-related network depicting protein–protein interactions in U87MG/T98G/U251MG/DBTRG-05MG. This network revealed processes distinctly regulated in DBTRG-05MG, e.g., microglia proliferation or neurofibrillary tangle assembly. POLE4 and HSF2BP were selected as DBTRG-discriminating driver genes based on the gene significance, module membership, and fold-change. Alongside WWOX, POLE4 and HSF2BP expression was used to stratify patients into cell lines-resembling groups that differed in, e.g., prognosis and treatment response. Some differences from a WWOX-related network were certified in patients, revealing genes that clarify clinical outcomes. Presumably, WWOX overexpression in DBTRG-05MG resulted in expression profile change resembling that of patients with inferior prognosis and drug response. Among these patients, WWOX may be inaccessible for its partners and does not manifest its anti-cancer activity, which was proposed in the literature but not regarding glioblastoma or concerning POLE4 and HSF2BP. Conclusion Cell lines data enabled the identification of patients among which, despite high expression of WWOX tumor suppressor, no advantageous outcomes were noted due to the cancer-promoting profile ensured by other genes.


Introduction
Glioblastoma (GBM) is one of the most aggressive and frequently occurring brain tumors.Despite extensive research efforts and advancements in treatment modalities, the prognosis for GBM patients remains dismal, with a median survival of only around 12 months (Blakstad et al., 2023).The molecular complexity of glioblastoma entails significant challenges in understanding its mechanisms and developing effective therapeutic strategies (Singh et al., 2004;Pang et al., 2019).Combined with virtually inevitable relapse that only worsens the prognosis, the significance of subsequent research on this tumor entity is justified (Birzu et al., 2020;Oronsky et al., 2020).Pathogenetic GBM-related mechanisms are still far from being unraveled, and the application of current anti-glioblastoma therapies brings only a marginal improvement of the patients' outcome due to the tumor location or acquired drug resistance (Stupp et al., 2009;Peignan et al., 2011;Delgado-Martin and Medina, 2020).One of the GBM traits hindering the effective treatment is the high heterogeneity that manifests on both intra-tumoral and inter-tumoral levels (Neftel et al., 2019;Dymova et al., 2021).Other clinical implications include the inconsistency between biomarkers identified in a single tumor biopsy or the need to consider clinically relevant subpopulations in planning of appropriate treatment regimen (Parker et al., 2015).This emphasizes a need for research identifying distinct molecular landscapes that may allow for a broader understanding of GBM complexity.
Recently, we investigated the phenotypic heterogeneity of four GBM cell lines (U87MG, T98G, U251MG, and DBTRG-05MG) and revealed some discrepancies in the results from in vitro assays following the ectopic overexpression of WW domain-containing oxidoreductase (WWOX) (Kaluzinska-Kolat et al., 2023;Varricchio et al., 2023).The literature on the role of this gene in brain development and pathology is generally of a reasonable extent and emphasizes its importance in nervous tissue (Kosla et al., 2020).However, it is acknowledged that less attention has been paid to determining the influence of WWOX on nervous system tumors than in other cancers (Winardi et al., 2013), even though WWOX affects the prognosis and treatment response of glioma (Liu et al., 2015).In GBM, the downregulation of this haploinsufficient tumor suppressor gene stems from the loss of heterozygosity and promoter methylation (Kosla et al., 2011).It has been shown that WWOX mitigates the glioblastoma's infiltrative potential, and its loss can promote the migration of cancer cells (Chiang et al., 2012(Chiang et al., , 2013)).Our previous research concluded that the anti-GBM activity of WWOX is mainly a consequence of reduced cell viability and invasion (Kaluzinska-Kolat et al., 2023).However, while most of the data indicated the antineoplastic activity of WWOX, discrepancies related to colony growth and formation were noticed in one of the above GBM cell lines, i.e., DBTRG-05MG.This suggests that DBTRG-05MG exhibited a distinct expression profile in comparison to U87MG/T98G/ U251MG, which prompted us to perform a High-throughput sequencing (HTS) on cellular variants from the previous study, with an intention to identify the differences that may help in stratifying GBM patients into groups of various clinical outcomes.The usage of high-throughput sequencing technologies over the years has played a considerable role in elucidating the molecular intricacies of glioblastoma (Jovcevska, 2020).To ensure the robustness of data in such studies, it is advisable to prepare biological replicates for RNA-sequencing (RNA-Seq), seeing that the use of only technical replicates is uncertain (Blainey et al., 2014).The current study considered this recommendation and aimed to identify molecular landscapes that distinguish DBTRG-05MG from U87MG/T98G/ U251MG and may clarify the phenotypic WWOX-related differences or translate them into the clinical context.

Cell culture, stable transduction, and confirmation of acquiring cellular variants
The cell culture, procedure of lentiviral transduction for WWOX upregulation, and its confirmation are described in our preceding study (Kaluzinska-Kolat et al., 2023).Briefly, four cell lines representing glioblastoma (U87MG, T98G, U251MG, DBTRG-05MG) were purchased from the European collection of cell cultures (ECACC), incubated at 37°C in a humidified atmosphere of 5% CO 2 , and cultured in MEM (T98G), EMEM (U87MG and U251MG) or RPMI-1640 (DBTRG-05MG) medium supplemented with L-glutamine, fetal bovine serum, and Antibiotic-Antimycotic.To overexpress WWOX, stable transduction was performed via the GIPZ™ lentiviral system in a polybrene-containing starvation medium.Following the exchange of viral medium to full medium, the clonal selection was performed using puromycin.In each glioblastoma cell line, the stable transductants represented either the "WWOX" cellular variant (treated with the GIPZ™ system) or the "Vec" cellular variant (treated with the control Puro-Blank system).After protein extraction, the efficiency of transduction was confirmed by Western Blot analysis.Additional characterization of four cell lines is summarized in Supplementary Table S1; data included the genetic alterations acquired from Cancer Cell Line Encyclopedia (Ghandi et al., 2019), as well as the putative disease subtypes obtained from the literature (Zhong et al., 2014;Kiseleva et al., 2016;Salvadores et al., 2020;Oliveira et al., 2022;Qiu et al., 2023;Varricchio et al., 2023).Total RNA was isolated using an Extracol reagent per the manufacturer's protocol (EURX, Gdansk, Poland).The quality of total RNA was assessed by Agilent 2,100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States) to confirm that RNA integrity number is over 8.0.Reverse transcription of RNA was performed using random primers (CAGE library preparation Kit; K.K. DNAFORM, Yokohama, Japan).The selection of RNA/cDNA hybrids was enabled by cap-trapping on streptavidin beads.Following the digestion of RNA using RNaseI/H, the linker ligation to 5′ and 3' cDNA ends enabled the construction of double-stranded cDNA libraries.Sequencing of libraries using the Cap analysis gene expression (CAGE) method was performed using 75 nt single-end reads on a NextSeq 2000 instrument (Illumina, San Diego, CA, United States).Data from this experiment are deposited in Gene expression omnibus (GEO) database (accession number GSE229210).The quality of the obtained data was evaluated via the FastQC tool (v0.11.9).Reads alignment and mapping to a human reference genome (hg38) were performed using the STAR method (Dobin et al., 2013).Counting tags at each side was performed using SAMtools in R (Rsamtools v2.2.3).In each cell line, CAGE-Seq was performed in biological triplicate for the "WWOX" and "Vec" cellular variants.

Bootstrapping and transforming CAGE-seq data, weighted gene co-expression network analysis (WGCNA), identification and functional annotation of driver genes
Read counts were subjected to the resampling bootstrap method (Al Seesi et al., 2014;Kulesa et al., 2015) to acquire a relevant quantity of samples for Weighted gene co-expression network analysis (WGCNA) (Qu et al., 2021).The latter was initially performed to obtain the consensus modules across cell lines (temporarily excluding DBTRG-05MG due to technical issues preventing analysis) via Biological network reconstruction omnibus (BioNERO) (Almeida-Silva and Venancio, 2022) and further using the classical WGCNA approach (Langfelder and Horvath, 2008) to investigate all four cell lines independently.At the BioNERO stage, bootstrapped counts were subjected to variance stabilizing transformation (Love et al., 2014) via the vstransform parameter set to "TRUE" and filtered by variance using variance_filter set to "TRUE" with n = 10,000 during the automatic data preprocessing.Consensus modules identified during BioNERO workflow (via consensus_modules() and consensus_trait_ cor() functions) were further investigated using the classical WGCNA approach.The transformation into an adjacency matrix was performed with power = 9 selected in virtue of the pickSoftThreshold() function with signed hybrid network type and with the scale-free topology fitting index (R2) > 0.8.Following the topological overlap matrix preparation, the hierarchical clustering was performed through the hclust() function with an average agglomeration method.Established genes-containing modules were correlated to a binary trait representing stable transductants (samples were denoted as "WWOX" or "Vec" according to the procedure they were subjected to during stable transduction).Expression profiles were visualized using a heatmap.2()function of gplots v3.1.1 package with Euclidean distance metric for row ordering.Identification of genes most correlated with both module and trait (henceforth referred to as "driver genes" or "drivers") was facilitated via the verboseScatterplot() function with threshold ≥0.7 for both Gene significance (GS) and Module membership (MM).Expression of driver genes was also visualized using gplots v3.1.1 package.Subsequently, the intramodular connectivity of top genes was established via exportNetworkToVisANT() function of WGCNA with further manual adjustment using VisANT 5.0 (Granger et al., 2016) and Cytoscape 3.9.0(Shannon et al., 2003).The top 50 genes within a module for each cell line were obtained except DBTRG-05MG, for which the top 25 genes were acquired from the module part consistent with other cell lines and the top 25 from the module part with the unique expression profile.Subsequently, driver genes for each cell line were independently subjected to the Gene-set enrichment analysis (GSEA) using the R-package WebGestaltR (Liao et al., 2019) with a functional annotation database set to "gene ontology (GO) -biological process." The top 20 annotations were retrieved for each list of genes (i.e., drivers combined from two modules of interest described in section 3.1).A minimum number of three query genes for a specific annotation was set.Genes were ranked on the basis of their log 2 foldchange (log 2 FC) values acquired from the limma-voom, which is described in the subsequent section.Lastly, an intersection analysis between lists of driver genes for all cell lines was performed using the UpSetR package in the R environment (Conway et al., 2017); for this step, the threshold of GS and MM was set to ≥0.65.

Differential expression analysis (cell lines data), building and annotating the network
Bootstrapped counts were also subjected to Differential expression analysis (DEA) using the limma-voom method (Law et al., 2014;Ritchie et al., 2015).Preprocessing data via calcNormFactors() and filtering lowly expressed genes (requirement of minimum ≥5 counts per million in ≥1 library) were followed by transformation via the voom() function.The model was fitted in limma using weighted least squares for each gene via lmFit(), and log 2 FC values for "WWOX" (case) versus "Vec" (control) in each cell line were acquired following the execution of makeContrasts() function with default parameters.Empirical Bayesian smoothing of standard errors preceded the acquisition of the top-ranked genes via the topTable() function.At this stage, it was intended to compare differences between cell lines and not only between cellular variants within a specific cell line.Thus, the typically reasonable log 2 FC threshold of |0.57| was reduced to |0.37| to facilitate investigating the distinct expression profile of DBTRG-05MG versus the other three cell lines.In the middle of the study workflow, this approach was considered justifiable since it allowed to identify inversely regulated genes whose expression was, e.g., gently increased in the "WWOX" variant compared to the "Vec" DBTRG-05MG, but at the same time evidently downregulated in other cell lines.Genes identified through this approach served as an input for Cytoscapeintegrated NDEx iQuery (Pillich et al., 2023) that enabled the construction of a protein-protein interaction network using data from the BioGRID repository (Oughtred et al., 2021).Moreover, the input list contained the driver genes simultaneously identified in at least three cell lines, as well as drivers that exhibited inverse properties in DBTRG-05MG versus U87MG/T98G/U251MG.These genes were identified using methods described in the previous section.If drivers presenting inverse regulation in DBTRG-05MG were not a part of an automatically built network, they were added manually via the Advanced network analysis tool (ANAT) (Signorini et al., 2021).Afterwards, the "STRINGify network" option of Cytoscape's plugin "stringApp" (Doncheva et al., 2019) was run onto the entire network to expand it with additional proteins that demonstrated a very high interaction score (confidence of ≥0.9).After the EntOpt Layout was applied to maximize the clarity of the network, log 2 FC values for all cell lines were imported from a separate text file and visualized inside nodes via the Cytoscape's plugin "enhancedGraphics" (Morris et al., 2014).Finally, network was subjected to the Over-representation analysis (ORA) via WebGestaltR with top 20 annotations retrieved for each list of genes (see section 3.3 for details) and minimum number of three query genes for a specific annotation.

Patients' data acquisition, designating groups, and survival analysis
The remaining methodological parts were intended to certify the cell lines-related observations using data from glioblastoma patients.RNA-Seq quantification files generated using STAR workflow (Dobin et al., 2013) were acquired from Genomic Data Commons (GDC) for glioblastoma patients of The Cancer Genome Atlas (TCGA-GBM) using the GDCRNATools v1.18.0 package in the R environment v4.2.3.Together with WWOX, the two most important driver genes discerning DBTRG-05MG from other cell lines (i.e., POLE4 and HSF2BP; see section 3.4 for details) were considered a geneset of which expression signature could reflect the expression profile change that occurred in cells following stable lentiviral transduction.Independent groups of patients representing (1) "Vec" of U87MG/T98G/U251MG, (2) "WWOX" of U87MG/T98G/ U251MG, (3) "Vec" of DBTRG-05MG, and (4) "WWOX" of DBTRG-05MG were determined on the basis of WWOX/POLE4/HSF2BP expression denoted as "high" (↑) or "low" (↓); specific cut-off values summarized in Supplementary Table S2.
Together with designating cellular variants-resembling groups, the acquisition of corresponding clinical data for TCGA-GBM patients from GDC enabled the analysis of Overall survival (OS) and Diseasefree survival (DFS).Prior to focusing on designated groups, TCGA-GBM patients were subjected to DFS analysis via Evaluate Cutpoints tool (Ogluszka et al., 2019) based on RNA-Seq expression data separately for WWOX, POLE4, and HSF2BP.Subsequently, patients with desired WWOX/POLE4/HSF2BP expression were selected via which() function in the R environment and further subjected to survival analysis using the surfvit(), of which results were visualized using ggsurvplot().For each patient from the designated groups, the supplementary clinical information and data on molecular subtypes were acquired from the literature (Brennan et al., 2013;Wang et al., 2017;Zakharova et al., 2022), whereas the drug sensitivity prediction for Temozolomide (TMZ) measured by Half-maximal inhibitory concentration (IC50) was obtained from the CancerRxTissue repository (Li et al., 2021).These additional data were visualized using pROC and ggplot2 R-packages (Robin et al., 2011).

Differential expression analysis (patients' data) and principal component analysis
After the groups of patients were established, they were subjected to DEA performed via limma-voom as in section 2.5, except for the log 2 FC threshold set to |0.57| at the current stage.Two comparisons were considered.To resemble the WWOX/POLE4/HSF2BP expression changes that were corresponding in U87MG/T98G/ U251MG cell lines, patients representing the "WWOX" cellular variant (WWOX↑ POLE4↓ HSF2BP↓) were compared to those representing the "Vec" variant (WWOX↓ POLE4↑ HSF2BP↑).In contrast, reflecting the expression changes in the "WWOX" variant versus the "Vec" variant of DBTRG-05MG required patients with WWOX↑ POLE4↑ HSF2BP↑ or WWOX↓ POLE4↓ HSF2BP↓ profile, respectively.Expression of genes regulated inversely between two comparisons was visualized on heatmap using gplots v3.1.1 as mentioned in section 2.3; the same genes were also functionally annotated using ORA as in section 2.4.At the same stage, subsidiary analysis of immune infiltration and gene expression implicated in the myeloid-related transcriptional states was performed among patients subjected to DFS analysis.Infiltration estimates were obtained via TIMER 2.0 (Li et al., 2020); the expression profile of patients was uploaded to the "Estimation" module of TIMER 2.0 with the "GBM" option under the "cancer type." The list of genes implicated in the myeloid-related transcriptional states was acquired from the literature (Gangoso et al., 2021;Ravi et al., 2022;Rajendran et al., 2023).Gene expression and infiltration estimates were visualized on the violin plots using the ggpubr v0.6.0 package in R environment v4.2.3.Afterwards, the results of patients-related DEA were explored for expression changes that reflect the fold-change data from cell lines implemented in the interaction network from section 2.4.Expression of identified genes was used for dimensional grouping performed via the Principal component analysis (PCA).These genes were treated as active variables, of which contribution to spatial partitioning across the first two principal components was visualized via FactoMineR and factoextra R-packages (Lê et al., 2008).The biplot of both quantitative and qualitative variables was visualized using the fviz_ pca_biplot() function with a confidence ellipse drawn for each group of patients.Additionally, gene expression differences between groups of patients were investigated for statistical significance and visualized as violin plots using the vioplot R-package.

Statistical analysis
GraphPad Prism v8 (GraphPad Software, San Diego, CA, United States) was employed for statistical analysis.The normality of distribution was determined by the Shapiro-Wilk test.The unpaired t-test or the Mann-Whitney test was used depending on the data distribution.Results with a value of p less than 0.05 were considered statistically significant.

Results and discussion
3.1.Gene modules that were consistent in U87MG, T98G, and U251MG, were not entirely coherent in DBTRG-05MG The initial step of the current analysis was intended to identify a common expression profile for U87MG, T98G, and U251MG cell lines owing to their coherence in results from the previous study following WWOX overexpression (Kaluzinska-Kolat et al., 2023).DBTRG-05MG was temporarily excluded from establishing consensus module-trait relationships since it caused technical issues preventing analysis.Nonetheless, the input gene list to this step was determined after CAGE-Seq data filtering to retain genes whose expression varied between "Vec" and "WWOX" samples for each cell line independently, including DBTRG-05MG.Up to 10,000 most variable genes per cell line were acquired, with 5,744 variable genes present in all four lists of genes (Supplementary Table S3).This geneset served as the input to the BioNERO pipeline that revealed two consensus modules, i.e., green (159 genes) and pink (525 genes), that distinguish "Vec" and "WWOX" cellular variants of U87MG, T98G, and U251MG cell lines (Figure 1A).These modules were visualized separately for each cell line (Figures 1B-D), including DBTRG-05MG (Figure 1E), of which the expression profile was not entirely consistent with U87MG, T98G, and U251MG.Approximately 20% of the green module (32 genes) and about 30% of the pink module (172 genes) tended to be inversely regulated in DBTRG-05MG compared to other cell lines.The distinct expression profile of DBTRG-05MG might coincide with the gene mutation profile that is much more altered in this cell line than in others.According to the Catalogue of somatic mutations in cancer (COSMIC), mutation quantity is doubled in DBTRG-05MG relative to U251MG and U87MG or is 1.6 times higher than in T98G (Forbes et al., 2011).
3.2.Cell lines were characterized by unique driver genes, but some of the overlapping ones exhibited inverse properties in DBTRG-05MG Insights into the WGCNA revealed that there are dozens to over a hundred cell line-specific driver genes, i.e., genes most correlated with both module and trait (Figures 2A,B).Herein, the threshold was set to ≥0.7 for both GS and MM.Some DBTRG-specific drivers met these requirements within unique parts of two gene modules for this cell line (DUSP4 and ZNF786 for green module; CARHSP1, FGGY, FOSL2, HSF2BP, LINC01547, NCBP2AS2, POLE4, RUSC1, SLCO4A1, THAP11, and UBTD1 for pink module).A short synopsis for those having the literature data in the context of GBM will be provided at the end of this section.Complete lists of driver genes and related GSEA are summarized in Supplementary Table S4.In brief, functional annotation revealed that T98G-specific drivers regulate dendrite development or actin filaments, while U87MG-specific drivers control development growth or vesicle-mediated transport in the synapse.Moreover, U251MG-specific drivers alter cell polarity or cellular drug response, whereas DBTRG-specific drivers orchestrate proliferation or axon development.Investigating intramodular connectivity of green or pink module for each cell line indicated some genes being top hubs, i.e., DYM, FADD, GNPDA1, TNIP1, MITD1, RPS2, EYA4, TMBIM6, PUF60, and FOSL2 (all were driver genes except for TMBIM6 that presented MM ≥ 0.7 but GS < 0.7).For some of these genes, their role in GBM is documented.Overexpression of Fas-associated via death domain (FADD) may sensitize GBM to cell death (Marin-Rubio et al., 2019).TNFAIP3-interacting protein 1 (TNIP1) mediates a signaling cascade that sustains GBM proliferation; a high level of TNIP1 was correlated with poor survival (Lei et al., 2020).Elevated expression of Microtubule interacting and trafficking domain-containing 1 (MITD1) forecasted unfavorable prognosis in patients with low-grade glioma and GBM.For the latter, MITD1 was positively correlated to homologous recombination deficiency (Dong et al., 2022).The gene Eyes absent homolog 4 (EYA4) plays a putative tumor-promoting role in nervous system tumors (Chong et al., 2023) and was found to promote cell proliferation in GBM (Li et al., 2018b).Transmembrane BAX inhibitor motif containing 6 (TMBIM6) is elevated in GBM but may be downregulated in other tumor entities; it may be associated with tumor growth and dependent signaling (Yi et al., 2023).Knockdown of the Poly(U) binding splicing factor 60 (PUF60) reduced GBM tumorigenicity and proliferation (Wang et al., 2022).
Subsequently, intersection analysis of the driver genes for all cell lines was performed to assess whether there are common genes for U87MG, T98G, U251MG, and DBTRG-05MG.At this stage, GS and MM threshold was set at ≥0.65 to facilitate the overlap.Four drivers (two in green and two in pink module) were found in all cell lines and presented congruent expression change between "WWOX" and "Vec" cellular variants; these were PRKRA and PNMA2 (green module) or GPS1 and ARF3 (pink module).The high expression ratio of the Protein activator of interferon-induced protein kinase (PRKRA) and its presumably oncogenic antisense counterpart CHROMR was linked to favorable survival of GBM patients (Sirvinskas et al., 2023).PNMA family member 2 (PNMA2) was downregulated in pediatric GBM cell line SJ-GBM2 following irradiation and might positively regulate the apoptotic process (Alhajala et al., 2018).No precise data on the role of GPS1 and ARF3 in GBM were found, although the latter gene might potentiate G 1 /S cell cycle transition, as noted in breast cancer (Casalou et al., 2020).
Of driver genes found in at least three cell lines, only those of DBTRG-05MG were inversely regulated in comparison to U87MG, T98G, or U251MG (Figures 2C,D).This suggested that particular attention should be drawn to driver genes of DBTRG-05MG that exhibit inverse relation to the "WWOX" variant relative to "Vec, " as estimated by GS or MM, as well as log 2 FC from DEA performed in parallel (Figure 2E).Reducing the GS and MM threshold to ≥0.65 enabled the acquisition of three drivers (i.e., PLCD3, SH3BP2, and SLC36A1) in addition to those mentioned at the beginning of this section.Some of these drivers have documented function in glioblastoma.Dual-specificity phosphatase (DUSP) family regulates GBM sensitivity to treatment (Prabhakar et al., 2014); specifically for DUSP4, its overexpression in glioblastoma significantly reduced proliferation and colony formation (Waha et al., 2010).Knockdown of Calcium regulated heat stable protein 1 (CARHSP1) alleviated glioblastoma radioresistance, presumably via the inflammatory signaling pathway.When receiving radiotherapy, individuals with higher CARHSP1 levels had unfavorable survival (Zhu et al., 2021).Fos-like 2 (FOSL2) regulates angiogenesis and plasticity of GBM, is involved in tumor development, and regulates the so-called natural evolution signature that contains genes most differentially expressed between new and old GBM lesions (Marques et al., 2021;Wan et al., 2021;Eisenbarth and Wang, 2023b).Furthermore, it was hailed as one of six master regulators of mesenchymal gene expression signature (Fedele et al., 2019).Phospholipase C delta 3 (PLCD3) was correlated to diacylglycerol and is converted to phosphatidic acids by diacylglycerol kinase α, of which inhibition was recently considered a promising strategy for GBM (Ohanian and Ohanian, 2001;Wang et al., 2021;Purow, 2022).Lastly, although Solute carrier family 36 member 1 (SLC36A1) expression presented high variability among GBM individuals, it was higher than in control samples (i.e., human cortex) (Schaffenrath et al., 2021).Interestingly, this can bring benefits to the clinic since SLC36A1 has been utilized as a drug-delivery platform for therapeutics (Thwaites and Anderson, 2011).Frontiers in Neuroscience 07 frontiersin.org3.3.WWOX, drivers, and related genes formed a network that was of high interconnectivity and was related to various biological processes, including those associated with the nervous system Following WGCNA and DEA (described in sections 2.3 and 2.4, respectively), it was decided to build a WWOX-related network (Figure 3).Included genes were the drivers identified in at least three cell lines (GS and MM ≥0.65) but also the "background" genes from DEA (log 2 FC ≥ |0.37|) and from the STRING database (genes with a very high interaction score, i.e., the confidence of ≥0.9).Details on network construction are provided in section 2.4.Such an approach to elaborating the network has enabled the identification of other genes that were not considered drivers but still exhibited the inverse expression profile in DBTRG-05MG relative to U87MG, T98G, and U251MG.These genes were ACSF2, ARFRP1, CLU, ENG, GET4, IGF2R, OGFR, PDF, S100A4, WDTC1, and YPEL3 (Figure 3A).
Acyl-CoA synthetase family member 2 (ACSF2) may regulate neuronal development, and even though it did not affect anchoragedependent growth in GBM cells, it is related to ACSVL3 that supports glioblastoma stem cell maintenance and tumorigenicity (Sun et al., 2014;Tomi-Andrino et al., 2022).The association between ADP ribosylation factor-related protein 1 (ARFRP1) and a higher risk of GBM-type glioma requires further investigation (Song et al., 2012) since available data are inconsistent (Atkins et al., 2019;Ali et al., 2021).Clusterin (CLU) levels were found elevated in glioblastoma stem cells which induced an anti-apoptotic state (Osuka et al., 2021); however, it is worthy of note that CLU can also promote apoptosis when localized in the nucleus (Kim and Choi, 2011) (the first mentioned reference investigated a secreted protein form).Endoglin (ENG) is an angiogenic biomarker determining the microvessel density of GBM (Afshar Moghaddam et al., 2015).Although its prognostic significance is inconsistent, it is overexpressed in actively proliferating endothelial cells and supports new vascular networks (Bastos et al., 2022).Insulin-like growth factor 2 receptor (IGF2R) is Investigation of driver genes for all cell lines with an emphasis on drivers that exhibit inverse properties in other cell lines relative to DBTRG-05MG.Separately for U87MG, T98G, U251MG, and DBTRG-05MG, the driver genes from (A) green module or (B) pink module were identified via WGCNA using GS and MM threshold of ≥0.7.The part of the green or pink module that was consistent in DBTRG-05MG and other cell lines was presented separately from the part with an inverse profile found exclusively in DBTRG-05MG.Intramodular connectivity for top genes was additionally included.Overlap of driver genes between cell lines was presented separately for (C) green and (D) pink modules.(E) Some driver genes in DBTRG-05MG were characterized by inverse GS and log 2 FC relative to other cell lines.For example, POLE4 was identified in the pink module, and the expression of this gene is higher in "Vec" variants of U87MG, T98G, and U251MG.In other words, POLE4 is negatively correlated to "WWOX" variants of these cell lines since the negative GS and log 2 FC are present.Regarding DBTRG-05MG, POLE4 expression was higher in the "WWOX" variant relative to "Vec"; thus, positive GS and log 2 FC values are visible.

FIGURE 3
The WWOX-related network containing genes from DEA, drivers from WGCNA, and top interactors from STRING.(A) EntOpt Layout was adopted to maximize the clarity of the network.The log 2 FC values ("WWOX" versus "Vec" comparison) for all cell lines were included to facilitate the interpretation; overexpressed in GBM samples compared with normal brain specimens, but no association between its expression and patients' outcomes was observed (Maris et al., 2015).The same authors underlined that the reports on the entire IGF family in glioblastoma are typically restricted to small series and have not yielded consistent results.Targeting the mitochondrial Peptide deformylase (PDF) was suggested as a promising approach for sensitization of GBM to chemotherapy; data were obtained from research that investigated an inhibitor of PDF, i.e., actinonin, that led to the activation of mitochondrial unfolded protein response, increased mitochondrial fission, as well as promoted the integrated stress response to increase apoptosis (Lan et al., 2022).S100 calcium-binding protein A4 (S100A4) was found to be necessary for maintaining glioblastoma stem cells self-renewal; moreover, ablation of S100A4-expressing cancer cells is sufficient to impede glioma growth in vivo (Chow et al., 2017).
From another point of view, the elaborated network can be grouped on the basis of fold-change values established between "WWOX" and "Vec" variants in U87MG, T98G, U251MG, or DBTRG-05MG (Figure 3B).Although most of the included genes had consistent profiles in all cell lines, the observations unique for DBTRG-05MG constituted more than 20% of the entire graph (26 genes out of 113).Functional annotation revealed that misfolded protein and cell death were similarly affected in all cell lines following WWOX overexpression.On the contrary, WWOX overexpression regulated microglia proliferation, neurofibrillary tangle assembly, and nervous system vasculogenesis differently in DBTRG-05MG than in other cell lines.A more detailed summary of gene ontology is collected in Supplementary Table 5.Some annotated processes can be certified using data from our previous research (Kaluzinska-Kolat et al., 2023).ORA of the entire network and the parts consistent for all cell lines has repeatedly indicated cell death and apoptotic signaling pathway regulation.Indeed, all four "WWOX" cellular variants were characterized by decreased cell viability, alongside the intensified apoptosis of T98G/DBTRG-05MG/U87MG cells following WWOX overexpression.Results for U251MG were statistically insignificant, although a consistent trend was visible (Kaluzinska-Kolat et al., 2023).Part of the network specific for DBTRG-05MG suggested that vascular tube formation might differ this cell line from the others; however, there was no in vitro assay investigating this process in our former study.Nevertheless, annotations indicating a unique regulation of proliferation align with previous data -proliferative potential after WWOX overexpression was increased only in DBTRG-05MG (Kaluzinska-Kolat et al., 2023).Annotated processes in each group of the network seem to be related; for example, the association between neural crest cells and vasculogenesis is known (Bergwerff et al., 1998;Wiszniak et al., 2015).Moreover, microglia cells are observed around Neurofibrillary tangles (NFTs), i.e., aggregates of hyperphosphorylated tau protein that are neuropathological indicators of Alzheimer's disease (AD) (Metaxas and Kempf, 2016).Although microglia cells can ingest NFTs, they do so inefficiently and may release tau aggregates into the environment, propagating disease pathology (Hopp et al., 2018;Perea et al., 2018).The functionality of tau protein in brain pathologies is not only limited to AD; its role in GBM is also documented (Pagano et al., 2021;Hedna et al., 2022).

Stratification by expression of WWOX and the two most relevant drivers explaining DBTRG-05MG otherness (POLE4 and HSF2BP) showed that patients differ in survival
Analysis was directed toward verifying data from cell lines using GBM patients.However, there was a need to reflect the dissimilarity observed in DBTRG-05MG relative to U87MG, T98G, and U251MG, the last three abbreviated as "UTU" from the first letters of cell line names.Thus, the two most relevant driver genes of DBTRG-05MG discerning this cell line from others were selected alongside WWOX to stratify patients and elaborate the groups resembling "Vec" and "WWOX" cellular variants.The most promising drivers from Figure 2E appeared to be POLE4 and HSF2BP, of which GS, MM, and log 2 FC values distinguished DBTRG-05MG from other cell lines.GS and MM values for all DBTRG-discerning drivers were visualized on a scatterplot emphasizing two selected genes; inverse fold-changes between DBTRG-05MG versus others were included (Figure 4A).After designating groups (see section 2.5 for details), patients resembling "WWOX" variants of "UTU" cell lines were compared to those representing "Vec" variants in the same cell lines.Likewise, individuals similar to the "WWOX" variant of DBTRG-05MG were compared to patients showing a "Vec"-like profile in this cell line.Only ~30% of patients from the TCGA-GBM cohort had an adequate expression profile of WWOX, POLE4, and HSF2BP.Survival analysis revealed that patients resembling the profile of the "WWOX" cellular variant in "UTU" cell lines had longer OS and DFS than those resembling the "Vec" variant.On the other hand, contrary results were observed for individuals representing the cellular variants of DBTRG-05MG.Overall, WWOX appears to be less significant in the presence of POLE4 and HSF2BP, although subsequent studies on larger datasets are required to evaluate their relation thoroughly.When analyzing each gene individually, DFS data from the TCGA-GBM cohort indicate a favorable effect of higher WWOX expression, whereas higher POLE4 or HSF2BP expression is somewhat unfavorable (Supplementary Figure S1).In further steps, we focused on patients representing DFS outcomes since the events caused by disease recurrence occur earlier than death (Birgisson et al., 2011;Han et al., 2014).All four groups of GBM patients were subjected to DEA (log 2 FC ≥ |0.57|), which revealed more than 400 genes that were inversely regulated among two comparisons (Supplementary Figure S2).Functional annotation of these genes suggested that patients differ in, e.g., response to axon injury, locomotory behavior, synapse organization, and myeloid cell differentiation (Supplementary Figure S2).Not only the relationship between neurodevelopmental and injury programs has recently emerged as a see the top-right corner for the legend.(B) The same network was visualized using Grid Layout to group nodes according to their fold-change.Next to each group, a number representing the quantity of nodes is displayed.Groups were independently subjected to gene ontology analysis; two annotations per group are provided in dashed rectangles.The top 20 annotations per group are summarized in Supplementary Table S5.crucial factor influencing treatment sensitivity, but also the injury responses may act as supplementary drivers of phenotypic heterogeneity (Brooks et al., 2022).Pathways associated with locomotory behavior were highly expressed in high-risk GBM patients (Wang et al., 2016).Synapse organization is related to a so-called synapse score showing the prominent ability to predict glioma prognosis, diagnosis, and grading (Ji et al., 2020).Various types of connectivity in glioma (e.g., neuron-to-glioblastoma or glioblastomato-glioblastoma) facilitate membrane depolarization followed by calcium transients in networks of tumor cells (Monje, 2020).The Justifying the selection of POLE4 and HSF2BP as the two most relevant drivers explaining the otherness of DBTRG-05MG.(A) Scatterplot presenting GS and MM of all DBTRG-discerning driver genes that were summarized in Figure 2E.POLE4 and HSF2BP are annotated to emphasize the relevance of their GS and MM values accumulating in the top-right corner for all cell lines.For this step, absolute GS values were used to facilitate the visualization of data from all cell lines at once.The log 2 FC values are also provided to illustrate the dissimilarity between the "WWOX" versus "Vec" comparison in DBTRG-05MG.Alongside WWOX, the expression of POLE4 and HSF2BP was used to stratify GBM patients in groups resembling cellular variants of U87MG, T98G, and U251MG cell lines (abbreviated as "UTU" from the first letters of their names), as well as separately the variants of DBTRG-05MG (abbreviated as "DBTRG").Data for patients presenting the differences in (B) OS and (C) DFS are summarized in Supplementary Table S2.10.3389/fnins.2023.1260409Frontiers in Neuroscience 11 frontiersin.orginclusion of myeloid cell differentiation in functional analysis suggests that groups of patients may differ in the context of myeloid-derived suppressor cells (MDSCs) that accumulate in late-stage glioblastoma and repress immune activation (Eisenbarth and Wang, 2023a).MDSCs-related annotations in gene ontology motivated us to perform subsidiary analysis.Although more sophisticated data are required to evaluate MDSCs level among established groups of patients, the estimations obtained via TIMER 2.0 (Supplementary Figure S3) suggested that the immune infiltration by natural killer (NK) cells and myeloid dendritic cells (DCs) was reduced in "Vec"-like "UTU" and "WWOX"-like "DBTRG" groups that include patients having unfavorable survival compared to other two groups.MDSCs are known to hinder the activity of NK cells, the latter being utilized in GBM immunotherapy owing to its ability to counteract tumor progression via cancer cell lysis (Bayik et al., 2020;Breznik et al., 2022;Hosseinalizadeh et al., 2022).Therefore, it is unsurprising that MDSCs are targeted to enhance NK-based treatment (Joshi and Sharabi, 2022).Regarding DCs, their reduction among cancer patients may be due to the preferential differentiation of MDSCs since these two cell types share a common progenitor (Ostrand-Rosenberg et al., 2012).Being a crucial cell type responsible for antigen processing and presenting, DCs initiate anti-tumor immune responses and are used as vaccines in cancer therapy (Srivastava et al., 2019;Zhou et al., 2023).Other immune components with various infiltration levels among established patient groups include the cancer-associated fibroblasts and CD4 + Th2 T-cells, summarized in Supplementary Figure S3 (detailed TIMER 2.0 estimations available in Supplementary Table S6).Additionally, genes implicated in the myeloid-associated transcriptional states or accompanying phenomena were collected from the literature (Gangoso et al., 2021;Ravi et al., 2022;Rajendran et al., 2023) and investigated for expression changes between cell lines-resembling groups.
Among differentially-expressed genes (Supplementary Figure S3), it is worth mentioning that APOC1, CCL2, CCL3, IL10, IRF8, NF1, and SPP1 significantly distinguished "UTU" groups of patients ("DBTRG" groups had various CCL3 expression and might differ in NF1 level).Except for NF1, a tumor suppressor gene that diminishes GBM angiogenesis and migration, the remaining genes were elevated in "Vec"-like "UTU" patients.Loss of NF1 was found to be correlated with increased immune infiltration within the tumor microenvironment and may prime cells for malignant transformation (Gangoso et al., 2021).CCL2 was increased in mesenchymal GBM and recruits MDSCs (Chang et al., 2016;Rajendran et al., 2023).Another member from the same chemokine family, CCL3, is probably responsible for chemotaxis and is expressed in activated microglia cells that stimulate glioma progression and development (Geribaldi-Doldan et al., 2020;Ravi et al., 2022;Rajendran et al., 2023).Furthermore, IRF8 contributes to immune evasion (Gangoso et al., 2021), APOC1 promotes GBM tumorigenesis and is involved in lipid metabolism or immunosuppression (Zheng et al., 2022;Rajendran et al., 2023), whereas SPP1 is a marker of glioma-associated macrophages that ensures immunosuppression and worsens survival (Rajendran et al., 2023).Lastly, IL10 is a driver force of tumor immune escape; this interleukin is released by a subset of myeloid cells localized in mesenchymal-like tumor regions (Ravi et al., 2022).Based on the presented subsidiary analysis of immune cellular components and myeloid-associated genes, it can be suggested that groups of patients having unfavorable DFS, especially "Vec"-like "UTU" patients, are affected by more immunosuppressive GBM tumors.Evading immunosurveillance or immunotherapy resistance are major obstacles to GBM treatment (Liu et al., 2021), and dynamic communication with the immune tumor microenvironment is inherent during GBM progression (Arrieta et al., 2023;Guo and Wang, 2023).
3.5.Some expression profile differences in the network established from cell lines data were confirmed in patients, revealing genes whose function can explain diverse clinical outcomes Ultimately, it was decided to utilize the results of differential expression analysis performed on patients and identify which genes demonstrate the cell lines-resembling expression profile change.Alongside WWOX/POLE4/HSF2BP, twelve genes from the network exhibited the expression profile that is consistent between patients and cell lines.These genes were DUSP4 (Prabhakar et al., 2014), SLC36A1 (Schaffenrath et al., 2021;Wei et al., 2022), EPB41L1 (Han et al., 2019), PNMA2 (Alhajala et al., 2018), YPEL3 (Phoa et al., 2020), SEC11C (Pernia et al., 2020), FGGY (Salzillo et al., 2016;Smith et al., 2021), LSM2 (Sun et al., 2022), IGFBP7 (Jiang et al., 2008), STOML2 (Qu et al., 2019;Ma et al., 2021), PKM (Mukherjee et al., 2013;El Atat et al., 2023), and SYNPO (Ning et al., 2021;Krishna et al., 2023).Some of them were shortly described in previous sections, but the possible roles of all genes are recapitulated in Figure 5A, where the results of PCA are also presented.At first glance, a separate cluster for each group of patients was visible, suggesting that individuals are characterized by distinct expression profile that appears to contribute to patients' survival.Namely, patients from both comparisons with more favorable DFS occupied the left side of the graph (red and purple groups are in the area with negative values on Dim1, i.e., on the left relative to Dim2).In contrast, individuals with shorter DFS were located on the right side of the graph (green and blue groups are in the area with positive values on Dim1, i.e., on the right relative to Dim2).To facilitate interpretation, groups were marked with the same colors as in Figure 4C.Moreover, patients occupying the top and bottom sides of the graph differed regarding WWOX expression (the top side corresponds to blue and red groups, whereas the bottom side corresponds to purple and green groups).Once these groups of patients are considered as "Vec"-like and "WWOX"-like equivalents of cellular variants, the difference between "Vec" and "WWOX" for "UTU"-resembling patients is directed diagonally to the left from top to bottom.In contrast, the difference between "Vec" and "WWOX" for DBTRG-resembling patients is directed diagonally to the right from top to bottom.This suggests that contrary to "UTU" cell lines, WWOX overexpression in DBTRG-05MG resulted in distinct expression profile change, resembling that of patients from WWOX↑ POLE4↑ HSF2BP↑ group characterized by inferior outcome presumably due to the events that support cancer development.The latter can be supported by the function of genes included in PCA and their contribution to grouping patients in the first two principal components.Differences in the expression of these genes are also summarized in Supplementary Figure S4 for "UTU"resembling patients and Supplementary Figure S5 for individuals resembling DBTRG-05MG cellular variants.Lastly, supplementary clinical and molecular data (Cerami et al., 2012;Brennan et al., 2013;10.3389/fnins.2023.1260409Frontiers in Neuroscience 12 frontiersin.orgWang et al., 2017;Li et al., 2021;Zakharova et al., 2022) have indicated that groups of patients occupying the right side of the graph (green and blue groups) were characterized by less prevalent Methylguanine methyltransferase (MGMT) promoter methylation or more frequent mesenchymal subtype of GBM (Figure 5B).It may explain shorter DFS in these groups since unmethylated MGMT promoter predicts an unfavorable outcome (Riemenschneider et al., 2010).Likewise, mesenchymal glioblastoma is the most malignant disease subtype with a signature that forecasts worse survival (Gangoso et al., 2021;Kim et al., 2021).Moreover, higher TMZ IC50 values were noted in the above-mentioned groups in comparison to the other two groups (Figure 5C), suggesting that patients occupying the right side of the graph (green and blue groups) are worse responders to chemotherapy.The groups of patients were not different regarding the Isocitrate dehydrogenase (IDH) status (only two samples were IDH-mutant) or chromosome 1p/19q co-deletion.However, there were some differences in copy number alterations (CNAs) or mutational status among the most frequently affected genes (Supplementary Table S7; Figure 5D).For example, patients characterized by worse outcomes had more frequent MUC16 mutation and the presence of NF1 Investigating the WWOX-dependent expression profile observed in cell lines using data from patients.(A) Principal component analysis indicated a separate cluster for each cell lines-resembling group of patients (see top-right corner for the legend).Incorporating data from survival analysis (groups are colored as in Figure 4C) revealed that patients with more favorable DFS occupied the left side of the graph (red and purple groups), whereas individuals with shorter DFS were located on the right side of the graph (green and blue groups).Moreover, patients occupying the top and bottom sides of the graph differed in WWOX expression (the top side corresponds to blue and red groups, whereas the bottom side corresponds to purple and green groups).(B) Supplementary data were provided for each group, indicating that, e.g., groups with longer DFS were characterized by more prevalent MGMT promoter methylation.This is complemented with (C) predicted temozolomide IC50 values acquired from the CancerRxTissue repository that indicate a statistically significant lower TMZ IC50 in groups with longer DFS, suggesting that these patients are better responders to chemotherapy than those with unfavorable DFS.(D) Mutation and copy number alteration (CNA) data revealed that patients characterized by worse outcomes had more frequent MUC16 and NF1 mutation but less frequent EGFR amplification and mutation.Unk, unknown; Mut, mutated; WT, wildtype. 10.3389/fnins.2023.1260409 Frontiers in Neuroscience 13 frontiersin.orgmutation.Alterations of MUC16 in GBM are associated with unfavorable prognosis (Ferrer, 2023), whereas the NF1 mutation is a marker of treatment-resistant gliomas (Basindwah et al., 2022).Even though CDKN2A/B deletion predicts worse survival (Hsu et al., 2022), CNA data among PCA-included patients are inconclusive.However, EGFR amplification and mutation are less frequent in green and blue groups, possibly contributing to their inferior outcomes (Higa et al., 2023) even though some literature data on its prognostic value are debatable (Hobbs et al., 2012;Li et al., 2018a;Alnahhas et al., 2021) or indicate the lack of influence on survival (Ma et al., 2020).The role of POLE4 or HSF2BP in glioblastoma is poorly understood.Even though some CpG sites functional in high-grade IDH wildtype glioma are located in the genomic location of POLE4 (Kessler et al., 2020), this gene was only mentioned as being implicated in nuclear excision repair (NER) (Bady et al., 2018).Interestingly, NER is a process by which a wide variety of structurally unrelated DNA lesions generated by chemotherapy and radiotherapy can be removed (Rominiyi and Collis, 2022), suggesting that POLE4 might control the treatment resistance of GBM.Outside the glioblastoma context, POLE4 was found to regulate proliferation, migration, and invasion of non-small cell lung cancer (Wang F. et al., 2021).The only GBM-related study mentioning HSF2BP investigated hypoxia signatures since the hypoxic interiors of glioblastoma hamper chemotherapy and drive angiogenesis, cell death inhibition, or tumor immune escape (Hu et al., 2012;Lathia et al., 2015;Yaghi et al., 2016;Monteiro et al., 2017).In the mentioned study, HSF2BP was upregulated in hypoxia-induced GBM, but no follow-up was performed for this gene (Bhushan et al., 2021).Nevertheless, a study on lung adenocarcinoma indicated that HSF2BP forecasts unfavorable patients' survival and promotes genomic instability (Huang et al., 2021).It is worth noting that such instability may be driven by tumor hypoxia (Tang et al., 2021), which could be considered a valuable hint for future research.
Summarizing all results from this section, it appears that even such a narrow set of genes is able to present differences among cell linesresembling groups of patients that conforms to discrepancies between DBTRG-05MG and other cell lines found in our previous study.Out of several biological processes investigated in vitro, changes in DBTRG-05MG concerned proliferation and tumor growth that were accelerated following WWOX overexpression (Kaluzinska-Kolat et al., 2023).Some of the PCA-included genes (DUSP4, EPB41L1, IGFBP7) orchestrate these processes in GBM (Jiang et al., 2008;Prabhakar et al., 2014;Han et al., 2019); it remains to be elucidated whether POLE4 also regulates glioblastoma proliferation as it does in lung cancer (Wang F. et al., 2021).Nonetheless, it should not be forgotten that the elaborated resemblance of findings from GBM cell lines using data from GBM patients is imperfect, and it cannot be assumed that POLE4 and HSF2BP levels were elevated in patients with unfavorable outcomes due to the higher levels of WWOX.However, it highlights that a specific group of patients (6-7% of the entire cohort) cannot benefit from relatively high expression of this tumor suppressor.A possible explanation sheds light on the accessibility of WWOX for its protein partners, which seems crucial for the manifestation of pleiotropic functions.Available literature data indicate that two cancer-promoting proteins, TMEM207 and VOPP1, can restrain the antineoplastic activity of WWOX by inhibiting its ability to associate with some partners (Bonin et al., 2018;Bunai et al., 2018).Taouis et al. admitted that this could explain the existence of cancerous cells expressing high levels of WWOX (Taouis et al., 2021).
Perhaps a similar situation is present in patients with WWOX↑ POLE4↑ HSF2BP↑ profile.Even though POLE4 and HSF2BP were not confirmed to be direct WWOX interactors, their overexpression might upregulate other proteins that hinder the accessibility of WWOX for its cooperators, which is crucial for the maintenance of anti-cancer activity.Further studies are highly encouraged, which should add new dimensions to our current understanding of WWOX multimodality.

Study limitations
Similar to other research, our study has some limitations.For instance, the conclusions drawn from in vitro experiments may not accurately reflect the conditions within an organism.It is advisable to conduct a further examination, especially since the data from TCGA-GBM patients included in the current study are encouraging.Moreover, our in vitro findings are based on "classical" commercially available cell lines instead of, e.g., glioma stem cells.Although cell lines established and maintained in serum-free neural stem cell media are increasingly used as experimental glioma models (Allen et al., 2016), they resemble but do not mimic original tumors (Garcia-Romero et al., 2016).Compared to conventional cell lines, they are still superior in retaining tumor-specific phenotypes but at the cost of throughput and overall complexity or standardization (Pollard et al., 2009;Boccellato and Rehm, 2022).Given the execution of CAGE-Seq in biological triplicate for each cellular variant, it was inevitable in this study to use more conventional models; however, this can facilitate the replication of the study due to the availability and inexpensiveness.It is also worth mentioning that one of the experimental models utilized in our research -U87MG -can be considered a "problematic" cell line.Historically, the original U87MG cell line was established in 1968 at the University of Uppsala (Ponten and Macintyre, 1968) and is currently denoted as "U87MG Uppsala" according to the Cellosaurus repository (accession number: CVCL_GP63).Herein, we used the cell line that is distributed by most biobanks such as American Type Culture Collection (ATCC) or ECACC (denoted as "U87MG" or "U87MG ATCC"; Cellosaurus accession number: CVCL_0022), which is a different cell line of unknown patient origin but with confirmed central nervous system origin (Allen et al., 2016).It has also been stated that studies utilizing this cell line still reflect brain cancer biology, and there is no need to relegate them (Dolgin, 2016), especially since U87MG has been used in thousands of publications, contributing to the scientific advancement in glioblastoma research (Pokorna et al., 2021).Focusing on discrepancies between DBTRG-05MG versus the other three cell lines can also be considered a limitation since there is a possibility that some valuable results from other comparison schemes were omitted.However, the findings from our previous research revealed that such a scheme is the most justified (Kaluzinska-Kolat et al., 2023).Another aspect is related to the TCGA-GBM cohort that, despite being one of the most considerable resources of glioblastoma patients' data, does not provide transcriptomic profiling for all individuals.This has complicated the stratification based on the WWOX/POLE4/HSF2BP level and limited the number of genes that were used to establish expression profiles.Lastly, provided explanations for some results (such as those related to immune cellular components and myeloidassociated genes) may require further verification using dedicated and expensive methods that are currently out of our reach but which will 10.3389/fnins.2023.1260409Frontiers in Neuroscience 14 frontiersin.orgbe undertaken in the future.Still, these explanations may aid the scientific community in choosing subsequent research directions.

Conclusion
The following inferences can be drawn from consecutive parts of this study.The expression profile of DBTRG-05MG is distinct from U87MG, T98G, and U251MG; some of the DBTRG-05MG driver genes are unique for this cell line (FGGY, NCBP2AS2, UBTD1) or may be present in other cell lines but are regulated inversely following WWOX overexpression (CARHSP1, DUSP4, FOSL2, HSF2BP, LINC01547, PLCD3, POLE4, RUSC1, SH3BP2, SLC36A1, SLCO4A1, THAP11, ZNF786).Drivers, top hubs, and other WWOX-dependent genes identified in this study regulate various processes, including those important for the nervous system.Genes that were not previously the subject of GBM-related research (e.g., FGGY, LINC01547, NCBP2AS2, DYM, GNPDA1, GPS1, GET4) warrant further investigation.Some driver genes with consistent expression change following WWOX overexpression were found in all cell lines (PRKRA, PNMA2, GPS1, ARF3), suggesting that a part of dependent processes might be regulated similarly in U87MG, T98G, U251MG, and DBTRG-05MG, which was the case in our previous study.Constructing the network has enabled the extension outside the drivers-related context, which showed other genes that may be inversely regulated in DBTRG-05MG relative to other cell lines.When split into relevant parts, functional annotation of the network suggested that processes such as apoptosis may be controlled congruously between cell lines, whereas microglia proliferation or neurofibrillary tangle assembly might distinguish DBTRG-05MG.Subsequently, it was possible to identify patients resembling cell line comparisons using WWOX, POLE4, and HSF2BP data, whose expression in various combinations indicated, e.g., distinct survival (OS and DFS), sensitivity to TMZ, or infiltration of immune cellular components.Some gene expression differences between cell lines were reflected in patients; genes manifesting equivalent profiles were DUSP4, SLC36A1, EPB41L1, PNMA2, YPEL3, SEC11C, FGGY, LSM2, IGFBP7, STOML2, PKM, and SYNPO.Their contribution to biological processes seems to clarify the cancer-related events that determined longer or shorter survival.The role of POLE4 or HSF2BP in GBM is poorly understood but may be related to nuclear excision repair that renders treatment resistance or might affect genomic instability, presumably in hypoxic interiors.It has been proposed that in patients with WWOX↑ POLE4↑ HSF2BP↑ expression profile (who account for ~6-7% of the entire cohort), WWOX protein may not be available for its cooperating partners, which is manifested by the absence of anti-cancer activity despite a high WWOX level.Such a circumstance was previously found in the literature but not in the case of glioblastoma.Whether POLE4 and HSF2BP bind WWOX directly is yet to be elucidated in future research, where particular attention should be drawn to their role and interaction networks in GBM.
All things considered, data from cell lines enabled the identification of the group of patients among which a high WWOX expression with simultaneously elevated expression of other genes (herein: POLE4 and HSF2BP) did not bring prognostic benefits and may be related to the more cancer-promoting profile.This group accounts for about 6-7% of the investigated cohort, which translates into a considerable percentage on a global scale and may complement the qualification of patients for personalized medicine in the future.

FIGURE 1
FIGURE 1Expression profile of green and pink consensus modules in all cell lines included in the study.(A) The consensus module-trait relationship calculated via BioNERO/WGCNA indicated that green and pink consensus modules were congruous in (B) U87MG, (C) T98G, or (D) U251MG.In contrast, the same modules were not entirely coherent in (E) DBTRG-05MG, for which an inverse profile was noticed in a portion of each consensus module.p < 0.05 (*), p < 0.001 (***).