MiRNA and Potential Prognostic Value in Non-Smoking Females with Lung Adenocarcinoma by High-Throughput Sequencing

Hao Wang,¹ Lijun Wang,² Gengyun Sun¹

¹Department of Respiratory and Critical Care Medicine, The First Affiliated Hospital of Anhui Medical University, Hefei, People’s Republic of China; ²Department of Respiratory Disease, Tongling People’s Hospital, Tongling, People’s Republic of China

Correspondence: Gengyun Sun, Department of Respiratory and Critical Care Medicine, The First Affiliated Hospital of Anhui Medical University, No. 218, Jixi Road, Hefei, Anhui, 230022, People’s Republic of China, Email [email protected]

Background: Non-smoking females with lung adenocarcinoma (LUAD) account for a unique disease entity and miRNA play critical roles in cancer development and progression. The purpose of this study is to explore prognosis-related differentially expressed miRNA (DEmiRNA) and establish a prognostic model for non-smoking females with LUAD.
Methods: Eight specimens were collected from thoracic surgery of non-smoking females with LUAD and implemented the miRNA sequencing. The intersection of our miRNA sequencing data and TCGA database were identified as common DEmiRNA. Then, we predicted the target genes of the common DEmiRNAs (DETGs) and explored the functional enrichment and prognosis of DETGs. A risk model by overall survival (OS)-related DEmiRNA was constructed based on multivariate Cox regression analyses.
Results: A total of 34 overlapping DEmiRNA were obtained. The DETGs were enriched in pathways including “Cell cycle” and “miRNAs in cancer”. The DETGs (KPNA2, CEP55, TRIP13, MYBL2) were risk factors, significantly related to OS, progression-free survival (PFS), and were also hub genes. ScRNA-seq data also validated the expression of the four DETGs. Hsa-mir-200a, hsa-mir-21, and hsa-mir-584 were significantly associated with OS. The prognostic prediction model constructed by the 3 DEmiRNA could effectively predict OS and can be used as an independent prognostic factor of non-smoking females with LUAD.
Conclusion: Hsa-mir-200a, hsa-mir-21, and hsa-mir-584 can serve as potential prognostic predictors in non-smoking females with LUAD. A novel prognostic model based on the three DEmiRNAs was also constructed to predict the survival of non-smoking females with LUAD and showed good performance. The result of our paper can be helpful for treatment and prognosis prediction for non-smoking females with LUAD.

Keywords: women, lung adenocarcinoma, miRNA, prognosis, bioinformatics

Introduction

Lung cancer is an issue of increasing global concern, with the most incidence and mortality worldwide.¹ About 80% of the lung cancers are non-small cell lung cancer (NSCLC) and lung adenocarcinoma (LUAD) is the most common subtype. Smoking is a major risk factor for lung carcinogenesis²; however, about a quarter of lung cancer cases are not related to smoking, and it is reported that more than half of female patients never smoke.³ Furthermore, the number of LUAD in female never smokers is gradually increasing.^4,5 Non-smoking females with LUAD account for a unique disease entity, showing different pathogenesis, epidemiological characteristics and course of disease.⁶ Therefore, it is crucial to explore the molecular mechanism that drives LUAD in female never smokers. Exploring prognostic factors and models can be helpful for treatment and prognosis prediction for non-smoking females with LUAD.

miRNAs are a highly conserved class of approximately 19–24 nucleotide-long non-protein coding RNA, which monitor cell proliferation, apoptosis, tumor growth and metastasis by negatively regulating gene expression.^7,8 Aberrant expressions of miRNAs have widely been reported to regulate NSCLC progression,⁹ including cell growth, stemness, drug resistance, and immune evasion.¹⁰ miRNAs could also act as biomarkers for NSCLC diagnosis, prediction of clinical outcomes.¹¹ Risk score model constructed by miRNAs can well predict lung cancer prognosis.¹² However, miRNAs in non-smoking females with LUAD have not been well explored.

In the present study, we investigated the differentially expressed miRNAs (DEmiRNAs) and explored the GO, KEGG, OS and PFS of their target genes. Additionally, we evaluated the prognostic value of the DEmiRNAs, and constructed three miRNA (hsa-mir-200a, hsa-mir-21 and hsa-mir-584) based risk model that could effectively predict overall survival (Figure 1).

Figure 1 Flow chart of the analysis process in our study.

Abbreviations: DEmiRNA, differentially expressed miRNA; DETGs, differentially expressed target genes; PPI, protein–protein interaction; OS, overall survival; PFS, progression-free survival; ROC, receiver operating characteristic curves; DETGs, target genes of the common differentially expressed miRNA.

Materials and Methods

Patients and Sample Collection

All tissue samples were collected from non-smoking females with LUAD who were treated by thoracic surgery in Anhui Provincial Cancer Hospital from June 2020 to June 2021. Patients treated with preoperative targeted therapy, chemotherapy, immunotherapy, and radiotherapy were excluded. The pathological type of lung adenocarcinoma was confirmed by pathologists. Tumor specimens and matched normal tissues were quickly frozen at −80 °C for further high-throughput sequencing. The present study was approved by The Human Investigation and Ethical Committees of Anhui Provincial Cancer Hospital (Ethics Number-41/2020) and performed according to the Declaration of Helsinki. All participants provided written informed consent.

Small RNA Sequencing and DEmiRNA Analysis

Total RNA was extracted from fresh tissues using TRIzol according to the manufacturer’s instructions (Thermo Fisher Scientific, Inc.). The purity and concentration of RNA were examined using NanoDrop 2000. The integrity and quantity of RNA were measured using the Agilent 2100/4200 system. Three micrograms of total RNA per sample was used to generate the small RNA library. Small RNA sequencing was performed using Illumina NextSeq500sequencing (BG, Beijing, China). Low quality reads were removed. Adaptor sequences and small RNAs which perfectly match in the human genome database (ftp://ncbi.nlm.nih.gov/genomes/Homo_sapiens) were considered valid for subsequent analysis. The “edgeR” R package was utilized to find the differentially expressed miRNA (DEmiRNA) with the threshold set at absolute log₂Fold Change >1 and P value <0.05.

Validation of DEmiRNA by TCGA

The miRNA expression data, mRNA expression data, single nucleotide polymorphism (SNP) data and corresponding clinical features of non-smoking females with lung adenocarcinoma were also downloaded from The Cancer Genome Atlas (TCGA). The DEmiRNA between tumor samples and normal tissues were determined using “edgeR” R package. The log₂ Fold Change and P value cutoff were set as 1 and 0.05. We applied “VennDiagram” R package to identify the common DEmiRNA in our miRNA-seq data and TCGA miRNA data. The common DEmiRNA was used for the following analysis. Tumor mutation burden (TMB) is a measure of the total number of somatic non-synonymous mutations per base pair. SNP data from TCGA was utilized to calculate TMB.

Common DEmiRNA and Their Target Genes

The common DEmiRNAs were predicted their target mRNA using StarBase (http://starbase.sysu.edu.cn/), miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/), miRDB (http://www.mirdb.org/), and TargetScan (http://www.targetscan.org/) datasets. The target genes were submitted to STRING v11.0 (https://www.string-db.org/) database, for the exploration of protein–protein interaction (PPI) network. Cytoscape software (version 3.8.0) was used to perform visualization of PPI network. Cytohubba, a plug-in of Cytoscape, was applied to find the key genes in the PPI network. Target genes were also inputted into DAVID Bioinformatics Resources 6.8 (https://david.ncifcrf.gov/) online database, to present Gene Ontology (GO) analysis.

Prognosis-Related Common DEmiRNA and Target Genes

Target genes’ expression level and corresponding clinical characteristics were obtained from TCGA dataset. We employed univariate Cox proportional-hazards regression analyses to evaluate the progression-free survival (PFS) and overall survival (OS) of common DEmiRNA and target genes. P value less than 0.05 was considered statistically significant.

The Expression of DETGs Was Validated by scRNA-Seq

The scRNA-seq data of 19 non-smoking women with lung adenocarcinoma patients (13 tumor tissues and 6 normal tissues) were also downloaded from the GEO database (GEO 131907). The package of Seurat was applied to process the scRNA-seq data. The data of the samples were batch corrected by the package of harmony. The FindVariableFeatures was used to find the top 2000 highly variant genes and performed principal component analysis. The first 30 principal components and resolution 0.6 were used to find clusters function to generate clusters. The FindAllMarkers were used to identify marker genes in each cluster.

Construction of the Risk Score by DEmiRNA

The univariate Cox proportional hazards regression analyses were employed to evaluate the association between OS and DEmiRNA. DEmiRNA with a p-value less than 0.05 was selected for construction the risk prognostic model. Then, we employed multivariate Cox proportional hazards regression analysis to produce coefficients which were acted as weights in the prognostic model. The risk scores equal to the sum of the product of each DEmiRNA and its weight. Based on the prognostic model, we divided patients into high-risk or low-risk groups. In order to test the robustness of the risk prognostic model, we randomly divided all samples of non-smoking females with LUAD into train-dataset and test-dataset. We conducted “survival” package of R to draw Kaplan–Meier (KM) curve, which was used to compare the OS between the high-risk group and low-risk group. We also applied “timeROC” package of R to draw the receiver operating characteristic (ROC) curves, which were used to validate the accuracy of the risk model in predicting the survival prognosis of the non-smoking females with lung adenocarcinoma.

Independent Prognostic Analysis

We conducted univariate and multivariate Cox proportional hazard regression analyses to assess whether the prognostic model we constructed, independent of the clinical characteristics (age, primary tumor, regional lymph node, stage, tumor mutation burden) in predicting OS of non-smoking females with lung adenocarcinoma.

Statistical Analysis

All statistical analyses were performed by statistical package R (version 4.1.0, https://www.r-project.org/). Fisher’s exact tests were used to evaluate the target genes in Gene Ontology. Univariate and multivariate Cox regression analyses were implied to assess whether the fact was related to prognosis. P value less than 0.05 was considered statistically significant.

Results

Differentially Expressed miRNA in Sequencing Data

We collect eight samples from thoracic surgery specimens (Supplementary Data 1). After completing sequencing and data quality control, we applied “edgeR” R package to identify differentially expressed miRNA between tumor and normal tissues in non-smoking females with LUAD. Finally, 120 differentially expressed miRNAs (DEmiRNA) were obtained, in which 66 DEmiRNAs upregulated, and 54 DEmiRNAs downregulated. The top 10 up- and downregulated DEmiRNAs were visualized by heat map, as shown in Figure 2.

Figure 2 The top 20 differentially expressed miRNA between tumor and corresponding normal tissues in sequencing data.

Identification of DEmiRNA in TCGA and Common DEmiRNA

Through screening of TCGA clinical data, we finally obtained 73 non-smoking females with lung adenocarcinoma and 3 normal samples. For the purpose of obtaining DEmiRNA, we conducted “edgeR” R package. At last, we got 196 DEmiRNA (70 upregulated, 126 downregulated) in TCGA data. We performed the top 10 up and down regulated DEmiRNAs by heat map, as shown in Figure 3A. In order to increase the credibility of DEmiRNA, we used the Venn diagram to select the intersection part of DEmiRNA in TCGA and our miRNA sequencing data, a total of 34 overlapping DEmiNRAs were obtained, as shown in Figure 3B.

Figure 3 The top 20 differentially expressed miRNA in TCGA and common DEmiRNA. (A) The top 20 differentially expressed miRNA (DEmiRNA) in TCGA. (B) Venn diagram shown the Common DEmiRNA.

PPI Network and Functional Enrichment of Differentially Expressed Target Genes (DETGs)

By predicting the target genes of 34 common miRNAs, 863 target genes were obtained. We used the “limma” package of R to explore the differentially expressed target genes (DETGs). Then, 123 DETGs (60 upregulated, 63 downregulated) in TCGA data set were acquired. This study performed the top 10 up and down regulated DETGs by heat map, as shown in Figure 4. We also explored the protein–protein interaction (PPI) network of DETGs. The PPI network was visualized by Cytoscape. There were 141 edges, 80 nodes, and 14 hub genes in the PPI network, as shown in Figure 5. In order to explore the influence of DEmiRNA on the pathogenesis in non-smoking females with lung adenocarcinoma, GO analysis on DETGs were performed. The 123 DETGs were grouped into BP (biological process), MF (molecular function), and CC (cellular component) categories. The most enriched GO terms in the biological process(BP) category was “signal transduction”, “positive regulation of RNA polymerase II promoter”. The most enriched GO terms in the cellular component (CC) category were “cytoplasm”, ‘nucleus’. The most enriched GO terms in the molecular function (MF) category were “protein binding” and “transcription factor activity, sequence−specific DNA binding” (Figure 6).

Figure 4 The heat map of the top 20 DETGs and 4 hub genes. The four hub genes (KPNA2, CEP55, TRIP13, MYBL2) were marked in red. CEP55 and MYBL2 are not only in the top 20 genes but also hub genes.

Figure 5 The protein–protein interaction network of DETGs. The small circle represents hub-genes. The red color represents up-regulated genes. The green color represents down-regulated genes.