Neutrophils correlate with hypoxia microenvironment and promote progression of non-small-cell lung cancer.

Hypoxia, a strong and selective pressure, has been involved in invasion, metastasis, and angiogenesis of tumor cells. Our study performed the transcriptome profiles of 666 non-small-cell lung cancer (NSCLC) patients. Various bioinformatic approaches were combined to evaluate the immune cell infiltration in the high hypoxia risk patients. In addition, in vitro experiments were performed to assess the effects of tumor-associated neutrophils (TANs) on NSCLC cells proliferation, migration and invasion and to reveal the underlying mechanisms. We divided NSCLC into two groups (Cluster1/2) based on the expression profiles of hypoxia-associated genes. Compared with the Cluster1 subgroup, the Cluster2 had a worse prognosis. Significant enrichment analysis revealed that PI3K/AKT/mTOR signaling pathway and TANs were highly related to hypoxia microenvironment. Eleven hypoxia-related genes (FBP1, NDST2, ADM, LDHA, DDIT4, EXT1, BCAN, IGFBP1, PDGFB, AKAP12, and CDKN3) were scored by LASSO COX regression to yield risk scores, and we revealed a significant difference in overall survival (OS) between the low- and high-risk groups. Mechanistically, CXCL6 in hypoxic cancer cells promoted the migration of TANs in vitro, and in turn promote NSCLC cells proliferation, migration and invasion. In summary, this study revealed a 11-hypoxia gene signature that predicted OS of NSCLC patients, and improved our understanding of the role of TANs in hypoxia microenvironment.

Introduction

Lung cancer is the leading cause of cancer death worldwide. Non-small cell lung cancer (NSCLC) comprises up to 80% of all lung cancer cases and has a poor prognosis [1]. Compared to the steadily increasing survival rates for most cancers, the survival rates for NSCLC are still unsatisfactory, mostly caused by resistance to anticancer treatment and metastasis [2]. Therefore, it is crucial to discover neo-markers to deliver greater prognostic value and understand the mechanisms of progression for NSCLC.

It has been demonstrated that tumor microenvironment (TME) plays a critical role in supporting the progression of NSCLC. The TME in NSCLC is composed of multiple components, including inflammatory cell infiltration, blood vessels, soluble factors, and hypoxic condition [3]. The hypoxic zone in solid tumors is invaded by a great many of immune-suppressant cells, such as tumor-associated neutrophils (TANs), myeloid derived suppressor cells (MDSCs), T-regulatory (Treg) cells and tumor-associated macrophages (TAMs) [4]. Emerging evidence suggests a strong link between the hypoxic TME and infiltrating immune cells. In hypoxic TME, lung cancer-derived exosomes reprogramed macrophages, which in turn promoted metastasis of lung cancer [5]. Hypoxia selectively upregulates PD-L1 on MDSCs via HIF-1α and blockade of PD-L1 abrogates MDSCs function [6]. As an important element of TME in NSCLC, TANs play crucial roles in supporting tumor growth. A pan-cancer study including NSCLC showed that TANs infiltration in the TME is a leading predictor of poor outcome [7]. Similarly, in NSCLC, a high ratio of neutrophils to lymphocytes is linked to a poor outcome [8]. However, the mechanisms involved in the association between TANs and hypoxic TME in the progression of cancer are still to a great extent unknown.

Hypoxia is a complex process that involves several hundred molecules [9]. Thus, models that integrate several hypoxia-related genes can improve the accuracy of prognostic predictions compared to a single gene. However, some of these characteristics were often uncertain and resulted from an unspecified genetic background. This study aimed to exploit machine learning methods to construct a robust prognostic gene panel. Based on the increased infiltration of TANs in high hypoxia risk tumors, we also investigated the specific mechanisms by which NSCLC cells recruit TANs under hypoxic conditions.

Materials and methods

Acquisition of lung cancer datasets

The set of sequence-based mRNA expression data (RNA-seq data) of NSCLC was downloaded from The Cancer Genome Atlas (TCGA) (https://cancergenome.nih.gov/). Another three gene expression arrays of human NSCLC datasets (GSE50081, GSE30219, and GSE37745) were downloaded from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo). Moreover, clinicopathological information and survival data from four cohorts (TCGA, GSE50081, GSE30219, and GSE37745) were also obtained for further prognostic analysis.

Identification of hypoxia status and assessing the immune infiltrates

To explore the function of hypoxia related genes in NSCLC, R package ‘ConsensusClusterPlus’ was applied to classified NSCLC patients into various clusters.

Infiltration levels of the abundances of immune cells in NSCLC were evaluated by CIBERSORT [10]. To avoid calculation errors, an immune infiltration estimation was performed using the ‘xCell’ R package [11].

Gene set variation analysis (GSVA)

The actin pathway, MYC targets V1, IL-6/JAK2/STAT3 pathway, PI3K/Akt/MTOR signaling, hypoxia, epithelial-mesenchymal transition (EMT), and angiogenesis gene set files obtained from the ‘Molecular Signatures Database’ were utilized for GSVA using ‘GSVA’ packages for R [12].

Establishment of prognostic model, receiver operating characteristic (ROC) curve and nomogram

A univariate Cox regression analysis was firstly performed on all of candidate genes to calculate the relationship between the expression of each gene and OS. The selected key genes were then further confirmed by the Lasso regression. The prognosis risk score was established, and a risk score based on the normalized expression data was generated. Patients were then divided into high- and low- risk groups according to the median cutoff of the prognosis risk score. The prognostic performance was evaluated by using time-dependent ROC curve analysis to estimate the predictive sensitivity and accuracy of our model. Using the survival and the rms package for R, age, gender, risk score, and stage were used to build a nomogram.

Cell lines, cell culture and hypoxic exposure

The NSCLC cells (H-1299 and A549) were acquired from the ATCC. For normoxia experiments, the cells were cultured at 37°C and 5% CO2. For hypoxia experiments, the cells were cultured in a hypoxic chamber under hypoxic conditions (1% O2).

Neutrophils from peripheral blood

Briefly, blood was obtained from NSCLC patients and erythrocytes were excised using the dextran sedimentation method. Using density centrifugation, neutrophils were clarified and re-suspended in RPMI 1640 containing 2% fetal bovine serum.

Neutrophil from NSCLC tissues

Briefly, fresh NSCLC tissues were sliced into small pieces and digested in RPMI-1640 supplemented at 37°C for 30 min with 0.002% DNase I and 0.05% collagenase IV. The leukocytes were harvested and CD66b+ neutrophils were isolated using the EasySep PE Selection Kit (Stemcell Technologies) according to the manufacturer’s instructions.

Cell proliferation and wound healing assay

To detect proliferation, 500 cells were plated on 6-well plates and maintained for 2 weeks. A 4% paraformaldehyde assay were used to stained with crystal violet.

Cells were loaded into 6-well plates for 48 h, and a wound was swabbed with a pipette tip. After incubation with serum-free medium for 48 h at 37°C, at the wound front migrating cells were pictured. The wound area was calculated with Image J software.

Assessment of mRNA expression

Total RNA from NSCLC cells was extracted according to the manufacturer’s instructions by the Trizol Reagent. Subsequently, complementary DNA (cDNA) was obtained via reverse transcription for RT-PCR according to the SYBR-Green PCR Mix.

Primers sequences used in our study were as follows:

CXCL1 forward 5′-CCCAAGAACATCCAAAGTGTG-3′,

reverse 5′-CATTCTTGAGTGTGGCTATGAC-3′;

CXCL2 forward 5′-CTCAAGAACATCCAAAGTGTG-3′,

reverse 5′-ATTCTTGAGTGTGGCTATGAC-3′;

CXCL5 forward 5′-GGAAGGAAATTTGTCTTGATCC-3′,

reverse 5′-TTTCCTTGTTTCCACCGTC-3′;

CXCL6 forward 5′-CGTTACGCTGAGAGTAAACC-3′,

reverse 5′-GTTCTTCAGGGAGGCTACC-3′;

CXCL8 forward 5′-ACTCCAAACCTTTCCACCC-3′,

reverse 5′-CAATAATTTCTGTGTTGGCGC-3′;

CCL3 forward 5′-CGGTGTCATCTTCCTAACC-3′,

reverse 5′-TCGCTGACATATTTCTGGAC-3′.

Plasmid constructs and transfections

For plasmid transfections, NSCLC cells were transfected with vectors of CXCR6-siRNA with the following target sequences: 5ʹ-ACGGATGTGTTCCTGGTGA −3ʹ or negative control RNA according to the manufacturer’s instructions using Lipofectamine 2000 (Roche, USA). Briefly, prior to treatment, Lipofectamine 2000/siRNA complexes were produced in reduced serum medium, OptiMEM (Invitrogen) at the recommended ratio of 20 pmol siRNA per 1ul Lipofectamine 2000. NSCLC cells were then treated with Lipofectamine 2000/siRNA complexes for 6 h before being replaced with RPMI 1640 containing 10% fetal bovine serum. NSCLC cells were also treated with Lipofectamine 2000/NC siRNA complexes as above descried.

Enzyme-linked immunosorbent assay (ELISA)

According to the manufacturer’s instructions, the level of CXCL6 was measured using the Human CXCL6 ELISA Kit (Abnova). A curve of absorbance versus the concentration of CXCL6 in the standard wells was generated.

Neutrophil chemotaxis assay

Neutrophil chemotaxis was assayed in a 3 μm Transwell system. Briefly, hypoxia or normoxia conditioned medium (CM) from NSCLC cells, CXCL6 in RPMI 1640 containing 2% FBS (Gibco) or hypoxia CM with CXCL6 neutralizing antibody was added to the lower wells. Neutrophils migrating into the lower chamber were gathered and numbered in the chambers.

Transwell invasion assay

According to the manufacturer’s instructions, cells were placed in Boyden chambers coated with Matrigel diluted 1:8 with DMEM (Gibco) on the submembrane surface. After 48 hours, the aforementioned cells were immobilized in 4% paraformaldehyde and dyed with 0.5% crystal violet, and were evaluated under microscope by counting cells in five different areas per each condition.

Statistical analysis

Intergroup comparisons were performed by one-way ANOVA or Student ‘s test. OS of the patients was determined by the Kaplan–Meier analysis and compared by the log-rank tests. The Cox proportional hazards regression model was used for univariate and multivariate survival analyses. LASSO analyses were applied to reduce the complexity of the model and select the significant hypoxia-related genes to construct the prognostic model. A value of P < 0.05 was regarded as statistically significant.

Results

Lung cancer patient characteristics and hypoxia-related gene identification

Three GEO databases with available OS data and clinical information (GSE50081, GSE30219, and GSE37745) were enrolled into one meta-cohort contained 666 lung cancer patients (Table S1). To classify the hypoxia status of lung cancer samples, we focused on 205 hypoxia-related genes that were obtained from previous research and Gene Set Enrichment Analysis (hallmark-hypoxia).

Then, we applied unsupervised clustering methods to classify 666 tumor patients into different molecular subgroups through 205 hypoxia related genes. Two distant subgroups, referred to as Cluster 1 and Cluster 2, were finalized (Figure 1a and 1B). The Kaplan-Meier curves showed significant differences in OS between the two clusters by log-rank test (Figure 1c). We further selected the GSE39582 datasets, which presents the most complete patient information, featuring the clinical discrepancies between these subgroups. Pathways involved in epithelial-mesenchymal transition (EMT), IL-6/JAK2/STAT3, PI3K/Akt/MTOR, Hypoxia, MYC targets V1 and Angiogenesis were activated in Cluster 2 (worse survival), whereas Cluster 1 (favorable survival), revealed enrichment of pathways correlated to tumor progression, including EMT, hypoxia, and angiogenesis (Figure. S1A). Further survival analysis showed that two gene sets (Hypoxia and PI3K/Akt/MTOR signaling) were associated worse survival (Figure 1d).

Figure 1.

Lung cancer patient characteristics and hypoxia-related gene identification. (A) Consensus matrices of NSCLC patients for k=2 based on 205-hypoxia-related genes in three GEO cohort; (B) NSCLC cases are divided into two subtypes based on unsupervised analysis; (C) Differences in patient overall survival with two clusters; (D) Cox proportional hazard regression of survival months and survival status were performed using gene set enrichment scores for the six hypoxia-associated gene signatures. Log-rank test was used to determine significant p values

The immune cell infiltration between high and low hypoxia risk NSCLC patients

Growing evidence shown that the hypoxic microenvironment may protect tumors from anti-tumor immune responses by promoting immune escape and suppressing anti-tumor immune cells. Using the CIBERSORT method, we assessed differences of 22 immune cell infiltration between Cluster 1 and Cluster 2 patients. Figure 2a and S2A summarizes the results obtained from 666 lung cancer patients in GEO databases. The proportion of immunosuppressive cells was significantly higher in Cluster 2 with high hypoxia risk (e.g., macrophages M2 and neutrophils), and NK cells resting (Figures 2B). Although Tregs cells did not differ between high and low hypoxic risk tumors, inactivated NK cells and immunosuppressive cells may drive the immunosuppressive microenvironment. Transcriptome file of three GEO databases was applied on the ESTIMATE method to estimate the infiltration of different immune cells. Cluster 2 with high hypoxia risk was significantly enriched in immunosuppressive immune cells, such as macrophages M2 and neutrophils as mentioned above (Figure 2c). We determined if any of the above immune cells with differential expression correlated with patient survival. This analysis showed that only neutrophils were associated with worse survival (Figure 2d). Therefore, targeting hypoxia may have significant clinical implications in improving immunotherapy.

Figure 2.

Immune Landscape Between Low and High Hypoxia Risk NSCLC Patients. (A) and (B) The 22 immune cells proportions obtained from 666 lung cancer patients in GEO databases; (C) Rows of the heatmap show expression of TME-infiltrating cell signatures calculated by xCell; (D) Kaplan-Meier plots of infiltrating immune cells with differential expression as described above

Validation of hypoxia-related gene set in the TCGA cohort

The hypoxia-related gene set was further validated by the TCGA cohort, and we condensed into two distinct hypoxia-related subtypes by applying consensus clustering analysis (Figure 3a). Patients in Cluster1 had a higher OS, while patients in Cluster2 yielded a poorer prognosis outcome (Figure 3b). All these genes were subjected to Cox regression analysis. A total of 60 genes were markedly correlated with the OS of TCGA, and LASSO Cox regression was used to filter out less relevant prognostic genes (Fig. S2B). The optimal gene signature consisting of 11 prognostic genes were identified (Fig. S2C). Using the corresponding coefficients of 11 prognostic genes, the individual-level risk score formula for each patient was calculated as following: risk score = – 0.24× FBP1 – 0.90× NDST2 + 0.13× ADM + 1.02× LDHA + 0.05× DDIT4 + 0.82× EXT1 + 0.06× BCAN + 0.12× IGFBP1 + 0.20× PDGFB + 0.22× AKAP12 + 0.26× CDKN3. As shown in Figure 3c, the heatmap revealed that the expression level of 11 genes (FBP1, NDST2, ADM, LDHA, DDIT4, EXT1, BCAN, IGFBP1, PDGFB, AKAP12, and CDKN3) varied accompanying with the higher risk scores. In addition, we calculated the relationship between risk score and cancer-related mortality, and the patients who had higher risk scores were supposed to have a shorter OS than lower risk group, suggesting that these 11 genes might be considered as risk genes in NSCLC (Figure 3d). Time-dependent ROC curve showed that our 11 prognostic gene model exhibited good performance and accuracy in prognosis prediction either for 1-year (AUC = 0.738), 3-year (AUC = 0.702) and 5-year (AUC = 0.666) OS of NSCLC patients (Figure 3e). We generated a nomogram with a concordance index of 0.733 to predict the probability of 5 and 10-year OS, by incorporating the 11-hypoxia gene signature, age, gender, and TNM stage (figure 3f). The calibration curve showed that the nomogram performed a high level of accuracy (Figure 3g).

Figure 3.

Validation of hypoxia-related gene set in the TCGA cohort. (A) Consensus matrices of NSCLC patients for k=2 based on 205-hypoxia-related genes in the TCGA cohort; (B) Differences in patient overall survival with two clusters; (C) Establishment of a prognosis-predictive model dividing patients into high and low risk groups; (D) Differences in patient overall survival with high and low risk groups; (E) Time-dependent ROC curves for the risk score in the TCGA dataset for predicting 1, 3, and 5-year OS; (F) Nomogram based on risk score, age, gender, and stage; (G) Calibration plots of the nomogram for predicting the probability of OS at 5 years in the TCGA dataset

CXCL6 is the key chemokine for the recruitment of TANs in hypoxic NSCLC cells

Based on the aforementioned important role played by TANs infiltration in hypoxic TME, we further examined the effect of conditioned media on neutrophils in hypoxia-treated NSCLC cells. As shown in Figure 4a, CM from hypoxia-induced A549 and H-1299 cells recruited more TANs derived from NSCLC tissues than CM from normoxic cells. To further screen the possible chemokines, candidate chemokines were identified through RT-qPCR. Hypoxia facilitated NSCLC cells to generate more neutrophil chemokines, and CXCL6 increased the most among all the increased chemokines (CXCL1, CXCL2, CXCL5, CXCL6, CXCL8, and CCL3) as detected by RT-PCR (Figure 4b). Consistently, ELISA revealed that the expression of CXCL6 was significantly increased under hypoxia condition (Figure 4c). Transwell migration assays were performed to determine that migration of TANs enhanced by hypoxic NSCLC cells was almost identical to that of 700 pg/ml CXCL6 (Figure 4d). In agreement, migration of TANs of hypoxic NSCLC was repressed by a CXCL6 antibody. Furthermore, obstruction of the CXCL6 receptor CXCR6 (siCXCR6) in TANs markedly compromised the migration of TANs. Similar results were obtained by the exploration of TANs derived from peripheral blood (Figure. S3A-D).

Figure 4.

CXCL6 is the critical chemokine induced by hypoxic NSCLC cell to recruit TANs derived from NSCLC tissues. (A) Quantification of neutrophil migration as assessed by transwell assays; (B) and (C) Expression of CXCL6 in hypoxic or normoxic NSCLC cells was examined by real-time PCR and ELISA; (D) Quantification of neutrophil migration as assessed by transwell assays

TANs promote NSCLC cells proliferation, migration and invasion

To examine the effect of neutrophils on NSCLC cells biology, we cocultured NSCLC cell lines with TANs derived from NSCLC tissues. Proliferation analysis demonstrated that TANs promoted NSCLC cells growth in vitro (Figure 5a). Moreover, compared with control, TANs significantly enhanced the migration and invasion abilities of A549 and H-1299 cells (Figure 5b and c). Similar results were obtained by the exploration of TANs derived from peripheral blood (Figure. S4A-C).

Figure 5.

TANs derived from NSCLC tissues promote NSCLC cells proliferation, migration and invasion. (A), (B) and (C): NSCLC cells cocultured with TANs or alone were subjected to colony formation, wound healing, and transwell invasion assays; (D) Schematic illustration of the crosstalk between CXCL6-overexpressing NSCLC cells and TANs in the TME

In this study, we verified hypoxia upregulated neutrophils recruitment and constructed a novel hypoxia associated gene signature that could improve the individualized prognosis prediction in NSCLC. Moreover, hypoxia increased CXCL16 production and attracting TANs. CXCL6-mediated TANs infiltration contributes to hypoxic-related cancer progression (Figure 5d).

DISCUSSION

Given the wide variation in prognostic outcomes of NSCLC, it is of great importance to build a robust classifier to segment patients with different risks and prognosis, which is essential to maximize the benefits of accurate assessment, individual therapy and timely long-term follow-up. Numerous data indicate that hypoxia and TME are involved in processes that promote the development of tumor cells and a more malignant phenotype [4,13]. Consistent with criteria prognostic parameters like nodal status, tumor stage and tumor grade, hypoxia and immune cells infiltrating TME have been proposed as prognostic factors for patient outcome [14,15]. Although methods such as PET imaging and immunohistochemical marker expression assays have been initially implemented to detect the degree of hypoxia in patients’ tumors, the specific identification methods are still fraught with many unknowns. In this study, the integrated mining of transcriptional profiles and microenvironment features was intended to establish a tool to solve this vital clinical issue.

In this study, we profiled the mRNA expression of 205 hypoxia-associated genes in three GEO databases. Two NSCLC subgroups Cluster 1/2 were identified by consistent clustering on the basis of hypoxia-related genes. The Cluster 1/2 subgroup not only affects disease prognosis and key signaling pathways, but it influences immune cells infiltration. Bioinformatic approach uncovered that hypoxia condition mainly associated with the PI3K/Akt/MTOR signaling pathway, the JAK2/STAT3 signaling pathway, the MYC targets V1 signaling pathway, and angiogenesis. Interestingly, a number of studies have found that hypoxia plays a role in the induction of tumor cell growth and metastasis by PI3K/Akt signaling and MYC target V1 signaling [16-18]. There is also much evidence to support the interaction between hypoxia and angiogenesis [19-21]. Moreover, we used Cox and LASSO regression to develop an 11-gene prognostic signature for TCGA cohort. The risk scores were derived by integrating mRNA expression levels and risk coefficients for 11 hypoxia-related genes, and significantly classified prognosis in NSCLC patients into low and high-risk clusters. There are a number of clinically available multigene based risk models that can predict the prognosis of NSCLC. Zhu and colleagues recently reported a 22-autophagy-related signature based on overall survival in patients with lung adenocarcinoma [22]. Recently, Wu et al. reported a seven long non‐coding RNAs prognostic model to predict OS in NSCLC patients. So far, there were limited prognostic models based on hypoxia-associated genes for patients with NSCLC. In an effort to exploit hypoxia-induced epithelial-mesenchymal transition gene signatures associated with clinical outcomes in NSCLC, Gao et al. demonstrated a 17 gene prognostic panel for NSCLC [23]. Our model shows good accuracy and stability in clinical outcome prediction either for 1-year (AUC = 0.738), 3-year (AUC = 0.702) and 5-year (AUC = 0.666) OS of NSCLC patients. Our data also revealed that all 11 genes in our model were associated with poor prognosis in NSCLC patients. These discovered hypoxia-associated genes have been repeatedly reported in cancer [24-27]. Previous studies reported that hypoxia-induced GBE1 expression promotes tumor progression by repressing FBP1 activities in lung adenocarcinoma [28]. In contrast, in the present study, we found a positive coefficient for FBP1, indicating that FBP1 was considered as risk gene in NSCLC. The underlying molecular mechanisms deserve further investigations.

Accumulating evidence indicates that hypoxia is an important feature of TME and that it can drive progression and metastasis by facilitating suppressive cells (TANs, Tregs, and TAMs), and producing immunosuppressive molecules. CIBERSORT revealed that patients in Cluster 2 with high hypoxia risk had significantly more abundant infiltrative TANs and M2 macrophages phenotype. Furthermore, TANs were associated with a poor prognosis, suggesting that our hypoxia model may predict the immune microenvironment. Hypoxia not only tightly regulates the production of specific chemokines, it also controls their action by regulating their receptors [29,30]. Chemokines play an important role in regulating tumor immunity [31]. Our findings therefore provide insight into the underlying mechanisms of recruitment-related chemokines in TANs infiltration, and found that the production of CXCL6 increased in NSCLC cells under hypoxic condition, and the blockade of CXCL6 almost halted TANs migration, suggesting that the recruitment of TANs is mediated mainly by CXCL6. CXCL6 induced proliferation and metastasis of lung cancer cell lines was confirmed in other studies, and the role of CXCL6 in the TME of NSCLC is unclear [32,33]. Our study broadened the understanding of CXCL6 in NSCLC progression. As previously mentioned, multiple studies have noted that TANs appear to mostly develop a pro-tumorigenic phenotype [34,35]. In this study, we revealed TANs enhanced the proliferative and invasive ability of NSCLC cells in vitro. Experimental studies are needed to further clarify the molecular mechanisms underlying the TANs related pro-tumorigenic phenotype in NSCLC.

Conclusions

In conclusion, we identified a hypoxia model based on 11-hypoxia gene signature, and investigated the association between hypoxic condition and infiltrating immune cells, especially TANs in NSCLC patients. Then, we confirmed that CXCL6 is a strong chemotactic cytokine for TANs infiltration in NSCLC and observed that the role of TANs on NSCLC proliferation, migration and invasion. Thus, our study provided more comprehensive insight into how hypoxia status influence on prognosis and the TME, and may help clinicians develop individualized hypoxia-targeted therapies in NSCLC.

Click here for additional data file.