Identification of circRNA-miRNA-mRNA Networks for Exploring the Fundamental Mechanism in Lung Adenocarcinoma

Introduction

Pulmonary malignancy is the most commonly diagnosed cancer and ranks first among all cancers according to their incidence and mortality.¹ Typically, lung cancer is classified into two major histologic types, non-small cell lung cancer (NSCLC) and small cell lung cancer (SCLC), each has several subtypes. As a common form of NSCLC, lung adenocarcinoma (LUAC) accounts for approximately 50% and is increasing year by year, especially in young adults, women and never-smokers.^2,3 Some serum tumor markers, such as squamous cell carcinoma antigen (SCCA) and carcinoembryonic antigen (CEA), and imaging modalities, like computed tomography (CT), are available for early-stage LUAC diagnosis. However, traditional tumor markers have relatively low specificity and sensitivity; meanwhile, CT might cause radiation damage on human body. Therefore, effective biomarkers are urgently needed for diagnosing early-stage LUAC. In recent years, with the development of molecular biology and the rise of translational medicine, novel markers and technologies are emerging for diagnosis and treatment of LUAC.^4–6

Circular RNAs (circRNAs) are categorized as a type of endogenous non-coding RNAs (ncRNAs) and have been regarded as transcription errors during the splicing process. Recently, accumulating evidences have revealed that circRNAs, which formed a covalent closed continuous loop without a polyadenylated tail and 5ʹ to 3ʹ polarity, are of vital significance in various cellular processions, like gene expression regulation, cell cycle progression, chromatin modification, etc.^7,8 Furthermore, circRNAs exhibit tissue-specific expression and incredible diversity, and they have stable structure and their abundance are associated with developmental age or stage. Moreover, many studies stated that aberrant expressions of circRNAs would play significant roles in the tumorigenesis and cancer progression.⁹ For instance, as shown in investigations, circPRKC1 was upregulated in lung adenocarcinoma and it could promote cell proliferation and migration by sponging miR545/589.¹⁰ Additionally, circRNA_102231, circFADS2, etc., have been proved to be candidate biomarkers and potential regulators for LUAC initiation and progression.^11,12 Researches into circular RNAs are becoming a hotspot in the field of RNA for their various functions and specific properties. Nowadays, microarray and RNA-sequencing (RNA-seq) were widely used for measuring transcripts, discovering aberrantly expressed gene and identifying new transcriptional units; meanwhile, bioinformatics methods are usually used for processing these data. The Gene Expression Omnibus (GEO) database provides microarray and RNA-seq data for researchers around the world. Thus, GSE101586,¹³ GSE101684¹⁴ and GSE104854¹⁵ were utilized for discovering lung cancer-related circRNAs.

In this study, public datasets were employed for exploring the aberrant circRNAs and their potential molecular function. Also, we aimed to predict the potential mechanisms of these circRNAs by constructing a ceRNA network. Furthermore, several experiments, such as agarose gel electrophoresis and qRT-PCR, were performed. Above all, we provided powerful evidence for the preliminary study of lung cancer-related circRNAs. Naturally, further experiments were supposed to be performed in the coming exploration.

Methods

Specimen of LUAD Tissues

Forty LAC and corresponding adjacent normal tissues were obtained from 20 patients in the Seventh Affiliated Hospital, Sun Yat-sen University. The tissues were collected frozen in liquid nitrogen immediately after surgical resection and stored at −80°C until RNA extraction. None of these patients received radiotherapy, chemotherapy or target therapy before the surgery. All the patients had been informed; and written informed consents were obtained from them. The Ethics Committee of the Seventh Affiliated Hospital, Sun Yat-sen University approved this study.

Cell Culture

Two LUAD cell lines, A549 and H1299, and a normal cell line, 16HBE, were used in this study (Chinese Academic of Science). The 16HBE cells were cultured in DMEM, and A549 and H1299 cells were cultured in RPMI-1640 medium (Gibco, Invitrogen, Carlsbad, CA) supplemented with 10% fetal bovine serum (FBS) and 1% antibiotics. All the cell lines were maintained at 37°C in a humidified atmosphere with 5% CO₂.

CircRNA Expression Profile Analysis

A flowchart of this study is shown in Figure 1. In considering the fewer samples explored in the single microarray, GEO database was used for searching any suitable datasets. Thus, GSE101586, GSE101684, and GSE104854 were selected (Details of these three datasets are shown in Table 1). Among them, GSE101586 and GSE101684 were processed in microarray, while GSE104854 was processed by high-throughput sequencing. We downloaded raw data and re-organized them as raw counts expression matrix. Then, limma R package was applied to normalize the raw counts (R code was provided in Table S2), plus, the fold-change and Student’s t testing were used for identifying differentially expressed circRNAs (DECs) in the two sets, respectively. The criterion for DECs was |fold-change|≧2 and P-value < 0.05. Furthermore, graphical heatmaps and RCircos plots were performed based on the circRNAs information.

Table 1 Information of Expression Profiling Microarray Datasets

Figure 1 A flow chart of this study.

Construction of circRNA-miRNA-mRNA Interaction Network

Cancer-Specific CircRNAs Database (CSCD) is an online tool for studying the structure model of cancer-specific circRNAs. In order to investigate the miRNAs sponged to DECs, CSCD was applied for searching the possible miRNA binding sites. Subsequently, the miRNA–mRNA interactions were obtained from MiRTarbase (http://miRTarBase.mbc.nctu.edu.tw/),¹⁶ TargetScan (http://www.targetscan.org/),¹⁷ and StarBase 3.0 (http://starbase.sysu.edu.cn/).¹⁸ Target mRNAs were intersected from the above.

Furthermore, the data of miRNAs and mRNAs, including their expression profile and clinical patient characteristics, were used for filtering the miRNAs and mRNAs which were significantly associated with patients’ survival. Thus, the key miRNAs and mRNAs were found by integrative analysis and the CircRNA-miRNA-mRNA network could be visualized by Cytoscape software.

Functional and Pathway Analysis

The Database for Annotation, Visualization, and Integrated Discovery (DAVID)¹⁹ (http://david.abcc.ncifcrf.gov/) provides systematic and integrative annotation tools for predicting biological function of genes. We used the DAVID database to investigate the potential function of the key genes from the previous step to identify the underlying function of circRNAs. Significant results were obtained with P value < 0.05.

In addition, we analyzed the key genes related pathways according to the pathway analysis of KOBAS,²⁰ an online database for functional gene set enrichment and functional annotation, and corrected P < 0.05 was considered statistically significant.

RNA Extraction and qRT-PCR

All the cell lines (16HBE, A549, and H1299) and 40 matched tissues were used to perform qRT-PCR. Total RNA of them was extracted by TRIzol Reagent (Invitrogen) and reversely transcribed with a reverse transcriptase kit from Takara Bio, Inc. In order to verify the circular form of these DECs, gel electrophoresis was performed and the PCR products were confirmed by Sanger sequencing. The glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as an internal control.

The relative expression levels of DECs were calculated using 2^−ΔΔt, and GraphPad Prism 7.0 was used to visualize the results. The conditions and the primers used for the qRT-PCR reactions are described in Tables S1.

Results

Identification of DECs

Three microarray datasets, GSE101516, GSE101684, and GSE104854 were included in this study. Using the limma R package, we studied the DECs in these 3 datasets, respectively (logFC > 2 and p < 0.05). Furtherly, we summarized the general characteristics of the DECs obtained from the previous step. First, we took intersection of these 3 datasets (Figure 2), and the heatmaps are shown in Figure 3. The circRNA hsa_circ_0072088 and hsa_circ_0082564 were found upregulated, while the hsa_circ_0008274, hsa_circ_0000519, and hsa_circ_0003528 were identified downregulated in the results. Also, as can be seen from Figure 3, hsa_circ_0072088, hsa_circ_0082564, hsa_circ_0008274, hsa_circ_0000519, and hsa_circ_0003528 were located in different chromosomes (the basic information of these 5 circRNAs is shown in Table 2).

Table 2 Basic Information of Differentially Expressed circRNAs

Figure 2 Identification of overlapping DECs in GSE101586, GSE101684 and GSE104854.

Figure 3 Heatmaps of differentially expressed circRNAs in 3 datasets (A–C) and the chromosomal location of 5 common circRNAs (D).

Finally, we used CSCD, a comprehensive cancer-specific circRNA database and plays a vital role in the research for the function and regulation of cancer-associated circRNAs, to study the structure patterns of DECs (Figure 4).

Figure 4 Structure patterns of hsa_circ_0072088 (A), hsa_circ_0082564 (B), hsa_circ_0008274 (C), hsa_circ_0000519 (D) and hsa_circ_0003528 (E) based on the CSCD database.

Identification of DECs by qRT-PCR

Considering the signature value of DECs, qRT-PCR was performed for validation of DECs expression. In this step, the cell lines (16HBE, A549, and H1299) and 20 matched tissues were used to perform qRT-PCR. In order to ensure the specialty and back-spliced junction of DECs, agarose gel electrophoresis and Sanger sequencing were performed (Figure S1). Then, qRT-PCR showed that hsa_circ_0082564 and hsa_circ_0000519 cannot be detected for their relatively low expression level in tumor tissue, and hsa_circ_0003528 presented the downregulated trend, with a P value > 0.05 (Figure S2). Finally, hsa_circ_0072088 and hsa_circ_0008274 were found differentially expressed in normal tissues compared to cancerous tissues with a P value < 0.05. Also, the circRNA expression validation was performed using cells with the same result (Figure 5).

Figure 5 Expression of hsa_circ_0072088 in tissues and cell lines (A, C). Expression of hsa_circ_008274 in tissues and cell lines (B, D). **P<0.01.*P<0.05.

Construction of a ceRNA Network

“miRNA sponge” has been largely acknowledged as a potential mechanism of circRNA working. Thus, a ceRNA network was constructed as per the DECs. CSCD database was used for circRNA-miRNA prediction, and the regulator relationships of miRNA-mRNA were obtained from MiRTarbase, TargetScan, and StarBase 3.0. The circRNA-miRNA-mRNA network was then constructed by Cytoscape (Figure 6). Interestingly, we found that miR-578 was predicted to be common downstream of hsa_circ_0072088 and hsa_circ_0008274, which indicated that there might be some functional connection between them.

Figure 6 Construction of circRNA-miRNA-mRNA network of hsa_circ_0072088 (A) and hsa_circ_008274 (B).

Functional Analysis of mRNAs

As a novel non-coding RNA, circRNA was still not well understood. Although the functional analysis was not directed on DECs, we considered that the mRNAs may play significant roles in ceRNA network and the analysis on mRNAs could reveal the potential mechanisms of DECs. Thus, DAVID and KOBAS online database were applied for further exploration. By submitting the mRNAs obtained above, GO and KEGG analysis were performed.

As can be seen from Figure 7, the GO analysis of hsa_circ_0072088 was mainly enriched in regulation of cellular response to growth factor stimulus, cellular response to transforming growth factor beta stimulus, response to transforming growth factor beta, etc., which revealed that the circRNA might play roles in tumor growth and apoptosis. Fewer cancer-related pathways were obtained from the KEGG analysis, and MAPK signaling pathway was worthy of special attention.

Figure 7 GO and KEGG analysis of hsa_circ_0072088 (A, C) and hsa_circ_008274 (B, D).

As for hsa_circ_0008274, we noticed that the GO analysis was mainly enriched in skeletal system development, mesenchyme development, muscle cell differentiation, etc., and it seemed like there was a faint relationship between them and lung cancer. Then, we noticed that epithelial to mesenchymal transition (EMT) was obtained, which had been proved to play important roles in cancer invasion and metastasis. The KEGG analysis showed promising results that some famous pathways, like mTOR signaling pathway, MAPK signaling pathway, and PI3K−Akt signaling pathway were obtained, which provides more possibility on the assumption that hsa_circ_0008274 may play significant roles in cancer initiation and progression.

Key circRNA-miRNA-mRNA Axis Analysis

In order to further understand the potential function of circRNAs and to provide valuable information for experiments, we reorganized the circRNA-miRNA-mRNA axis according to the results of functional analysis.

From the previous steps, multiple GO terms and KEGG pathways were found, and mRNAs were clustered in any different functions. Next, we reorganized the mRNAs which were significantly associated with cancer progression and linked their upstream miRNA and circRNA (Table 3). We thought these axes could provide more information and more appropriate for specific experiments.

Table 3 Cancer-Related circRNA/miRNA/mRNA Axis

Discussion

CircRNAs, which lacked 5ʹ-caps and 3ʹ-polyadenylated tails, was first put forward by Sanger and had been ignored in traditional transcriptome studies for decades.²¹ Particularly, with next-generation sequencing and computational biology methods developed, a great deal of circRNAs was identified and found to play significant roles in multiple tissues and cells.²² Nowadays, exploring the underlying mechanisms and functions of circRNAs has attracted more and more attention. Among the studies, many researchers have indicated that circRNAs play significant roles in different malignant tumors. For instance, hsa_circ_0013958 was identified to be overexpressed in the plasma and tissue of non-small cell lung cancer, and its expression level was meaningfully related with the tumor-node-metastasis (TNM) and lymphatic metastasis.²³ Also, hsa_circ_0058246 was found to be upregulated in gastric cancer (GC) and associated with tumor recurrence.²⁴ In addition, increasing studies have shown that circRNAs may serve as important regulators in various cancers.^25–27 However, the novel circRNAs and the function they played were still needed to explore.

In this study, the public database, GEO, was used to collect the circRNA expression. GSE101586, GSE101684, and GSE104854 were chosen for analyzing in lung adenocarcinoma. Considering for removing any other factor, we downloaded the raw data of these three datasets, and reorganized them into matrix file. Limma R package was used to explore the differentially expressed circRNAs, respectively. Also, in order to focus on the most significant circRNAs, we took the intersection of 3 datasets. Three heatmaps and a Circos plot were drawn to show relative expression and the location of DECs. In the result, hsa_circ_0072088 and hsa_circ_0082564 were found upregulated, while the hsa_circ_0008274, hsa_circ_0000519, and hsa_circ_0003528 were found downregulated in lung adenocarcinoma. Furthermore, qRT-PCR and agarose gel electrophoresis were performed, and the results were verified by Sanger sequencing. Hsa_circ_00720881 and hsa_circ_0008274 were then identified indeed different expressed in tissues and cells. As conserved endogenous RNAs, circRNAs could sponge the miRNAs and play significant roles in regulating miRNA target genes.^28,29 In order to explore the aforementioned two circRNAs that sponged miRNAs, an online tool, CSCD, was used, and it could predict the MREs within 50bp downstream and upstream of junction point of circRNA.³⁰ Then, miRNA-targeted mRNAs were predicted by using MiRTarbase, TargetScan, and StarBase 3.0; thus the ceRNA network was constructed.

Furthermore, the genes were input following the collection of then, and then GO and KEGG analysis were performed. As can be seen, the mRNAs, which were predicted the downstream of hsa_circ_0008274, were enriched in various functions. We highlighted the functions related to lung cancer, like “regulation of cellular response to growth factor stimulus”, “cellular response to TGF-β stimulus”, “response to transforming growth factor beta”, and “MAPK signaling pathway”. Also, we noticed the “EMT”, which was proved to be associated with cancer invasion and metastasis, were obtained in the hsa_circ_0008274. The pathways, such as mTOR signaling pathway, MAPK signaling pathway, and PI3K−Akt signaling pathway, give more enlightenment on researches into tumor growth and progression. Thus, we further reorganized the ceRNA network and selected the most interested circRNA/miRNA/mRNA axes, like “EMT”, “Cellular response to TGF-β stimulus”, “MAPK signaling pathway”, and “PI3K-AKT signaling pathway”, with details noted in Table 3.

Nowadays, biological big data was emerging, and bioinformatics was gaining more and more attention. Some public databases, like GEO and the cancer genome atlas (TCGA), were widely used for data analysis. Since the TCGA has no public circRNA data, we just chose the GEO database to study and validate with our own tissues and cell lines. As is known, the bioinformatics analysis provided guidance and suggestions for follow-up experiments, saved money for researches, and improved efficiency. Thus, a series of analysis were performed in this study. Besides, Sun et al and Guan et al also performed bioinformatical analysis for circRNAs in hepatocellular carcinoma and gastric cancer, which used different methods and studied the ceRNA in different angles.^31,32 Thus, it would be reasonable to believe that bioinformatics provide guidance for our study. Certainly, further rigorous experiments are necessary to validate these findings. Given that the results are based on computational biology, further in-depth studies on circRNAs are indispensable to verify their possible roles in lung adenocarcinoma.