RNA Sequencing Identifies Novel Signaling Pathways and Potential Drug Target Genes Induced by FOSL1 in Glioma Progression and Stemness

Introduction

FOSL1 (also known as FOS-related antigen 1, Fra-1) is a component of AP-1 transcription factor complexes, which are formed by the dimerization of JUN family members (c-Jun, JunB, and JunD) with FOS family members (c-Fos, FosB, Fra-1, and Fra-2).¹ FOSL1 is overexpressed in most solid tumors^2–9 including GBM,^10,11 and its overexpression is correlated with tumor progression and worsened patient survival.^12–15 Thus, FOSL1 has been a popular focus for therapeutic intervention of tumors.¹⁶ The FOSL1 protein is mostly located in the nucleus and is involved in various cellular functions,¹⁷ such as in cancer cell proliferation, invasion/metastasis,¹⁸ epithelial-to-mesenchymal transition (EMT),^7,19,20 cancer stemness,^21,22 chemo-resistance,²¹ as well as modulation of tumor immunity.²³ Notably, FOSL1 is implicated in promoting the transition of proneural (PN) to mesenchymal (MES) states in GBM.^11,24 The regulation of FOSL1 is governed by many mechanisms, including transcriptional control by STAT3 and posttranslational changes including phosphorylation and deacetylation. These findings have been found not only in our research^25,26 but also in other research groups.^17,22,25 Phosphorylation of FOSL1 at the C-terminal residues T223, T230, S252, and S265 increases its stability and prevents proteasome degradation independent of ubiquitin. In contrast, phosphorylation at T217, T222, and T231 enhances FOSL1’s transactivation.¹ Notably, phosphorylation at T231 is essential for FOSL1 transactivation by 12-O-tetradecanoylphorbol-13-acetate (TPA) via extracellular signal-regulated kinase (ERK) in tumor cells.^27,28 Deacetylation at K116 in the FOSL1 DNA-binding domain (DBD) by histone deacetylase (HDAC) like HDAC6 enhances DNA binding, thereby increasing its transactivation and promoting colorectal cancer (CRC) progression and increases CRC stemness.^1,22

The PN subtype transition to MES (PMT) of glioma stem cells (GSCs) confers resistance to chemo-radiation therapy in glioblastoma (GBM) patients and is strongly associated with poor patient outcomes;^29–32 FOSL1 is a crucial regulator of the transition from PN to MES, associated with glioma development.^24,32 FOSL1 is essential for preserving the stem cell characteristics of the MES subtype of GBM, which makes it a promising candidate for therapeutic targeting. While it may be difficult to pharmacologically suppress FOSL1 due to its cancer-causing properties, targeted therapy by methods like CRISPR/Cas9 or PROTAC that specifically target the abnormal expression of FOSL1 could be a potentially effective strategy for treating patients with mesenchymal GBM.^24,33 Our findings indicate that FOSL1 is linked to the advancement of cancer, specifically in increasing the growth and spread of glioma, as well as enhancing the stem cell-like properties of glioma cells,^18,25 but the specific mechanisms behind FOSL1-induced PMT and treatment resistance have not been fully understood. Our observations indicated that the PN subtype PDX-L14/X456 exhibited a markedly decreased level of FOSL1 in comparison to the other glioma cell lines (see Figure 4D in reference).¹⁸ To elucidate the role of FOSL1 in the progression and development of GBM, we injected FOSL1 into PDX-L14 cells, which exhibit limited FOSL1 expression, to enhance the expression of this gene. Subsequently, we conducted RNA-seq and transcriptome studies. In this study, we reported crucial gene profiles and signaling pathways that are induced by FOSL1 in glioblastoma. Our analytical workflow is briefly outlined in Figure 1.

Figure 1 Study Flowchart.

Material and Methods

Cell Culture, Plasmids, and Transfection

Human glioblastoma cell lines: A172 (RRID: CVCL_0131) and U87MG (HTB-14, RRID: CVCL_0022) were obtained from ATCC. The cells were cultured in DMEM enriched with 10% FBS, 50 units/mL penicillin, and 50 µg/mL streptomycin at 37°C. The PDX-L14 line, a patient-derived xenograft (PDX), was generated by implanting PDX-tumor tissue cubes subcutaneously into the flanks of male or female 6–8 weeks old nude mice under anesthesia (ketamine/Xylazin 90/6 mg/kg BW) in accordance with the approach detailed in our prior article.²⁵ IACUC (Institutional Animal Care and Use Committee) NO. 23–18 (20–14 renewal). Following this, a suspension of individual cells was cultured in DMEM/F-12 media enriched with 10% FBS, 50 units/mL penicillin, and 50 µg/mL streptomycin. All investigations utilized cells devoid of mycoplasma. The FOSL1 human GFP-tagged ORF clone (RG202104) and the control vector PS100001 were acquired from OriGene.

Transfection of DNA Constructs

Transfection of DNA constructs was carried out when glioma cells reached approximately 50–75% confluency. FOSL1 tagged with GFP, along with corresponding controls, were transfected using lipofectamine 3000 transfection reagent (Invitrogen, Carlsbad, CA), according to the manufacturer’s instructions. The extent of overexpression was assessed 48 or 72 hours post-transfection via Western blotting or immunofluorescence imaging, if required. All studies were done in triplicates.

RNA Isolation, cDNA Synthesis, and qPCR

The procedures for isolating total RNA, synthesizing cDNA, and amplifying by PCR were conducted according to the previously described methods.³⁴ The cells were subjected to RNA extraction using TRIzol, and the resulting RNA was quantified using a Nanodrop. The iScript cDNA Synthesis Kit was employed to reverse transcribe 0.75 μg of purified RNA. The process of reverse transcription was conducted using random hexamers at a temperature of 25°C for a duration of 5 minutes, followed by a temperature of 42°C for 30 minutes, and finally a temperature of 85°C for 5 minutes. The cDNA was diluted by a factor of ten and subjected to amplification using iQ SYBR Green Supermix (Bio-Rad Laboratories, Inc) following the manufacturer’s instructions. The amplification was carried out under the following conditions: an initial denaturation step at 95°C for 3 minutes, followed by 40 cycles of denaturation at 95°C for 10 seconds, and annealing/extension at 61°C for 45 seconds. The CFX Connect Real-Time PCR Detection System was employed for quantitative gene analysis. Each experiment was conducted in three biological replicates, and GAPDH was utilized as the reference gene.³⁴ The relative expression for each gene was determined using the 2^−ΔΔCt technique. The primers for FOSL1, ITGA5, SDC1, PHLDB2, TNFRSF8, ADAM8, TLR7, STEAP3, POU3F2, IFT1, FBXO16, ARL3, BEX1, and GAPDH were provided in Table 1.

Table 1 List of Primers Used in the Study

Library Preparation and RNA-Sequencing

PDX-L14 cells were transfected with GFP-tagged FOSL1 and corresponding controls. At 72 h, the total RNA of each sample was extracted from the cells using the TRIzol reagent. Whole transcriptome sequencing RNAseq) was conducted by Novogene (Sacramento, CA) using 3 replicates of samples for each group. Messenger RNA (mRNA) was isolated from total RNA utilizing poly-T oligo-conjugated magnetic beads. The mRNA was fragmented, and the first strand of cDNA was subsequently generated with random hexamer primers. The second strand of cDNA was subsequently generated utilizing either dUTP for a directional library or dTTP for a non-directional library. The quality of the library was evaluated by Qubit and real-time PCR for quantification, and a bioanalyzer for size distribution analysis. Quantified libraries were aggregated and sequenced on the Illumina platform according to effective library concentration and data specifications.

RNA-Seq Data and Computational Analysis

Original image data from high-throughput sequencing on Illumina platforms are converted to sequenced reads (Raw data/Raw reads) using CASAVA base calling. The raw data are stored in FASTQ (fq) format files, which include sequences of reads along with corresponding base quality scores. Differential gene expression analysis was performed using the DESeq2 software with normalization achieved through the DESeq method. The criteria for identifying differentially expressed genes were set at |log2(Fold Change)| ≥ 2 and p<0.05. Genes that were differentially expressed with a significant false discovery rate (FDR<0.05)³⁵ were selected for further analysis. FDR was calculated using the Benjamini-Hochberg (BH) method.

Venn Diagram

The Differential Expression Venn diagram was used to illustrate the differential genes of the compared treatment group and the control group. The sum of all the numbers in the circle represents the total number of genes in the compared groups, and the overlapping area indicates the number of differential genes shared between the groups.

Cluster Analysis

The differential gene set was created by combining all the genes that were expressed differently in the comparison group. For experiments with more than two groups, cluster analysis was performed on different gene sets, allowing genes with similar expression patterns to be clustered together. The FPKM values of genes were clustered using mainstream hierarchical clustering, and the rows were homogenized using Z-score. The heat map facilitated the grouping of genes or samples that exhibited comparable expression patterns. The color in each grid represents the result achieved after standardizing the expression data rows (often ranging from −2 to 2), rather than the actual gene expression value. Therefore, the colors in the heat map are not only compared horizontally (the expression of the same gene in different samples) but also compared vertically (the same sample). Clustering occurs not only across groups but also among individual samples.

Volcano Plots

Volcano plots were utilized to deduce the overall distribution of genes that were expressed differentially. The clusterProfiler³⁶ software was utilized for enrichment analysis, encompassing GO Enrichment, DO Enrichment, KEGG, and Reactome database Enrichment. The gene set enrichment analysis (GSEA) software^37,38 was utilized to conduct gene enrichment analysis and examine the signaling pathways controlled by FOSL1.

To address the concern about gene set size bias, we acknowledge that clusterProfiler may be influenced by potential gene set size biases and pathway redundancies. To mitigate this, we utilized the default settings in clusterProfiler, which incorporated statistical adjustments, including the hypergeometric test and multiple testing correction (eg, Benjamini-Hochberg), to minimize false positives and account for variability in gene set size. Regarding pathway redundancies, we acknowledge that overlapping results can arise from pathway redundancies. To address this, we employed the simplify function in clusterProfiler, which clusters similar Gene Ontology terms based on semantic similarity, thereby reducing redundancy in enriched pathways.

To minimize batch and off-target effects, we used randomization of sample processing and sequencing across different experimental runs to ensure that any technical variation was evenly distributed across conditions. We applied rigorous quality control measures to assess the integrity of RNA-seq data, including visual inspection of data distribution (through PCA or clustering) to ensure that the experimental conditions were well represented and that no confounding factors from technical artifacts were present.

Immunohistochemistry (IHC)

IHC staining was performed on 5-µm thick microarray slides. The slides underwent baking at a temperature range of 55–60°C for a duration of 30 to 60 minutes. Subsequently, they were subjected to pretreatment in a citrate buffer with a pH of 6 for 10 minutes at a temperature of 100°C. The samples were subsequently placed in a 4°C incubator overnight with primary antibodies (diluted at a ratio of 1:100), followed by a 60-minute incubation at room temperature with secondary antibodies (diluted at a ratio of 1:200). The detection method employed was a streptavidin-biotin immunoperoxidase approach utilizing diaminobenzidine (DAB). The slides were stained with Mayer’s hematoxylin for a duration of 2 minutes and observed using a light microscope.

HSORE Determination

The staining intensity in the TMA was evaluated as either negative or positive in three separate bright regions, with a minimum of 100 cells per field. The HSCORE for antigens STEAP3 and ARL3 was calculated using the formula HSCORE = Ʃpi (i + 1), where “i” denotes the intensity levels (0, 1, 2, or 3), and “pi” indicates the percentage of cells at each intensity.^18,39,40 Two impartial evaluators conducted a blind evaluation and assessment of the slides.

Tissue Microarray (TMA)

Glioma tissue arrays obtained from glioma patients were procured from BioCoreUSA Corporation (https://biocoreusa.com/default.aspx). The tissue arrays are commercially available and do not need ethics committee approval. The biopsy features encompassed information such as the individual’s age, gender, the specific organ or anatomical place affected, the degree of severity, and the pathological diagnosis based on sections stained with hematoxylin and eosin (HE). The slides (product no. GL1001b) comprised 75 glioma examples: grade II included 51 cases (47 astrocytomas, 2 oligodendrogliomas, and 2 oligoastrocytomas); grade III included 12 cases (anaplastic astrocytoma); grade IV included 12 cases (glioblastoma); and there were 10 cases of normal brain tissue.

Immunofluorescence Staining

The glioma cells were immobilized using a 4% solution of paraformaldehyde (PFA), made permeable with a 0.02% solution of Triton-X-100, and then treated with a 10% solution of goat serum to prevent non-specific binding. The cells were subsequently cultured at a temperature of 4°C overnight with a primary antibody (diluted at a ratio of 1:100) in a solution that prevents non-specific binding. This was followed by incubation with a secondary antibody conjugated with either FITC or Alexa Fluor 568 for a duration of 1 hour at room temperature. The cultures were rinsed with PBS and then placed on slides with VECTASHIELD mounting media containing DAPI. Samples were imaged using a Zeiss AxioImager or Zeiss Confocal Microscope LSM700. Total nuclear immunofluorescence intensity with FOSL1 was determined using the Fiji software.

Results

1. Differentially expressed mRNA in glioma PDXL-14 overexpressing FOSL1. To identify the underlying mechanisms of FOSL1 in glioblastoma development, we utilized the proneural subtype PDX-L14, which expresses a very low level of FOSL1.¹⁸ After PDX-L14 cells were transfected with GFP-tagged FOSL1 and corresponding control, we performed RNA-sequencing (RNA-seq) and transcriptome analyses of PDX-L14 cells after FOSL1 overexpression and corresponding control (3 replicates for each group). High transfection efficiency was demonstrated by significantly increased FOSL1 mRNA levels measured by qPCR (Figure 2A, top) and elevated FOSL1 protein levels detected by Western blot (Figure 2A, bottom). Additionally, we confirmed that the nuclear localization of FOSL1 was not changed by the GFP tag (Figure 2B).

Figure 2 High transfection efficiency of FOSL1 overexpression in PDX-L14 cells. (A) Increased expression of FOSL1 mRNA (top) detected by qPCR and protein levels by Western blot (bottom) after PDX-L14 cells were transfected with GFP-tagged FOSL1 and corresponding control. (B) PDX-L14 cells transfected with a GFP-tagged FOSL1 and control for 72 h were stained for FOSL1. Magnification 40x.

2. FOSL1 overexpression in PDXL-14 resulted in differentially expressed mRNA of genes in glioma. Pair-wise differential expression analyses showed a total of 493 DEGs (fold change >0.5; adjusted p<0.05) identified in the FOSL1 group when compared with the control (Figure 3A). As illustrated in the Venn diagram (Figure 3B), 9 genes consistently exhibited differential expressions across all three pairs, namely FOSL1, RGPD6, IGFN1, OAS2, IFI44L, KRT75, STEAP3, ITGA5, and ARL3. Additionally, 32 genes showed differential expressions in the comparison of PDX-L14 control 1 vs PDX-L14 FOSL1 overexpression 1, as compared to the pairs of PDX-L14 control 2 vs PDX-L14 FOSL1 overexpression 2. Furthermore, 10 genes were found to be differentially expressed in the comparison of PDX-L14 control 1 vs PDX-L14 FOSL1 overexpression 1, as opposed to the pairs of PDX-L14 control 3 vs PDX-L14 FOSL1 overexpression 3. Similarly, 11 genes exhibited differential expressions in the comparison of PDX-L14 control 2 vs PDX-L14 FOSL1 overexpression 2, as compared to the pairs of PDX-L14 control 3 vs PDX-L14 FOSL1 overexpression 3. Among 493 DEGs, we identified 152 genes that were upregulated and 341 genes that were downregulated from triplicate samples of PDX-L14 cells under both control conditions and those overexpressing FOSL1, as shown in Figure 3C. Expression heatmap of differentially expressed mRNA in triplicate samples of PDX-L14 cells under control conditions and in cells overexpressing FOSL1 was shown in Figure 3D. These results indicate that FOSL1 overexpression in PDXL-14 resulted in differentially expressed mRNA of genes in glioma.

Figure 3 Transcriptome analyses of human PDX-L14 cells with FOSL1 overexpression. (A) Volcano plot of differentially expressed genes from RNAseq data. The x-axis represents the log2 fold change in gene expression between samples, while the y-axis shows the statistical significance of the differences, either as -log10 p-value or log10 adjusted p-value. The blue dashed line indicates the threshold for differential gene screening. Green dots represent downregulated genes (n=341), and red dots represent upregulated genes (n=152). (B) Venn diagram illustrating the overlap of differentially downregulated and upregulated mRNA in triplicate samples of PDX-L14 cells under control conditions and in cells overexpressing FOSL1. (C) Of the 493 DEGs, 152 genes were identified as upregulated, while 341 as downregulated in triplicate samples of PDX-L14 cells under both control conditions and FOSL1 overexpression. (D) Cluster analysis of FPKM values: The overall results of FPKM cluster analysis were visualized using the log2 (FPKM+1) value. Red indicates genes with high expression levels, while blue indicates genes with low expression levels. The gradient from red to blue reflects the range of log2 (FPKM+1) values from high to low.

3. Analysis of common differential FOSL1 target genes and their functional enrichment. By analyzing differentially expressed genes, we aimed to identify significant associations with biological functions or pathways related to the overexpression of FOSL1. The enrichment analysis was performed using the clusterProfiler software³⁶ encompassing GO enrichment, KEGG enrichment, Reactome enrichment, DO enrichment, and DisGeNET enrichment.

In Figure 4A, the left panel highlights the top 30 most significant GO Terms, including nuclear-transcribed mRNA catabolic process, mitochondrial protein complex, and cytochrome-c oxidase activity. GO terms related to CSC, RNA processing, NF-κB pathway, and cell proliferation were specifically displayed in the right panel of Figure 4A. For KEGG enrichment results (Figure 4B), the top 23 KEGG most significant pathways were selected, including oxidative phosphorylation, NOD-like receptor signaling pathway, proteasome, IL-17 signaling pathway, antigen processing and presentation, and ferroptosis. One of the FOSL1 target genes identified in this study, STEAP3, was shown to have high expression at both mRNA and protein levels. This gene has been previously reported to promote ferroptosis.⁴¹

Figure 4 Functional enrichment results of FOSL1 target genes. (A) GO enrichment results: Left: highlights the top 30 most significant GO Terms. Right: Displays GO terms related to CSC, RNA processing, and cell proliferation. The abscissa represents the ratio of the number of differential genes associated with the GO term to the total number of differential genes, and the ordinate represents the GO Term. The size of a point indicates the number of genes annotated to a specific GO term, and the color gradient from red to green represents the significance level of enrichment (padj). (B). KEGG enrichment results: this displays the top 23 most significant KEGG pathways. The abscissa represents the ratio of the number of differential genes associated with the KEGG pathway to the total number of differential genes, and the ordinate represents the KEGG pathway. The point size indicates the number of genes annotated to a specific KEGG pathway, and the color gradient from red to green represents the significance level of enrichment.

Figure 5A, the left panel displays the top 20 most significant Reactome pathways, such as peptide chain elongation, eukaryotic translation initiation, and elongation. The right panel of Figure 5A specifically highlights Reactome pathways associated with CSC signaling, cell proliferation, cell cycle, and RNA processing. In the context of DO Enrichment Analysis, the Human Disease Ontology served as a community-driven, standards-based ontology that provides a disease interface between data resources. Terms with p < 0.05 were deemed significant enrichment. From the DO enrichment results, the top 20 most significant DO Terms were selected for display (Figure 5B), with malignant glioma identified as a target of the FOSL1 gene. Additionally, DisGeNET enrichment analysis results (Figure 5C) showcased the 20 most significant DisGeNET pathways, including multiple malignant cancers.

Figure 5 Functional enrichment results of FOSL1 target genes. (A) Reactome enrichment analysis results: Left: This displays the top 20 most significant Reactome pathways. Right: This highlights Reactome pathways associated with CSC signaling, cell proliferation, and RNA processing. The abscissa represents the ratio of the number of differential genes associated with the Reactome pathway to the total number of differential genes, and the ordinate represents Reactome Pathway. The size of a point indicates the number of genes annotated to a specific Reactome pathway, and the color gradient from red to purple represents the significance level of the enrichment. (B) DO enrichment results: This displays the top 20 most significant DO terms. Malignant glioma is identified as a target of the FOSL1 gene. The abscissa represents the DO Term, and the ordinate represents the significance level of the pathway enrichment, with higher values corresponding to higher significance. (C) DisGeNET enrichment analysis results: This displays the top 20 most significant DisGeNET pathways, including multiple malignant cancers. The abscissa represents the ratio of the number of differential genes associated with the DisGeNET pathway to the total number of differential genes, and the ordinate represents DisGeNET Pathway. The size of a point indicates the number of genes annotated to a specific DisGeNET pathway, and the color gradient from red to purple represents the significance level of the enrichment.

These analyses offer valuable insights into the functional enrichment of differential FOSL1 target genes and their associations with various biological functions, pathways, and diseases.

4. Directed Acyclic Graph (DAG) in glioma PDXL-14 overexpressing FOSL1. In PDX-L14 glioma cells overexpressing FOSL1, DAG was utilized to visualize the enriched Gene Ontology (GO) terms of differentially expressed genes and their hierarchy (Figure 6). The top 5 enriched phrases for the biological process (GO bp) are SRP-dependent cotranslational protein targeting to membrane, cotranslational protein targeting to membrane, protein targeting to endoplasmic reticulum (ER), establishment of protein localization to ER, and translational initiation (Figure 6A). For cellular components (GOcc), the top 5 enriched terms are cytosolic ribosome, ribosomal subunit, ribosome, inner mitochondrial membrane protein complex, and mitochondrial protein complex (Figure 6B). The five most enriched terms for molecular functions (GOmf) are: structural constituent of ribosome, electron transfer activity, cytochrome-c oxidase activity, heme-copper terminal oxidase activity, and oxidoreductase activity, which acts on a heme group of donors with oxygen as the acceptor (Figure 6C). These results indicate that in glioblastoma, the FOSL1 gene is connected to and hierarchically influences multiple organelle functions, such as those in mitochondria.

Figure 6 Gene Ontology (GO) enrichment analysis represented as directed acyclic graph (DAGs) in glioma PDXL-14 overexpressing FOSL1. In this project, the DAG of biological processes, molecular functions, and cellular components were sequentially drawn to provide a comprehensive visualization of the enriched GO terms and their relationships. (A) Biological Process (GObp): DAG was employed to visualize the enriched GO term of differentially expressed genes and their hierarchy. In this graph, each branch represents the hierarchical relation, with functional specificity increasing from the top to bottom. The top 5 results of GO enrichment analysis are selected as primary nodes (depicted by boxes) in the DAG, with related GO terms connected hierarchically. The degree of enrichment is represented by varying shades of color, with darker shades indicating higher enrichment levels. (B) Cellular Component (GOcc): DAG depicting the enriched GO terms associated with cellular components. (C) Molecular Function (GOmf): DAG visualizing the enriched GO terms related to molecular functions.

5. Gene set enrichment analysis (GSEA) of RNA-seq data in glioma PDX-L14 overexpressing FOSL1. We further explored the signaling pathways modulated by FOSL1 in glioma PDX-L14 cells. The GSEA results demonstrated that FOSL1 activates IκB (Figure 7A) and NF-κB signaling pathways (Figure 7B), as well as GSC-related signaling pathways, including BMP (Figure 7C) and WNT (Figure 7D). FOSL1 also promoted pathways associated with cell proliferation, such as PI3K/AKT (Figure 7E) and VEGF (Figure 7F), while enhancing cell-to-cell communication (Figure 7G). Additionally, FOSL1 upregulated TRP channel expression/activity (Figure 7H) and influenced axon guidance in glioma (Figure 7I). Furthermore, FOSL1 affected RNA processing and splicing (Figure 8A–D), miRNA splicing and capping (Figure 8E and F), proteasome and ubiquitination signaling (Figure 8G–I), as well as cell cycle progression and cell death (Figure 8J and K), and TNF receptor-associated factors (TRAF6) signaling (Figure 8L). These results indicated that the FOSL1 gene is related to multiple biological functions, including cell proliferation, cancer cell stemness, angiogenesis, and ion channel activities. These data are consistent with our previous findings on FOSL1 interaction with NF-κB⁴² and TRPM7 channels.²⁵

Figure 7 Signaling pathway undergoing positive changes in response to FOSL1 overexpression in glioblastoma. Gene set enrichment analysis of RNA-seq data shows that FOSL1 activates the IκB (A) and NF-κB signaling pathways (B), as well as GSC-related signaling pathways, including BMP (C) and WNT (D). FOSL1 also promoted pathways associated with cell proliferation, such as PI3K/AKT (E) and VEGF (F), while enhancing cell-to-cell communication (G). Additionally, FOSL1 upregulated TRP channel expression/activity (H) and influenced axon guidance in glioma (I).

Figure 8 Signaling pathways exhibiting negative changes in response to FOSL1 overexpression in glioblastoma. FOSL1 exhibited effects on RNA processing and splicing (A–D), mRNA splicing and capping (E and F), proteasome and ubiquitination signaling (G–I), cell cycle progression and cell death (J and K), and TRAF6 signaling (L).

6. The prognostic significance of a risk signature composed of FOSL1 target genes was evaluated across various stratified cohorts of glioma cases. By examining genes with differential expression in FOSL1 overexpression and identifying in at least two pairs comparison as shown in Figure 3B, we found 8 genes that were upregulated and 4 genes that were downregulated. To assess the prognosis value of these genes, survival analysis was conducted using TCGA data via GlioVis software. The Kaplan-Meier illustrated that glioma patients exhibiting elevated expression levels of genes, including Integrin subunit alpha 5 (ITGA5), Syndecan 1 (SDC1), Pleckstrin Homology Like Domain Family B Member 2 (PHLDB2), TNFRSF8 (CD30), ADAM8, TLR7, STEAP3, and POU3F2 (Figure 9A–H), all upregulated in response to FOSL1 overexpression, and these patients had shorter overall survival compared to those with lower expression levels of these genes. Conversely, patients with glioma displaying increased expression levels of genes such as IFIT1, FBXO16, ARL3, and the Brain-Expressed X-linked (BEX1) (Figure 9I–L), downregulated following FOSL1 overexpression, exhibited better prognosis compared to those with low expression levels. These data support the notion that FOSL1 serves as a prognostic marker in glioblastoma.

Figure 9 Survival analysis assessing the prognosis significance of FOSL1 target genes using the TCGA data through the GlioVis software. Kaplan-Meier curves revealed in panels (A)–(H) showed that glioma patients with elevated expression levels of genes including ITGA5, SDC1, PHLDB2, TNFRSF8 (CD30), ADAM8, TLR7, STEAP3, and POU3F2, which are upregulated in response to FOSL1 overexpression, exhibited poor outcomes. Conversely, panels (I–L) illustrate that patients with glioma who have increased expression levels of genes such as IFIT1, FBXO16, ARL3, and BEX1, downregulated following FOSL1 overexpression, demonstrated improved prognosis.

7. Validation of FOSL1 target genes was performed by assessing mRNA expression through qPCR and protein expression through IHC. To validate the RNA-seq findings regarding the upregulation and downregulation of genes induced by FOSL1 overexpression, an FOSL1 plasmid was transfected into A172, U87MG, and PDX-L14 cells. Following this, qPCR was used to measure the mRNA levels of these genes. As shown in Figure 10A, for upregulated genes such as ITGA5, SDC1, PHLDB2, TNFRSF8, ADAM8, TLR7, STEAP3, and POU3F2, mRNA expression increased following FOSL1 introduction into the three glioma cell lines. Specifically, the fold changes were as follows: ITGA5: 1.72 to 2.1, SDC1: 2.12 to 5.97, PHLDB2: 2.01 to 10, TNFRSF 8: 2 to 3.43, ADAM8: 2.19 to 3, TLR7: 1.9 to 3.45, STEAP3: 2.17 to 3.21, and POU3F2: 2.17 to 3.2. Conversely, genes such as IFIT1, FBXO16, ARL3, and BEX1 were consistently downregulated across all three cell lines upon FOSL1 overexpression, with fold changes for IFIT1 ranging from 0.44 to 0.62, FBXO16: 0.29 to 0.52, ARL3 0.38 to 0.61, and BEX1 0.32 to 0.47 (Figure 10B). IHC results confirmed the expression of STEAP3 in both the cytoplasm and nucleus, showing STEAP3 upregulation in grade III and IV glioma patients (Figure 10C and E, left and middle panes). In contrast, ARL3 was primarily expressed in the cytoplasm and downregulated in grade III and IV glioma patients (Figure 10D and E, right panel). This data from both mRNA and protein levels effectively validated and confirmed the FOSL1 targeted genes identified in the RNA-seq analysis.

Figure 10 Validation of FOSL1 target genes was carried out by evaluating mRNA expression levels using qPCR and protein levels by IHC. FOSL1 plasmid was transfected into A172, U87MG, and PDX-L14 cells, followed by qPCR analysis to assess the mRNA levels of these genes. (A) The mRNA expression of upregulated genes, such as ITGA5, SDC1, PHLDB2, TNFRSF8, ADAM8, TLR7, STEAP3, and POU3F2, was confirmed to be elevated. (B) Genes including IFIT1, FBXO16, ARL3, and BEX1 were consistently downregulated across the three glioma cell lines following FOSL1 introduction. (C) The IHC results verified the expression of STEAP3 in both the cytoplasm and nucleus, indicating its upregulation in grade III and IV glioma patients. (D) ARL3 was predominately expressed in the cytoplasm and exhibited downregulation in grade III and IV glioma patients. (E) Statistical analysis of STEAP3 (left and middle panels) and ARL3 (right panel). In Panels A and B, § indicates p<0.05 compared to the control in A172 cells. A t-test was used to analyze each gene expression compared to control. * indicates p<0.05 compared to control in U87MG cells, and # indicates p<0.05 compared to control in PDX-L14 cells. In Panel E, One-way ANOVA was used to analyze each gene expression between normal and glioma patients and among different grades of glioma. * indicates p<0.05, **indicates p<0.01, and *** indicates p<0.001.

Discussion

Previous studies demonstrated that FOSL1 exerts its cancer-causing effects by controlling multiple signaling pathways, such as MAPK, PI3K/Akt, Wnt/β-catenin, and TGF-β, which have a role in cell growth, survival, cell death, and the spread of cancer cells (see our recent review.³³) In addition, recent findings have also shown that noncoding RNAs (ncRNAs) have significant involvement in the control of FOSL1 expression and function. FOSL1 has been discovered to regulate the expression of several cytokines, chemokines, and extracellular matrix (ECM) proteins that play important roles in the tumor microenvironment (TME).³³ In the current study, we identified a number of novel signaling pathways in which FOSL1 is involved, including nuclear-transcribed mRNA catabolic process, mitochondrial protein complex, cytochrome-c oxidase activity, oxidative phosphorylation, NOD-like receptor signaling pathway, proteasome, IL-17 signaling pathway, antigen processing and presentation, and ferroptosis. This finding indicates that FOSL1 may have more diverse functions in different biological processes. FOSL1’s role is widespread in angiogenesis,⁴³ senesce,⁴⁴ and neural development.⁴⁵ Our most recent findings show that FOSL1 promotes glioma stemness and tumorigenesis by inhibiting p27 ubiquitination and promoting FOSL1 transcription (unpublished data).

Through performing an integrated bioinformatics analysis,^46,47 in this study, transcriptome data following FOSL1 overexpression in glioma cells identified 8 upregulated genes ITGA5, SDC1, PHLDB2, TNFRSF8 (CD30), ADAM8, TLR7, STEAP3, and POU3F2 and 4 downregulated genes IFIT1, FBXO16, ARL3, and BEX1 with prognostic values in glioblastoma, based on the TCGA data. These 8 upregulated and 4 downregulated genes were also validated using qPCR. Furthermore, IHC analysis of TMA samples from glioma patients revealed that STEAP3 was upregulated and correlated with glioma grade, whereas ARL3 was downregulated and also associated with glioma grade.

STEAP3, a six-transmembrane epithelial antigen of prostate 3, has been implicated in regulating the cell cycle and immune responses in glioma, as shown from both TCGA and CGGA data.⁴⁸ Our findings demonstrated that there is a high expression of STEAP3 in glioma at both the mRNA and protein levels (Figure 10). STEAP3, along with five other copper-binding proteins,⁴⁹ including angiogenin (ANG), coagulation factor V (F5), interleukin 1α (IL1A), lysyl oxidase L1(LOX1), LOX2, forms a six-gene risk-score model that accurately predicts unfavorable prognosis in glioma,^48–50 which is confirmed by our findings. STEAP3 also acts as an iron reductase involved in intracellular iron homeostasis, a process linked to ferroptosis.⁴¹ Ferroptosis is evident in the enriched KEGG pathways following FOSL1 overexpression (Figure 4B), suggesting that STEAP3 likely contributes to this mechanism. There is not substantial direct literature providing STEAP3 as a therapeutic target for glioma patients, its biological function and the context of iron metabolism in cancer suggest that it is a candidate worthy of further exploration in glioma research.

ARL3, ADP-ribosylation factor-like 3, is a member of the ADP-ribosylation factor (ARF) family.⁵¹ Members of this family are activated through GDP to GTP exchange.⁵² ARL3 primarily localizes to the mitochondria, endosomes, lysosomes, and proteasome, playing a crucial role in autophagy.⁵³ It has been recently identified as a binding partner of STAT3, enhancing its phosphorylation and nuclear accumulation.⁵² ARL3 plays a significant role in angiogenesis and immune cell infiltration within the GBM immune microenvironment.⁵¹ It influences the infiltration of immune cells in the GBM environment. Decreased ARL3 expression is linked to an increased percentage of dendritic cells (DCs), macrophages, NK cells, and CD4 T cells while showing a decreased percentage of CD8 T cells. Moreover, ARL3 acts as an inhibitor of angiogenesis by targeting VEGFA, making it a promising candidate for antiangiogenic treatment in GBM.⁵¹ Our study reveals that ARL3 is decreased in glioma patients and exhibits a negative correlation with glioma grade (Figure 10). As of the current knowledge, there is insufficient literature to claim that ARL3 is a well-established therapeutic target for glioma patients, but its biological significance suggests it could be a potential area for future glioma research.

Several other genes, which play pivotal roles in GBM tumorigenesis and progression, were also validated by qPCR. ITGA5 encodes a protein that is a part of the integrin alpha chain family, crucial for cell surface adhesion and signaling. Data from TCGA and CGGA demonstrates that high expression of ITGA5 is positively associated with aggressive clinicopathological features and poor survival outcomes in glioma patients. Additionally, ITAG5 is implicated in remodeling glioma immune infiltration and the immune microenvironment, leading to increased immune cell infiltration when ITGA5 expression is elevated.⁵⁴ Moreover, there is a positive correlation between ITGA5 expression and ECM factors. ITGA5-induced vascular mimicry (VM) has been observed to promote resistance to temozolomide (TMZ) and bevacizumab by altering glioma vascularization. Methylation patterns of ITGA5 have also shown an association with its expression levels and tolerance to TMZ.⁴⁸ Targeting ITGA5 in animal models has been shown to have anti-tumor effects such as reduced tumor growth and angiogenesis.⁵⁵

SDC1 has been identified to be associated with BRD4 expression according to the CGGA database. Downregulation of SDC1 expression in in vitro U87MG and U251 glioma cell cultures significantly suppresses GBM cell proliferation and inhibits DNA replication.⁵⁶ Additionally, the BRD4 inhibitor MZ1 exerts its anti-cancer effects by suppressing SDC1 in GBM.⁵⁶ In the context of radioresistant GBM, SDC1 plays a role in transporting transglutaminase (TGM2) from the cell membrane into the cytoplasm. Once in the cytoplasm, SDC1 binds to flotillin 1 (FLOT1) and BHMT, facilitating the transportation of TGM2 to lysosomes. The formation of a polymer involving these four molecules leads to the fusion of lysosomes and autophagosomes, ultimately enhancing the radioresistance of GBM.⁵⁷

PHLDB2 plays a pivotal role in EMT and serves as a significant mediator of cell migration, metastasis, and maintenance of cancer stem cells in gastric and colorectal cancers.^58,59 In glioma, m6A-modified circEPHB4 enhances the stability of SOX2, leading to the transactivation of PHLDB2. This process contributes to glioma growth, metastasis, and the maintenance of stemness.⁶⁰ In colorectal cancer treatment, PHLDB2 could be a potential target for overcoming cetuximab resistance.⁶¹

TNFRSF8 (CD30), ADAM8, and TLR7 are three genes that play crucial roles in the signaling pathway for the positive regulation of NF-κB. CD30 is a cell surface protein receptor belonging to the tumor necrosis factor receptor superfamily. It is prominently expressed in immune cells such as activated T cells, B cells, and natural killer (NK) cells, with higher expression observed in malignant cells compared to normal cells.⁶² CD30, a transmembrane protein of 120 kDa, interacts with TRAF to activate the NF-κB and extracellular signal-regulated kinase (ERK) pathways. Due to its significance, CD30 is regarded as a promising therapeutic target for aggressive lymphoid malignancies.⁶² In solid tumors like triple-negative breast cancer (TNBC), cytokine-related pathways, including CD30, along with other cytokines like CCL25, CXCL13, IL12RB2, IL21, TNFRSF13C, CCL7, and GDF5, have been utilized to predict prognosis for TNBC patients.⁶³ ADAM8, one of the 21 members of the ADAM family of type I transmembrane proteins,^64,65 regulates tumorigenesis, cell adhesion, invasiveness, and proliferation.^66–68 Chronic CD30 signaling in B cells can lead to lymphoma by promoting the expansion of plasmablasts and B1 cells, indicating CD30 could be a potential target for lymphoma treatment.⁶⁹ CD30 expression on CD4 T cells is required in the pathogenesis of experimental autoimmune encephalomyelitis. Blocking CD30 may be a way to interfere with the development of autoimmune diseases of the central nervous system.⁷⁰

Similar to matrix metalloproteinase (MMPs), ADAM proteins have been implicated in regulating tumorigenesis, invasiveness, and proliferation.^68,71 Specifically, ADAM8 is involved in mediating chemoresistance induced by TMZ and enhanced invasion of GBM cells. It does so by modulating pERK1/2 and pAkt pathways and by shedding invasion-relevant substrate CD44 and c-met (HGF). Targeted inhibition of ADAM8 in future therapeutic approaches could potentially optimize TMZ chemotherapy and hinder the formation of recurrent GBM.⁷² High expression of ADAM8 is correlated with poor clinical outcomes in pancreatic ductal adenocarcinoma (PDAC). A peptidomimetic ADAM8 inhibitor, BK-1361, can prevent ADAM multimerization, reduce the invasiveness of PDAC cells, and decrease tumor burden and metastasis in mice.⁷³

TLR7 is not only expressed in immune cells, such as monocytes, macrophages, dendritic cells, T lymphocytes, B lymphocytes, and natural killer cells, which are involved in the inflammatory response and immune regulation but also in non-immune cells, such as tumor cells.⁷⁴ T cell-independent eradication of experimental glioma by intravenous TLR7 agonist-loaded nanoparticles⁷⁵ is achieved not through direct antitumor effects but through a proinflammatory shift in macrophages and increased ROS production.⁷⁵ Overreaction of TLRs can lead to high levels of IFN and other cytokines, resulting in chronic inflammation. Consequently, small-molecule TLR antagonists have been investigated as potential cancer treatments.⁷⁴ TLR7/8 agonists as potential therapeutics for tumor-targeted immunotherapy have been developed in various cancers.⁷⁶

POU3F2 acts as a master regulator of neuronal differentiation, functioning as a POU-domain transcription factor in developmental biology.⁷⁷ Studies have shown that inhibiting POU3F2 expression leads to a significant reduction in proliferation, migration, and invasion in small cell lung cancer (SCLC), prostate cancer, and melanoma.^78–81 Additionally, in glioma, Cui et al reported that tumor suppressor miR-146a directly targets POU3F and SMARCA5, two transcription factors that mutually regulate each other. This targeting significantly impairs GBM stemness and enhances resistance to TMZ.⁸² POU3F2 has shown certain biological characteristics and functions that suggest its potential as a cancer drug target. As of now, there is indeed a lack of widely available and clinically approved direct inhibitors or antibodies specifically targeting POU3F2 in the context of cancer. Research in this area is still in progress, and no such targeted therapies have reached the mainstream of cancer treatment.

IFIT1, an interferon-induced protein containing tetracopeptide repeats, is pivotal in suppressing growth and promoting apoptosis in cancer cells. Interferon is known to upregulate IFIT1 expression while inhibiting the expression of O-6-methylguanine-DNA methyltransferase (MGMT).⁸³ Research by Zhang et al suggests that IFIT1 may serve as a predictive biomarker for favorable clinical outcomes. Additionally, their findings indicate a negative correlation between IFIT1 and MGMT expression induced by interferon. Combining IFIT1 with MGMT expression provides a more accurate prediction of prognosis in newly diagnosed glioblastoma cases.⁸³

It was found that IFIT1 is highly expressed in the central nervous system and is significantly positively expressed in the lesion tissues of central nervous system tuberculosis patients. It is suggested that IFIT1 is an important substance involved in central nervous system tuberculosis, indicating its potential as a biomarker and possible drug target for related diseases.⁸⁴

FBXO16, a component of the SCF E3 ubiquitin ligase complex, interacts with β-catenin to facilitate its degradation.⁸⁵ Khan et al demonstrated that FBXO16 functions as a tumor suppressor in GBM. In glioma, reduced FBXO16 expression fails to target β-catenin in tumor cells, leading to its nuclear accumulation and subsequent activation of Wnt signaling. This promotes glioma cell proliferation and malignancy.⁸⁶ Under normal growth conditions, FBXO16 proteasomally degrades β-catenin in a GSK-3β-independent manner. The low expression of FBXO16 in glioma cells leads to nuclear accumulation of β-catenin, resulting in the activation of the Wnt signaling pathway, which makes glioma cells highly proliferative and malignant. This suggests that FBXO16 is closely related to the occurrence and development of glioblastoma and can be considered as a potential drug target for glioblastoma treatment.⁸⁶ Under normal growth conditions, FBXO16 proteasomally degrades β-catenin in a GSK-3β-independent manner. The low expression of FBXO16 in glioma cells leads to nuclear accumulation of β-catenin, resulting in the activation of the Wnt signaling pathway, which makes glioma cells highly proliferative and malignant. This suggests that FBXO16 is closely related to the occurrence and development of glioblastoma and can be considered as a potential drug target for glioblastoma treatment.⁸⁶ Additionally, in ovarian cancer, FBXO16 interacts with the RPM3 domain of hnRNPL through its C-terminal region, inducing hnRNPL’s proteasomal degradation. Inadequate degradation of hnRNPL results in increased ovarian cancer cell proliferation in vitro and tumor growth in vivo, mimicking the effects observed with FBXO16 deficiency.⁸⁷ The dysregulated expression of FBXO16, along with P3D1, DCAF10, LRFN4, PTPN2, SAYSD1, and ZNF426, is significantly associated with ovarian cancer progression, and these seven genes can serve as independent prognostic markers for the disease.⁸⁸

The BEX family proteins encompass five human proteins, including BEX1, BEX2, BEX3, BEX4, and BEX5. These proteins are widely expressed across various tissues and are involved in neuronal development. While BEX1 and BEX3 have been shown to exhibit tumor-suppressive properties, BEX2 appears to function as an oncogene.⁸⁹ In specific patient subgroups with acute myeloid leukemia and chronic myeloid leukemia, BEX1 expression is notably diminished.⁸⁹ In glioma patients, Foltz et al identified BEX1 and BEX2 as genes that undergo silencing, accompanied by significant promoter hypermethylation. These epigenetically silenced genes may serve as tumor suppressors, biomarkers for disease diagnosis and detection, and therapeutically reversible modulators of critical regulatory pathways involved in glioma pathogenesis.⁹⁰ Our study further supports the tumor-suppressive role of BEX1 in glioma. BEX1 engages in an interaction with RUNX3, and this interaction serves to impede RUNX3 from suppressing the transcription of β-catenin, As a consequence, the Wnt/β-catenin signaling pathway gets activated, which in turn is crucial for maintaining the stemness in hepatoblastoma (HB) and CSC-hepatocellular carcinoma (HCC). In the case of non-CSC-HCC, when the expression of BEX1 is downregulated, it results in the inhibition of the activation of the Wnt/β-catenin signaling pathway. This indicates that BEX1 holds the potential to be a valuable therapeutic target for the treatment of HB and CSC-HCC.⁹¹

GO enrichment analysis revealed that RNA processing and ubiquitin-like protein ligase binding are associated with the overexpression of FOSL1. Our unpublished immunoprecipitation data indicates that FOSL1 physically interacts with TRIM21, a ubiquitin E3 ligase, leading to the degradation of the tumor suppressor gene p27 in glioblastoma cell lines and contributing to glioblastoma tumorigenesis. Additionally, FOSL1 interacts with hnRNPE1 and hnRNPH1, further supporting its role in RNA processing. FOSL1 is an important drug target for glioma patients. Targeting FOSL1 will probably downregulate the genes positively upregulated by FOSL1 and upregulate the genes negatively regulated by FOSL1 in glioma patients.

Limitations of This Study

Due to the heterogeneity of GBM and GSC, we selected three cell lines, A172, U87MG, and PDX-14, representing three typical GBM subtypes: the neural, mesenchymal, and proneural subtypes, as previously discussed.²⁵ However, this study has a limitation. Specifically, it does not include more PDX lines and animal studies. Another limitation is that we did not validate all the genes upregulated or downregulated by FOSL1. Instead, we validated 8 upregulated and 4 downregulated genes with |log2(Fold Change)| ≥ 2 and p<0.05 under FOSL1 overexpression conditions. All these genes showed prognostic value in glioma. Finally, In the current study, we used transcriptome to identify the FOSL1-related signaling molecules and pathways. However, multiple omics, including genomics, epigenomics, proteomics, and metabolomics, are worth integrating with this transcriptome data to discover more targets for GBM patients.^92,93

Conclusion

This transcriptome study elucidates that FOSL1 plays a multifaceted role in various biological processes and signaling pathways. Specifically, FOSL1 is participates in the regulation of key signaling pathways, such as the NF-κB, Wnt, and BMP pathways, which are critical for maintaining glioma stemness. Additionally, FOSL1 is involved in RNA processing and ubiquitination regulation. Notably, we identified 8 upregulated and 4 downregulated genes associated with glioblastoma prognosis. After this study, future efforts should be on targeting these genes to reverse the proneural to mesenchymal transition.