Prognostic value and immunological role of FOXM1 in human solid tumors

FOXM1 gene expression in human cancers

According to our findings, FOXM1 expression was significantly upregulated in 22 types of solid tumors from the TCGA database (Figure 1A). Further investigation of FOXM1 protein expression was undertaken in the HPA cohort. FOXM1 protein expression in BRCA, CESC, LUAD, STAD and THCA tumor tissues was significantly higher (Figure 1B). Next, our study validated FOXM1 mRNA expression in 7 matched LUAD normal and tissues by quantitative real-time PCR (qRT-PCR). LUAD tissues showed significantly increased FOXM1 expression levels compared with normal tissues (Figure 1C). Immunohistochemistry was further used to assess FOXM1 protein expression levels, demonstrating that LUAD had higher levels of FOXM1 protein expression than normal tissues (Figure 1D).

Figure 1. Expression levels of FOXM1 in pan-cancer. (A) Boxplot of the mRNA expression of FOXM1 in 31 types of solid tumor and normal tissues from TCGA database. (B) Representative IHC staining of FOXM1 in BRCA, CESC, LUAD, STAD and THCA normal and tumor tissues in HPA. (C) Expression of FOXM1 in tumor group were upregulated than that in normal group. (D) IHC analysis of FOXM1 in LUAD tissues. Representative images are shown. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001.

We next evaluated FOXM1 expression according to tumor stage, gender and age. Higher FOXM1 expression was dramatically related to advanced tumor stage in ACC, BRCA, ESCA, KICH, KIRC, KIRP, LIHC, LUAD, LUSC, TGCT (Figure 2). Interestingly, ACC, ESCA, KICH, KIRC, KIRP, LUAD patients with stage IV tumors expressed more FOXM1 than patients with stage I tumors. A comparison of FOXM1 expression in tumors according to gender and age is shown in Supplementary Figure 1.

Figure 2. The association of FOXM1 expression with the pathological stages of cancers. p < 0.05 was considered significant.

Pan-cancer prognostic value of FOXM1

To investigate how FOXM1 expression relates to prognosis, survival analyses were performed by the TCGA cohort. Patients with higher FOXM1 expression had poor overall survival (OS) in ACC, KICH, KIRC, KIRP, LGG, LIHC, LUAD, MESO, PAAD, SARC, SKCM, UVM based on Kaplan-Meier analysis. Conversely, patients with higher FOXM1 expression had better OS in THYM. In order to avoid the bias resulting from non-tumor related death, the relationship of FOXM1 expression with a poor disease-specific survival (DSS) was further evaluated. The results, much like that of the OS analysis, demonstrated that higher FOXM1 expression significantly predicted a poor DSS. There was a negative correlation between disease-free interval (DFI) and FOXM1 expression in KIRP, LUAD, PAAD, SARC, and THCA. Our study finally analyzed the association between FOXM1 expression and progression-free interval (PFI). Patients with higher FOXM1 expression had poor PFI in ACC, KIRC, KIRP, LGG, LIHC, LUAD, MESO, PAAD, PRAD, SARC, SKCM, THCA, UVM (Figure 3).

Figure 3. Kaplan-Meier survival curves comparing the expression of FOXM1 in pan-cancer. (A) Kaplan–Meier OS curves for the low expression group and high expression group of FOXM1 in 13 types of cancers; (B) Kaplan–Meier DSS curves for the low expression group and high expression group of FOXM1 in 11 types of cancers; (C) Kaplan–Meier DFI curves for the low expression group and high expression group of FOXM1 in 5 types of cancers; (D) Kaplan–Meier PFI curves for the low expression group and high expression group of FOXM1 in 13 types of cancers. p < 0.05 was considered significant.

According to the forest plots, FOXM1 could significantly affect the OS of 13 types of solid tumors. The expression of FOXM1 represented a high-risk indicator in ACC, BLCA, COAD, KICH, KIRC, KIRP, LGG, LIHC, LUAD, MESO, SKCM, UCEC, UVM, but a low-risk indicator in THYM. For DSS, similar results were observed. Regarding FOXM1 and DFI, Cox regression analysis revealed that FOXM1 expression impacted the DFI of KIRP, LUAD, LUSC, PRAD, SARC, THCA. In terms of PFS, forest plot revealed that the hazard ratios for FOXM1 were significant for 18 cancer types (Figure 4).

Figure 4. Forest plots of Cox regression analyses in pan-cancer. (A) Association of FOXM1 expression with OS; (B) Association of FOXM1 expression with DSS; (C) Association of FOXM1 expression with DFI; (D) Association of FOXM1 expression with PFI. p < 0.05 was considered significant.

Association of FOXM1 with tumor immune microenvironment

In order to understand how FOXM1 expression correlates with immune infiltration, we used the Estimating Relative Subsets of RNA Transcripts (CIBERSORT) method. Different subpopulations of invasive macrophages showed a correlation with FOXM1 expression. There was a positive association of FOXM1 expression with macrophage M1 in 9 types of cancers. BRCA, KIRC, LUAD, and STAD revealed a positive association of FOXM1 expression with macrophage M0, whereas KIRP and THYM showed a negative correlation. A positive correlation was found between FOXM1 expression and macrophage M2 infiltration in TGCT and THYM, but a negative correlation was found in KIRP, LIHC, and THCA (Figure 5). Moreover, there was a correlation between FOXM1 expression and infiltration levels of CD4+T cells, B cells, CD8+T cells, mast cells, dendritic cells, T cells follicular helper, NK cells and Tregs (Supplementary Figure 2).

Figure 5. Association of FOXM1 expression with macrophages. p < 0.0001 was considered significant.

Next, it was explored whether the expression of FOXM1 and immune scores correlated using the ESTIMATE algorithm. FOXM1 expression was positively associated with immune and stromal scores of KIRC and THCA, but negatively associated with GBM, LUAD, LUSC, PAAD, SKCM, STAD, and UCEC. Additionally, a negative relationship was found between FOXM1 expression and the ESCA and TGCT immune scores, and a negative relationship was found between FOXM1 expression and the BRCA, COAD, LIHC, and SARC stromal scores (Figure 6).

Figure 6. Association of FOXM1 expression with the immune score and stomal score. p < 0.05 was considered significant. p < 0.001 was considered significant.

Association of FOXM1 expression with TMB and MSI

We further explored the association of TMB and MSI with FOXM1 expression. A positive correlation was found between FOXM1 and TMB in 19 cancer types. While in ESCA and THYM, negative relationships were found (Figure 7A). Additionally, FOXM1 and MSI had a positive correlation in 13 cancer types, including ACC, CESC, CHOL, COAD, GBM, LIHC, LUSC, OV, PAAD, SARC, STAD, TGCT, UCEC, UVM (Figure 7B).

Figure 7. Associations of FOXM1 expression with TMB and MSI in cancers. (A) Radar map of association of FOXM1 expression with TMB, (B) Radar map of association of FOXM1 expression with MSI. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001.

Correlations of FOXM1 with immune checkpoint-associated genes

A gene coexpression analysis was conducted to investigate how FOXM1 is related to immune checkpoint genes. More than 40 immune checkpoint genes were assessed. In LIHC, SKCM, and THCA, FOXM1 expression correlated significantly with the major checkpoint genes. Interestingly, CD276 was positively associated with FOXM1 expression in most cancers, with the exception of CESC, CHOL, COAD, GBM, READ, SARC and UCS. Forty-one immune checkpoint markers had significant associations with FOXM1 expression in TGCT. Furthermore, CESC, CHOL, COAD, READ, USC, and UVM showed relatively small correlations (Figure 8).

Figure 8. Heatmap of the association of FOXM1 expression with 47 common immune checkpoints gene levels in 31 types of cancers. For each pair, the top left triangle represents the P-value, and the bottom right triangle represents the correlation coefficient. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001.

Immunotherapeutic response prediction value of FOXM1

As FOXM1 was closely related with tumor immune microenvironment, we sought to confirm whether it could predict response to immune checkpoint inhibitor treatment based on IMvigor 210 cohort. In patients with high FOXM1 levels, OS was significantly longer than in those with low levels (Figure 9A). Furthermore, higher FOXM1 was tested in responders than that in non-responders (Figure 9B). Anti-PD-L1 therapy was significantly more effective when FOXM1 levels were higher. (Figure 9C).

Figure 9. Validation of the immunotherapeutic predictive value of FOXM1. (A) Kaplan–Meier OS curves for patients with FOXM1 high and low expression subgroups in IMvigor 210 cohort. (B) FOXM1 expression was higher in responders than that in non-responders. (C) The proportion of patients with response to anti- PD-L1 therapy in FOXM1 high and low expression groups. p < 0.05 was considered significant.

Gene set enrichment analysis

Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) pathways enrichment was analyzed to investigate the biological significance of FOXM1 expression in various cancers (Supplementary Figures 3 and 4). In KEGG analysis, FOXM1 was related to “cell cycle” and “DNA replication” pathways in multiple tumors. Besides, FOXM1 could modulate some pathways, such as “cytokine receptor interactio”, “complement and coagulation cascades”. In GO analysis, FOXM1 could modulate some pathways, such as “antigen binding”, “B cell mediated immunity”, “T cell receptor complex” in CHOL; “adaptive immune response based on somatic recombination”, “immunoglobulin production” in KIRC; “antigen binding”, “Fc receptor mediated stimulatory”, “humoral immune response circulating immunoglobulin”, “negative regulation of immune response”, “positive T cell selection” in TGCT; “antigen binding”, “humoral immune response”, “immunoglobulin complex” in UCS. FOXM1 also regulated many other pathways associated with “cell cycle” and “chromosome segregation”. These results demonstrated that FOXM1 widely regulated cell proliferation and immunity signaling pathways.

Functional analysis of FOXM1 in lung cancer cells

We additionally selected a lung cancer cell line A549 as an example and performed functional experiments to further verify our bioinformatics findings. We designed and transfected siRNA of FOXM1 to knock down endogenous FOXM1 expression. The results of CCK8 assays identified that knockdown of FOXM1 with siRNA resulted in the decrease of the proliferation ability of A549 cells (Figure 10A). According to wound healing test, FOXM1 downregulation markedly diminished cell migration of A549 cells (Figure 10B). Similarly, transwell migration analysis confirmed that FOXM1 downregulation could inhibit A549 cells migration. Transwell invasion analysis also showed that after knockdown of FOXM1, A549 cells revealed a decrease in invasion capacity (Figure 10C).

Figure 10. FOXM1 knockdown inhibits lung cancer cell proliferation, migration and invasion. (A) The proliferation curves of A549 cells transfected with si-NC and si- FOXM1. (B) The wound healing assay results for A549 cells transfected with si-NC and si-FOXM1 at 0, 36 h, respectively. (C) The transwell migration and invasion assay for A549 cells transfected with si-NC and si-FOXM1. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001.

Data collection

RNA sequencing datasets, clinicopathological information and survival information for 31 solid tumors were acquired from the TCGA database (https://portal.gdc.cancer.gov/) on March 20, 2021. Table 1 provides a general overview of pan-cancer. Patient transcriptomic and clinical data for BLCA receiving anti-PD-L1 therapy (atezolizumab) were obtained from the IMvigor 210 trial data [43]. R version 3.6.3 was used for all data analysis.

FOXM1 expression analysis

Using the Wilcoxon test, we compared FOXM1 expression levels in tumor and normal tissue. Then, a box plot was displayed by the ‘ggpubr’ R package. From the HPA database, immunohistochemistry images of the FOXM1 protein were obtained for both normal and cancerous tissues. We also explored differences in FOXM1 expression across different stages, genders and ages using chi-square or Fisher’s exact tests. It is statistically significant if p < 0.05.

Isolation of RNA and qRT-PCR analysis

The First Affiliated Hospital of Nanjing Medical University provided the LUAD tissue. TRIzol reagent (Invitrogen, USA) was applied to extract total RNA, and PrimeScript RT reagent Kit (Takara Bio) was used to reversely transcribe into cDNA. For the detection of FOXM1 mRNA expression levels, qRT-PCR was conducted with an ABI StepOnePlus system (Applied Biosystems) and Sybrgreen-based qRT-PCR assays. GAPDH was used for normalizing the samples. Using 2^−ΔΔCt method, relative expression levels of mRNA were determined and each sample was performed in triplicate. The primers: FOXM1 (forward: 5′-TCTGCCAATGGCAAGGTCTCCT-3′ and reverse: 5′-CTGGATTCGGTCGTTTCTGCTG-3′).

Immunohistochemistry

The First Affiliated Hospital of Nanjing Medical University provided Archival formalin fixed paraffin-embedded (FFPE) specimens of LUAD tissues and corresponding para-tumor tissues. All of the tumor samples were again independently confirmed by two pathologists. Thin sections (5 mm) were cut from FFPE biopsy specimens. Following deparaffinization, the sample was rehydrated with gradient alcohol. The antigens were repaired by microwaving with citrate buffer followed by incubation with 3% H₂O₂. A primary antibody for FOXM1 (catalog no., bs-21487R; 1:200 dilution; USA) was incubated overnight at 4°C, followed by a secondary antibody at 32°C for 30 minutes. Afterwards, diaminobenzidine chromogen was used to stain. In the final step, hematoxylin counterstaining was followed by dehydration and cover-slipping with permanent media.

Survival analysis

Survival data from the TCGA was downloaded to evaluate the prognostic value of FOXM1. Based on the median FOXM1 expression, we divided each tumor sample into two groups. A Kaplan-Meier method was performed by the R packages ‘survival’ and ‘survminer’ to determine the correlation of FOXM1 with survival. Hazard ratios with 95% confidence intervals were collected by Cox proportional hazard regression analysis in 31 cancer types.

Immune correlation analysis

CIBERSORT helps distinguish 22 tumor infiltrating immune cells from other cells in tissues [44]. The associations of FOXM1 expression with 22 immune cell subtypes were estimated by the CIBERSORT algorithm.

ESTIMATE is a tool that predicts tumor purity [45]. Stromal and immune scores were obtained using R packages ‘estimate’ and ‘limma’. The correlation between FOXM1 and these scores was determined using Spearman correlation analysis.

TMB and MSI have been found to have close links with the immune response in recent studies [20, 25]. The association of TMB/MSI with FOXM1 expression was evaluated using Spearman correlation analysis. The results were shown by the R ‘fmsb’ package in radar chart. The association of FOXM1 expression with common immune checkpoint genes was further evaluated using Pearson’s correlation analysis. The results were shown by ‘Reshape 2’ and ‘R Color Brewer’ packages in heatmaps.

Gene set enrichment analysis

GSEA was conducted with ‘cluster-Profiler’ R package to explore the relevant signaling pathways between FOXM1 high and low expression groups in GO and KEGG. The enrichment significance criteria were |NES|>1, NOM p < 0.05, together with FDR q < 0.25.

Cell culture and siRNA transfection

The National Institute of Cells (Shanghai, China) provided lung cancer cell line A549 for this study. RPMI1640 medium with 10% fetal calf serum (Gibco) was used to culture the cells, which were maintained at 37°C with 5% CO₂ in a humidified atmosphere. The siRNA targeting FOXM1 were acquired from ShanghaGenePharma Co., Ltd., (Shanghai, China): Sense, 5′-GCCAAUCGUUCUCUGACAGAATT-3′ and antisense, 5′-UUCUGUCAGAGAACGAUUGGCTT-3′. A negative control siRNA (si-NC) was also obtained: Sense, 5′-UUCUCCGAACGUGUCACGUTT-3′ and antisense, 5′-ACGUGACACGUUCGGAGAATT-3′. FOXM1 siRNA or si-NC was transfected into A549 cells with Invitrogen™ Lipofectamine 2000 (Thermo Fisher Scientific, Inc., USA) as directed by the manufacturer.

Cell proliferation assay

The proliferation of A549 cells was examined by the CCK8 assay kit, as instructed by the manufacturer. In 96-well plates, transfected A549 cells were planted at a density of 2000 cells/well. Absorbance value of each well at 450 nm was detected.

Wound healing assay

We cultured transfected cells in a monolayer until confluence was achieved. A 200 μl pipette tip was applied to make an incision-like gap. The wound area was acquired with a microscope at 0 and 36 h after wounding. Analyzing the cell migration data was done with Image-Pro Plus 7.0 software.

Transwell migration and invasion assay

Using a transwell chamber system, migration and invasion were assessed. For the migration assay, 24 hours after transfection, 10,000 cells containing medium supplemented with 2% serum were planted onto the upper chamber. For the invasion assay, a ratio of 1:4 was used for mixing matrigel and serum-free medium, followed by seeding the mixture into an upper chamber. Next, 10,000 cells were placed in the upper chamber. Each chamber was placed into a complete medium in both assays. When 24 hours had passed, using 95% ethanol, we fixed the invaded cells on the lower chamber surface for 30 minutes. Next, crystal violet at a concentration of 0.1% was applied to stain the cells for about 15 minutes. Under a microscope, three random fields of stained cells were observed and counted.

Statistical analysis

Statistical analysis of the in vitro A549 cell assays was conducted by SPSS 20.0 software. The mean ± standard deviation (SD) of at least three experiments were presented. We defined statistically significant differences as two-sided p < 0.05 for Student’s t-tests.