Introduction

Lung cancer is the most common cause of cancer deaths worldwide, and approximately 80 % of the patients have non-small cell lung cancer (NSCLC). Current mode of lung cancer therapy mainly depends on the TNM disease staging and tumor histological classification, and the treatment-of-choice for early stage NSCLCs constitutes radical surgery.

However, about 30 to 50 % of early stage NSCLC patients (in stages I and II) die within 5 years after radical surgical treatment as a consequence of disease recurrence or metastasis. Results of recent studies have shown some benefit of adjuvant chemotherapy in surgically treated NSCLC patients with 5-year survival advantage of 4–15 % [1]. However, the factors that might help to prospectively identify patients who will most likely benefit from any or specific type of postoperative chemotherapy remain unknown, especially among individuals with stages I–IIA disease.

It has been suggested that evaluation of molecular biomarkers for prediction and prognostics could improve the management of patients with NSCLC [2]. Recent advances in genomics and proteomics of lung cancer have established new candidate markers with potential clinical value. Based on the results of high-throughput microarray and real-time PCR data, several prognostic gene expression signatures have been created to stratify patients for appropriate treatment; however, the proposed gene sets differ significantly from one another. Thus, further studies on gene expression signature as well as verification of previously reported data have to be performed [3].

The quantitative real-time PCR (qRT-PCR) seems to be the most appropriate method for both the microarray-based prognostic gene expression signatures validation and their further practical application in the clinics [4]. Recently, the capabilities of qRT-PCR have been expanded with the development of low density gene expression arrays that allow simultaneous examination of multiple, user-defined genes [5].

In the current study, we evaluated the expression of 23 genes reported previously as putative prognostic markers in NSCLC. The genes used for evaluation were defined according to a number of microarray-based studies (three genes: AGFG1/HRB, STC1, SLC2A1 that overlap between NSCLC prognostic signatures), comprehensive real-time PCR analyses, and/or protein-evaluating assays [6–17]. The selected genes are listed in Table 1. We used low density microarray approach to simultaneously evaluate the expression of all the analyzed genes and to compare gene expression level in tumor and corresponding non-tumor lung tissues. We also analyzed the association between the expression level of the particular genes and the overall survival (OS) and disease-free survival (DFS).

Table 1 List of the genes analyzed in the study
Full size table

Materials and Methods

Patients and samples

The study was performed on 109 pairs of tumor and matched unaffected lung tissue specimens obtained from early stage NSCLC (stages I and II) patients who underwent a curative surgery at the Bialystok Medical University Hospital. All patients gave the written informed consent for specimen collection and clinicopathological data processing. The study was approved by the ethics committee of the university.

All patients underwent complete resection of the tumor and mediastinal lymph nodes and were followed up for at least 3 years or until death. None of them received chemo- or radiotherapy before or after the surgery. The OS was estimated as the time from the date of the surgery to the date of death (event) or of the last control visit (censoring). The DFS was defined as the time from the date of surgery to the date of disease relapse or death (events), whichever occurred first, or to the date of the last visit (censoring). Patients’ clinicopathological characteristics (age, gender, histological type, and pathological TNM staging) are summarized in the “Results” section.

Tissue samples were collected intraoperatively and processed immediately after surgical removal; after the macroscopic visual assessment, the pieces of tumor tissue and unaffected lung tissue from the same lobe or lung of the patient were frozen in liquid nitrogen followed by storage at −80 °C. Prior to processing, the cryo-sections of frozen tissue specimens were stained with hematoxylin–eozyn and evaluated by an experienced pathologists (LC) to confirm the suitability of tumor cell content. Only the tumor samples which contained at least 50 % of tumor cells on a microscopic section as well as unaffected lung tissue samples without malignant cells were used for further processing.

RNA extraction

Total RNA was isolated from tissue specimens using magnetic extraction method. Briefly, about 40–50 mg of frozen tissue was disrupted in 500 μl of Lysis Buffer (Biomerieux, France) with a TissueRupter (Qiagen, Germany) and incubated with Proteinase K (Sigma-Aldrich, Germany) for 2 h in 56 °C. Nucleic acids from deproteinated cell lysates were extracted automatically on EasyMag machine (Biomerieux, France) according to the producer’s protocol. The 100-μl resulting RNA extracts were stored in −80 °C before further processing.

The concentration of the RNA in the extracts was determined on the NanoDrop 2000c spectrophotometer (Thermo Scientific, UK) and about 500 ng of the RNA was transcripted into cDNA in a reaction with High Capacity RNA-to-cDNA Master Mix (Applied Biosystems, Foster City, CA) according the producer’s recommendations.

mRNA expression level

An mRNA expression level of 23 genes was evaluated in the tumor and unaffected lung tissues with the comparative real-time PCR (RT-PCR) method using a TaqMan low density array analysis [18]. Ribosomal 18S RNA gene with a relatively low level of the expression variability in lung tissue [9] was used to normalize for the differences in the input cDNA concentration. For each sample, the amplification of all the transcripts was performed simultaneously in the MicroFluid Cards (Applied Biosystems) that contained manufactory loaded and dried commercially available primers/probe sets for gene expression examination (Assays-on-Demand, Applied Biosystems). The Assay-on-Demand accession numbers are summarized in Supplementary Table 1. Each channel of a card was loaded with 100 μl of the reaction mixture containing 50 μl 2XTaqMan Gene Expression Master Mix (Applied Biosystems) and 20 μl of a cDNA solution (corresponding to 100 ng of total RNA). The amplification was performed on a ABI PRISM 7900HT Sequence Detection System equipped with the SDS v.2.4 software for baseline and Ct calculations. The cycling conditions were as follows: 50 °C for 2 min followed by 95 °C for 10 min hold and 40 cycles of 95 °C for 15 s and 60 °C for 60 s. Each sample was analyzed in duplicate. To minimize run-to-run variations, the same card was used to evaluate gene expression in tumor and matched normal lung tissues from a particular patient. The raw Ct data from each run were exported into the Excel software for further processing. All gene expression measurements were analyzed on the logarithmic scale (log-expression). For each of the 23 analyzed genes, the two measurements were averaged and the ratio, i.e., fold change of the averages for the tumor and normal tissues was computed.

Statistical analysis

The difference in the log-expression between the tumor and normal tissue (log-fold-change) for each gene was analyzed by using a linear mixed effects model (with a random tissue effect) taking into account the correlation between the repeated expression measurements obtained for a single patient [19]. OS and DFS probabilities were estimated by using the Kaplan-Meier method. The median follow-up time was estimated by using the “reverse” Kaplan-Meier method [20]. The influence of clinical factors and gene expression on OS and DFS was analyzed by using the proportional hazards model. First, a model including age, sex, TNM stage, and histological type was fitted to the data. The model was then simplified by retaining only the factors showing a statistically significant effect. Then, for each gene, the model was extended by including the log-fold-change as a covariate. The effect of the covariate was modeled by using fractional polynomials of order two, which allow flexible modeling of the effect [21]. The statistical significance of the effect was then tested by using the likelihood ratio test.

All the applied tests were two-sided. Given that a number of genes were tested, gene-specific p values were corrected for multiple testing by using Hommel’s method [22] with the aim to preserve the overall 5 % level of significance. The analyses were performed by using the SAS 9.2 and STATA 11 statistical software.

Results

Patient characteristics

A total of 109 early stage NSCLC patients aged from 39.8 to 78.1 years (mean 61.8, standard deviation 8.4 years) were included in the study. The majority of the patients (77.1 %) were males. Among the patients, 36 (33.0 %) had lung adenocarcinoma (AdC), 48 (44.0 %) had squamous cell carcinoma (SqCC), and the remaining 25 (23.0 %) had large cell carcinoma (LCC) of the lung. Fifty-three (48.6 %) patients had TNM stage I disease: 18 (16.5 %) with stage IA and 35 (32.1 %) with stage IB. Fifty-six patients (51.4 %) had stage II disease: 25 (22.9 %) with stage IIA and 31 (28.5 %) with stage IIB. The median follow-up time was equal to 48.8 months. During the follow-up, 49 patients had disease recurrence and 36 of them had died.

Differential gene expression between tumor and non-tumor lung tissues

The expression of 23 different genes was analyzed simultaneously in each sample using a TaqMan low density array analysis from Applied Biosystems. We were unable to successfully amplify a number of gene transcripts in some cases possibly due to their tissue level below the detection limit of the assay. As a result, for none of the genes, the log-fold-change could be defined for all the patients. Twenty out of 23 analyzed genes showed statistically significant differences in expression between tumor and non-tumor tissues (Table 2). For 13 genes (ITGB1, ITGB3, CXCL1, CXCL9, CXCL10, CXCL12, CXCR3, CXCR4, TNFA, CHKA, AGFG1, SLC2A1, and CTC1), the expression in the tumor was lower than in the normal sample, and for 7 genes (ITGA5, IL8, IL6, CXCL2, CXCL3, CXCL5, and CXCL11), an increase in expression in tumor was observed (Table 2).

Table 2 Log-fold-changes in gene expression level between tumor and non-tumor lung tissue (the All NSCLC column) and the differences in the log-fold-change between SqCC and LCC (the SqCC vs. LCC column) and between AdC and LCC (the AdC vs. LCC column)
Full size table

The analysis of the effect of age, gender, TNM stage, or histological type on gene expression revealed that only histological type had substantial influence on the log-fold-change. In particular, for all the genes except of CXCL5, CXCL6, CXCR2, IL10, and SLC2A1, the effect of histological type was statistically significant at the 5 % significance level (after adjusting the gene-specific p-values for multiple testing). In general, for the genes for which histological type had a significant effect, the log-fold-change in SqCC was statistically significantly higher than in AdC or LCC. For the latter two types, the log-fold-change was not statistically significantly different, except of genes CXCR3 and TNF, for which the change was higher for AdC than for LCC. The results are summarized in Table 2.

Analysis of overall survival

In the group of 109 analyzed patients, 36 deaths were observed. The median follow-up time was equal to 48.8 months. The left hand side panel of Fig. 1 presents the OS curve for the whole group of analyzed patients.

Fig. 1
figure 1

Left: The estimated overall survival probability. Right: The estimated disease-free survival probability

Full size image

In the proportional hazards model including patients’ age, gender, TNM stage of the disease (I vs. II), and histological type of tumor, only TNM stage was found to be statistically significantly influencing OS (Table 3). In the model simplified by retaining only TNM, the hazard ratio (HR) for TNM II vs. I was estimated to be equal to 7.5 (95 % confidence interval, CI: [3.09, 18.29]; p < 0.001). When the model containing TNM was extended by including the gene-specific log-fold-change for each gene separately as a covariate, only the effect of CXCL5 was found to be statistically significant (multiple-testing-corrected p = 0.04). Five of the genes (CXCL6, CXCR2, CXCR3, IL10, and TNF) were excluded from the survival analyses due to the lack of expression data for more than 50 % of the patients.

Table 3 Multivariable analysis of the prognostic effect of patients’ clinicopathological characteristics on overall survival and disease-free survival (proportional hazards model)
Full size table

The association between the logarithmically transformed OS hazard ratio and log-fold-change in the CXCL5 expression between the tumor and normal lung tissue is shown in Fig. 2 (the panel of the left-hand side). The estimated form of the dependence of the log-HR on the log-fold-change of CXCL5 is presented as a solid line. To provide information about the individual patient values of the log-fold-change, deviations (indicated by the closed circles) from the predicted log-HR curve for individual patients were added to the plot. The sudden drop of the curve at the left hand side is mainly due to a single observation (not shown on the plot) for a TNM IIA patient with the log-fold-change of −0.36, who was alive at the last follow-up visit at around 36 months postsurgery. If we disregard this boundary effect, the plot suggests that the hazard of death continuously decreases with the increasing fold change for the gene. The plot does not suggest any particular threshold value for the log-fold-change.

Fig. 2
figure 2

Left: The estimated form (solid line) of the dependence of the OS log-HR on the log-fold-change of CXCL5. The shaded region indicates the 95 % CI limits. The closed circles indicate deviations from the predicted log-HR curve for individual patients and provide an information about the individual patient values of the log-fold-change. Right: The estimated OS curves for two groups of patients with the log-fold-change below or above the observed median (0.0794). The plots are based on data for 74 patients with available information about the expression of the gene

Full size image

For illustrative purposes, in the right hand side panel of Fig. 2, we present the OS curves for two groups of patients created arbitrarily by considering the log-fold-change below or above the observed median (0.0794). The curves indicate the higher risk of death for the patients with the log-fold-change below 0.0794, i.e., with the gene expression in the tumor tissue smaller than exp(0.0794) = 1.08 times the expression in the normal tissue.

Analysis of disease-free survival

In the group of 109 analyzed patients, there were 49 deaths and disease recurrences. Figure 1 (the panel on the right-hand side) presents the DFS curve for the whole group of patients.

In the proportional hazards model including patients’ age, gender, TNM stage (I vs. II) of the disease, and histological type of tumor, only TNM stage was found to statistically significantly influence DFS (Table 3). In the simplified model, HR for TNM II vs. I was estimated to be equal to 2.68 (95 % CI: [1.47, 4.90]; p = 0.001).

Similar to the results of the OS, a statistically significant effect of the log-fold-change, adjusted for the effect of TNM, was found only for gene CXCL5 (multiple-testing-corrected p = 0.004). The left-hand-side panel of Fig. 3 presents the estimated form (solid line) of the dependence of the log-HR on the log-fold-change for CXCL5. The shape of the curve is very similar to the one shown in Fig. 3 for OS and suggests that the hazard of disease recurrence continuously decreases with the increasing fold change in the gene expression. For illustrative purposes, in the right hand side panel of Fig. 3, we present the DFS curves for two groups of patients created by considering the log-fold-change below or above the observed median (0.0794). The curves indicate the higher risk of disease recurrence for the patients with the log-fold-change in the CXCL5 expression below the median value.

Fig. 3
figure 3

Left: The estimated form (solid line) of the dependence of the DFS log-HR on the log-fold-change of CXCL5. The shaded region indicates the 95 % CI limits. The closed circles indicate deviations from the predicted log-HR curve for individual patients and provide an information about the individual patient values of the log-fold-change. Right: The estimated DFS curves for two groups of patients with the log-fold-change below or above the observed median (0.0794). The plots are based on data for 74 patients with available information about the expression of the gene

Full size image

Discussion

The aim of the study was to evaluate the expression of 23 genes as putative prognostic markers in early stage NSCLC. Twenty-three genes had been selected for the examination. These included three genes evaluated by microarray approach (AGFG1/HRB, STC1, and SLC2A1) that overlap between NSCLC prognostic signatures, CHKA suggested as a new prognostic predictor in early stage NSCLC: three genes encoded for intergins α5, β3, and β1 which had been prognostic in the comprehensive analysis of 15 integrin genes by Dingemans et al. [17]: and 16 genes for cytokines or their receptors that have been previously implicated in lung cancer biology [11, 23, 24] and/or demonstrated to be a part of a prognostic signature [12, 15]. A relative quantification RT-PCR method was applied as a research method. The method allows estimation of the fold changes in gene expression between tumor and non-tumor lung tissues from the same patient and is widely used for microarray data verification [4, 25]. To minimize technical, run-to-run variations, a TaqMan low density array analysis that allows a simultaneous examination of the expression of dozens of different genes in each sample was used [5].

In our study, the expression of the majority of the analyzed genes (20 out of 23) in tumor and matched normal lung tissues differed statistically significantly after applying a multiple-testing correction. Among the differentially expressed genes, there were integrin genes ITGA5, ITGB3, and ITGB1; genes for cytokines or their receptors CXCL1, CXCL2, CXCL3, CXCL5, CXCL9, CXCL10, CXCL11, CXCL12, CXCR3, CXCR4, IL6, IL8, and TNF; a choline kinase alpha gene CHKA; as well as genes AGFG1 and CTC1 encoded for nucleoporin-like protein and stanniocalcin-1 (SLC2A1), respectively. On the other hand, tumor-associated alterations of the expression level of 3 genes (CXCL6, CXCR2, and IL10) did not reach statistical significance.

Changes in gene expression were more pronounced in SqCC as compared with AdC or LCC. On the other hand, no significant differences in gene expression in tumor related to matched non-tumor lung tissue were found between AdC and LCC samples, except for the activity of TNF and CXCR3. Based on these results, the analyzed genes seem to be involved in the process of NSCLC evolution, especially in the case of SqCC of the lung, but further examinations are necessary to clarify their contribution to tumor growth and progression.

However, we failed to evaluate the expression of particular genes in some cases possibly due to their tissue level below the detection limit of the assay. Similar results were reported by Seike et al. who had been unable to detect TNF-α and interleukin (IL)-10 expression in a significant percentage of examined tumor and non-tumor tissues [15]. Of interest, the lack of IL-10 expression in specimens from stage I NSCLC patients was found to be associated with significantly worse outcome [16].

One of the major goals of our study was to evaluate the prognostic significance of the tumor-associated alterations in the expression level of the 23 selected genes in stages I and II NSCLCs. In the analyses of OS and DFS, five genes (CXCL6, CXCR2, CXCR3, IL10, and TNF) had to be excluded from the evaluations due to lack of the expression data for more than 50 % of the patients. For the remaining 18 genes, only the expression of one gene (CXCL5) was statistically significantly associated with patients’ overall and disease-free survival after applying a multiple-testing correction (p = 0.04 and p = 0.004, respectively). Upregulation of CXCL5 expression in tumor was found to be a favorable prognostic factor for both OS and DFS, and none of the analyzed clinicopathologic characteristics (patients’ age, gender, stage, or tumor histological type) influenced these associations. Gene product CXCL5 (also known as ENA-78/SCYB-5) is a CXC-type chemo-attractive cytokine that was primary described as a powerful attractant for granulocytes and an inflammatory mediator [26]. More recently, CXCL5 contribution to cancer growth and metastasis has been demonstrated [24, 27]. Specifically, CXCL5 has been found to be important in mediating tumor-associated angiogenesis [12, 28] as well as in promoting growth, migration, and invasion of tumorigenic cell lines derived from NSCLC [29] and other cancers [30–32]. The CXCL5 expression has been shown to be elevated at both mRNA and protein levels and affect patients’ survival in a number of human tumors, including in pancreatic [31, 33], colorectal [34], breast [35], ovarian [36], or prostate [30] carcinomas. However, survival data are limited and inconsistent [31, 37–39]. In the study of White et al., risk of the recurrence after lung cancer resection in NSCLC patients’ was found to be associated with enhanced expression of the sum angiogenic CXC chemokines [12].

In contrast to CXCL5, none of the remaining genes examined in our study was found to be associated with patients’ postoperative outcome. It has been repeatedly emphasized that the direct comparison of gene expression data is difficult because potential variations in applied analytical methods and patient selection algorithms may greatly compromise the reproducibility of the results [3, 9, 40, 41]. The analysis of gene expression signatures reported for a defined tumor phenotype or patients’ prognosis reveals very minimal overlapping between them in terms of gene identity. Although a particular signature validation on several independent sample sets generally demonstrated high degree of reproducibility [42], cross-examination of public available gene expression data by different statistical methods generated different results [40, 43]. RT-PCR-based gene expression studies usually used for microarray-based data validation often failed to confirm the examined effect of the gene. Consequently, in the recent critical review of 16 published studies reporting the development of the gene expression-based prognostic signatures in NSCLC, Subramanian and Simon found little evidence that any of the signatures are ready for clinical application [3]. Besides, when effects of multiple genes are investigated, high probability of false-positive finding exists if multiple-testing correction methods are not applied. To illustrate, in the study by Lau et al., who examined 158 putative prognostic genes in a cohort of 147 surgically treated NSCLC patients with RT-PCR method, 24 genes influenced patients’ survival, but the effect of only five of them remained significant after the false discovery rate adjustment for multiple testing [44]. Thus, multiple-testing corrections are very important to properly control the overall probability of a false-positive finding.

In conclusion, in our study, 18 out of the 23 genes previously found to be relevant for NSCLC formation and/or early staged NSCLC patients’ postoperative outcome demonstrated differential expression in tumor and matched unaffected lung tissues. Changes in expression were more pronounced in squamous carcinomas as compared to tumors of non-squamous histology, pointing out the possible contribution of the genes to SqCC carcinogenesis. Among the 23 genes previously suggested to be relevant for early staged NSCLC patients’ postoperative outcome, only CXCL5 showed a statistically significant prognostic effect.