MicroRNA and mRNA Interaction Network Regulates the Malignant Transformation of Human Bronchial Epithelial Cells Induced by Cigarette Smoke.

This study analyzes the correlation and interaction of miRNAs and mRNAs and their biological function in the malignant transformation of BEAS-2B cells induced by cigarette smoke (CS). Normal human bronchial epithelial cells (BEAS-2B) were continuously exposed to CS for 30 passages (S30) to establish an in vitro cell model of malignant transformation. The transformed cells were validated by scratch wound healing assay, transwell migration assay, colony formation and tumorigenicity assay. The miRNA and mRNA sequencing analysis were performed to identify differentially expressed miRNAs (DEMs) and differentially expressed genes (DEGs) between normal BEAS-2B and S30 cells. The miRNA-seq data of lung cancer with corresponding clinical data obtained from TCGA was used to further identify lung cancer-related DEMs and their correlations with smoking history. The target genes of these DEMs were predicted using the miRDB database, and their functions were analyzed using the online tool "Metascape." It was found that the migration ability, colony formation rate and tumorigenicity of S30 cells enhanced. A total of 42 miRNAs and 753 mRNAs were dysregulated in S30 cells. The change of expression of top five DEGs and DEMs were consistent with our sequencing results. Among these DEMs, eight miRNAs were found dysregulated in lung cancer tissues based on TCGA data. In these eight miRNAs, six of them including miR-96-5p, miR-93-5p, miR-106-5p, miR-190a-5p, miR-195-5p, and miR-1-3p, were found to be associated with smoking history. Several DEGs, including THBS1, FN1, PIK3R1, CSF1, CORO2B, and PREX1, were involved in many biological processes by enrichment analysis of miRNA and mRNA interaction. We identified the negatively regulated miRNA-mRNA pairs in the CS-induced lung cancer, which were implicated in several cancer-related (especially EMT-related) biological process and KEGG pathways in the malignant transformation progress of lung cells induced by CS. Our result demonstrated the dysregulation of miRNA-mRNA profiles in cigarette smoke-induced malignant transformed cells, suggesting that these miRNAs might contribute to cigarette smoke-induced lung cancer. These genes may serve as biomarkers for predicting lung cancer pathogenesis and progression. They can also be targets of novel anticancer drug development.

Introduction

Lung cancer is one of the most common carcinomas in men and women around the world. It is the first and second leading cause of cancer-related deaths in men and women, respectively (1, 2). There were 2.09 million new lung cancer cases and 1.76 million lung cancer deaths, which accounts for about 18.4% of all cancer deaths around the world in 2018 (3). The incidence of lung cancer is closely associated with cigarette smoking (2, 4). The risk of developing lung cancer in smokers is nearly ten times higher than that in non-smokers (5, 6). However, it is still not clear how normal lung epithelial cells become cancerous change in cigarette smokers.

It is well-known that the initiation and development of lung cancer are associated with abnormal expression of oncogenes and tumor suppressor genes. Numerous evidence suggested that the change in gene expression, which affects the occurrence and progression of tumors is closely related to epigenetic modification (7). Epigenetic modification could be DNA methylation, microRNAs (miRNAs), histone modifications, and nucleosome remodeling. These modifications are independent but could interact with each other to regulate gene expression (8). Epigenetic disruptions could promote the acquisition of a cancerous phenotype and aggressive behavior in lung cancer cells as well as primary or acquired resistance to treatment (9).

MiRNAs are highly conserved non-coding RNAs and consist of 18–24 nucleotides (nt) that are involved in the post-transcriptional regulation of gene (10). An individual miRNA is able to regulate many different transcripts. It is also believed that miRNAs can regulate more than one in three coding RNAs in the genome (11). MiRNAs participate in many vital biological processes through pairing with target mRNAs and regulating their expression (12, 13). The imbalance of miRNAs is usually associated with the progression and suppression of cancer, suggesting that miRNAs may play important roles as oncogenes or tumor suppressors (14).

The rapid development of high-throughput next-generation sequencing technology made it possible to identify changes in single bases in coding sequences of specific genes during lung tumorigenesis. A meticulous and thorough analysis of these data identified various vital genes and signaling pathways related to the tumor resulting in a better understanding of the mechanism of occurrence, development, and prognosis of lung cancer. Using novel technology and bioinformatics analysis, The Cancer Genome Atlas (TCGA, http://cancergenome.nih.gov/) project has previously identified panels of genetic mutations contributed to or were associated with the cause of a variety of cancers (15). Recently, the TCGA had shown studies on lung adenocarcinoma (LUAD) and squamous cell carcinomas (LUSC) at the molecular level (16, 17).

The aim of this study is to analyze the correlation and regulating mechanism of the regulatory network of miRNAs and mRNAs during carcinogenesis. An in vitro cell model of malignant transformation was established by exposing normal lung epithelial cells BEAS-2B to cigarette smoke (CS). Using high-throughput sequencing analysis, we analyzed the miRNA and mRNA expression profile in BEAS-2B cells with or without CS exposure. The differential expression miRNAs (DEMs) and differentially expressed genes (DEGs) were selected, and the integrative miRNA-mRNA network was analyzed. We identified some critical genes involved in carcinogenesis. This study will provide potential target candidates for novel drug development.

Methods

Preparation of Malignantly Transformed Cells

The CS-exposed malignant transformed cell model was established as described previously. Briefly, exponentially growing BEAS-2B cells were plated onto transwell membrane (Corning, USA) with 1 × 105 in a single well (18). CS was produced using an automatic smoking machine, and the CS was pumped into an inhalation exposure chamber. Cells were directly exposed to CS for 10 min every day at a smoke concentration of 20%. This procedure was continued until the cells reached 30 passages and named S30 cells (18).

Scratch Wound Healing Assay

The normal BEAS-2B cells and S30 cells (2 × 105) were seeded into 6 well plates and cultured at 37 °C. Cells were allowed to grow up to 100% confluence and a scratch was made in the plate using a P10 pipette. The cells were cultured in fresh serum-free DMEM medium. Images were collected at 0 and 24 h under an inverted microscope (Olympus, Germany) and quantitatively analyzed with NIH ImageJ software.

Transwell Migration Assay

The normal BEAS-2B cells and S30 cells (2 × 105) were seeded in the upper chambers (pore size, 8 μm) of the 6 well plate (Corning, USA) in 1 ml serum-free medium. The lower chambers were filled with 2 ml complete medium with 10% FBS, and the plate was incubated under standard conditions for 24 h. At the end of incubation, after removing the cells in the upper surface of the membrane with a cotton swab, cells in the lower chamber were fixed with methanol and stained with 0.5% crystal violet solution. The images were taken with an inverted microscope (Olympus, Germany) and analyzed using NIH ImageJ software.

Colony Formation Assay

1 × 103 normal BEAS-2B cells and S30 cells were plated in 0.35% agarose on top of a 0.7% agarose base supplemented with complete medium. The medium was renewed every 2–3 days. The colonies were stained with 0.5% crystal violet (Sigma, USA) for 20 min at room temperature. The colony formation rate was calculated by the following equation: colony formation rate = the number of colonies/number of seeded cells × 100%.

Tumorigenicity Assay

Five-week-old male BALB/c-nude mice of SPF grade were purchased from Beijing Vital River Laboratory Animal Technology Company Limited (Beijing, China). All nude mice were housed in the Laboratory Animal Center Soochow University. The animal experiment protocol was approved by the Laboratory Animal Ethics Committee of Experiment Animal Center of the Soochow University (Suzhou, China). Approximately 5 × 106 normal BEAS-2B cells or S30 cells were injected subcutaneously into the right flank of male BALB/c-nude mice (5 mice were used for BEAS-2B cells injection and 10 mice for S30 cells injection). Animals were euthanized 45 days after injection, and tumors were collected and photographed.

RNA Extraction and Sequencing

Total cellular RNA was extracted from S30 and normal BEAS-2B cells using TRIzol reagent (Invitrogen, USA) according to the manufacturer's protocol. Small RNA sequencing was performed on the Illumina Hiseq 2500 platform (Illumina, San Diego, CA). NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB, USA.) was used to prepare the small RNA sequencing library. To determine the known and novel miRNAs, unique clustered reads were aligned with the reference genome and database obtained from miRBase 20.0 (http://www.mirbase.org/). The miRDeep2 algorithm was used to predict novel miRNA precursors. The expression levels were estimated by transcript per million (TPM) and mRNA sequencing was performed on the Illumina HiSeq 4000 platform. The Illumina TruSeq RNA kit was used for preparing the mRNA sequencing library. The mRNAs with expression profiles that differed between the samples were normalized as fragments per kilobase of transcript per million mapped reads (FKPM). The DEGSeq package was used to analyze the differential expressed miRNAs (DEMs) or mRNAs (DEGs). P-value < 0.05 and |log2 (foldchange)| ≥ 1 were regarded as the threshold of significantly differential expression.

Data Source and Processing

The NSCLC (LUAD and LUSC) miRNA-Seq datasets and related clinicopathology information were obtained from the Xena (https://xena.ucsc.edu). The LUAD miRNA expression data included a total of 493 samples consisting of 448 LUAD samples and 45 normal adjacent samples. The LUSC miRNA expression data included a total of 380 samples comprising 336 LUSC samples and 44 normal adjacent samples. The limma package was used to identify the differential expressed miRNAs in LUAD and LUSC when compared with corresponding normal adjacent samples. The differentially expressed miRNAs were defined by a threshold of p-value < 0.05 and |log2 (foldchange)| ≥ 1.

Real-Time Quantitative PCR

A total of 1.5 μg RNA isolated from each sample was reversely transcribed into complementary DNA (cDNA) using Revert Aid First Strand Complementary DNA Synthesis Kit (for mRNA detecting, Thermo, USA) or Mir-X™ miRNA First-Strand Synthesis Kit (for miRNA detecting, Clontech Laboratories, USA) according to the manufacturer's instructions. Quantitative PCR (qPCR) was performed using NovoScript® SYBR Two-Step qRT-PCR Kit (novoprotein, China) with a QuantStudioTM 6 Flex qRT-PCR system (USA). The internal control for the quantitive analysis of miRNA and mRNA were U6 and GAPDH, respectively. The primer used for qPCR were listed in Table 1.

Table 1

Primers used in this study.

Symbol		Sequence
miR-106b-5p	Forward	5′-TAAAGTGCTGACAGTGCAGAT-′3
miR-589-5p	Forward	5′-TGAGAACCACGTCTGCTCTGAG-′3
miR-96-5p	Forward	5′-TTTGGCACTAGCACATTTTTGCT-′3
miR-181a-5p	Forward	5′-AACATTCAACGCTGTCGGTGAGT-′3
miR-361-3p	Forward	5′-TCCCCCAGGTGTGATTCTGATTT-′3
IGFBP3	Forward	5′-AGAGCACAGATACCCAGAACT-′3
	Reverse	5′-GGTGATTCAGTGTGTCTTCCATT-′3
KRT17	Forward	5′-GCCGCATCCTCAACGAGAT-′3
	Reverse	5′-CGCGGTTCAGTTCCTCTGTC-′3
FAM129A	Forward	5′-GCCTGGAAGGAACGATCCG-′3
	Reverse	5′-GGCCACCATCGCTTTGATCTT-′3
FLNC	Forward	5′-GCTCGTGTCCATAGACAGCAA-′3
	Reverse	5′-CTGGGGCACCTTGTTCTGG-′3
TIE1	Forward	5′-ACGACCATGACGGCGAATG-′3

Reverse primers of miRNAs and U6 primers are provided by Mir-X™ miRNA First-Strand Synthesis Kit.

Analysis of Gene Expression and Smoking History

To validate the correlation between expression of miRNAs and patients' smoking history, all valid LUAD samples in the TCGA database were divided into four groups according to smoking history, including (1) lifelong non-smoker (n = 66); (2) current smokers (n = 104); (3) Current reformed smoker for >15 years (n = 116); (4) current reformed smoker for ≤15 years (n = 144). The expression of miRNAs in lung adenocarcinoma tissues of each group was compared.

miRNA-mRNA Integrative Network

For identification of the candidate miRNA-mRNA network in smoking-induced malignant transformed cells, two separate steps were carried out. First, the target mRNAs of interest miRNAs were predicted through the miRDB database (http://mirdb.org/miRDB/). Second, the intersection of differently expressed genes and target genes was taken to screen the potential target genes of miRNAs in smoking-induced malignant transformed cells. These different expression target genes and miRNAs were used to construct the miRNA-mRNA regulation network through the Cytoscape software (V3.7.1, https://cytoscape.org).

Enrichment Analyses

Metascape (http://metascape.org/gp/index.html) is an effective and efficient tool for experimental biologists to comprehensively analyze and interpret OMICs-based studies in the big data era (19). The database was used to perform the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, which is used to predict the potential biological functions of the overlapping genes of the DEGs and target genes.

Statistical Analysis

SPSS 22.0 was used for statistical analysis. Values were presented as mean ± standard deviation (SD). Difference analysis between two groups was performed by using student t-test. A p < 0.05 was considered statistically significant. Correlation between differentially expressed miRNAs and predicted target mRNAs were calculated by Pearson's correlation. A p < 0.05 was regarded as statistically significant.

Results

CS-Induced Malignant Transformation in BEAS-2B Cells

To validate the malignant change of S30 cells, the normal BEAS-2B cells and S30 cells were seeded in soft agar. As shown in Figure 1A, cells formed significantly more and bigger colonies in S30 cells compared to the normal BEAS-2B cells. Besides, colony formation rate in S30 cells was remarkably higher than in the normal BEAS-2B cells (Figure 1B). Moreover, the normal BEAS-2B cells and S30 cells were used to generate xenograft tumors in nude mice. Among the ten mice injected with S30 cells, 3 developed tumor tissue (Figures 1C,D). While no tumor was found in the five mice injected with normal BEAS-2B cells.

Figure 1

CS-induced malignant transformation in BEAS-2B cells in vitro and in vivo. (A) Representative photographs of the colony formation assay of the normal BEAS-2B cells and S30 cells. (B) Graph of soft agar colony forming rate of normal BEAS-2B cells and S30 cells. **p < 0.01 vs. BEAS-2B. (C) Photographs of tumors excised 45 days after injection of normal BEAS-2B cells and S30 cells into nude mice. (D) Representative HE staining histopathologic image of tumor tissues excised 45 days after injection of S30 cells into nude mice.

CS Promoted the Migratory Ability of BEAS-2B Cells

To investigate the effect of CS in cell migration, we performed scratch wound healing and transwell cell migration assays. Scratch wound healing assay indicated that the migratory ability was significantly increased in S30 cells compared to the normal BEAS-2B cells (Figures 2A,B). As shown in Figures 2C,D, further transwell migratory assay demonstrated that the migrated cells were significantly increased in S30 cells compared to the normal BEAS-2B cells. These results suggested that long-term exposure to CS could promote the migratory ability of BEAS-2B cells.

Figure 2

CS promoted the migratory ability of BEAS-2B cells. (A) Representative images of scratch wound healing assay of normal BEAS-2B cells and CS-induced malignant transformed cells (S30). (B) Graph of wound closure rate in scratch wound healing assay of normal BEAS-2B cells and S30 cells. **p < 0.01 vs. BEAS-2B. (C) Representative images of transwell assay of normal BEAS-2B cells and S30 cells. (D) Graph of migrated cells in transwell assay of normal BEAS-2B cells and S30 cells. **p < 0.01 vs. BEAS-2B.

Differentially Expressed miRNAs Between S30 Cells and Normal BEAS-2B Cells

To test whether CS affects the miRNA-mRNA regulatory network in BEAS-2B cells, the miRNAs in normal BEAS-2B cells and S30 cells were quantitatively analyzed using small RNA sequencing. Compared with the normal BEAS-2B cells, the S30 cells showed dysregulation of 42 miRNAs that had significantly different expression levels with 2-fold change as a cut off (Figures 3A,B, Supplementary Table 1). Of these 42 miRNAs, 28 were upregulated (67%), and 14 were downregulated (33%) in the S30 cells (Figure 3C). The top five most significantly aberrant expression miRNAs are marked in the scatter plot; miR-106b-5p, miR-589-5p, and miR-96-5p were upregulated, and miR-181a-5p and miR-361-3p were downregulated (Figure 3B). The qPCR results of the top five miRNAs showed the increased miR-106b-5p, miR-589-5p, and miR-96-5p and decreased miR-181a-5p and miR-361-3p expression in S30 cells compared to normal BEAS-2B cells (Figure 4).

Figure 3

Graphical representation of the 42 miRNAs differentially expressed between S30 cells and normal BEAS-2B cells. (A) Heatmap of the 42 differentially expressed miRNAs (DEMs) between the CS-induced malignant transformed cells (S30) and normal BEAS-2B cells. The colors in the heatmap represent the normalized expression values, with low expression values being colored in shades of green and high expression values in shades of red. (B) Volcano plots were generated to visualize the distribution of DEMs between normal BEAS-2B and S30 cells. The top five most significantly dysregulated miRNAs are marked. (C) Counts of miRNAs upregulated or downregulated in S30 cells.

Figure 4

qPCR analysis of the top five most significantly dysregulated miRNAs in CS-exposed cells. Data are presented as relative fold induction compared with normal BEAS-2B cells. *p < 0.05, **p < 0.01 vs. BEAS-2B.

Differentially Expressed mRNAs Between S30 Cells and Normal BEAS-2B Cells

Next, we investigated the expression patterns of mRNAs using transcriptome resequencing. Compared with the normal BEAS-2B cells, the S30 cells showed dysregulation of 753 mRNAs that had significantly different expression levels with 2-fold change as a cut off (Figures 5A,B). Of these 753 mRNAs, 273 were upregulated (36%), and 480 were downregulated (64%) in the S30 cells (Figure 5C). The top five most significantly dysregulated mRNAs are marked in the scatter plot; IGFBP3 and KRT17 were upregulated, and FAM129A, FLNC, and TIE1 were downregulated (Figure 5B). The qPCR results of the top five mRNAs validated the increased IGFBP3 and KRT17 and decreased FAM129A, FLNC, and TIE1 expression in S30 cells compared to normal BEAS-2B cells (Figure 6).

Figure 5

Graphical representation of the 753 mRNAs differentially expressed between S30 cells and normal BEAS-2B cells. (A) Heatmap of the 753 differentially expressed genes (DEGs) between the CS-induced malignant transformed cells (S30) and normal BEAS-2B cells. The colors in the heatmap represent the normalized expression values, with low expression values being colored in green and high expression values in red. (B) Volcano plots were generated to visualize the distribution of DEGs between normal BEAS-2B and S30 cells. The top five most significantly dysregulated genes are marked. (C) Counts of genes upregulated or downregulated in S30 cells.

Figure 6

qPCR analysis of the top five most significantly dysregulated genes in CS-exposed cells. Data are presented as relative fold induction compared with normal BEAS-2B cells. *p < 0.05, **p < 0.01 vs. BEAS-2B.

Integrated Analysis of the DEMs in S30 Cells and Lung Cancer Samples

To explore whether the DEMs' expression is altered in lung cancer tissues, we analyzed the miRNAs sequencing data of lung cancer, including lung adenocarcinomas (LUAD) and squamous cell lung carcinomas (LUSC), in the TCGA database. A total of 8 miRNAs were found dysregulation in S30 cells, LUAD and LUSC samples with a similar tendency of change. Among these 8 miRNAs, 5 were upregulated (miR-96-5p, miR-93-5p, miR-589-5p, miR-4661-5p, and miR-106b-5p) and 3 were downregulated (miR-190a-5p, miR-195-5p, and miR-1-3p) in the three datasets (Figure 7).

Figure 7

Identification of differential expressed miRNAs in S30 cells and lung cancer samples. DEM_UP/DEM_DOWN: up-regulated/ down-regulated miRNAs in S30 cells; LUAD_UP/LUAD_DOWN: up-regulated/ down-regulated miRNAs in LUAD samples; LUSC_UP/LUSC_DOWN: up-regulated/ down-regulated miRNAs in LUSC samples.

Association of miRNA Expression With Smoking History

Among the five up-regulated miRNAs, three miRNAs, including miR-96-5p, miR-93-5p, and miR-106-5p, showed a higher expression in current smoking LUAD patients when compared with the lifelong non-smokers (Table 2). Three of the screened down-regulated miRNAs, including miR-190a-5p, miR-195-5p, and miR-1-3p, showed lower expression in current smoking LUAD patients when compared with the lifelong non-smokers (Table 2). Moreover, miR-96-5p and miR-106b-5p are overexpressed in the current reformed smoker for >15 years, while miR-190a-5p has lower expression in the current reformed smoker for >15 years when compared with the lifelong non-smoker (Table 2).

Table 2

The expression of miRNAs in the LUAD patients with different smoking history.

miRNAs	Significant	1 (n = 66)	2 (n = 104)	3 (n = 116)	4 (n = 144)
miR-96-5p	UP	4.25 ± 0.74	4.58 ± 1.07*	4.39 ± 1.12	4.55 ± 1.04*
miR-93-5p	UP	11.70 ± 0.87	12.09 ± 0.87**	11.65 ± 0.86	11.79 ± 0.97
miR-589-5p	UP	6.66 ± 0.65	6.72 ± 0.79	6.26 ± 0.71**	6.49 ± 0.77
miR-4661-5p	UP	2.03 ± 0.71	2.14 ± 1.05	1.94 ± 0.84	2.10 ± 1.02
miR-106b-5p	UP	7.72 ± 0.65	8.21 ± 0.72**	7.75 ± 0.68	8.03 ± 0.72**
miR-190a-5p	DOWN	1.87 ± 0.77	1.52 ± 0.74**	1.74 ± 0.65	1.57 ± 0.56**
miR-195-5p	DOWN	5.15 ± 0.98	4.73 ± 0.84**	5.11 ± 0.82	5.02 ± 0.87
miR-1-3p	DOWN	3.37 ± 1.41	2.44 ± 1.25**	3.41 ± 1.34	3.04 ± 1.28

Lifelong non-smoker (<100 cigarettes smoked in Lifetime) = 1; Current smoker (includes daily smokers and non-daily smokers or occasional smokers) = 2; Current reformed smoker for > 15 years (> 15 years) = 3; Current reformed smoker for ≤ 15 years (≤15 years) = 4 (*p < 0.05 vs. non-smoker; **p < 0.01 vs. non-smoker).

Integrative Analysis of Correlation of miRNA and mRNA in S30 Cells

To understand the potential functions of the smoking-related differentially expressed miRNAs, and to explore miRNA-mRNA interaction in S30 cells, we predicted the target genes of miRNAs and performed an intersection analysis with the gene expression data to identify genes that were inversely co-expressed with miRNAs. A total of 2,477 target genes of low-expressed miRNAs and 2,295 target genes of high-expressed miRNAs were screened by searching miRDB database. Consequently, 25 up-regulated genes and 70 down-regulated genes were found to have at least one negatively regulated miRNA-mRNA pair for smoking-related differentially expressed miRNAs (Figure 8A, Supplementary Table 2). The miRNAs-DEGs network was generated by Cytoscape software, as showed in Figure 8B.

Figure 8

Integrative analysis of miRNA-mRNA regulatory network in S30 cells. (A) Venn Diagram depicting the distribution of negatively correlated miRNA-mRNA pairs within the four datasets. (B) Regulatory network for 95 negatively correlated miRNA-mRNA, including 25 DEG_UP/DEM_DOWN_Targets pairs and 70 DEG_DOWN/DEM_UP_Targets pairs. DEG_UP/DEG_DOWN: up-regulated/down-regulated genes in S30 cells; DEM_UP/DEM_DOWN: up-regulated/ down-regulated miRNAs in S30 cells.

Enrichment Analysis of Correlation of miRNA and mRNA in S30 Cells

To further explore the function of the negatively correlated miRNA-mRNA pairs, 95 up-regulated or down-regulated target genes in S30 cells were selected for mapping into the Metascape database and subjected to functional enrichment analysis. As shown in Figure 9A, GO analysis demonstrated that these target genes are associated with several cancer-related, especially tumor migration related GO terms, including “positive regulation of cell motility,” “regulation of cell adhesion,” “mononuclear cell migration,” “cell junction organization,” “extracellular structure organization” and so on. Among these enriched DEGs, several DEGs, including THBS1, FN1, PIK3R1, CSF1, CORO2B, and PREX1, were involved in many biologic processes which derived from enrichment analysis of negative miRNA-mRNA correlations (Figure 9B). Moreover, the KEGG pathway enrichment analysis suggested that these target genes were significantly correlated with “TNF signaling pathway,” “Small cell lung cancer,” “Rap1 signaling pathway,” “PI3K-Akt signaling pathway,” “mTOR signaling pathway,” “FoxO signaling pathway,” “Focal adhesion,” “ECM-receptor interaction,” and some other cancer-related pathways (Figure 10A). In particular, “Focal adhesion” and “ECM-receptor interaction” are closely related to cell migration. In addition, THBS1, FN1, PIK3R1, and IRS1, were involved in many KEGG pathways which derived from enrichment analysis of negative miRNA-mRNA correlations (Figure 10B).

Figure 9

Gene Ontology (GO) enrichment analysis for negatively correlated miRNA-mRNA. (A) The top enriched GO terms are shown in the bubble chart. The Y-axis represents the biologic process GO terms, and the X-axis represents the rich factor (rich factor = amount of differentially expressed genes enriched in the GO term/amount of all genes in background gene set). The color and size of the bubble represent enrichment significance and the number of genes enriched in the GO term, respectively. (B) Network diagram of top enriched GO terms for the negatively correlated miRNA-mRNA. DEG_UP: up-regulated differentially expressed genes; DEG_Down: down-regulated differential expression genes; GO_BP: Gene Ontology of biological process.

Figure 10

KEGG pathway enrichment analysis for negatively correlated miRNA-mRNA. (A) The top enriched KEGG pathways are shown in the bubble chart. The Y-axis represents the KEGG pathways, and the X-axis represents the rich factor (rich factor = amount of differential expressed genes enriched in the pathway/amount of all genes in the background gene set). The color and size of the bubble represent enrichment significance and the number of genes enriched in the pathway, respectively. (B) Network diagram of top enriched KEGG pathways for the negatively correlated miRNA-mRNA. DEG_UP: up-regulated differential expressed genes; DEG_Down: down-regulated differential expression genes.

Discussion

There are nearly 1.3 billion cigarette smokers in the world, which leads to 5 million cancer related deaths every year (20). Cigarette-smoking is a notable risk factor for multiple pathologies (21–23), among them lung cancer takes the lead with smokers having a much higher risk than non-smokers. Our previous studies have suggested that chronic exposure to CS could induce malignant transformation of the human bronchial epithelial cell line (BEAS-2B) and tumorigenesis (18, 24). In recent years, studies have indicated that miRNAs play an essential role in tumor initiation, development, and metastasis as well as the cellular response to stress by modulating the expression of their target genes (25–27). In this study, we investigate the effect of chronic exposure to CS on the expression of miRNA and mRNA in BEAS-2B cells and further examined the interaction of miRNA and mRNAs.

Based on our high throughput sequencing data and the TCGA database analysis, we found significant dysregulation of 6 smoking-related miRNAs in S30 cells compared with normal BEAS-2B cells. Among these miRNAs, miR-190a is found downregulated in aggressive neuroblastoma (NBL). Overexpression of miR-190a contributed to the inhibition of tumor growth and prolonged the dormancy period of fast-growing tumors (28). A recent study showed that miR-190a could inhibit the metastasis of breast tumor by involving in estrogen receptor (ERα) signaling (29). miR-195 usually serves as a tumor suppressor in several cancer types, such as gastric cancer (30), renal cancer (31), cervical cancer (32), liver cancer (33), and osteosarcoma (34), and its downregulation was related to lymph node metastasis and advanced clinical stage (32). Similarly, miR-1 can regulate multiple behavior of the tumor cells, such as proliferation (35, 36), migration (37), apoptosis (38, 39), and metabolism (40). In addition, miR-106b and miR-93 are both the members of miR-106b~25 cluster, which have regarded as significant oncogenic drivers as well as potential biomarkers and therapeutic targets in various tumors (41–44). Moreover, Several studies have demonstrated that miR-96 could act as an oncogene (45–47) or tumor suppressor (48, 49) depending on the different types of cancer. Although these miRNAs have extensively been reported to be associated with many other kinds of cancer, their roles in lung cancer has yet been demonstrated.

Numerous studies have established the regulatory relationships between miRNA and mRNA expression (50, 51). CS-induced DEMs can bind to 3′UTR regions of several genes and down-regulate their expression, indicating that these miRNAs may contribute to the pathogenesis of smoking-related diseases. It has been reported that negatively regulated miRNA-mRNA pairs are significantly contributed to the initialization and development of different kinds of cancer (52–54). In order to further comprehend the roles of the miRNA-mRNA pairs in CS-induced lung cancer, we selected 95 dysregulated target mRNAs of the 6 CS-related miRNAs and found that they are involved in several cancer-related signaling pathways including TNF signaling pathway, Rap1 signaling pathway, PI3K-Akt signaling pathway, mTOR signaling pathway, FoxO signaling pathway, ECM-receptor interaction, and so on. Meanwhile, the GO enrichment analysis results indicated that these target genes were participated in a series of cell adhesion and migration biological processes, suggesting these miRNA-mRNA pairs related to the process of epithelial-mesenchymal transformation (EMT). EMT is considered to be an important regulator of metastasis by promoting the invasion and spread of tumor cells to distant organs (55). Among these enriched DEGs, IRS1, PIK3R1, THBS1, and FN1 are related to more than 4 KEGG pathways. As an adaptor of the insulin growth factor-1 receptor, IRS1 plays an essential role in cell growth and proliferation, primarily via the Akt pathway, and it was reported to be regulated by several miRNAs through direct or indirect action (56–59). Moreover, studies have demonstrated that PIK3R1 was directly targeted by miR-127 (60), miR-21 (61), miR-155 (62), and some other miRNAs in different kinds of cancers. It's reported that the activity of phosphoinositide 3- kinase (PI3K) is activated by many oncogenes and the PI3K family members are involved in a serious of biological processes and the genesis and progression of various tumors (63). Thrombospondin 1 (THBS1) is a secreted protein with multiple biological functions (64), including a potent anti-angiogenic activity and activation of latent transforming growth factor beta (TGF-β) (65, 66). A recent study suggested that the expression of THBS1 was modulated by multiple miRNAs (67). Moreover, it's reported that fibronectin 1 (FN1) is crucial to many cellular processes, including cell proliferation, adhesion, migration and differentiation (68, 69), and the expression of FN1 was regulated by miR-1271 (70), miR-9 (71), and miR-206 (72). Similar to previous studies, we identified the negatively regulated miRNA-mRNA pairs in the CS-induced lung cancer, which were implicated in several cancer-related (especially EMT-related) biological process and KEGG pathways in the malignant transformation progress of lung cells induced by CS. Further study will be needed to explore the targeting relationships of these miRNAs and their target mRNAs and their possible roles on cancer-related molecular mechanisms for the development of novel targeted therapy for CS-induced lung cancer.

In conclusion, our study demonstrated that the expression profiles of miRNA and mRNA were significantly dysregulated in BEAS-2B cells with long-term exposure to CS. Smoking induced miRNAs are associated with EMT and carcinogenesis.

Data Availability Statement

The datasets generated for this study can be found in the Sequence Read Archive (SRA) database (https://trace.ncbi.nlm.nih.gov/Traces/sra/) with identifier SRP182926 and SRP181756. The LUAD and LUSC datasets analyzed for this study can be obtained from UCSC Xena (https://xenabrowser.net/datapages/).

Ethics Statement

The animal study was reviewed and approved by The Laboratory Animal Ethics Committee of Experiment Animal Center of the Soochow University. Written informed consent was obtained from the owners for the participation of their animals in this study.

Author Contributions

JL conceived and designed the study. JW performed and analyzed the experiments. XY, NO, SZ, HY, and XG assisted to collect and analyze the data. JW wrote the manuscript. JT and TC were of immense help in the modification of the manuscript. All authors read and approved the final manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.