Identification of microRNAs as potential biomarkers for lung adenocarcinoma using integrating genomics analysis

Lung adenocarcinoma (LUAD) is the most common histological subtype of non-small cell lung cancer, but novel biomarkers for early diagnosis are lacking. Extensive effort has been exerted to identify miRNA biomarkers in LUAD. Unfortunately, high inter-lab variability and small sample sizes have produced inconsistent conclusions in this field. To resolve the above-mentioned limitations, we performed a comprehensive analysis based on LUAD miRNome profiling studies using the robust rank aggregation (RRA) method. Moreover, miRNA-gene interaction network, pathway enrichment analysis and Kaplan-Meier survival curves were used to investigate the clinical values and biological functions of the identified miRNAs. A total of six common differentially expressed miRNAs (DEMs) were identified in LUAD. An independent cohort further confirmed that four miRNAs (miR-21-5p, miR-210-3p, miR-182-5p and miR-183-5p) were up-regulated and two miRNAs (miR-126-3p and miR-218-5p) were down-regulated in LUAD tissues. Pathway enrichment analysis also suggested that the above-listed six DEMs may affect LUAD progression via the estrogen signaling pathway. Survival analysis based on the TCGA dataset revealed the potential prognostic values of six DEMs in patients with LUAD (P-value<0.01). In conclusion, we identified a panel of six miRNAs from LUAD using miRNome profiling studies. Our results provide evidence for the use of these six DEMs as novel diagnostic and prognostic biomarkers for LUAD patients.


INTRODUCTION
The incidence and mortality of lung cancer have increased steadily worldwide in the past several decades. Lung cancer is a fatal disease that is tightly associated with cigarette smoking. Apart from other histological subtypes, lung adenocarcinoma (LUAD) is the most common histological subtype among never-smokers [1]. Despite many achievements made in anti-cancer therapy over years, the survival of LUAD is far from satisfactory. Due to a lack of specific clinical symptoms, most patients with LUAD are diagnosed at a late stage, leaving little chance for effective treatment. Thus, exploring novel cancer-specific biomarkers for LUAD patients will help monitor tumor progression and guide clinical treatment.
MiRNAs are a group of small molecules that facilitate tumor progression via regulation of the function of protein-coding mRNAs [2]. MiRNAs are stable molecules that exist in tissue and fluid samples and exhibit potential as biomarkers for diagnostic, prognostic and therapeutic applications in multiple malignancies [3]. The application of high-throughput miRNA profiling methods during the past decade, such as RNA sequencing and microarrays, has enabled researchers to identify a group of miRNAs as biomarkers in cancer diagnosis [4]. The majority of previous studies were primarily based on small specimen sizes and various technological platforms, and no consistent conclusion has ever been made. To overcome these limitations, Vosa et al developed a robust rank aggregation (RRA) method for comparing several ranked gene lists and identifying the most commonly overlapping miRNAs [5]. By applying the RRA method, Vosa et al identified 15 aberrantly expressed miRNAs in lung cancer [6]. However, there has been no attempt to investigate LUAD-specific miRNAs using the RRA method.
In the present study, we first performed RRA analysis based on high-throughput miRNome profiling studies and identified candidate miRNAs in LUAD tissues. We next validated the aberrant miRNAs identified in the first step in an independent cohort. Finally, we predicted the biological functions of the above-identified DEMs and investigated the clinical values of these miRNAs in The Cancer Genome Atlas (TCGA) cohort.

Functional prediction of DEMs
The main function of miRNA is the silencing of gene expression via binding to specific target sites. Therefore, to elucidate the biological function of DEMs, we performed miRNA-target interaction and pathway enrichment of predicted targets. All predicted targets of each miRNA were obtained from TargetScan (http:// www.targetscan.org/vert_71/), and only those targets with a cumulative weighted context score less than -0.4 were selected for further analysis. Figure 5 shows the miRNA-target interaction of each candidate miRNA using Cytoscape software (version 3.4.0). The number of predicted targets of a single miRNA ranged from 13 to 169. Seventeen genes (ARHGAP6, BRMS1L, FAM175B, FAM217B, FRS2, KCNJ6, L3MBTL3, MBNL1, MITF, PDCD4, RABGAP1L, RASA2, RGS17, SAMD12, SERP1, SHC4 and TPD52) were targeted by more than one miRNA. We also performed pathway enrichment analysis of the predicted targets. Figure 6A shows that five KEGG pathways, including the "MAPK signaling pathway", "Neurotrophin signaling pathway", "Pathway in cancer", "Regulation of actin cytoskeleton" and "Axon guidance pathway", were associated with more than three miRNAs. Another algorithm (DIANA-miRPath v3.0) predicted the pathways affected by multiple miRNAs as "Pathway in cancer", "Regulation of actin cytoskeleton", "MAPK signaling pathway" and "Axon guidance pathway", which were also common pathways affected by the largest number of miRNAs (Table 2). Notably, we found that the "Estrogen signaling pathway" was the most significant pathway affected by six miRNAs ( Figure  6B). Activation of the estrogen signaling pathway was a reported promoter for LUAD [17], and our findings support critical roles for these six miRNAs in LUAD development.

The association between DEMs and clinical outcome in LUAD
To explore the prognostic values of DEMs in LUAD patients, we obtained the LUAD TCGA dataset from SurvMicro (http://bioinformatica.mty.itesm.mx/ SurvMicro), which is a web-based tool for assessing miRNA-based prognostic signatures [18]. We evaluated the association between miRNAs and patients' overall survival using two panels, including four up-regulated miRNAs (miR-21-5p, miR-210-3p, miR-182-5p and miR-183-5p) and two down-regulated miRNAs (miR-126-3p and miR-218-5p). Patients in each panel were divided into two groups, high-risk and low-risk groups, according to the prognostic index calculated by multivariate survival analysis using the SurvMicro web tool. Figure 7A shows that the panel of four up-regulated miRNAs was significantly associated with overall survival of LUAD, with a P-value of 0.0035. Patients in the high-risk group exhibited lower cumulative survival rates than the low-risk group (

DISCUSSION
In the present study, we combined nine miRome profiling studies and identified LUAD-specific miRNAs from a total of 595 LUAD and 168 non-cancerous tissue samples using the RRA method. A panel of four upregulated miRNAs (miR-21-5p, miR-210-3p, miR-182-5p and miR-183-5p) and two down-regulated miRNAs (miR-126-3p and miR-218-5p) were identified as commonly aberrantly expressed miRNAs in LUAD. Functional analysis revealed that these six miRNAs may be involved in LUAD development via modulation of the estrogen signaling pathway. Our clinical investigation further supports the prognostic value of these six miRNAs in LUAD patients.
In recent years, high-throughput profiling methods have been widely used to identify tumor-specific miRNAs as

37
Endocrine and other factor-regulated calcium reabsorption 3 www.impactjournals.com/oncotarget biomarkers for cancer diagnostic, prognostic and therapeutic applications [19,20]. However, due to the use of different methods and platforms, the conclusions varied among those profiling studies. The RRA method was specifically designed to identify the most commonly overlapping factors and make these studies comparable [5]. An increasing number of studies have attempted to examine tumor specific-miRNAs using the RRA method. Researchers used this method in multiple malignancies, including pancreatic, liver, renal, breast, colon, ovarian, bladder, and nasopharyngeal tract cancer [16,[21][22][23][24][25][26][27]. Võsa U et al [6] first identified a panel of 15 aberrantly expressed miRNAs in lung cancer, including seven up-regulated miRNAs (miR-21, miR-210, miR-182, miR-31, miR-200b, miR-205 and miR-183) and eight down- regulated miRNAs (miR-126-3p, miR-30a, miR-30d, miR-486-5p, miR-451a, miR-126-5p, miR-143 and miR-145). Their results provided solid evidence for the use of these miRNAs as promising diagnostic and prognostic biomarkers in lung cancer. However, the authors mixed different histological types of lung cancer together and did not analyze miRNA signatures separately. Therefore, their findings based on mixed samples may have limited value in the clinic. Therefore, this study only included miRNome profiling studies of LUAD. Our findings improve our understanding of the molecular mechanism and identify more specific tumor biomarkers for LUAD [28].
In non-small cell lung cancer, miR-218 was reported to be significantly down-regulated and affect cell migration and invasion [36,37]. Patients with low miR-218 expression exhibited a worse prognosis than patients with high miR-218 expression [38]. Using LUAD cell lines, Sher YP et al [39] found that the expression of miR-218 was significantly inhibited and involved in the process of brain metastasis. In our study, using an independent cohort, we also found that miR-218 was significantly down-regulated in LUAD tissues (fold change=0.29; P<0.0001). These reports support the need for further experimental and clinical exploration of miR-218 and other miRNAs in LUAD.
To further explore the underlying molecular association between DEMs and LUAD, in this study, we also performed miRNA-target interaction network and pathway enrichment analyses. We found that several KEGG pathways, such as "Regulation of actin cytoskeleton" and "MAPK signaling pathway", were affected by nearly all six miRNAs. A very recent study demonstrated that "Regulation of actin cytoskeleton" and "MAPK signaling pathway" contribute to disease progression and drug resistance in LUAD cells [40], which suggests the critical role of DEMs in LUAD. Notably, we also found that the "Estrogen signaling pathway" was the most significant pathway affected by the six miRNAs, with a P-value of 4.44E-05. For decades, epidemiological studies have demonstrated an increase in LUAD in women. A mouse model demonstrated that estrogen is an essential factor during LUAD carcinogenesis [17]. Ikeda K et al [41] observed that the estrogen concentration in LUAD was significantly higher than controls. Moreover, a clinical investigation revealed that postmenopausal women who received estrogen plus progestin treatment exhibit an increased number of deaths from lung cancer [42]. The estrogen signaling pathway was also involved in cellular cytoskeleton organization [43] and cross-talk with the MAPK signaling pathway [44]. Taken together, biological function analyses of these six miRNAs strongly suggest critical roles in LUAD, which may provide a new strategy for LUAD diagnosis and prognosis.
In conclusion, we identified a panel of six miRNAs (miR-21-5p, miR-210-3p, miR-182-5p and miR-183-5p, miR-126-3p and miR-218-5p) in LUAD and evaluated their expression and prognostic values in an independent cohort. Our biological informatics analysis provided evidence for the important roles of these miRNAs in LUAD via modulation of the estrogen signaling pathway. Our findings warrant further clinical investigations to validate these six miRNAs as biomarkers for the early diagnosis and prognosis prediction of LUAD.

Selection strategy for miRNome profiling studies
MiRNome profiling studies that used highthroughput methods, including second-generation sequencing, miRNA microarray and polymerase chain reaction (PCR) methods, which were designed for the parallel quantification of large numbers of miRNAs (96-or 384-well microplates) in LUAD tissue were included. MiRNome profiling studies based on only cell lines and non-LUAD subtypes were excluded. Rank lists of aberrantly expressed miRNAs from each study were mapped according to miRBase Version 21. Non-human miRNA, including viral and other species' miRNAs, were excluded. The miRNA lists were separately recorded as up-or down-regulated miRNAs.

Robust rank aggregation analysis
RRA analysis was performed using R software and an R package named "Robust Rank Aggreg" (freely available at comprehensive R Archive Network website, http://cran.r-project.org/). All miRNA lists were re-ranked according to the P-value assigned for each miRNA. Bonferroni corrections of P-values were calculated to avoid false positive results. Only corrected P-values less than 0.05 were considered statistically significant.

MiRNA-gene interaction network analysis
The TargetScan web tool (www.targetscan.org/) was used to predict the gene target of each DEM. Only gene targets with a cumulative weighted context score<-0.4 were selected for inclusion in the miRNA-gene interaction network. The miRNA-gene interaction network was visualized using Cytoscape software (version 3.4.0).

KEGG pathway enrichment analysis
KEGG pathway enrichment analysis for the predicted gene targets of each DEM was performed using the DAVID web tool (https://david.ncifcrf.gov/). The false discovery rate (FDR) of each KEGG pathway was transformed by -log10 and visualized as a heatmap. The combinatorial effects of multiple miRNAs in KEGG pathways were calculated using the miRPath algorithm in the DIANA web tool (http://www.microrna.gr/ miRPathv2).

Survival analysis
LUAD patients from The Cancer Genome Atlas (TCGA) cohort were used to investigate the prognosis values of DEMs in LUAD. The prognostic data and Kaplan-Meier curve were calculated on SurvMicro web tools (http://bioinformatica.mty.itesm.mx/SurvMicro).

Statistical analysis
All data in the present study were analyzed using GraphPad Prism 6.0 software. The expression of each miRNA is expressed as a mean ± standard error of measurement (SD). Comparisons of each miRNA between cancer and normal groups were performed using Student's t-test. A P-value less than 0.05 was considered statistically significant.

Author contributions
Zhuo Peng worked for study search, quality check, data extraction, and analysis as a principal investigator. Longfei Pan, Zequn Niu, Wei Li and Xiaoyan Dang contributed for study conception. Lin Wan and Rui Zhang worked as methodologist. Shuanying Yang critically revised the manuscript and coordinated the study. All authors reviewed and finally approved the manuscript.