MicroRNA-31 Predicts the Presence of Lymph Node Metastases and Survival in Patients with Lung Adenocarcinoma

This article is featured in Highlights of This Issue, p. 5259

Lung cancer is among the deadliest type of cancer in men and women. In 2012, there were an estimated 226,160 new cases of lung cancer and 160,340 deaths in the United States alone (1). Non–small cell lung cancer (NSCLC) accounts for approximately 80% of all lung cancer and is further classified into subtypes, including adenocarcinoma (ADC), squamous cell carcinoma (SCC), and large cell lung carcinoma. Lung ADC is the most common lung cancer subtype comprising approximately 40% of all lung cancer. Survival has not improved significantly over several decades despite enormous research efforts.

MicroRNAs (miRNA) are 20 to 22 nucleotide long, single-stranded noncoding RNAs, which regulate gene expression by transcriptional suppression (2). It is now well known that miRNAs are involved in many different biologic processes such as cell proliferation, cell death, and fat metabolism (3). Several studies have used miRNA expression profiling in primary tumors as a potential tool to identify novel signatures associated with diagnosis, staging, progression, prognosis, as well as response to treatment (4, 5). The potential of miRNAs as prognostic or diagnostic biomarkers has been intensely evaluated and has shown great promise in certain cancer types, such as chronic lymphocytic leukemia (6), breast cancer (6), and lung cancer (7, 8). In lung ADC, one study reported that highly expressed pre-miR-155 and low expression of let-7 were correlated with poor prognosis (9). Another study identified miR-34a as a potential marker of relapse in surgically resected NSCLC (10), emphasizing the significance of miRNAs in NSCLC biology.

Recent advances in sequencing-based techniques allow the capture of a more complete set of transcription factor-binding information, genomic interactions, high-resolution DNA methylation, as well as mRNA transcripts in cells. The technique has further been used in studies to classify human cancers (11) and distinguish between induced pluripotent (iPS) and cancer stem cells (CSC; ref. 12). However, most of these studies have been conducted on in vitro cultured cell lines. Thus, we sought to apply this sequencing-based approach to primary lung ADC tumor samples in order to identify novel associations between miRNAs and clinical variables, including survival, to better understand the role of miRNAs in lung ADC pathobiology.

Patients with lung ADC with mediastinal lymph node biopsy establishing early-stage disease have been shown to be associated with lower local recurrence rates and improved overall survival (13). High-dose radiotherapy using stereotactic ablative body radiation (SBAR) or stereotactic body radiation therapy (SBRT) in patients that are unfit for surgical treatment has been shown to be a successful approach (14). Currently, there are several ongoing clinical trials and several that are expected to soon report survival outcomes to evaluate the efficacy of SBRT for operable patients (i.e., ACOSOG Z4099/RTOG 1021, a phase III trial comparing sublobar resection with SBRT for high-risk patients with stage I NSCLC). SBRT is typically not preceded by surgical lymph node staging and patients routinely undergo positron emission tomography (PET)/computed tomography (CT) for staging. However, the sensitivity to detect metastatic tumor in lymph nodes smaller than one centimeter on PET/CT is low (32.4–40%; refs. 15, 16). Molecular markers would, thus, be of great use as they would allow for a more accurate risk assessment of nodal metastasis in addition to PET scans, which ultimately could improve patient selection for SBRT beyond classical staging. In addition, a better understanding of the degree of invasiveness of an individual tumor could guide the field of radiation oncology in the future to individualize radiation field expansions beyond the gross tumor to include invasive microscopic disease and, at the same time, maximize normal tissue sparing. For these reasons, we sought to identify molecular signatures of primary lung ADCs associated with the presence of lymph node metastases.

In this study, we conducted genome-wide miRNA-seq in samples from patients of lung ADC and used bioinformatic analyses to identify novel and annotated miRNAs that are differentially expressed between patients with and patients without lymph node metastasis, as well as examine their potential prognostic value. We selected one miRNA that was highly upregulated in patients of lung ADC with lymph node metastases for further characterization, including in silico mining of publicly available data generated by The Cancer Genome Atlas (TCGA) consortium and in vitro assays to characterize the function of this miRNA in lung ADC cell lines. Our study applied miRNA-seq, for the first time, to study microRNomes in samples from patients of lung ADC for the discovery of novel miRNA signatures that are potential biomarkers for outcomes in lung ADC.

Sixty-four frozen tissue specimens from patients of lung ADC collected between 2003 and 2012 were obtained through the OSUCCC Tissue Procurement Shared Resource, supported, in part, by the NCI Cancer Center Support Grant (CCSG-P30) and the Cooperative Human Tissue Network (CHTN/NCI), Midwestern Division, based on an institutional review board (IRB)-approved research protocol for this retrospective study. Other investigators may have received specimens from the same subjects. The remaining ADC tissues were flash frozen within two hours of surgical resection and stored at −80°C until analysis. Samples were received with a coded final pathology report and a quality assessment report verifying collection of tumor and/or normal adjacent tissue.

MiRNA samples from fresh-frozen tissue were extracted using PureLink microRNA Isolation Kit (Life Technologies Corporation). Briefly, 5 mg fresh-frozen tissue samples were mixed with 300 μL Binding Buffer (L3) and homogenized using a tissue homogenizer. Homogenized samples were centrifuged, and the supernatants were mixed with 300 μL 70% ethanol. Solutions containing RNA were purified by two rounds of spin column cleanup, and RNA was eluted with 50 μL sterile RNase-free water.

Total RNA samples were processed by the flashPAGE fractionator (Ambion) and flashPAGE Clean-Up Kit (Ambion). Enriched small RNA was then processed according to the SOLiD Small RNA Expression Kit protocol (Applied Biosystems), and purified small RNAs were ligated with 5′ and 3′ adapter Mix using RNA ligase. Ligated products (40–60 bases in length) were reverse transcribed and purified on Novex 10% TBE–Urea gel. Subsequently, 15 to 18 cycles of PCR were carried out by amplifying the purified cDNA with barcoded PCR primer sets provided with the kit, which differed by a unique six-nucleotide sequence. Amplified products were loaded on Novex 6% TBE gel (Invitrogen) and the gel bands containing 110 to 130bp fragments were excised. Amplified products were purified from the excised gel band and then loaded on the Applied Biosystems SOLiD next generation high-throughput sequencing system for data acquisition. The quality of the samples and libraries were verified on the Agilent Bioanalyzer (17).

To validate the identified targets from SOLiD system, quantitative real-time (qRT) PCR was used on fresh-frozen tissue samples. The miR-31 (miR-31-5p) and the small nuclear RNA U6 were analyzed in triplicates by TaqMan MicroRNA Assay kits (Applied Biosystems). A 30 ng sample of RNA was processed by the TaqMan MicroRNA Reverse Transcription kit (Applied Biosystems). In brief, 5 μL RNA was mixed with 1 mmol/L of each deoxyribonucleotide triphosphate, 50 units of Multiscribe Reverse Transcriptase, 5 × reaction buffers, four units RNase inhibitor, and 5× gene-specific RT primers mix in a final reaction volume of 15 μL. Reactions were then incubated at 16°C for 30 minutes, 42°C for 30 minutes, and 85°C for 15 minutes with a final hold at 4°C.

Afterward, 15 μL of the cDNA solution was diluted by nuclease-free water to a final volume of 100 μL. qPCR was run on a 7900 HT PCR machine with denaturation step at 95°C for 10 minutes, followed by 40 cycles of a denaturation step at 95°C for 15 seconds, and an annealing/elongation step at 60°C for 60 seconds. Fold-change for miRNA in tissue samples was then calculated using the equation 2^−ΔCt test/2^−ΔCt control. Paired Student t test was conducted to identify if miR-31 is differentially expressed between cancer and matched normal adjacent tissues (P ≤ 0.05). Welch t test was conducted to test if miR-31 expression is different between patients with and without lymph node metastasis to account for unequal variances. For sensitivity and specificity analysis, receiver operator characteristic (ROC) as well as area under the curve (AUC) was determined using CART (Classification and Regression Trees).

From among a total of 64 tissue samples, 10 samples (6 with N₀ stage and 4 with N₁₊ stage) were sequenced using the SOLiD platform. First, all reads were trimmed by removing the 3′ adaptor sequence using the cutadapt program (18), and the length distribution of reads for each sample has been plotted (Supplementary Fig. S1). Subsequently, we used a “Sequential Trim Alignment (SeqTrimAlign)” strategy to map all reads to the hg19 reference genome in color space (19), in which the last one color base of those unmapped reads in the previous alignment round would be trimmed, and resubmitted to the next round of alignment using Bowtie (20) until a minimum read length was reached–in our case, 19 color bases. This strategy is believed to fit specifically for solid data and could promote a high mapping ratio without sacrificing much performance. Detection of novel miRNAs was done by miRDeep2 using the above alignment, which uses the position and frequency of reads uniquely aligned to the genome (“signature”) with respect to a putative RNA hairpin and scores the miRNA candidate employing a probabilistic model based on miRNA biogenesis (17). The more positive the score, the more reliable the prediction would be. We selected the predicted miRNAs as novel candidates if they had a score of 2.0 or more and were present in at least three samples. Applying these selection criteria, 65 novel miRNAs were identified in our data set, and these were then included in the subsequent expression analysis.

Expression analysis of miRNA-Seq data was conducted using the R/Bioconductor package EdgeR (21), which is designed for use with digital gene expression data. First, we counted the number of reads uniquely mapped to miRNA regions according to the reference database miRBase (22) and novel miRNAs identified as described earlier; then, the normalization and differential expression between the N0 group and N1+ group was assessed in EdgeR by calculating an exact test P value analogous to the Fisher exact test.

Target prediction was done for miR-31, which was identified to be expressed higher in ADC N₁₊ patients than in N₀ patients. Currently, there is no single bioinformatics tool with accurate miRNA target prediction. However, integration of various computational methods is a common approach to improve prediction accuracy and to create an optimal framework for deciphering biologic functions of miRNAs. We used the miRWalk database (23), which reports predicted miRNA–mRNA interactions on the 3′ UTRs of known genes calculated by several established target prediction programs including TargetScan (24), miRanda (25), and RNAhybrid (26). Statistically significant miRNA–mRNA relationships were extracted from the results using two criteria: P value 0.05 or less and identification by at least five of the six selected target prediction programs.

The TCGA miRNA-seq data set with clinical information was downloaded on 18 July, 2012. The set included 341 level-3 lung ADC miRNA data sets. Among these, 233 out of the 341 patients from the TCGA lung ADC cohort had outcome as well as staging information available and were nonmetastatic. EdgeR normalized data of these 233 cases were used for correlative analysis. For analysis of the TCGA data set, Welch t test was conducted to test if miR-31 expression is different between patients with and without lymph node metastasis to account for unequal variances.

To assess if miR-31 controls the epithelial–mesenchymal transition (EMT) program, we conducted correlative analysis of miR-31 expression and genes known to be up- or downregulated during EMT (27). A total of 174 lung adenocarcinoma samples containing both gene expression and miR-31 expression were extracted from the 233-patient set with available miRNA expression data. As only level-3 data was analyzed, no additional normalization was done. Correlation coefficients were calculated using log2(EMT gene expression) and log2(miR-31 expression) data. The zero miR-31 expression values or EMT gene expression values were replaced by the minimum nonzero expression value of miR-31 or EMT gene.

Study of dependency between miR-31 expression levels and miR-31 promoter methylation status was analyzed using Illumina Infinium Human DNA Methylation 450 beadchip data. A total of 177 lung adenocarcinoma samples containing both methylation and miR-31 expression were extracted from the 233-patient set. Spearman correlation coefficient was calculated using β-values of methylation probe cg05146756 (chromosome 9: 21548871) and log2(microRNA31 expression). The zero miR-31 expression values or methylation β-values were replaced by the minimum nonzero expression value of miR-31 or β-values of the cg05146756 probe.

Two hundred and thirty-three, out of a total of 341 patients, from the TCGA lung ADC cohort had clinical information available and were used for correlative analysis. TCGA cohort staging was done according to either American Joint Committee on Cancer (AJCC) sixth or seventh edition, and tumor size was not available. The decision to classify M_x patients as M₀ for the purpose of this study was made on the clinical experience that pathologic synoptic reports commonly report M_x for clinically nonmetastatic tumor patients and that surgical primary tumor removal is typically only carried out in these patients. MiR-31 expression was log2 transformed to construct normal distribution. The associations between overall survival of all patients with the expression levels of miR-31, as well as the survival of T₂N₀ patients with the expression level of miR-31 were analyzed by using the Kaplan–Meier method, log-rank test, and Cox proportional hazard regression models. Statistical significance was accepted for P values less than 0.05.

Graphs were generated with GraphPad Prism 5 for Windows.

Analysis of biologic functions of target genes by miR-31 was done using the Ingenuity Pathways Analysis software. Target genes of miR-31 were uploaded to IPA to identify relevant cellular and metabolic functions and molecular networks. Pathways with P values less than 0.015 were determined to be significant and were included in further analysis.

Human lung ADC cell lines H23, H1573, and H2228 were purchased from the American Type Culture Collection (ATCC) and passaged for less than six months in RPMI1640 medium containing 10% FBS. The ATCC authenticates its cell lines using short tandem repeat (STR) profiling. The antibodies against β-Actin, Vimentin, Twist, and SNAI1 were purchased from Santa Cruz Biotechnology. Antibodies against phospo-AKT, AKT, phospho-ERK1/2, ERK1/2, and Actin were purchased from Cell Signaling Technology. Cells were preincubated with 10 μmol/L AZD6244 from AstraZeneca for 24 hours before conducting the assays as outlined further.

H23, H1573, and H2228 cells were stably infected with the human pre-miRNA expression construct Lenti-miR-31 vector and human miRZip-31 anti-miR-31 miRNA construct, respectively (System biosciences). The empty vectors were used as a control. The pre-miR-31, miRZip-31, anti-miR-31, and control constructs were packaged with pPACKH1 Lentivector Packaging Plasmid mix (System Biosciences) in a 293TN packaging cell line. The Transdux reagent (System Bioscience) was used for virus transduction, and infected cells were selected by fluorescence-activated cell sorting (FACS) analysis (FACSCalibular, BD Bioscience).

In vitro cell migration and invasion assay were conducted using Boyden chambers (BD bioscience) that use 8-μm micropore membranes without Matrigel (for migration assay) or with Matrigel (for invasion assay). Both assays were carried out according to the manufacturer's instructions. The cells were resuspended in 0.1% BSA in RPMI1640 medium and seeded in the upper chamber at a concentration of 5 × 10⁴/500 μL. The chambers were incubated in the wells containing RPMI-1640 medium with 10% FBS for 24 or 48 hours. Filters were fixed with 5% glutaraldehyde and stained with 1% crystal violet. The cells on the upper surface of the filters were removed by swabbing with a cotton swab, and the cells that had migrated to the reverse side were counted in 10 random fields under a microscope at ×100 or ×200 magnification.

Cultured cells were placed into 24-well plates at 2 × 10⁴ cells/well in triplicate. After three days, the numbers of cells were quantified using a Cell Counting Kit 8 (Dojindo) as described previously (28).

The cells were lysed with radioimmunoprecipitation assay (RIPA) buffer containing Protease/Phosphatase Inhibitor Cocktail (Cell Signaling Technology) and 10 or 50 μg protein lysates per lane were separated on 4 to 20% SDS–PAGE and, then, electrotransferred to nitrocellulose or polyvinylidene difluoride (PVDF) membranes. Membranes were blocked with blocking solution 1% or 5% BSA in Tris-Buffered Saline and Tween 20 (TBST) buffer and incubated with primary antibody in 3% BSA in TBST, followed by incubation with appropriate horseradish peroxidase (HRP)-conjugated secondary antibody. Specific proteins were detected using the enhanced chemiluminescence system (GE Healthcare).

We selected 43 cases with 64 tissue specimens available from patients of lung ADC (Supplementary Table S1) and conducted genome-wide miRNA-seq on tissues from 10 cases (Table 1, and analysis flow chart Supplementary Fig. S2). In four of those 10 cases, lymph node metastasis was present (pN₁ and pN₂, referred to as N₁₊) and, in six cases, mediastinal lymph node dissection was carried out but no lymph node metastasis was identified (pN₀). Using the EdgeR program (21) with a fold-change cutoff of 2 and a P value less than 0.05, we identified 32 annotated and three novel miRNAs that were differentially expressed between the N₁₊ and N₀ groups. Of 35 differentially expressed miRs, 29 were upregulated in N₁₊ cases and six were downregulated (Supplementary Fig. S3). Among these, miR-9, miR-31, the miR-34 family, and miR-224 have been previously reported to be deregulated in human cancers (29–32). Interestingly, we found that miR-31was upregulated in N₁₊ patients in our study cohort (Fig. 1A), consistent with miR-31 previously reported to be upregulated in metastatic colon cancer (33) but contrary to the reported high expression of miR-31 associated with reduced frequency of metastases in breast cancer (34).

We then conducted RT-qPCR of miR-31, which has not yet been reported to be associated with lymph node or distant metastasis in lung ADCs, on the specimens from the original 10 patients as well as on an additional 54 tissue specimens (clinical information of patients in Supplementary Table S1). There was good concordance between miR-31-5p qPCR (ΔCt) and miR seq (RPM normalized reads) data with a Pearson correlation coefficient of 0.554 (n = 10). Consistent with previous reports (35), miR-31 expression was higher in lung ADC cancer tissues compared with their matched normal adjacent tissues (paired t test, P = 2.86E-5, n = 21; Fig. 1B). Comparing miR expression by qPCR in a validation cohort of cancer tissues from patients with lymph node metastasis (N₁₊, n = 10) and patients without lymph node metastasis (N₀, n = 23), there was a higher expression of miR-31 in N₁₊ patients than N₀ patients (4.23-fold increase in N₁₊ patients, P = 0.009; Fig. 1C). Sensitivity/specificity analysis revealed a receiver operator characteristic area under the curve of AUC of 0.79 in the combined set (cohorts 1 + 2, n = 43; Supplementary Fig. S4).

We then tested our miR-31 lung ADC lymph node metastasis signature in an external validation set. The publicly available TCGA miR-seq data set with limited clinical information was downloaded on 18 July, 2012. Of the evaluable 249 patients with staging and survival information available, 140 patients had N₀M₀ disease, 93 patients had positive lymph nodes but no distant metastasis (N₁₊M₀), and 16 patients had distant metastasis (M₁; clinical information of 233 M₀ patients in Supplementary Table S2). Among the 93 cases of the N₁₊M₀ patient group, 49 patients had N₁ disease, 43 patients had N₂ disease, and one patient had N₃ disease. MiR-31 was 2.5-fold upregulated in patients of primary lung ADC with lymph node metastasis (N₁₊M₀) compared with those without (N₀M₀; t test, P = 0.031; Fig. 1D). The 16 patients with M₁ disease showed a statistically nonsignificant, but numerically higher, expression of miR-31 compared with M₀ patients. Similar results were obtained if the 43 patients with prior cancer diagnosis were excluded from the analysis (2.96-fold, t test, P = 0.02 for miR-31 for N₀M₀ vs. N₁₊M₀ comparison).

After establishing that miR-31 is associated with the presence of lymph node metastasis and given the knowledge that lymph node metastasis is a well-established prognostic factor, we examined the prognostic value of miR-31 in the publicly available TCGA data set. Factors associated with increased hazard of death on univariate analysis included miR-31 expression (log-rank test, P = 8.50E-04), N stage (log-rank test, P = 1.27E-4 for N₁ stage, P = 2.59E-6 for N₂₊ stage), and T₃ stage (log-rank test, P = 1.61E-03).

To formally evaluate if miR-31 expression is significantly associated with survival independent from lymph node involvement, we conducted Cox proportional hazards modeling on 233 lung ADC cases. In a multivariate Cox proportional hazard regression model including miR-31 expression, age, sex, lymph node stage, tumor stage, race, prior cancer diagnosis, and smoking history, the following factors were significant predictors of survival: miR-31 expression, nodal stage, and prior cancer diagnosis (Table 2).

After understanding that miR-31 expression is of prognostic value in lung ADC, independent of lymph node status, we conducted an exploratory correlative analysis of miR-31 in the T₂N₀ patient subset, given the controversy related to possible benefits of adjuvant chemotherapy in these patients and the need for novel biomarkers in this patient group (36, 37). We found that high miR-31 expression is associated with adverse outcome in patients with T₂N₀ lung ADC (log-rank, P = 0.0194) and, remarkably, there was no death observed for almost two years in the group with low miR-31 expression in the TCGA cohort (Supplementary Fig. S5).

Many studies have shown that downstream targets of miRNAs often regulate cancer initiation and progression (38, 39). Therefore, we wanted to further examine the targets of miR-31 and the pathways those targets are involved in to better understand the function of this miRNA. We applied the TargetScan software designed for in silico prediction of miRNA targets to identify the targets of miR-31. We found that 143 genes are predicted to be regulated by miR-31 when requiring prediction by at least five of the six selected target prediction algorithms using a P value cutoff of 0.05. We conducted IPA analysis on the set of miR-31 targets and found that the top five canonical pathways for miR-31 targeted genes are CDK5, PTEN, p70S6K, extracellular signal–regulated kinase (ERK)/mitogen-activated protein kinase (MAPK), and phosphoinositide 3 kinase (PI3K)/AKT signaling (Supplementary Fig. S6). The top five categories of associated biologic functions are cell death and survival, cell cycle, cell morphology, cellular growth and proliferation, as well as cellular function and maintenance (Supplementary Table S3).

As described earlier, miR-31 was significantly upregulated in patients of lung ADC with lymph node metastasis (N₁₊M₀) compared with those without (N₀M₀). To functionally characterize miR-31, we tested migration, invasion, and proliferation in three lung ADC cell lines (H23, H2228, and H1573; Supplementary Table S4). H23 and H2228 cells were transduced with a Lenti-miR vector containing miR-31 precursor to overexpress miR-31 in this cell line. H1573 cells were transduced with a vector containing a miRZip-31 anti-miR-31 miRNA construct to knockdown miR-31. The expression levels of miR-31, following treatment, were confirmed by qRT-PCR (Fig. 2A; Supplementary Figs. S7A and 8A). In response to overexpression of miR-31, cell invasion into Matrigel-coated transwell membranes was markedly increased (Fig. 2B and Supplementary Fig. S7B). In addition, overexpression of miR-31 significantly increased migratory and proliferative ability of H23 and H2228 cells (Fig. 2C and D, Supplementary Figs. S7C and D). Conversely, knockdown of miR-31 in H1573 cells significantly reduced cell invasion (Supplementary Fig. S8B), migration (Supplementary Fig. S8C), and proliferation (Supplementary Fig. S8D). We then tested if the in silico-identified miR-31-targeted pathways were indeed altered in our cell culture models. Pathway analysis identified ERK/MAPK and PI3K/AKT signaling as two of the top pathways related to genes that are predicted targets of miR-31. Identified pathways were analyzed for changes in activation by Western blot after miR-31 overexpression or knockdown. We observed no major signaling alterations in PI3K/AKT; however, we observed a reduction in ERK1/2 signaling in our H1573 miR-31 knockdown cell line and an increase in ERK1/2 activation in our H23 and H2228 miR-31 overexpression cell lines (Fig. 3A–C). To test if ERK1/2 signaling is necessary for miR-31 to induce an invasive phenotype, we assessed proliferation and migration in the H23 and H2228 cell lines overexpressing miR-31 after treatment with the MEK inhibitor AZD6244 (10 μmol/L). ERK1/2 signaling inhibition, indeed, suppressed the previously observed increases in proliferation or migration in miR-31–overexpressing cell lines compared with control cells (Fig. 4A–C; Supplementary Fig. S9A–C), suggesting that ERK1/2 pathway activation is a main mediator of the miR-31–induced invasive phenotype. ERK1/2 is known to be an important regulator of proliferation and migration (40). This observed association of miR-31 expression and ERK pathway activation is consistent with our observed phenotype in migration, invasion, and proliferation in cells that have miR-31 knockdown and overexpression. We then tested if the association between miR-31 expression and ERK phosphorylation can also be observed in tissue samples. There was no correlation identified (data not shown). To explore possible reasons for this observation, we assessed if miR-31 expression levels are similar in cell lines and in the tissues profiled (Supplementary Fig. S10). MiR-31 expression levels were, overall, lower in tissues than in cell lines, which might contribute to the observed lack of correlation of pERK and miR-31 in tissues.

Epithelial–mesenchymal transition (EMT) has been reported as a key process in tumor invasion, metastasis, and tumorigenicity. Several EMT-activating markers including Vimentin (VIM), ZEB1, SNAI1, and TWIST1 have been identified to be involved in this process (27, 41). We investigated the relationship between expression of these EMT-related markers and miR-31 expression in our lung ADC cell lines. As shown in Supplementary Fig. S11, the expression of Vimentin, TWIST1, and SNAI1 were increased after overexpression of miR-31 in H23 cells. A difference in expression in those proteins was not observed in H1573 miR-31 knockdown cells (data not shown). Given these observations, we sought to systematically assess if miR-31 could be a regulator of EMT. To this end, we analyzed associations of miR-31 expression and genes known to be up- or downregulated during EMT using the TCGA data set (Supplementary Table S5). Notably, the three EMT markers Vimentin, TWIST1, and SNAI1, which had been shown to be upregulated in the H23 miR-31–overexpressing cell line, were also positively correlated in the TCGA dataset. Overall, there were numerous EMT markers significantly correlated with miR-31 expression, indicating that miR-31 might contribute to the control of EMT. It remains to be determined if ERK signaling activation and EMT are co-incident events or if one leads to the other. However, we present data indicating that ERK signaling is necessary for the miR-31–induced invasive phenotype.

To understand if miR-31 is controlled by promoter methylation, the TCGA dataset was examined for correlation of miR-31 promoter methylation and miR-31 expression. There was a strong anticorrelation (correlation coefficient = −0.61, P = 4.4e-13), which suggests that DNA promoter methylation leads to repression of miR-31 expression.

In this study, we applied miRNA-seq, a next-generation sequencing (NGS) technique, to identify miRNAs associated with lymph node metastasis in tissue samples from patients of lung ADC. We identified and validated miR-31 to be associated with the presence of lymph node metastasis and patient survival. We further showed that miR-31 modulates the migratory, invasive, and proliferative behavior of lung ADC cell lines in culture.

Among the 32 annotated and three novel miRNAs identified with our comprehensive approach, there were several miRNAs, such as miR-9, miR-31, the miR-34 family, and miR-224, that have been previously reported to be deregulated in human cancers (29–32). MiR-31 was previously reported to be upregulated in colon cancer tissues compared with benign tissues, and was shown to increase proliferation and migration in intestinal cell lines (42), suggesting that miR-31 could be associated with tumor progression and have prometastatic functions. In squamous cell cancer of lung, high expression of miR-31 was associated with poor patient survival (43). In our study, we found miR-31 to be upregulated in lung ADC tissues from patients that have lymph node metastasis of adenocarcinoma, and further show, in vitro, that miR-31 leads to ERK1/2 pathway activation and promotes migration, invasion, and proliferation (44). Although we show that ERK signaling is necessary for the miR–31-mediated invasive phenotype and while there is evidence that ERK signaling can promote EMT (45, 46), it remains to be determined if miR–31-mediated ERK signaling activation and EMT are coincident events or if one leads to another. High expression of miR-31 in breast differs and has been reported to be associated with reduced frequency of metastasis (34). MiR-31 likely suppresses metastasis in breast cancer by targeting the metastasis-promoting genes ITGA5, RDX, and RhoA (47), and overexpression of miR-31 in established breast cancer metastases can lead to regression of those metastases (48). In prostate cancer, high levels of miR-31 are associated with better outcomes (49, 50). Previously, it has been shown that miR-31 has tumorigenic effects in lung cancer cell lines via repression of LATS2 and PPP2R2A, genes with known tumor-suppressive function (35). In smooth muscle cells, miR-31 expression appears to be controlled by ERK1/2 signaling and lead to LATS2–degradation-mediated increase in proliferation (51). However, a more recent report suggests DICER1, also a gene with known tumor-suppressive function, and not LATS2 or PPP2R2A is the functionally relevant miR-31 target in lung squamous cell carcinoma (43). These findings suggest that miR-31 plays an important role in tumor progression and metastasis formation in multiple cancers. However, the function of mir-31 appears to depend on the genetic background of the tumor.

We initially sought out to identify novel markers of lymph node metastasis. Given the increased use and success of SBRT for early-stage lung cancer and the low sensitivity of PET/CT for the detection of tumor positive small lymph nodes (<1 cm), markers to more accurately assess the risk of nodal involvement beyond imaging studies are very useful. We identified miR-31, a miRNA known to be deregulated in many cancers, to be upregulated in lung ADCs of patients with involved lymph nodes. MiRNAs appear to have potential for further development as markers to assess the risk of lymph node involvement using primary tumor biopsy tissues. Furthermore, we showed that miR-31 is an independent prognostic marker for lung ADC even when controlling for known prognostic factors such as tumor stage.

Given the controversy over patients with T₂N₀ lung ADC and possible benefits from adjuvant chemotherapy (36, 47, 52), we conducted an exploratory survival analysis of T₂N₀ patients in the TCGA data set to assess if the miR-31 signature for lymph node metastasis is of prognostic value in this patient subpopulation. We found that high miR-31 expression is associated with adverse outcome in patients with T₂N₀ lung ADC and, remarkably, there were no deaths observed for nearly two years in the low miR-31 expression group in the TCGA cohort. However, given that the TCGA data at this time includes only 70 death events (∼30% of included cases), this observation needs to be interpreted with caution. Our findings, if validated, may provide a rationale to evaluate miRNAs as biomarkers to determine which T₂N₀ patients might benefit from additional cancer therapy.

In summary, our findings have several implications: (i) miRNA-seq of tumor tissue is a feasible approach for biomarker discovery; (ii) miRNAs could predict not only the presence of lymph node metastasis but also clinical outcome; (iii) the mechanisms of how miR-31 controls invasion and metastasis appear to be complex and context specific as miR-31 plays a prometastatic role in lung ADC and lung squamous cell carcinoma, but plays an antimetastatic role in breast cancer and is a good prognostic marker in prostate cancer.

No potential conflicts of interest were disclosed.

The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Center For Advancing Translational Sciences or the NIH.

Conception and design: W. Meng, R. Cui, J. Perry, H. Nakanishi, C. Croce, A. Chakravarti, V. Jin, T. Lautenschlaeger

Development of methodology: W. Meng, S. Volinia, H. Nakanishi, V. Jin, T. Lautenschlaeger

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): W. Meng, R. Cui, J. Perry, V. Dedousi-Huebner, S-S. Suh, P. Ross, T. Lautenschlaeger

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): W. Meng, Z. Ye, R. Cui, J. Perry, A. Huebner, Y. Wang, S. Volinia, A. Chakravarti, V. Jin, T. Lautenschlaeger

Writing, review, and/or revision of the manuscript: W. Meng, R. Cui, J. Perry, V. Dedousi-Huebner, A. Huebner, B. Li, H. Nakanishi, L.W. Ayers, C. Croce, A. Chakravarti, V. Jin, T. Lautenschlaeger

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): W. Meng, V. Dedousi-Huebner, A. Huebner, B. Li, H. Nakanishi, T. Kim, L.W. Ayers, V. Jin, T. Lautenschlaeger

Study supervision: C. Croce, A. Chakravarti, V. Jin, T. Lautenschlaeger

This work was financially supported by the Department of Radiation Oncology grant (to T. Lautenschlaeger), the Comprehensive Cancer Center grant (to A. Chakravarti), and the Department of Biomedical Informatics grant (to V. Jin), The Ohio State University. The OSUCCC Tissue Procurement Shared Resource is supported, in part, by the Ohio State University NCI Cancer Center Support Grant (No. CCSG-P30). The project described was supported by Award Number Grant 8UL1TR000090-05 from the National Center For Advancing Translational Sciences.

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.