Identification of JAK2 and FOXM1 expression as novel candidate biomarkers for predicting the benefit of immunotherapy in lung squamous cell carcinoma.

BACKGROUND: Lung squamous cell carcinoma (LUSC) accounts for about 30% of all non-small cell lung cancers (NSCLC). However, only a small percentage of LUSC patients gain benefit from immune checkpoint inhibitors (ICIs).
METHODS: This study analyzed LUSC patients from The Cancer Genome Atlas (TCGA), which were divided into 2 groups: PD-L1 high-expression/TMB-high (TPH) and PD-L1 low-expression/TMB-low (TPL) group based on programmed death-ligand 1 (PD-L1) expression and tumor mutational burden (TMB) status. The differences in tumor-infiltrating immune cells were estimated between the 2 groups. The overlap of differentially expressed genes and proteins (DEGs and DEPs) between 2 groups were used as candidate biomarkers. Kaplan-Meier curves were used to evaluate the association between risk score and overall survival (OS).
RESULTS: More abundant immune infiltration fractions were found in TPH group. Janus kinase 2 (JAK2) and forkhead box protein M1 (FOXM1) were identified as DEGs between the TPH and TPL groups. Subsequently, we developed a risk score that combined the expression of JAK2 and FOXM1 in an effort to accurately determine the survival risk of LUSC patients. Patients with high-risk [hazard ratio (HR), median OS, 43.1 months 1.924; 95% confidence interval (CI): 1.256 to 2.945; P=0.002) had shorter survival than those with low-risk (median OS, 70.0 months). External data verification found that JAK2 and FOXM1 were significantly expressed at a higher level in the responders receiving immunotherapy (P=0.038 and P=0.009, respectively).
CONCLUSIONS: The expressions of JAK2 and FOXM1 can be used as novel candidate biomarkers for predicting the benefit of immunotherapy in LUSC. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

Non-small cell lung cancer (NSCLC) is currently the leading cause of cancer-related deaths worldwide and lung squamous cell carcinoma (LUSC) is among the most common histological subtypes of NSCLC (1-3). Unlike lung adenocarcinoma (LUAD) with oncogenic driver variations, treatment for LUSC is limited and conventional chemotherapy has endured as the standard therapy (4-6). Campbell et al. found that NFE2L2, KDM6A, RASA1, NOTCH1 and HRAS were significantly mutated only in LUSC,including over 47% of LUAD and 53% of LUSC samples had more than five predicted neoepitopes, which indicated that immunotherapy strategy had great potential, the change of neoepitopes antigen includes the site p.Glu79Gln of NFE2L2 (7).

Recently, immune checkpoint inhibitors (ICIs) targeting programmed cell death 1 (PD-1), its ligand (PD-L1), and CTLA-4 have revolutionized the treatment of both LUAD and LUSC. Some studies found express of PD-L1 was more frequent in central type, tumor location that could predict expression status of PD-L1, and could potentially serve as clinical response to immunotherapy. To date, several anti-PD-1/PD-L1 or CTLA-4 antibodies including nivolumab, pembrolizumab, atezolizumab, and ipilimumab combination with or without chemotherapy have been approved as first-line treatment for patients with advanced LUSC. However, only a small subset of patients with NSCLC (20–25%) can benefit from ICIs (8). Increasing evidence indicates the prognostic value of immune-related signatures in the tumor microenvironment. Qu et al. found that genes related to prognosis in LUAD include: ICAM, MS4A1 and IL-16. LUAD had more memory B cells and CD8+T cells, while M0 macrophages were less than to early N stage. Meanwhile, for LUSC, only GSTA1 and HAS1 genes are associated with prognostic factors related to immune-related. It was found that the high invasiveness of naive CD4+ T cells was related to the more advanced T stage, while the CD8+T cells and M1 macrophages were higher in the infiltration advanced N stage in LUSC (9). Thus, the identification of predictive biomarkers of response to ICIs is emergent in LUSC. Emerging evidence has suggested that positive PD-L1 expression and high tumor mutation burden (TMB) could predict the response to ICIs in NSCLC (10,11); Tumor can evade immune surveillance by inhibiting the activation of immune cells with up-regulation of PD-L1. Compared with traditional therapies that directly target cancer cells, that anti-PD-1/PD-L1 antibodies reactivate the patient’s immune system to eradicate tumors, to improve patients anti-tumor immunity with different tumor types, including lung cancer Anti-tumor immunity (12). Some studies showed that when atezolizumab is as single-agent for PD-L1 high expression patients, the ORR, PFS and overall survival (OS) in treatment was significantly better than that of PD-L1 negative patients (13). However, these biomarkers still have some limitations. Firstly, not all patients with positive expression of PD-L1 benefit from immunotherapy. Secondly, the immunotherapy can occasionally result in hyperprogressive disease. Thirdly, the use of PD-L1 expression is controversial because of its varied definition, cut-off value, and spatial as well as temporal heterogeneity (14). Finally, TMB alone cannot fully predict response to ICIs; the combination of TMB with other biomarkers may improve the predictive ability. Thus, it is an urgent need to find novel effective biomarkers for precise immunotherapy.

With the progress of next-generation sequencing (NGS) technology, a growing number of high-throughput sequencing data on cancer have been published in public databases of the Gene Expression Omnibus (GEO), providing a novel strategy to solve the abovementioned dilemma (15).

Previous studies have confirmed the use of gene mutational signature as biomarkers for predicting the therapeutic efficacy of immunotherapy. However, the efficacy of immunotherapy in LUSC cannot yet be predicted completely and accurately. Due to the complexity of the immune response and immunotherapy, we hypothesized that in addition to genomic studies, transcriptomic and proteomic investigation might also be essential for accurately predicting the clinical benefit.

Janus kinases (JAKs) are a family of non-receptor tyrosine kinases, which are related to autoimmune diseases and malignant tumors. As a member of the Janus kinase family, JAK2 is an important part of the signal transducer and activator of transcription (JAK-STAT) signaling pathway, and plays an important role in promoting tumor phenotypes such as tumorigenesis, invasion, metastasis, proliferation, survival, angiogenesis, anti-apoptosis, and immune evasion (16-18).

Forkhead box transcription factor M1 (FOXM1) is a well-known proliferation-associated transcription factor which belongs to the forkhead box (FOX) transcription factor family, and is characterized by a conserved DNA-binding domain referred to as the FOX. It is involved in a variety of biological functions, including cell proliferation, cell cycle regulation, angiogenesis, cell migration, tumor invasion, aging, DNA damage repair, stem cell expansion, renewal, and immunoregulation (19). The loss of Cyclin-dependent kinase (CDK) family function can cause the expression of DNA damage repair genes to be silenced, which may increase the sensitivity of tumor cells to platinum-based chemotherapy (20).

In this study, we explored the differences in gene expression profiles between PD-L1 high-expression/TMB-high (TPH) and PD-L1 low-expression/TMB-low (TPL) groups of patients, and identified JAK2 and FOXM1 as novel candidate biomarkers to predict the clinical benefit of ICIs in LUSC. The results of this investigation are of great clinical significance for personalized immunotherapy in LUSC.

We present the following article in accordance with the REMARK reporting checklist (available at https://dx.doi.org/10.21037/atm-21-2186).

Methods

Participants

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). Messenger RNA (mRNA) expression profiles and relative clinical information of patients with LUSC were downloaded from The Cancer Genome Atlas (TCGA) data portal. We processed and normalized TCGA level 3 RNA-seq data, and the fragments per kilobase of exon model per million mapped fragments (FPKM) values were used for the gene expression levels. Reverse phase protein array (RPPA) data were also downloaded at the gene level. In total 16,383 expressed genes used from 502 LUSC samples and 237 expressed proteins across 325 LUSC samples.

The clinical follow-up information and gene expression profiles of NSCLC tissue in GSE126044 (21) and GSE136961 (22) were downloaded from the GEO (http://www.ncbi.nlm.nih.gov/geo/) database (23). The microarray data of GSE126044 and GSE136961 were based on GPL16791 and GPL24014 platforms, respectively, and the former included 16 patients (Submission date: Feb 04, 2019), while the latter included 21 patients (Submission date: Sep 05, 2019) with NSCLC receiving anti-PD-1 treatment.

Identification of DEGs

A total of 317 patients with LUSC were divided into 2 groups based on PD-L1 expression and TMB status: the top 25% of patients with high TMB and high level of PD-L1 expression were defined as the TPH group (23 participants), and the bottom 25% of patients with low TMB and low level of PD-L1 expression were defined as the TPL group (21 participants). The differences in clinicopathological characteristics between TPH and TPL groups were compared. The ‘limma’ package (http://bioinf.wehi.edu.au/limma) was used to identify the differentially expressed genes (DEGs) in the TPH and TPL groups and P<0.05 and |logFC|>1 were used as the cut-off criterions to select DEGs for mRNA and protein expressions. Finally, the overlapping DEGs between the 2 omics were detected ().

Figure 1

Study design and workflow overview.

Immune infiltration in LUSC

The relative percent of 22 immune cells in each LUSC sample was calculated using the CIBERSORT algorithm, which included gene expression of 22 leukocyte subtypes (24). At the same time, the gene expression profile data were used to quantify the infiltration of immune cells in tumor tissues by single-sample gene set enrichment analysis (ssGSEA) (25), and the enrichment score of 29 immune cells was obtained. Then, the infiltration level of 22 immune cells and the enrichment score of 29 immune cells were compared between TPH group and TPL group by Wilcoxon ranked-sum test.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis

GO terms analysis of JAK2 and FOXM1 was performed using the R package ‘clusterProfiler’ (https://rdocumentation.org/packages/clusterProfiler/). The bar-chart of GO results for cellular component (CC), biological process (BP), and molecular function (MF) were ranked by P value and exhibited. The KEGG pathways of DEGs at the significant level (P<0.05) were employed.

Protein-protein interaction (PPI) network construction and analysis

PPI networks were integrated from the Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org/cgi/input.pl; version: 11.0) (26). Cytoscape software (version: 3.6.1; https://cytoscape.org) (27) was used to establish and visualize the interaction network of DEGs.

Risk score construction and survival analysis

In order to determine an effective prognostic feature, we derived the risk score for survival analysis with JAK2 and FOXM1 gene- and protein-expressions, which was accessed on the basis of the following equation, where .m is the gene expression and .p is the protein expression. The median risk score of participants was used as cut-off value; afterwards, samples were divided into high and low-risk groups. Kaplan-Meier survival analysis was performed to compare the prognosis between 2 groups.

Validation for the performance of risk score

The performance of the risk score was verified through 2 external datasets, which included 16 and 21 patients with NSCLC receiving immune checkpoint blockade (ICB), respectively. Receiver operating characteristic (ROC) curves were applied to assess the sensitivity and specificity of ICB-response based on the risk score, and the R package ‘pROC’ (28) was used to quantify the area under the curve (AUC). The association between clinical outcome and the level of JAK2 or FOXM1 expression was evaluated through Wilcoxon ranked-sum test and shown by boxplot.

Statistical analysis

The ‘limma’ package (http://bioinf.wehi.edu.au/limma) was used to identify the differentially expressed genes (DEGs) in the TPH and TPL groups with P<0.05 and |logFC| >1. The enrichment analysis was performed using the R package ‘clusterProfiler’. PPI networks were visualized by Cytoscape software. Univariate and multivariate Cox regression analyses were used to evaluate the association of variables with OS in the ‘survival’ R package, as well as the log-rank test, was applied to compare the prognostic value. In this study, we used area under the curve (AUC) as the performance measurement method for predictive models, which was plotted using the ‘survivalROC’ R package, and all statistical tests were performed using R-3.6.3.

Results

Different immune infiltration and survival analysis between TPH and TPL

The clinical information of TPH and TPL LUSC patients is listed in . To investigate the association between PD-L1 expression, TMB, and immune infiltration in LUSC, we firstly calculated the percentages of 22 leukocyte cells of TPH and TPL participants using the CIBERSORT algorithm and then compared immune cell fractions. We found that the TPH group had more abundant macrophages M1 (P<0.01), resting natural killer (NK) cells (P<0.05), activated memory CD4 T cells (P<0.001), and CD8 T cells (P<0.05), but fewer memory B cells (P<0.05), activated dendritic cells (P<0.05), and neutrophils (P<0.05). It manifested the patients with high expression of PD-L1 and high TMB harbored more abundant immune cell infiltration, which might predict good adaptive immune response and survival when receiving immunotherapy ().

Table 1

Clinical and demographic characteristics of LUSC patients

Clinical characteristic	TPH (n=23)	TPL (n=21)	P value
Age (years)
Mean [SD]	64 [12]	68 [8.6]	0.12
Gender			0.41
Female	35% [8]	19% [4]
Male	65% [15]	81% [17]
Smoking history			0.66
1	4.3% [1]	14% [3]
2	30% [7]	33% [7]
3	8.7% [2]	4.8% [1]
4	52% [12]	48% [10]
No	4.3% [1]	0% [0]
Stage			0.39
I	43% [10]	62% [13]
II	35% [8]	29% [6]
III	22% [5]	9.5% [2]
KRAS			1.0
No	96% [22]	100% [21]
Yes	4.3% [1]	0% [0]
EGFR mutation			1.0
No	96% [22]	100% [21]
Yes	4.3% [1]	0% [0]
EML4-ALK fusion			1.0
No	96% [22]	100% [21]
Yes	4.3% [1]	0% [0]
Pulmonary function			0.11
Normal	48% [11]	76% [16]
Abnormal	52% [12]	24% [5]
Radiotherapy			0.9
No	87% [20]	81% [17]
Yes	13% [3]	19% [4]

Smoking history: 1, Lifelong non-smoker (<100 cigarettes smoked in lifetime); 2, current smoker (includes daily smokers non-daily/occasional smokers); 3, current reformed smoker for >15 years; 4, current reformed smoker for ≤15 years. TPH, PD-L1 high-expression/TMB-high; TPL, PD-L1 low-expression/TMB-low; LUSC, lung squamous cell carcinoma.

Figure 2

Comparison of immune microenvironment, DEGs, proteins, and enrichment signaling pathways and survival between TPH and TPL groups. (A) Comparison of immune cell fractions between TPH and TPL groups by CIBERSORT algorithm. TPH group had more abundant macrophages M1 (P<0.01), resting NK cells (P<0.05), activated memory CD4 T cells (P<0.001), and CD8 T cells (P<0.05), but fewer memory B cells (P<0.05), activated dendritic cells (P<0.05), and neutrophils (P<0.05). (B) Comparison of enrichment levels of 29 immune cells, immune-related pathways, and activity of immune-related function in the 2 groups by ssGSEA. Compared with the TPL group, the TPH group harbored more immune cells, immune-related pathways, and activities of immune-related function. (C) DEGs between TPH and TPL groups. (D) DEPs between TPH and TPL groups. (E) Barplot of enrichment GO terms between TPH and TPL groups. (F) Barplot of enrichment KEGG signaling pathways between TPH and TPL groups. DEGs, differentially expressed genes; TPH, PD-L1 high-expression/TMB-high; TPL, PD-L1 low-expression/TMB-low; NK, natural killer; ssGSEA, single-sample gene set enrichment analysis; DEPs, differentially expressed proteins; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes. Each boxplot is labelled with asterisks indicating the P values (*P<0.05, **P<0.01, ***P<0.001).

We then quantified the enrichment levels of 29 immune cells, immune-related pathways, and activity of immune-related function in the two groups by ssGSEA score. We discovered that compared with the TPL group, the TPH group harbored more immune cells, immune-related pathways, and activities of immune-related function, which might further predict favorable adaptive immune response and clinical outcome when receiving immunotherapy (). Both median OS (mOS) and median progression-free survival (mPFS) were worse in the TPL group than in the TPH group. The mOS in the TPH group was 61.9 months, while that in the TPL group was 19.6 months. However, TPH and TPL both still have some limitations for prognosis (Figure S1).

DEGs and differentially expressed proteins (DEPs) between the 2 groups

The DEGs and DEPs between the TPH and TPL groups are shown in , respectively.

GO and KEGG enrichment analysis of signaling pathways between the 2 groups

Compared with the TPL group, GO analysis indicated that TPH mainly enriched signaling pathways such as regulation of reactive oxygen species (ROS) metabolic process, intrinsic apoptotic signaling pathway, aging, and response to metal ion (BP) (). The KEGG analysis indicated that TPH mainly enriched signaling pathways such as cellular senescence, FoxO signaling, and proteoglycans in cancer ().

Identification of JAK2 and FOXM1 as differentially expressed genes and proteins between the two groups

The ‘Limma’ package with multiple testing correction methods was used to identify DEGs between TPH and TPL groups (), respectively. Finally, JAK2 and FOXM1 were detected as the overlapping DEGs between the 2 omics.

Enrichment-pathway and PPI network analysis of JAK2 and FOXM1

The GO and KEGG biological pathway enrichment analyses demonstrated that signaling pathways JAK2 was involved in included ErbB, VEGF, Rap1, TNF, MAPK, JAK-STAT, and the AGE-RAGE signaling pathway in diabetic complications, EGFR tyrosine kinase inhibitor (TKI) resistance, platinum drug resistance, PI3K-AKT, PD-L1 expression and PD-1 checkpoint pathway in cancer, and FoxO and HIF-1 signaling pathway. In addition to some of the same signaling pathways of JAK2 such as JAK-STAT, AGE-RAGE signaling pathway in diabetic complications, EGFR TKI resistance, platinum drug resistance, PI3K-AKT, PD-L1 expression and PD-1 checkpoint pathway in cancer, and FoxO and HIF-1 signaling pathway, FOXM participation also included the p53 and cell cycle signaling pathway (). Based on the STRING database, about 20 proteins were predicted to directly interact with JAK2 or FOXM1 (). The JAK2 gene is associated with many genes in the STAT family, such as STAT3 and STAT2. Specifically, STAT3 increases tumor cell proliferation (29), survival and invasion and activates tumor-promoting inflammation, but also suppresses anti-tumor immune responses. Therefore, STAT3 is a promising target for cancer therapy. The FOXM1 gene is related to cell-cycle pathways (), which control cell proliferation, and cancer is a disease of inappropriate cell proliferation (30).

Figure 3

PPI network and pathway enrichment analysis for JAK2 and FOXM1. (A) PPI networks of related proteins of JAK2 and FOXM1. (B,C) GO terms and KEGG pathways of JAK2 and FOXM1. PPI, protein-protein interaction; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes

Validation the robustness of risk score as prognostic factors for immunotherapy through internal and external datasets

To verify the robustness of JAK2 and FOXM1 gene expression, we proposed a new risk score by integrating the gene and protein expression of JAK2 and FOXM1 through LUSC patients. Through univariate and multivariate Cox regression analysis it was revealed that the risk score was significantly associated with OS (). The Kaplan-Meier curve demonstrated the survival rate of the high-risk group was significantly lower than that of the low-risk group. Thus, JAK2 and FOXM1 expression might serve as potential prognostic factors to distinguish the OS of patients with LUSC (). The risk score prognosis model achieved AUC =0.709 and showed better prediction of survival than PD-L1 and TMB in regard to overall survival (). The distribution of patient risk scores and survival status were shown in .

Table 2

Univariable and multivariable Cox regression analysis

Variables	Univariable Cox regression analysis					Multivariable Cox regression analysis
	HR	95% CI of HR		P value		HR	95% CI of HR		P value
		Lower	Higher				Lower	Higher
TMB (high vs. low)	0.803	0.600	1.074	0.140		1.047	0.968	1.133	0.247
Gender (male vs. female)	1.095	0.669	1.790	0.717		2.102	0.511	8.645	0.303
Cultural background (white vs. other)	1.031	0.491	2.165	0.935		0.973	0.136	6.936	0.978
Smoking history (yes vs. no)	0.854	0.685	1.065	0.162		0.715	0.382	1.339	0.296
Pulmonary function (normal vs. abnormal)	0.892	0.413	1.925	0.771		0.318	0.070	1.428	0.134
Radiations (yes vs. none)	2.347	0.904	6.091	0.079		0.744	0.148	3.721	0.719
Stage (I/II vs. III/IV)	1.499	0.932	2.409	0.094		3.495	0.856	14.262	0.016
PD-L1 expression (high vs. low)	1.196	0.751	1.906	0.449		1.917	0.662	5.553	0.230
Risk score (high vs. low)	1.924	1.256	2.945	0.002**		6.935	1.579	30.450	0.010*

*, P<0.05; **, P<0.01. CI, confidence interval; HR, hazard ratio; TMB, tumor mutational burden; PD-L1, programmed death-ligand 1.

Figure 4

Validation the robustness of JAK2 and FOXM1 gene expression and protein expression as prognostic factors with internal database (TCGA). (A) The association of OS with different risk groups based on risk score (P=0.037). Vertical hash marks indicate censored data. (B) The association of PFS with different risk groups based on risk score (P=0.027). Vertical hash marks indicate censored data. (C) Comparison of the sensitivity and specificity for the prediction of OS based on the risk score and other clinical parameters. (D) Risk score and the OS status and time of patients. The dotted line in the middle of divides participants into low-risk and high-risk groups. TCGA, The Cancer Genome Atlas; OS, overall survival; PFS, progression-free survival.

Validation of the association of risk score with immunotherapeutic efficacy through external datasets GSE126044 and GSE136961

We used the external datasets GSE126044 and GSE136961 to validate the robustness of risk score in participants who received immunotherapy. The risk score prognosis model achieved AUCs of 0.727 and 0.67 for predicting the response of NSCLC patients treated with anti-PD-1 and anti-PDL1 antibodies, respectively (). In addition, we used the external dataset GSE136961 to validate the robustness of FOXM1 gene expression in participants who received immunotherapy. As shown in , the level of FOXM1 expression was much higher in participants who had durable clinical benefit than those who had not received such benefits with immunotherapy (P=0.043). In GSE126044, the results revealed that the level of both JAK2 and FOXM1 expression was significantly higher in participants who responded to immunotherapy than in those who did not (P=0.038 and P=0.009, respectively) ().

Figure 5

Validation the robustness of JAK2 and FOXM1 gene expression as prognostic factors for immunotherapy with external datasets GSE126044 and GSE136961. (A) The performance of the risk score in predicting ICB response in GSE126044 and GSE136961. (B) Validation with external datasets GSE136961 revealed the level of FOXM1 expression was much higher in patients who had durable clinical benefit than those who did not with immunotherapy (P=0.043). (C) Validation with external datasets GSE126044 suggested the level of JAK2 expression was significantly higher in the patients who responded to immunotherapy than those who did not (P=0.038). (D) Validation with external datasets GSE126044 indicated the level of FOXM1 expression was significantly higher in patients who responded to immunotherapy than those who did not (P=0.009). ICB, immune checkpoint blockade.

Discussion

LUSC accounts for approximately 45% of primary lung cancers in men and 25% in women, which can be divided into central type and peripheral type. Some studies found that PD-L1 expression was more frequent in central type, tumor location that could predict expression status of PD-L1, and could potentially serve as clinical response to immunotherapy. With the advance of ICIs in LUSC, it is urgent to explore effective biomarkers for LUSC patients. PD-1 can be used as a predictor of the efficacy of multiple cancers treated in tumor tissues by immunotherapy, The proportion of PD-L1 antigen-presenting cells had relationship with Clinical response (31). Although PD-L1 and TMB are the most common predictors for the patients with LUSC receiving immunotherapy, they could not meet patient needs due to their individual limitations. We recommend that a preoperative biopsy to detect driver gene mutations is necessary to guide the decision of neoadjuvant therapy.

Studies have revealed that patients with NSCLC who harbor both positive expression of PD-L1 and high TMB demonstrate significantly favorable objective response rate (ORR) and PFS than those with only 1 of or no such variables (14). Thus, we divided the LUSC patients into two participant groups (TPH and TPL) according to their expression of PD-L1 and TMB status. Through the DEG profiles, we identified JAK2 and FOXM1 as the DEGs between these two groups.

Results of GO and KEGG pathway analyses revealed the signaling pathways in which JAK2 or FOXM1 are involved, which was consistent with the role that JAK2 or FOXM1 played in the development of cancer, such as proliferation [JAK-STAT (32-34) and ErbB/PI3K-AKT signaling pathway (35), PI3K-AKT/FoxO (36)], anti-apoptosis (p53 signaling pathway), angiogenesis (VEGF signaling pathway), and immune response [JAK-STAT (37)]. The JAK2 component is an important part of JAK-STAT signaling pathway. Prior studies reported that JAK-STAT signaling mediates almost all immune regulatory processes (38). It has also been shown that FOXM1 could promote the activation, proliferation, and survival of immune cells (39,40). We discovered that levels of JAK2 and FOXM1 expression were correlated with immune infiltration in LUSC, which indicated that JAK2 and FOXM1 not only played a key role in the tumorigenesis, but were also involved in the regulation of immune response. Thus, we speculated that patients with higher levels of JAK2 or FOXM1 expression might gain durable clinical benefits when receiving treatment with ICIs. Further analysis with external datasets GSE126044 and GSE136961 supported our hypothesis. Some of our findings were consistent with those of previous investigations (41,42). For example, Liu et al. (41) revealed that a higher level of JAK2 expression was correlated with the tumor infiltrating lymphocytes, especially infiltrating T cells, in breast cancer. The difference was that the study by Liu et al. showed that expression of JAK2 was an independent indicator of favorable prognosis in patients with breast cancer. Further investigation is required to clarify that particular finding.

Survival analysis was conducted to validate the potential clinical benefits of assessing the risk score to guide ICI strategies, which demonstrated that patients with high risk scores were associated with worse OS and PFS (P=0.0018 and P=0.015, respectively). Some previous studies explained that JAK2 and FOXM1 were involved in promoting the pathogenesis of malignancy (32-39). In addition, the predictive performance of the risk score model on ICB-response indicated that it may be a potential immune-related biomarker for LUSC. Although the risk score can effectively aid in predicting the prognosis of LUSC patients, there is still a lack of clinical trials. If more comprehensive LUSC immunotherapy data emerges in the future, we believe it will be possible to better verify the performance of the risk score. Overall, this was the first study to identify that JAK2 and FOXM1 expression could independently predict the clinical outcome of LUSC patients receiving ICIs. Wu et al. found tumor profile of advanced NSCLC in May 2021, comprehensively characterize the disease characteristics of patients with advanced non-small cell lung cancer from multiple aspects. The results of innovative drug camrelizumab combined with chemotherapy for the first-line treatment of stage III Chinese LUSC patients, which confirmed the excellent efficacy of camrelizumab in the first-line treatment of lung squamous cell carcinoma (43). With the discovery of new drug targets and the continuous emergence of new combination treatment options, our investigation may provide insight into identifying the correct candidates for precision immunotherapy.

The article’s supplementary files as