Evaluation of MicroRNA Expression Profiles That May Predict Recurrence of Localized Stage I Non–Small Cell Lung Cancer after Surgical Resection

Lung cancer is the leading cause of death among those afflicted with malignancies in the United States. Surgical resection is the primary treatment for the most common (70%; ref. 1), non–small cell variant of lung cancer (NSCLC), but the cure rate remains disappointing. Even for those in the earliest stage of the disease (stage IA; ref. 2), the 5-year survival rate is only 73% (2). Although most studies have not shown any significant benefit from adjuvant therapies for early-stage lung cancer (3), molecular biological investigations are under way to supplement clinical and pathologic observations identify cases with a bad prognosis so as to administer therapy only to patients at a high risk of recurrence (4).

MicroRNAs (miRNA) are small noncoding RNAs that regulate protein expression. The miRBase miRNA registry (13.0; ref. 5) includes 885 human miRNAs. Alterations in miRNA expression have been seen in several forms of cancer, including those of the lung (6). Three previous studies have suggested that miRNA profiling can be used for prognostication in lung cancer (6–8).

The expression of 752 miRNAs in cancer tissue resected from 77 cases of stage I NSCLC was examined in this study. A comparison between the groups of cases that had recurrence after surgery and those that did not shows that the two case groups have distinct miRNA expression profiles. That miRNA expression can be used to predict recurrence is shown using two different classifier methods in leave-one-out and Monte Carlo cross-validation analyses.

This study was approved by an institutional review board. Eighty pathologic stage I NSCLC cases resected during 2000 to 2008 at Roswell Park Cancer Institute, with approximately half known to have a recurrence, were chosen for this study. Clinical data was acquired from Roswell Park Cancer Institute's tumor registry and through chart reviews. Tissue cores with >70% tumor content was obtained from formalin-fixed and paraffin-embedded (FFPE) tissue.

Small RNA was extracted from deparaffinized and proteinase K–treated FFPE core tissue (20–40 mg) using the High Pure miRNA Isolation Kit (Roche). Concentration and quality of RNA was assessed by absorbance spectrometry and electrophoresis on NanoDrop (Thermo Fisher Scientific, Inc.) and Bioanalyzer 2100 (Agilent Technologies) instruments. For a third of the cases, RNA preparations were of poor quality. Consequently, RNA was extracted again from FFPE tissue.

Levels of mature miRNAs miR-21, miR-30a (3p), miR-103, miR-106a, miR-200b, and miR-623 and U6 small nuclear 2 (U6B) RNA were measured using TaqMan miRNA assays (Applied Biosystems). Twenty nanograms of RNA was reverse-transcribed in 15 μL at 42°C for 30 min using a miRNA- or U6B-specific oligonucleotide. MiRNA-specific primers and 1.33 μL of reverse transcribed reactions were used in triplicate 45-cycle PCR of 20 μL volume in a 7900HT thermocycler (Applied Biosystems); the denaturation step at 95°C was done for 15 s, and the annealing and extension step was done at 60°C for 1 min. Fluorescence from the binding of miRNA-specific, carboxy-fluorescein dye–conjugated probes to amplification products was measured during PCR, and SDS software (2.3; Applied Biosystems) was used to identify cycle threshold (C_t) values. Statistical analysis and graphical plotting was done using Prism 5.0b software (GraphPad).

This work was performed by Exiqon in two experiments using 20 (10 each from the two groups) or 59 of the total 79 cases that yielded RNA of good quality. Reference RNA was generated by pooling equal amounts of RNA from all samples used in an experiment. miRCURY miRNA Power Labeling Kit (Exiqon) was used to 3′- or 5′-end label 0.5-μg samples or reference RNA with Hy3 or Hy5 (Exiqon) dye, respectively, before they were cohybridized to miRCURY locked nucleic acid (9) microarrays (10.0; Exiqon) on an HS4800 hybridization station (Tecan) for 18 h at 56°C with agitation. Slides were scanned on a G2565BA microarray scanner (Agilent Technologies) for analysis using ImaGene 8.0 software (BioDiscovery). Arrays had probes for multiple small RNAs of various species printed in quadruplicate on randomly distributed spots of 90 μm diameter. Human-specific probes covered 722 known and 30 proprietary (miRPlus, Exiqon) miRNAs, and 12 non-miRNA small RNAs. Raw data is available online in the ArrayExpress database (accession E-TABM-727).³

Hy3 and Hy5 signal intensities from 78 hybridizations of good quality were processed together as follows using the limma package (2.17; ref. 10) for Bioconductor (2.3) in R (2.8.1). After background correction using the convolution model–based normexp method with an offset of 10 (11), values were normalized first by the global Loess regression (12), with a smoothing value of 1/3, and then by limma's “rquantile” method to achieve identical distributions of Hy5 intensities across all arrays. Summarized intensity and intensity ratio (expression relative to reference) values for a probe were the mean of values from its four spots when the maximum was <150% of the minimum; the median was used otherwise. Only data for human RNAs (764) was considered from this point onwards. Two hundred and ninety RNAs, with Hy3 intensities >50% of the within-array interquartile range in >50% of samples, were considered expressed. One sample was excluded as a reassessment of clinical records indicated that the patient did not have stage I NSCLC. Expression values relative to reference were used to identify differentially expressed genes; log 2–transformed values were used for classification and hierarchical clustering using uncentered Pearson correlation and average linkage measures in R or TM4 (13) MultiExperiment Viewer (4.3).

The CMA R package (1.3; ref. 14) was used for classification using support vector machines (SVM; linear kernel). For leave-one-out cross-validation, training sets consisting of all but the single test samples were used to identify the top six differentially expressed miRNAs by the Wilcoxon rank sum test. As is recommended (15), a nested cross-validation loop was used to choose the best value for the C SVM kernel parameter from 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, and 50. Alternately, top-scoring gene pairs (TSP) were used for classification using the tspair R package (1.2; ref. 16). For Monte Carlo cross-validation (17) based on similar methodologies, 1,000 randomly selected 60-sample training sets were used to classify the 17 remaining samples. Model construction in multivariate analysis using Cox proportional hazards regression for clinical variables gender, age, pathologic stage of cancer (IA or IB), tumor histology, smoking history (ever or never), and presence or absence of TSP- or SVM-based classifier signature, with recurrence as the primary outcome, was performed by stepwise elimination of nonsignificant variables (P > 0.1) using SPSS software (15.0; SPSS).

A Microsoft Excel add-in (0.953)⁴

Among the 77 surgically treated pathologic stage I NSCLC cases for which miRNA expression profiles were examined, clinical recurrence of the cancer was documented in 37 with a median time of 13 months (range, 4–57) after surgery, with 36 having it within 32 months. The time of detection of recurrence was not known for four, and the time to death was used instead. The 40 recurrence-free cases were followed for at least 32 months, with 17 followed for at least 5 years. The characteristics of the two case groups are summarized in Table 1, and case-specific information is provided in Supplementary Table S1. For univariate analysis, categorical (gender, race, pathologic stage, tumor histology, type of surgery, history of smoking) and continuous (age) clinical variables were compared between the two case groups using χ² and two-tailed t tests, respectively, and no statistically significant association with recurrence was found.

RNA of <100 nucleotides in length was extracted from resected and FFPE-preserved cancer tissue. Median preservation ages for groups for recurrence and no recurrence were 2.5 and 1.5 years, respectively. For RNA samples of good quality, prepared from 79 tissues, the mean of ratio of absorbances at 260 and 280 nm was 1.92 (range, 1.68–2.05; SD, 0.09), of ratio of absorbances at 260 and 230 nm, 1.65 (1.23–2.13; 0.21), and of RNA integrity numbers (18), 2.4 (1.7–2.6; 0.2).

A two-color, locked nucleic acid (9) microarray platform with probes for 752 human miRNAs and 12 human non-miRNA small RNAs was used to quantitate gene expression relative to a reference generated by pooling all RNA samples. Signal intensities were corrected by background subtraction, and global Loess and quantile normalizations. Two hundred and ninety (38%) RNAs (279 miRNAs), with intensities >50% of the within-array interquartile range (corresponding to 1.1–2.5 times the median values) in >50% samples, were deemed as being expressed (Supplementary Table S2).

Levels of mature miRNAs miR-21, miR-30a, miR-103, miR-106a, miR-200b, and miR-623, chosen randomly, were assayed in eight arbitrary RNA samples by reverse transcription-PCR (RT-PCR) to validate expression data obtained by microarray hybridization. Significant correlation was observed between log 2–transformed array Hy3 signal intensities and RT-PCR C_t values for all six miRNAs (Fig. 1). The relative expression for each of the miRNAs with respect to that of the small nuclear RNA U6B, quantified as the difference in C_t or log 2–transformed signal intensity values, also correlated significantly in both Pearson and Spearman tests (Supplementary Table S3). For miR-30a, miR-103 and miR-106a, compared with RT-PCR, array-based quantitation reduced intersample differences by 1.5-fold, 4.6-fold, and 2.3-fold, respectively; differences were exaggerated 1.3-fold for miR-623 and unaffected for miR-21 and miR-200b (Fig. 1).

Using the nonparametric Wilcoxon rank sum (Mann-Whitney) test at a critical value of P = 0.01, 130 (47%) of the 279 expressed miRNAs were found differentially expressed (Supplementary Table S2). In the case group for recurrence, 59 were underexpressed; Table 2 lists 10 each of the most underexpressed and overexpressed miRNAs. Eight of 10 (80%) and 35 of 40 (88%) miRNAs were common to the most differentially expressed miRNAs identified by two other group comparison tests: t test assuming equal group variances (101 miRNAs with P < 0.01), and Welch-Satterthwaite t test for unequal group variances with t statistics based on 10,000 grouping permutations (105 miRNAs with P < 0.01; Supplementary Table S2).

Unsupervised hierarchical clustering of the 77 cases based on the levels of all 279 expressed miRNAs shows that cases that had recurrence cluster separately from those without recurrence (Fig. 2A). Figure 2B shows a heat-map illustrating the different miRNA expression profiles of the two groups; 56 miRNAs were differentially expressed between the two groups after multi-test correction for a false discovery rate of 0.5%.

Based on the differences between the two groups for the expression of many miRNAs, signatures comprised of multiple miRNAs that could predict recurrence were sought. For leave-one-out cross-validation, a linear SVM classifier based on the six most differentially expressed miRNAs was built using a training set of 76 samples and tested on the one remaining test sample. The linear SVM classifier method has the advantage of having only one adjustable kernel parameter, cost (C). To avoid biasing the estimation of the error rate by adjusting this parameter on the test set (15), an internal cross-validation was first performed to select the optimal value for C. Four miRNAs, miR-30c-1*, miR-200b*, miR-510, and miR-630, were present in all 77 six-miRNA classifiers, whereas miR-657 and miR-146b-3p were present in 59 and 28, respectively. Figure 3A shows the predictions against a Kaplan-Meier time-to-recurrence plot [hazard ratio (HR), 3.6; 95% confidence interval (95% CI), 1.8–7.1; P = 0.0002]. HR was 2.9 (95% CI, 1.4–6.4; P = 0.006) for overall survival (Fig. 3C). In a Monte Carlo cross-validation analysis using 1,000 random splits into 60-sample training and 17-sample test sets, such six-miRNA classifiers had a mean prediction accuracy of 69% (range, 35–94%; 95% CI, 68–70), and mean (95% CI) sensitivity, specificity, and positive and negative predictive values, respectively, of 67% (66–68%), 72% (71–73%), 70% (69–71%), and 69% (68–70%). Only 22 (7.9%) of the 279 expressed miRNAs were present in 98% of the 1,000 six-miRNA classifiers, the six most common being miR-200b* (100%), miR-30c-1* (77%), miR-510 (68%), miR-630 (67%), miR-657 (29%), and miR-146b-3p (27%). The intercase and intergroup variabilities in the expression of miRNAs is shown in Fig. 4A.

In leave-one-out cross-validation using the TSP classifier method (19), which considers even nondifferentially expressed genes, prediction accuracy was 83%, with miR-124*, which is not differentially expressed between the two case groups, and miR-585 being the paired genes in all 77 classifiers. Figure 3B shows the predictions against a Kaplan-Meier time-to-recurrence plot (HR, 9.0; 95% CI, 4.4–18.2; P < 0.0001). HR was 5.2 (95% CI, 2.4–11.4; P < 0.0001) for overall survival (Fig. 3D). In Monte Carlo cross-validation, the classifiers had a mean prediction accuracy of 72% (range, 24–100%; 95% CI, 71–72), and mean (and 95% CI) sensitivity, specificity, and positive and negative prediction values, respectively, of 68% (67–69%), 75% (74–76%), 73% (72–74%), and 70% (69–71%). Ten (4%) of the 279 expressed miRNAs were present in 98% of the classifiers, the three most common being miR-124* (60%), miR-585 (51%), and miR-708 (17%). The intercase and intergroup variabilities in expression of miR-124* and miR-585 is shown as an inset in Fig. 4A.

As expected, as the training set size in Monte Carlo cross-validation was increased from 20 to 75, prediction accuracy improved from 60% to 70% and 61% to 81% for the SVM and TSP methods, respectively (Supplementary Fig. S1).

Multivariate analysis using Cox proportional hazards regression led to a model consisting of only the male gender (P < 0.05), and the presence of SVM (P < 0.05) and TSP (P < 0.01) classifier signatures. Other clinical variables did not add to the prediction of recurrence in a statistically significant fashion.

Of 12 non-miRNA small RNAs quantitated by microarray hybridizations, 11 were considered expressed. Except for SNORD6, SNORD13 (U13), and SNORD14B (U14B) small nucleolar RNAs, all others, including U6 RNAs were differentially expressed between the two case groups (Supplementary Table S2). For use for normalization when quantitating miRNAs in future NSCLC studies, candidate reference miRNAs were therefore sought. Analysis of expression values relative to the reference of all nondifferentially expressed miRNAs across the 77 samples showed that 25 miRNAs had a mean between 0.9 and 1.1, coefficient of variation <25%, and a ratio of intergroup means between 0.9 and 1.1. The 25 miRNAs were ranked for gene expression stability using NormFinder (20), which uses a model-based variance estimation approach and attempts to minimize estimated intragroup and intergroup variations, and geNorm (21), which estimates stability using a pairwise comparison method (Supplementary Table S4). Three miRNAs, miR-488, miR-503, and miR-647, were among the five best reference miRNAs identified using either algorithm. The intercase and intergroup variabilities in expression of the three miRNAs, and of U6B, SNORD2, and SNORD6 are shown in Fig. 4B.

Pathway enrichment analysis was performed using the NIH DAVID Tools (22) or miRGator (23) online software on the gene targets predicted in the miRBase Targets database⁶

Prognostication of early stage lung cancer to identify cases likely to recur can help guide the administration of potentially harmful adjuvant therapies to those most likely to benefit from it. With this goal, the expression of 752 miRNAs in 77 stage I NSCLC cases were quantitated using microarrays. This examination of more than four-fifths of known miRNAs (5) using more specific and sensitive locked nucleic acid probes (9) is larger than that described in any published work on lung cancer. Formalin fixation does not adversely affect miRNA quantitation (24), and FFPE tissue, deliberately chosen for easy availability for future applications, was successfully used for extracting RNA for hybridizations of good quality. Based on a stringent criterion for detection of miRNA expression for more robust analysis, 279 miRNAs were considered expressed (Supplementary Table S2). They included more than half of the 43 miRNAs differentially expressed in lung cancer compared with normal lung as per a stage-oblivious study (6). About half of the miRNAs were differentially expressed between cases with recurrence and those without (Supplementary Table S2), although artifacts (25) may have arisen because of a slight intergroup difference for the ages of preserved tissues (Supplementary Table S2).

Consistent with previously published results (6–8), this study shows that miRNA profiling can be used to predict the recurrence of stage I NSCLC after resection. Yanaihara and colleagues (6), identified the association with prognosis of lung adenocarcinoma for 8 of 147 examined miRNAs; 7 of them were also differentially expressed between this study's case groups. MiRNAs let-7a-2 and miR-155 were both associated with poor outcome with a risk ratio of 3. Yu and colleagues (7), evaluated the expression of 157 miRNAs and found a five-miRNA signature to be an independent predictor of overall survival in NSCLC, although for stage I NSCLC (38% of training set), the signature was weakly associated (P = 0.03) with relapse-free survival. Two of the five miRNAs (let-7a and miR-221) were differentially expressed between this study's case groups. By quantitating 328 miRNAs, Raponi and colleagues (8), associated miR-146b and miR-155 with overall survival in squamous cell NSCLC (HR, 5.9 and 2.3, respectively); the former is also identified in this study as a marker for recurrence. Supplementary Figs. S2 and S3 show that the significant miRNAs identified in these studies could predict recurrence among this study's cases as well. This also suggests that miRNA profiling studies on NSCLC, unlike studies that have examined mRNA expression (26, 27), could have more consistent results. Because the previous studies did not quantitate the expression of most of the classifier-constituting miRNAs identified here, it is not possible to validate the SVM- or TSP-based classifiers against the data from those studies.

Although these three studies showed that miRNA profiles can be prognostic, they assessed the expression of <40% of currently known miRNAs using less powerful platforms, required fresh-frozen tissue, included cases with tumors of all stages, or looked for association with overall survival. Patients with lung cancer beyond stage I receive therapy apart from surgery, and this confounds the evaluation of naïve biological aggressiveness as adjuvant therapy affects clinical end points. Patients may also have died from causes other than the cancer. The studies did not detail the extent of surgical resection. This is important because recurrence after limited (wedge) resection is significantly more likely than after lobectomy (28).

Leave-one-out and Monte Carlo cross-validation analyses show that miRNA profiles obtained in this study could predict the recurrence of stage I NSCLC. Misclassification was decreased in the former, suggesting an improvement in classifier performance by increasing training set size, as is shown in Supplementary Fig. S3. Compared with SVM, the simpler and normalization-dispensable TSP classifier method was more accurate at predicting recurrence (83% versus 70%). Prediction accuracy of a miRNA-based classifier has not been reported thus far. The HRs associated with the classifiers, 3.6 (SVM) and 9.0 (TSP), are better than the only one (2.4) reported for a miRNA-based classifier for recurrence in stage I NSCLC (ref. 7; Fig. 3A and B). This may be because of the higher number of miRNAs examined in this study. The prediction accuracies and HRs for the miRNA-based classifiers identified here are similar to those described using mRNA-based classifiers (4, 26, 29–32). Multivariate analysis supports the conclusion that the miRNA-based classifiers are independent of clinical variables.

As has been observed before (33, 34), expression of small “housekeeping” RNAs varied widely among the cases; seven, including U6B which has been used for normalization when quantitating miRNAs in previous studies on NSCLC (6, 7), are differentially expressed between the two case groups (Fig. 4B; Supplementary Table S2). However, this could have resulted from the differential expression of genes responsible for miRNA processing as has been shown to occur as lung cancer evolves (35, 36). Among the miRNAs that were not differentially expressed, levels of miR-488, miR-503, and miR-647 varied the least among cases as well as the two case groups (Fig. 4B). These miRNAs may be the best possible reference miRNAs for use in future studies on FFPE tissue from early stage NSCLC.

S. Knudsen focuses on developing and commercializing microarray-based technology for cancer treatment. The other authors disclosed no potential conflicts of interest.

We thank Jason North and Wiam Bshara, Department of Pathology, Roswell Park Cancer Institute, for assisting with the selection of tissue cores.

Grant Support: Roswell Park Cancer Institute, Buffalo, NY and Medical Prognosis Institute, Hørsholm, Denmark.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.