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 posttranscriptional 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.

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 × 10 5 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 × 10 5 ) 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 serumfree 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 × 10 5 ) 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 × 10 3 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 × 10 6 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 R Multiplex Small RNA Library Prep Set for Illumina R (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

Symbol
Sequence Reverse primers of miRNAs and U6 primers are provided by Mir-X TM miRNA First-Strand Synthesis Kit.
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 TM miRNA First-Strand Synthesis Kit (for miRNA detecting, Clontech Laboratories, USA) according to the manufacturer's instructions. Quantitative PCR (qPCR) was performed using NovoScript R 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.

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.  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.

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.

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  (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).

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).

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).

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 upregulated 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.

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 cancerrelated, 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).

DISCUSSION
There are nearly 1.3 billion cigarette smokers in the world, which leads to 5 million cancer related deaths every year (20). Cigarettesmoking 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)(26)(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 fastgrowing 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)(42)(43)(44). Moreover, Several studies have demonstrated that miR-96 could act as an oncogene (45)(46)(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)(53)(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)(57)(58)(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 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. 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.

FUNDING
This study was supported by the National Natural Science Foundation of China (81573178 and 81172707) and the Suzhou Science and Technology Project (SS201832). The study was also supported by Jiangsu Key Laboratory of Preventive and Translational Medicine for Geriatric Diseases as well as the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD).