MiRNA expression profiling in adenocarcinoma and squamous cell lung carcinoma reveals both common and specific deregulated microRNAs.

The current study investigated the expression signatures of miRNAs in lung adenocarcinoma (LUAD) and squamous cell lung carcinoma (LUSC). miRNA profiling was performed using microarray in 12 LUAD and 12 LUSC samples and adjacent normal tissues. In LUAD, 107 miRNAs were significantly deregulated, whereas 235 miRNAs were deregulated in LUSC. Twenty-six miRNAs were common between the 2 cancer subtypes and 8 were prioritized for validation, in addition to 6 subtype-specific miRNAs. The RT-qPCR validation samples included 50 LUAD, 50 LUSC, and adjacent normal tissues. Eight miRNAs were validated in LUAD: 3 upregulated - miR-7-5p, miR-375-5p, miR-6785-3p, and 5 downregulated - miR-101-3p, miR-139-5p, miR-140-3p, miR-144-3p, miR-195-5p. Ten miRNAs were validated in the LUSC group: 3 upregulated - miR-7-5p, miR-21-3p, miR-650, and 7 downregulated - miR-95-5p, miR-140-3p, miR-144-3p, miR-195-5p, miR-375, miR-744-3p, and miR-4689-3p. Reactome pathway analysis revealed that the target genes of the deregulated miRNAs in LUAD were significantly enriched in cell cycle, membrane trafficking, gene expression processes, and EGFR signaling, while in LUSC, they were enriched in the immune system, transcriptional regulation by TP53, and FGFR signaling. This study identified distinct miRNA profiles in LUSC and LUAD, which are common and specific miRNAs that could be further investigated as biomarkers for diagnosis and prognosis.

1. Introduction

Lung cancer (LC) is the leading cause of death across the world.[ nonsmall-cell lung cancer (NSCLC) accounts for approximately 85% of all lung cancers, with adenocarcinomas (LUAD) and squamous cell carcinomas (LUSC) being the main histological subtypes.[ They can be further subdivided based on their driver mutations, expression, and methylation profiles.[

The role of miRNAs in regulating gene expression signatures in 15 epithelial cancer types from the Cancer Genome Atlas was recently shown, providing evidence for a complex map of interactions underlying the relationship between miRNA regulation and the hallmarks of cancer.[ The crucial role of miRNA deregulation through various mechanisms is exerted on many targets, including oncogenes and tumor suppressor genes. This determines the complexity of the investigation of miRNA function, which may also vary between tumor types and stages. A future challenge would be to identify and validate the specific miRNA expression profiles and critical targets of the miRNAs involved in each cancer type and establish their contribution to malignant transformation.[

Many previous studies have examined miRNA signatures in patients with NSCLC, either by microarray or NGS analysis, to identify reliable biomarkers for diagnostic and prognostic purposes. Dysregulated expression of many miRNAs in NSCLC has been reported.[ However, the validation of tissue-based tumor miRNA biomarkers has proven difficult, leading to inconsistent results in independent studies that could not be replicated. The reasons for this might be various, related to differences in the sample composition and size, quality, inter- (stage, age, sex, ethnicity), intratumor heterogeneity of the included samples, and differences in the methodologies used in the preanalytic phase for collection and storage of materials. In addition, further variability and difficulty for between-study comparison come from the variable platforms and methods used for miRNA expression studies, endogenous controls, and statistical analysis of the results, containing different numbers of miRNAs. Moreover, based on their initial studies, groups usually choose only a small number of differentially expressed miRNAs to replicate in further invalidation samples.[

In parallel, others have taken different approaches, such as data mining and bioinformatics analysis of available public datasets, to identify and validate differentially expressed miRNAs.[ Due to the considerable difficulties in gathering high-quality tissue material for expression studies, there are few published papers, and many investigators have explored circulating miRNAs as biomarkers using serum or plasma samples.[ However, these were not the focus of our study.

Several recent reviews have summarized studies exploring miRNAs as NSCLC subtype classifiers and diagnostic markers.[ The miRNAs that have been repeatedly confirmed in various studies to be differentially expressed and to distinguish between LUAD and LUSC so far include miR-375 and miR-205, while the sets of miRNAs identified to distinguish LUAD and LUSC from normal tissue vary greatly between studies with little overlap.[

The current study aimed to determine the molecular profile of 2 major subtypes of NSCLC: adenocarcinoma (LUAD) and squamous cell lung cancer (LUSC). Comprehensive microarray miRNA expression analysis in a set of fresh-frozen tissues and their adjacent normal tissues was used to define the specific profiles of deregulated miRNAs. We validated the most significant findings in independent samples of paired fresh-frozen tumors and normal tissues using RT-qPCR. Correlation analysis, followed by prediction of target analysis, gene ontology, and pathway analysis, was performed to better understand the gene targets of the deregulated miRNAs, affected molecular pathways, and mechanisms and processes related to tumorigenesis in the 2 subtypes of lung cancer. The workflow of the research is illustrated in Figure 1.

Figure 1.

Workflow of the current study.

2. Materials and Methods

2.1. Tissue specimens

In this study, 124 tissue samples were obtained from patients who underwent thoracotomy and various pulmonary resections for lung cancer. Samples were collected at the University Hospital for Pulmonary Diseases (“St. Sofia”) between 2009 and 2017. Tumors and adjacent nontumor tissues were quickly frozen and stored at –80°C until RNA extraction. A pathologist reviewed the samples and diagnosed them as adenocarcinoma or squamous cell lung carcinoma (n = 124; 62 LUAD, 62 LUSC). Clinicopathological characteristics of the patients are shown in Table 1. The Ethics Committee of the Medical University of Sofia approved this study, and all participants provided written informed consent.

Table 1

Clinicopathological characteristics of the patients with NSCLC included in the study.

Characteristics	Lung adenocarcinoma (n = 62)	Squamous cell lung cancer (n = 62)
Females	22	12
Males	40	50
Average age female	62	63
Average age males	61	62
Tumor stage (T)
T1	9	13
T2	36	24
T3	14	23
T4	3	2
Lymph metastases (N)
N0	41	46
N1	8	7
N2	13	12
Distantmetastases (M)
M0	42	42
M1	20	20

2.2. RNA isolation

Total RNA was extracted from fresh frozen tissues using an miRNeasy mini kit (Qiagen, Germany) according to the manufacturer instructions. The quality and purity of the samples were quantified using a NanoDrop spectrophotometer (Thermo Fisher Scientific, USA). Accurate measurements were performed using a Qubit instrument (Thermo Fisher Scientific, USA). The RNA integrity number (RIN) was evaluated using an Agilent 2200 TapeStation System (Agilent Technology, USA). Only nondegraded samples with RIN > 7.0 were used for further microarray analysis.

2.3. Microarray analysis of MiRNAs

Total RNA (100 ng per sample) of 12 LUAD, 12 LUSC, and their adjusted nontumor tissues was evaluated using the Agilent Technologies miRNA profiling system. For miRNA profiling, we used the Agilent G4870C Human miRNA Microarray Kit Release 21.0, 8 × 60 K, following the manufacturer instructions for labeling and hybridization. Agilent Feature Extraction software (version 3.4) was used to analyze the acquired array images from Agilent Scanner G2505C. The 90th percentile normalization of the raw microarray data and the Mann–Whitney U test were performed using Agilent GeneSpring GX 7.3.1 Software. Differentially expressed (DE) miRNAs were selected through fold change (FC) and volcano plot filtering, and visualized by hierarchical clustering. The criterion for DE miRNAs was FC > 2 for upregulated and downregulated miRNAs, and the Benjamini-Hochberg FDR (BH-FDR) adjusted value ≤0.05 was considered statistically significant.

2.4. RT-qPCR for validation and association with clinicopathological characteristics of the patients

Based on the results of the microarray analysis, 14 DE miRNAs in lung tumors were selected. Furthermore, they were validated in a group of newly collected tissue samples from 50 LUAD and 50 LUSC patients and their adjacent nontumor tissues. Quantitative real-time reverse transcription PCR (RT-qPCR) was performed using SYBR technology (Qiagen, Germany), the miScript SYBR Green PCR Kit, and the miScript Primer Assay (Qiagen, Germany) on a 7900 HT Fast Real-Time PCR system (Applied Biosystems). All tests were performed in triplicate following the manufacturer protocol. RNU6B was used as the endogenous control. Relative quantification of data was performed using the 2-∆∆Ct method. We assessed the outliers with a Grubbs test and removed it. The expression data for all miRNAs were tested for normal distribution using the Shapiro–Wilk test. The Mann–Whitney U test was performed to determine the significance of expression, and a value ≤ 0.05 was considered statistically significant. Statistical analysis was performed using SPSSv20, and the Seaborn library in Python was used for visualization of the data.[

2.5. Correlation-based networks

To analyze the biological relationships between miRNAs, correlation-based networks (CBNs) were constructed. nonparametric Spearman rank correlation was used to calculate the correlation coefficients and p-values. Correlation matrices were generated using SciPy, a Python-based ecosystem.[ Only statistically significant correlations (P ≤ .05) were illustrated via correlation heatmaps with Matplotlib and Seaborn libraries in Python.[ A correlation coefficient of ≥ 0.5 and Benjamini-Hochberg FDR adjusted P-value ≤ 0.05 were taken for the design of the CBN threshold. The networks were visualized using NetworkX.[

2.6. Prediction of miRNA-target gene interactions

The gene targets of the miRNAs were predicted by the “multiMiR” package in the R software, which includes fourteen different databases divided into 3 groups. The first group contained 3 validated miRNAtarget databases (TarBase, miRTarBase, and miRecords), the second group contained 8 predicted miRNA–target databases (TargetScan, DIANA-microT, miRanda, MicroCosm, PITA, ElMMo, miRDB, and PicTar), and the third group included 3 disease- and drug-related miRNA databases (PhenomiR, miR2Disease, and Pharmaco-miR).[ miRNA target data were obtained from TarBase, miRTarBase, and miRecords.

2.7. Gene ontology (GO) annotation and pathway analysis

GO describes the biological functions of genes in 3 aspects: biological processes (BP), cellular components (CC), and molecular functions (MF). The Reactome Pathway Database was used to perform enrichment analysis of target genes to further understand the potential functions of DE miRNAs in LUAD and LUSC. All enrichment analyses were performed using the GSEApy Python package.[ Adjusted P-value cutoff of < 0.05 was used to indicate significant GO enrichment and Reactome pathway enrichment of miRNA targets were determined based on a P-value cutoff of 0.05.

3. Results

3.1. Differentially expressed mirnas in LUAD and LUSC

Using microarray analysis, we assessed the expression levels of 2549 human mature miRNAs and found 107 miRNAs to be significantly differentially expressed (FC > 2.0 and BH-FDR adjusted P-value ≤ 0.05) in LUAD compared with the adjacent normal tissues. Thirty-nine miRNAs were significantly downregulated and sixty-eight significantly upregulated (Fig. 2A). Two hundred thirty-five miRNAs were significantly differentially expressed (FC > 2.0 and BH-FDR adjusted P-value ≤ 0.05) in LUSC compared to adjacent normal tissues; 140 of them were significantly upregulated and 95 were significantly downregulated (Fig. 2B). Red and blue colors reflect high and low expression levels, respectively, based on the log2-transformed scale. Volcano plots illustrating the distribution of miRNA expression. The horizontal y-axis represents the log10 corrected P-value (tumor tissue/normal tissue), the red squares indicate upregulated miRNAs, and the blue squares indicate downregulated miRNAs. The Venn diagram represents common and type-specific DE miRNAs (Fig. 2C).[ The expression of DE miRNAs, which displayed significant statistical differences, was evaluated using a hierarchical clustering dendrogram and heatmap (Supplementary Figure 1 and 2, http://links.lww.com/MD/H37).

Figure 2.

Identification of differentially expressed miRNAs in NSCLC. (A) Volcano plot of DE miRNAs in LUAD compared with the adjacent normal; (B) Volcano plot in LUSC compared with the adjacent normal; The horizontal x-axis represents the log2 of the fold change (tumor tissue/normal tissue), the red squares indicate upregulated miRNAs and the blue squares indicate downregulated miRNAs. The y-axis represents the decimal logarithms of corrected P-value. Criteria: log2FC > 1.0 (fold change > 2) and BH-FDR ≤ 0.05; (C) The Venn diagram represents the number of unique DE miRNAs identified within and among 2 subtypes of NSCLC. LUAD—lung adenocarcinoma, LUSC—lung squamous cell lung cancer, DE—differently expressed.

3.2. Validation of miRNA microarray results by RT-qPCR in LUAD and LUSC

Of the 26 common DE miRNAs, 8 were prioritized and 6 subtype-specific miRNAs were selected: 3 from DE miRNAs in LUAD tissues and 3 from LUSC. All DE miRNAs were selected on the basis of statistically significant differences in expression in the current analysis and available data from the literature.

To validate the results of the miRNA microarray, we evaluated the expression levels of 14 miRNAs in an independent group of 50 LUAD and 50 LUSC tissue samples and their adjacent normal tissues by RT-qPCR. Those miRNA were miR-4689, miR-101-3p, miR-140-3p, miR-144-3p, miR-744-3p, miR-95-3p, miR-650, miR-6785-5p, miR-7-5p, miR-21-3p, miR-375, miR-139-5p, miR-195-5p and miR-195-3p. For box plot preparation, which showed miRNA expression in LUAD and LUSC, we used log2RQ values. Comparison between the levels of expression in the initial (microarray) and replication sets (RT-qPCR) was performed using the Kolmogorov-Smirnov test. The test indicated differences in the distribution of miRNA expression between the 2 datasets in both LUAD and LUSC groups. Figure 3A shows the expression levels of the selected miRNAs in LUAD and LUSC. According to the microarray results, red and blue boxplots represent the expression levels of miRNAs, and green and purple boxplots represent the expression data from RT-qPCR. We validated 8 DE miRNAs in the LUAD cohort by RT-qPCR. Three of them were upregulated - miR-7-5p, miR-375-5p and miR-6785-3p, and 5 downregulated - miR-101-3p, miR-139-5p, miR-140-3p, miR-144-3p and miR-195-5p. Ten DE miRNAs in the LUSC group were validated using RT-qPCR. Three were upregulated miR-7-5p, miR-21-3p, and miR-650, and 7 were downregulated miR-95-5p, miR-140-3p, miR-144-3p, miR-195-5p, miR-375, miR-744-3p, and miR-4689-3p. Five miRNAs were common between LUAD and LUSC.

Figure 3.

(A) Comparison between the levels of expression in the microarray and the RT-qPCR. According to microarray results, red and blue boxplots represented expression levels, and green and purple boxplots the expression data from RT-qPCR. The Kolmogorov-Smirnov test was applied, and P ≤ .05 indicated the difference in the distribution of miRNA expression between the 2 datasets. Eight differentially expressed (DE) miRNAs in LUAD cohort were validated with RT-qPCR - miR-7-5p, miR-375-5p, miR-6785-3p, miR-101-3p, miR-139-5p, miR-140-3p, miR-144-3p, and miR-195-5р. Ten DE miRNAs were validated with RT-qPCR in LUSC group: miR-7-5p, miR-21-3p, miR-650, miR-95-5p, miR-140-3p, miR-144-3p, miR-195-5p, miR-375, miR-744-3p, and miR-4689-3p; (B) Box plot and p values of comparison between the levels of expression in LUAD and adjacent nontumor tissues by RT-qPCR; (C) Box plot and p values of comparison between the levels of expression in LUSC and adjacent nontumor tissues by RT-qPCR; (D) Box plot and p values of comparison between the levels of expression in LUAD and LUSC by RT-qPCR. The y-axis shows the log2RQ values. Mann–Whitney U test was performed to evaluate the significance of expression level differences, P ≥ .05 consider as nonsignificant. LUAD = lung adenocarcinoma, LUSC = lung squamous cell lung cancer, MA = microarray.

The Mann–Whitney U test was performed to evaluate the significance of the differences in expression levels. The results demonstrated that the expression levels of 12 miRNAs differed significantly between LUAD and normal tissues (ranging from P ≤ .05, P = .0001). Only the expression levels of miR-650 and miR-195-5p were not significantly different between the LUAD and normal tissues (Fig. 3B). In LUSC, the expression levels of 13 miRNAs differed significantly between tumor and normal tissues, while only miR-95-3p was not significantly different from normal tissue (Fig. 3C). MiR-101-3p, miR-139-5p, miR-140-3p, and miR-144-3p were significantly downregulated in both LUAD and LUSC (Fig. 3 B–D). MiR-375 was significantly downregulated in LUSC tumors (Fig. 3C) and upregulated in LUAD tumors (Fig. 3B). In miRNA biogenesis, the stem-loop structure of miR-195 is further processed by AGO and gives rise to the guide strand miR-195-5p and the passenger strand miR-195-3p. This miRNA arm switching preference varies according to tissue type and development stage.[ The expression of miR-195-5p is significantly downregulated in LUSC tumors (Fig. 3C).

3.3. Association between miRNA expression levels and clinicopathological characteristics

Statistically significant associations between the expression levels of the studied miRNAs are presented in Supplementary figures 3 and 4, http://links.lww.com/MD/H38. The results of RT-qPCR indicated a statistically significant association between their expression levels and different clinical features in LUAD and LUSC. In LUAD, miR-101-3p and miR-195-5p showed statistically significant associations with tumor stage (P = .030 and P = .044, respectively) and metastasis (P = .009 and P = .035, respectively). The expression of miR-144-3p (P = .017) and miR-6785 (P = .050) was associated with tumor stage in LUAD. MiR-650, miR-7-5p, miR-375, miR-139-5p, and miR-195-3p were associated with both nodal stage and distant metastasis status and mir-744-3p (P = .014) and miR-95-3p (P = .048) were associated with only metastasis status (Supplementary figure 3, http://links.lww.com/MD/H38). In LUSC association with tumor stage were miR-4689 (P = .050), miR-101-3p (P = .015), miR-140-3p (P = .024), miR-144-3p (P = .032), miR-6785 (P = .023), miR-21-3p (P = .049), miR-375 (P = .040), miR-139-5p (P = .046) and miR-195-3p (P = .022). The expression of miR-744-3p (P = .048) was decreased in LUSC patients positive for metastasis, and we found a significant association between miR-95-3p expression and nodal and distant metastases (Supplementary figure 4, http://links.lww.com/MD/H38).

3.4. Correlation-based networks

To analyze the biological and functional relationships between miRNAs expressed in LUAD and LUSC, we conducted Spearman rank correlation to calculate the correlation coefficients and p-values. As illustrated in Figure 4A, in LUAD tumor tissue, there were strong correlations between miR-4689 and miR-195-3p, miR-101-3p and miR-144-3p, miR-6785-5p, and miR-195-3p, whereas in LUSC, miR-101-3p correlated with miR-21-3p and miR-144-3p, miR-140-3p with miR-21-3p, and miR-6785-5p with miR-195-3p (Fig. 4B). Dark pink and blue reflect high and low correlation coefficients, respectively. Only the significant correlations are shown.

Figure 4.

Correlation heatmap of miRNA expression in NSCLC. (A) Correlation heatmap of miRNA expression in LUAD; (B) Correlation heat map of miRNA expression in LUSC; Dark pink and blue colors reflect high and low correlation coefficient. Only significant correlations (P-value ≤ 0.05) are shown.

To further investigate the correlation between miRNAs, we designed correlation-based networks with a correlation coefficient threshold of ≥ 0.5 and a Benjamini-Hochberg FDR adjusted P-value ≤ 0.05. The networks are presented as graphs, where nodes (vertices) are connected by edges (links). The nodes illustrate miRNAs, and the edges are based on correlation coefficients.

In LUAD tumors, there were 2 defined clusters of miRNAs (Fig. 5A), whereas in LUSC tumors, there was no clear separation between clusters (Fig. 5B). The miRNAs with the highest number of significant interactions in LUAD were miR-101-3p (with 8 links), miR-140-3p, miR-95-3p, miR-195-5p, and miR-195-3p with 6 links, whereas for LUSC, the miRNAs with the most interactions were miR-744-3p and miR-21-3p with 8 edges, and miR-140-3p with 7 edges.

Figure 5.

Correlation-based networks in NSCLC. (A) Correlation-based networks in LUAD. (B) Correlation-based networks in LUSC. The threshold of the graph of CBN was a correlation coefficient of ≥ 0.5 and BH-FDR adjusted P-value ≤ 0.05.

3.5. Gene ontology (GO) annotation and pathway analysis

To understand the different biological functions of the replicated miRNAs in LUAD and LUSC, we analyzed their target genes validated by different functional analyses. In the current study, the gene set in LUAD comprised 8455 genes and that in LUSC comprised 7926 genes. The GO annotation of the target genes in LUAD showed that they are involved in biological pathways such as regulation of gene expression and transcription, apoptotic processes, cell cycle and migration, and macromolecule biosynthetic processes (Fig. 6A). The molecular functions of these genes revealed enrichment in RNA binding, cadherin binding, protein kinase activity, ubiquitin activity, and transcription regulatory region DNA binding (Fig. 6C). As for cell components, these genes were related to focal adhesion, nuclear body, centrosome, Golgi subcompartment, microtubule organizing, spindle microtubule, etc (Fig. 6E).

Figure 6.

Gene ontology classification of the predicted targets for differentially expressed miRNAs validated by RT-qPCR. (A) Top 15 results of GO BP enrichment in LUAD; (B) Top 15 results of GO BP enrichment in LUSC; (C) Top 15 results of GO CC enrichment in LUAD; (D) Top 15 results of GO CC enrichment in LUSC; (E) Top 15 results of GO MF enrichment in LUAD; (F) Top 15 results of GO MF enrichment in LUSC. The x-axis represents the decimal logarithm of the adjusted P-value.

The target gene set of LUSC showed genes mainly involved in BP, such as regulation of the apoptotic process, regulation of gene expression and transcription, protein ubiquitination and regulation of intracellular signal transduction, and protein modification by small protein conjugation (Fig. 6B). MF indicates RNA binding, protein kinase activity, cadherin binding, protein kinase binding, ubiquitin ligase activity, etc (Fig. 6D). The cell component showed enrichment of these genes in focal adhesion, nuclear body, Golgi subcompartment, chromatin, and cytoskeleton.

For the Reactome analysis, the LUAD dataset was significantly enriched in cell cycle, NGF signaling, membrane trafficking, and EGFR signaling (Fig. 7A). In contrast, the LUSC target gene set was significantly enriched in cell cycle pathways, human immune response, transcriptional regulation by TP53, and FGFR signaling (Fig. 7B).

Figure 7.

Reactome pathway analyses of target genes. (A) Top 20 results of Reactome pathway analysis in LUAD; (B) Top 20 results of Reactome pathway analysis in LUSC. The x-axis represents the decimal logarithm of the adjusted P-value.

4. Discussion

In the current study, we identified specific miRNA expression profiles in the 2 major subtypes of NSCLC and determined both specific and common miRNAs for each subtype. In the microarray samples for the LUAD subtype, 107 miRNAs were differentially expressed, of which 68 miRNAs were upregulated and 39 miRNAs were downregulated (Fig. 2A). In LUSC samples, there were 140 miRNAs with over 2-fold increased expression and 95 miRNAs with reduced expression (Fig. 2B). 26 miRNAs were common to both subtypes (Fig. 2C). Based on these results, we prioritized and selected 14 miRNAs for further validation by RT-qPCR. Concerning the intertumor heterogeneity characterized by stage, age, sex and ethnicity, we carefully selected a stratified validation cohort of patients from the same population, with the same proportions of tumors in each subtype as observed in the microarray dataset. Special care was taken to standardize the preanalytic phase in order to get high quality biological samples that will allow us to detect robust and novel findings relevant for the Bulgarian population and compare them with the existing literature data.

The validation with RT-qPCR confirmed the statistically significant change in expression levels for 8 of them in LUAD, with 3 of them upregulated (miR-7-5p, miR-375-5p and miR-6785-3p) and 5 downregulated (miR-101-3p, miR-139-5p, miR-140-3p, miR-144-3p and miR-195-5p) (Fig. 3A). In LUSC, the analysis confirmed the direction of change for 10 of them, with 3 upregulated (miR-7-5p, miR-21-3p, and miR-650) and 7 downregulated (miR-95-5p miR-140-3p, miR-144-3p, miR-195-5p, miR-375, miR-4689-3p, and miR-744-3p) (Fig. 3A).

We identified 5 overlapping miRNAs that were significantly dysregulated in both the subtypes (Fig. 3D). Four of them showed a common direction of change: miR-7-5p was upregulated; miR-140-3p, miR-144-3p and miR-195-5p were down-regulated. MiR-375 has been repeatedly shown to be upregulated in LUAD and downregulated in LUSC samples. Three miRNAs were significantly deregulated only in LUAD: miR-6785-3p, miR-101-3p, miR-139-5p, and 5 miRNAs in LUSC - miR-21-3p, miR-650, miR-95-5p, miR-4689-3p, and miR-744-3p. The expression data for miR-6785-3p and miR-4689-3p have not been previously reported.

The levels of miR-140-3p and miR-101-3p were downregulated in both LUAD and LUSC samples in the microarray and validation analyses, although the levels of miR-101-3p in LUAD were lower than those in LUSC (Fig. 3A). MiR-101-3p was significant associated with tumor stage in both NSCLC subtypes, and LUAD tumors were associated with distant metastasis. The expression of miR-140-3p was associated with tumor stage in LUSC (Supplementary figure 3 and 4, http://links.lww.com/MD/H38). Reactome enrichment of all selected target genes was associated with the cell cycle, gene expression, and related signaling pathways in both LUAD and LUSC, as well as in immune system-related genes for LUSC (Fig. 7A,B). Halverson et al analyzed miRNA expression specific for each subtype in an extensive study of 241 NSCLC cases from Norway, LUAD, and LUSC, and found that miR-140-3p and miR-101-3p were also among the downregulated miRNAs. The predicted targets of miR-140-3p and miR-101-3p were highly correlated with specific gene sets, and a distinct pattern of biological processes with high immune activity for the LUAD proximal inflammatory and LUSC secretory subtypes, and upregulation of cell cycle-related processes in the LUAD proximal proliferative subtype, LUSC primitive, and LUSC classical subtypes were found. Low levels of miR-140-3p and miR-101-3p most likely led to the upregulation of the G2M checkpoint. However, since miR-101-3p is expressed in both immune cells and lung tissue, the immune microenvironment could also contribute to these suppressed expression levels. Furthermore, since miR-140-3p is highly expressed by B cells, low levels of miR-140-3p may indicate the absence of B cell infiltration.[ JAK-STAT inhibition has been recently shown to impair K-RAS-driven lung adenocarcinoma progression and is a promising therapy for K-RAS-driven lung LUAC.[ Finally, miR-140-3p overexpression was found to attenuate cisplatin resistance and reduce stem cell-like properties in LUAD cells by suppressing Wnt/β-catenin signaling.[

In our samples, the expression of miR-144-5p was downregulated in both LUAD (Fig. 3B) and LUSC patients (Fig. 3C), with significantly lower expression observed in LUAD patients (Fig. 3D). The expression of miR-144-5p was associated with tumor stage in both LUAD and LUSC tumors (Supplementary figure 3 and 4, http://links.lww.com/MD/H38). Both strands of miR-144 (miR-144-5p, passenger strand; miR-144-3p, guide strand) showed decreased expression in cancer tissues. Low miR-144-5p expression significantly predicted poor prognosis in patients with LUSC based on TCGA database analysis.[ Downregulated miR-144-3p contributes to the progression of lung adenocarcinoma by elevating the expression of the histone methyltransferase EZH2.[ EZH2 expression is upregulated in LUAD cells and promotes LUAD cell invasiveness and metastasis.[ MiR-144-3p effectively suppressed TGF-β1-induced aggressive lung cancer cell invasion and adhesion by regulating the Src–Akt–Erk pathway.[ In LUAD but not in LUSC patients, the serum IL-1β level was correlated with miR-144-3p, suggesting that it might affect miR-144-3p at the transcriptional level. Both were independent risk factors for LUAD prognosis. Therefore, although IL-1β and miR-144-3p might mediate inflammation-promoted tumorigenesis in patients with LUAD, the exact mechanism and their interaction remain unclear.[ MiR-144-3p regulates cisplatin resistance in lung cancer cells by targeting Nrf2 Furthermore, overexpression of Nrf2 can regulate the expression of miR-144-3p by binding to ARE in the miR-144-3p promoter region.[ A newly identified lncRNA, GAS6-AS2, upregulated MAPK6 by sponging miR-144-3p in NSCLC tissues and cells.[ Thus, GAS6-AS2 is an effective therapeutic target for NSCLC treatment. Upregulation of the lncRNA X-Inactive Specific Transcript (XIST) is associated with cisplatin resistance in NSCLC by downregulating miR-144-3p.[

The expression of miR-195-5p was significantly downregulated in both LUAD and LUSC and was associated with the presence of distant metastasis in both NSCLC and in LUAD subtypes with tumor stage (Supplementary figure 3, http://links.lww.com/MD/H38). MiR-195-5p has recently been reported to act as a tumor suppressor in various cancers, including NSCLC. It has been shown that miR-195-5p has lower expression in lung cancer and may act as a diagnostic biomarker.[ MiRNA-195-5p may regulate lung cancer growth and metastasis by regulating its direct target, FOXK1 (Forkhead box k1).[ Overexpression of FOXK1 and its interaction with the FHL2 and FOXK1-CCDC43 axis in colorectal cancer cells are associated with EMT, invasion, and metastasis.[ Functional assays showed that miR-195-5p overexpression inhibits the invasion, proliferation, and migration of NSCLC-derived cells by directly targeting CPNE1, a calcium-dependent phospholipid-binding protein. Pathway analysis revealed decreased expression of p-AKT, p-ERK, and Snail after transfection with miR-195-5p mimics in lung adenocarcinoma and squamous cell lines.[ However, the mechanisms by which miR-195 represses the tumorigenesis of NSCLC cells are not fully understood, and data are gradually accumulating. MiR-195-5p regulates cell cycle progression, migration, and invasion of NSCLC cells by targeting TrxR2, CEP55, VEGF, MYB, CHEK1, HDGF, or IGF1R.[ By targeting cyclin D3, miR-195-5p causes cell cycle arrest at the G1 phase and survivin induces apoptosis and senescence in NSCLC cells. Interestingly, BIRC5, which codes for survivin, is upregulated in both adenocarcinoma and squamous cell carcinoma tissues, and high expression of BIRC5 is associated with poor survival in LUAD but not LUSC.[

The miRNA in our analysis revealed that the most evident and consistent upregulation in LUAD and downregulation in LUSC samples, both in the training and validation samples, was miR-375 (Supplementary Figure 1 and 2, http://links.lww.com/MD/H37; Fig. 3D). In the microarray experiment, the expression in LUAD was upregulated, with 3,23 FC, while in LUSC, it was downregulated with 11,7 FC. Jin et al found that miR-375 expression was significantly upregulated in LUAD and small cell lung carcinoma but downregulated in LUSC. In our study, the expression of miR-375 was also upregulated in LUAD, although not statistically significant, and was downregulated in LUSC. The difference between the expression levels of miR-375 in both tumor tissues quantified by RT-qPCR was in the same direction as the microarray results (Fig. 3A) and was highly significant (P = .0001) (Fig. 3D). The expression of miR-375 was associated with nodal stage and distant metastases in LUAD. In LUSC, miR-375 was associated with tumor stage (Supplementary figure 4, http://links.lww.com/MD/H38).

Many previous studies have reported that miR-375 expression is significantly upregulated in adenocarcinoma and small cell lung carcinoma, but downregulated in squamous cell carcinoma.[ Further investigation revealed that miR-375 directly targets ITPKB mRNA and promotes cell growth.[ In addition, other targets were identified, such as crucial members of established signaling pathways in lung cancer, including calcium, insulin, Jak-STAT, PPAR, MAPK, mTOR, TGF-β, and Wnt pathways. Reduced miR-375 is correlated with tumor progression and poor prognosis of oral squamous cell carcinoma patients, and functional analysis showed that miR-375 binds to the 3’-UTR of the insulin-like growth factor 1 receptor.[ Extensive bioinformatics analysis of public datasets, such as TCGA and the Gene Expression Omnibus database, followed by GO enrichment and target prediction analysis, revealed critical pathways related to carbohydrate metabolism. Suggestions were made regarding miR-375 associated genes that might participate in LUAD by network analysis, such as FGF2, PAX6, and RHOJ.[ The lncRNA ROR1-AS1 enhances lung adenocarcinoma metastasis and induces epithelial-mesenchymal transition by sponging miR-375.[

We found that miR-7-5p was significantly upregulated in LUAD and LUSC tumor tissues. MiR-7-5p connects different signaling pathways through feedback and feedforward loops and is involved in differentiation, regulation of proliferation, apoptosis, migration, chromatin remodeling, and multidrug resistance.[ MiR-7 acts mainly as a tumor suppressor in various cancer types[; however, evidence suggests that it may also act as an oncogene in lung cancer.[ Elevated miR-7 expression levels in lung tumors compared to normal tissues and metastases have been found in previous studies.[

In a study by Chou et al, EGFR increased the expression of miR-7 through the Ras/ERK/Myc pathway, which suppressed the transcriptional repressor ERF. Interestingly, EGFR is a direct target of miR-7, and abnormalities in this fine-tuned regulatory mechanism may contribute to its oncogenic ability. The PI3K/AKT pathway also regulates miR-7 expression; however, the exact mechanism must be explored further.[ LINC00240 suppresses invasion and migration in nonsmall cell lung cancer by sponging miR-7-5p.[ In our study, the expression of miR-7 was higher in patients without lymph node metastasis and was associated with the presence of distant metastasis in LUAD (Supplementary figure 3, http://links.lww.com/MD/H38). Cheng et al observed that low miR-7 expression levels were associated with poorer tumor differentiation, advanced tumor stage, higher incidence of lymph node metastasis, advanced p-TNM stage and shorter 5-years overall survival rate.[ In our group of LUAD and LUSC patients, mainly tumor stage I-II were included, and most of them were without lymph nodes and distant metastases. Tumor heterogeneity may also explain the observed upregulation of miR-7 in LUAD and LUSC samples in the current study.

Hong et al identified the chromatin remodeling protein SMARCD1 as a target of miR-7 in lung cancer cells, leading to chemoresistance in lung cancer cells by affecting the association of SMARCD1 with p53.[ Poly ADP-ribose polymerase 1 (PARP1) is a direct target of miR-7-5p and PARP1 expression is downregulated by miR-7-5p. MiR-7-5p impeded Dox-induced HR repair by inhibiting the expression of HR repair factors (Rad51 and BRCA1), which resulted in the resensitization of SCLC cells to doxorubicin. MiR-7-5p targets PARP1 to exert its suppressive effects on HR repair, indicating that alteration of the expression of miR-7-5p may be a promising strategy for overcoming chemoresistance in SCLC therapy.[

MiR-21-3p and -5p were generated from the same precursor (premiR-21). However, the role and biological function of the passenger strand miR-21-3p are not clearly understood. According to previous research, miR-21 is one of the most frequently overexpressed miRNAs in cervical cancer, laryngeal squamous cell carcinoma, breast cancer, colon cancer,[ and lung cancer tissues.[

Our microarray results showed increased expression of both strands. MiR-21-3p was selected for the validation study, as it was significantly elevated in LUAD and LUSC samples in the microarray experiment (4,2 FC in LUAD and 2,34 FC in LUSC) (Supplementary figure 1,2, http://links.lww.com/MD/H37). The expression levels of miR-21-3p were increased in both LUAD and LUSC samples compared to normal tissue samples in the validation study, although the differences were not as pronounced for LUSC samples (Fig. 3 A-D).

GO annotation of the target genes in the LUSC sample showed that the top 15 BPs were related to the regulation of gene expression and apoptosis. Several enriched terms in the reactome were associated with transcriptional regulation by TP53 and the immune system, such as cytokine signaling and interferon signaling. Overexpression of miR-21 affects signaling pathways related to apoptosis, cell proliferation, survival, angiogenesis, inhibition of nuclear factor kB (NF-kB), and the PI3K/Akt/NF-κB signaling pathway in NSCLC in vitro and in vivo.[ Direct targets of miR-21 include PDCD4, Smad7, Bcl2, EGFR, Cas8, TGF-b.[ miR-21 expression level is associated with smoking, lymph node metastasis, and poor prognosis in NSCLC.[

Moreover, dysfunctional miR-21-3p has been reported to promote cancer cell proliferation and invasive ability, as well as to increase cisplatin resistance by targeting RNA-binding protein with multiple splicing (RBPMS), a regulator of chromosome condensation and POZ domain-containing protein 1 (RCBTB1), and zinc finger protein 608 (ZNF608) in ovarian cancer cells.[ The processes affected by these genes include the regulation of transcriptional activity through SMAD proteins (RBPMS), cell cycle and transcriptional control through chromatin remodelling (RCBTB1), and transcriptional regulation (ZNF608). Recent studies have reported that miR-21-3p overexpression can induce microsatellite instability by targeting MSH2, a key mismatch repair gene.[

In our sample set, miR-139-5p was consistently downregulated in LUAD and LUSC compared with that in normal adjacent tissues. In the validation sample, the expression levels were in the same direction, but were significantly different between the LUAD and LUSC samples (Fig. 3D). MiR-139-5p is a common tumor-associated miRNA that may act as a tumor suppressor gene during the occurrence and development of NSCLC. MiR-139-5p downregulation was associated with lymph node metastasis and distant metastasis in LUAD and tumor stage in LUSC (Supplementary Figures 3 and 4, http://links.lww.com/MD/H38). Our results further supported the involvement of miR-139-5p in the development of this malignancy. The mechanisms of action of miR-139-5p include inhibiting lung adenocarcinoma cell proliferation, migration, and invasion by targeting MAD2L1 and HDGF in NSCLC and CXCR4 in laryngeal squamous carcinoma cells.[ In another study, miR-139-5p was shown to play a pivotal role in lung cancer by inhibiting cell proliferation and metastasis and promoting apoptosis by targeting oncogenic c-Met.[ Differential expression of miR-139-5p plays a significant role in tumorigenesis, metastasis, and recurrence, suggesting that it may be used as a promising biomarker for cancer diagnosis, prognosis, and therapy.[ Recently, the role of miR-139-5p in cisplatin (DDP) sensitivity in nonsmall cell lung cancer (NSCLC) cells was investigated. It has been shown to enhance cisplatin sensitivity by inhibiting cell proliferation and promoting apoptosis by targeting the homeobox protein HoxB2.[ Monitoring the expression levels of miR-139-5p and its effects could be a novel approach to reverse DDP resistance and increase chemosensitivity in the treatment of NSCLC.

Compared with adjacent nontumor tissues, miR-650 was significantly upregulated in LUAD and LUSC, and was significantly higher in LUAD than in LUSC (Fig. 3B-D). MiR-650 deregulation has been reported in gastric cancer, glioma, chronic lymphocytic leukemia, and hepatocellular carcinoma tissues and is considered an miRNA with prognostic and predictive value.[ In a study by Huang et al, LUAD patients with low miR-650 expression had significantly better prognosis.[ In addition, the expression of miR-650 was negatively correlated with the response to docetaxel by regulating Bcl-2/Bax expression by targeting ING4 and stimulating cell proliferation and invasion in the Wnt-1/β-catenin pathway.[

In our LUSC patients, we observed significantly decreased expression of miR-744 (P = .017) (Fig. 3C) and increased expression in LUAD samples compared to that in normal tissues (P = .017) (Fig. 3B). This finding suggests a different role of miR-744 in the carcinogenesis of both types of NSCLC. Altered expression of miR-744 has been reported in numerous cancers. In prostate cancer, miR-744 promotes cancer progression by activating Wnt/β-catenin signaling through binding to the negative regulators SFRP1, GSK3β, TLE3, and NKD1[ and by regulating LKB1 expression.[ This suggests an oncogenic role for this miRNA. In glioblastoma, colorectal cancer, hepatocellular carcinoma, and ovarian cancer, miR-744 serves as a tumor suppressor by inhibiting proliferation and invasion and inducing cell death.[Few studies have investigated the role of miR-744-5p in NSCLC, particularly in cell lines.

Chen et al found significantly decreased expression of miR-744-5p in NSCLC cell lines, and in vitro studies showed that overexpression of miRNA could inhibit cell proliferation, colony formation, and cell invasion.[ Wang et al illustrated the carcinogenic role of XIST in promoting tumorigenesis and the progression of NSCLC by targeting miR-744 and increasing RING1 expression by activating the Wnt/β-catenin signaling pathway.[ Sui et al found that MAFG-AS1 and MAFG sponge miR-744-5p, and their upregulation could promote cell growth in LUAD samples.[ In our study, the expression of miR-744-5p was significantly decreased in patients with metastases (Supplementary figure 3 and 4, http://links.lww.com/MD/H38). Zhu et al analyzed the data of 456 primary NSCLC samples from TCGA (validation set) and 87 formalin-fixed paraffin-embedded NSCLC tissues (training set), and concluded that upregulated miR-744 was associated with poor clinical outcomes. MiR-744 increased the invasive ability of lung cancer cell lines, whereas knockdown of endogenous miR-744 decreased the proliferation and migration of the cells. Furthermore, miR-744 could directly and indirectly upregulate c-FOS expression via the MAPK pathway and contribute to carcinogenesis. They explored the in vitro effect of miR-744 in boosting the chemoresistance and radioresistance of lung cancer cells.[

To our knowledge, this is the first study to show significant changes in the expression of miR-744, miR-4689-3p and miR-6785-3p in fresh frozen LUAD and LUSC tissue samples. The expression of the miR-6785-3p was higher in T3-T4 of LUAD tumors and the levels of both miR-6785-3p and miR-4689-3p were associated with advanced LUSC tumors (Supplementary figure 3 and 4, http://links.lww.com/MD/H38). It is worth mentioning that the correlation-based network analysis in LUSC revealed that miR-744-3p was among the miRNAs with the most interactions, together with miR-21-3p and miR-140-3 (Fig. 5B). Therefore, we suspected that the expression levels of miR-744-3p and its signaling pathways play a vital role in the carcinogenesis of LUSC. However, further functional analysis in larger samples and validation is necessary to better ascertain the biological function and prognostic potential of miR-744 in LUAD and LUSC.

Correlation-based network analysis revealed 2 defined clusters of miRNAs in LUAD tumors (Fig. 5A). Analysis of miRNAs with the highest number of significant interactions revealed several key hub miRNAs in LUAD. Intriguingly, the first smaller cluster was centered around miR-195-3p, while the second contained miR-195-5p and the other hub miRNAs, miR-101-3p, miR-140-3p, miR-95-3p, and miR-375, which were used to differentiate between LUAD and LUSC. In LUSC on the other hand, key hub miRNAs are again miR-195-5p and miR-744-3p to form 1 big cluster (Fig. 5B). Although there are several overlapping miRNAs, the differences between the clusters revealed in LUAD and LUSC further point to the complex nature of molecular pathway deregulation and specific interactions, even when comparing the 2 lung cancer subtypes. In both NSCLC subtypes, miR-195-5p role has been previously reported, however this is the first time that the role of miR-195-3p is revealed in LUAD, based on the defined clusters following network analysis. The study of interactions between multiple miRNAs provides an opportunity to better understand the subtype differences.

However, our correlation analysis has the limitations of all in silico methods. Correlation-based network analysis draws our attention to key miRNAs whose expression, regulated pathways, and complex network interactions might help disentangle the similarities and differences between the 2 major subtypes of lung cancer. However, bioinformatics analyses need to be further supported by real functional or expression data to estimate the precise role of the prioritized miRNAs in lung carcinogenesis and to determine their potential role as valuable biomarkers.

5. Conclusions

In conclusion, this study identified and validated distinct miRNA profiles for LUSC and LUAD, including both common and subtype-specific miRNAs. We found 5 commonly deregulated miRNAs, with miR-7-5p being upregulated, while miR-140-3p, miR-144-3p and miR-195-5p were down-regulated. Our study adds to the accumulated evidence of the differential expression of miR-375 in the 2 lung cancer subtypes, which was upregulated in LUAD and downregulated in LUSC. In addition, 3 miRNAs (miR-6785-3p, miR-101-3p, and miR-139-5p) were identified with specific deregulation in LUAD only. Five miRNAs were found to be significantly deregulated in LUSC only: miR-21-3p, miR-650, miR-95-5p, miR-4689-3p, and miR-744-3p.

Based on the analysis of the target genes of the significantly deregulated miRNAs, both common and specific molecular pathways and biological processes in LUAD and LUSC were identified, adding new knowledge to our understanding of the complex mechanisms of carcinogenesis in lung cancer. This could serve as the basis for future biomarker studies for precise diagnosis and prognosis, contribute to the understanding of resistance to current therapy, and help in the search for new and more efficient drug targets.

Author contributions

All authors participated in conceiving the concept and writing the manuscript.