Construction and Analysis of a circRNA-Mediated ceRNA Network in Lung Adenocarcinoma

Introduction

Lung cancer is the most common malignancy and the leading cause of cancer-related death worldwide.¹ Lung adenocarcinoma (LUAD) is the most common form of lung cancer, accounting for approximately 40% of lung malignancies.^2,3 Difficulties with early diagnosis and the lack of effective treatments are associated with the poor prognosis of this malignancy.⁴ It is therefore important to identify molecular markers associated with patient survival, which may improve early diagnosis and contribute to the development of gene-targeted therapies.

Circular RNAs (circRNAs) are a new class of noncoding RNA molecules that show a prevalent, stable, conserved, and specific expression pattern.⁵ CircRNAs are abundant in many eukaryotic cells and human tissues, and certain circRNAs are conserved among different species.⁶ Moreover, circRNAs can be detected in tissue samples, saliva, or plasma and show a cell-specific or stage-specific expression pattern.⁷ CircRNAs act as competitive inhibitors by binding to miRNAs (referred to as “miRNA sponges”), or they can act as target mimics to block the activity of specific microRNAs (miRNAs).^8,9 CircRNAs provide binding sites for the miRNA either in the non-coding transcript or in the 3'-untranslated region of a specific gene; circRNAs are miRNA sponges and can be used as a new class of biomarker.¹⁰ CircRNAs enter the circulation and can be detected in the serum.¹¹ Because circRNAs downregulate the expression of miRNAs, they have been investigated in many human cancers, such as colorectal cancer, lung, and stomach cancers.¹² Wan et al showed that the expression level of the circRNA circ-ITCH is lower in lung cancer tissues than in adjacent normal tissues, and in the lung cancer cells lines A549 and NCI-H460, circ-ITCH inhibits miR-214 and miR-7 by acting as a sponge.¹³ CircRNA_100876 (circ-CER) is significantly upregulated in non-small cell lung cancer (NSCLC) tissues compared with paired adjacent non-tumorous tissues.

In this study, we analyzed the circRNA-miRNA-mRNA-competing endogenous RNA (ceRNA) regulatory network using GEO datasets. Differentially expressed circRNAs (DECs), miRNAs (DEMs), and mRNAs (DEGs) were identified and compared between LUAD samples and controls. The interaction between DEMs and predicted miRNAs, and the association between the expression of DEGs and LUAD patient survival were analyzed using The Cancer Genome Atlas (TCGA) samples. The present bioinformatics analysis may help identify novel biomarkers for LUAD.

Materials and Methods

Source of Microarray Data

The public GEO database is a functional genomic database found at http://www.ncbi.nlm.nih.gov/geo. Four datasets (GSE101684, GSE101586, GSE135918, and GSE130779) were downloaded. The circRNA dataset GSE101684 included four LUAD tissues and four adjacent normal lung tissues, and the samples were tested by the GPL21825 074301 Arraystar Human CircRNA microarray V2 platform. The circRNA dataset GSE101586 included four LUAD tissues and four lung tissues, and the samples were tested by the GPL19978 Agilent-069978 Arraystar Human CircRNA microarray V1 platform. The miRNA dataset GSE135918 included five LUAD tissues and five adjacent normal lung tissues, and the samples were tested by the GPL18058 Exiqon miRCURY LNA microRNA array, 7th generation [miRBase v18, condensed Probe_ID version] platform. The mRNA dataset GSE130779 included eight LUAD tissues and eight adjacent non-tumor tissues; the platform used for the microarray was GPL20115 Agilent-067406 Human CBC lncRNA + mRNA microarray V4.0 (Probe name version).

Identification of DECs, DEMs, and DEGs

Differential analyses of GSE101684, GSE101586, GSE135918, and GSE130779 were performed to identify DECs, DEMs, and DEGs between LUAD tissues and control normal tissues using the online tool GEO2R.¹⁴ The criterion for DECs was P <0.05 and |logFC| ≥2. The criterion for DEMs was P <0.05 and logFC ≤-1, and the criterion for DEGs was P <0.05 and logFC ≥1. A volcano plot of the DECs was constructed using the R software gplots package.

Construction of the ceRNA Network

Circinteractome is a web tool for exploring circRNAs and their interacting proteins and miRNAs.¹⁵ Circinteractome was used to predict the potential miRNAs (https://circinteractome.nia.nih.gov/). The predicted potential miRNAs intersecting with DEMs from GEO in the online website Venny (http://bioinformatics.psb.ugent.be/webtools/Venn/) were regarded as potential miRNAs of circRNAs in LUAD. StarBase 3.0 (http://starbase.sysu.edu.cn/index.php), an open-source platform for studying RNA interactomes,¹⁶ was used to predict mRNAs. The intersecting genes that commonly appear in predicted mRNAs and DEGs were regarded as potential mRNAs in LUAD. Finally, a circRNA-miRNA-mRNA network was constructed, and the data visualized using Cytoscape software 3.7.1.

Survival Analysis

The association between the mRNAs selected by the ceRNA network and survival was analyzed using UALCAN (http://ualcan.path.uab.edu/index.html) and GEPIA (http://gepia.cancer-pku.cn/). First, we performed the relative level of mRNAs in LUAD. The association between DEGs and LUAD patient survival and the expression of DEGs based on TCGA samples were analyzed using UALCAN. Pathological stage analysis and the prognostic value of mRNAs patients in the disease overall survival curve in LUAD were performed for further verification using GEPIA. P < 0.05 was considered statistically significant.

Tissue Specimens

A total of 5 pairs of LUAD tissues and adjacent normal tissues were collected from LUAD patients who had surgery in the Seventh Affiliated Hospital, Sun Yat-sen University. None of these patients received preoperative chemotherapy, radiotherapy or targeted therapy. The tissue specimens frozen and stored at liquid nitrogen until further analysis.

RNA Extraction and qRT-PCR

RNA from tissues was extracted by using TRIzol reagent (Abbion) under the manufacturer’s protocol. RNA was quantified using a NanoDrop-ND1000 (thermo) system, qPCR used Biorad CFX96 touch system. SYBR master mix (hamQ) and cDNA Synthesis Kit (TAKARA) were used in this study. The specific primers were shown as below: hsa_circ_0072088 (Forward, 5'- GGATTCACCCACAGCGG-3'. Reverse, 5'-CTACACCCGAAGCTGGCT-3'). 2−ΔΔCt was used for quantitative analysis.

Statistics

Data are displayed as mean ± SD, t-test was used to assess the differences. P < 0.05 was considered as statistically significant.

Results

Identification of DECs

The flowchart of the validation and construction of the ceRNA network is shown below (Figure 1). To identify DECs in LUAD, GSE101684 of four LUAD patients and GSE101586 of five LUAD patients were selected (Table 1). Using |logFC| ≥2 and P < 0.05 as the cut-off criteria, 34 DECs were identified in paired lung adenocarcinoma tissues and adjacent non-tumor tissues (5 in GSE101586 and 29 in GSE101684; Tables S1 and S2). A volcano plot of the DECs was plotted using the R software gplots package (Figure 2A and B). The intersection of the two datasets was determined using the online website Venny, and hsa_circRNA_103809 (hsa_circ_0072088 using the circBase database transformation) was upregulated in the results.

Table 1 Basic Information of the Four Microarray Datasets from GEO

Figure 1 Flow chart of the approach utilized in this study.

Figure 2 Identification of potential miRNAs and mRNAs that bind to circRNAs in LUAD. (A) The volcano plot of DECs in LUAD from GSE101586 dataset. (B) The volcano plot of DECs in LUAD from GSE101684 dataset. (C) The intersection analysis of predicted miRNAs and downregulated DEMs. (D) The intersection analysis of predicted mRNAs and upregulated DEGs.

Identification of Potential DEMs

To identify DEMs in LUAD, GSE135918 of five LUAD patients was selected (Table 1). Using logFC≤-1 and P < 0.05 as the cut-off criteria, 431 downregulated DEMs were identified in lung cancer tissues and non-tumor lung tissues (Table S3). MiRNAs that might bind to hsa_circ_0072088 were predicted by Circinteractome; 24 miRNAs were identified (Table S4). The intersection of DEMs and predicted miRNAs was determined, hsa-miR-532-3p and hsa-miR-942 were downregulated (Figure 2C).

Identification of mRNAs of DEMs

MiRNAs negatively regulate the expression of mRNAs. To identify mRNAs regulated by hsa-miR-532-3p and hsa-miR-942, mRNAs were predicted using the StarBase database. The StarBase database can only predict the target genes of hsa-mir-942-5p, and 6674 mRNAs (3299 mRNAs in hsa-miR-532-3p and 3375 mRNAs in hsa-miR-942-5p) were predicted (Table S5). After removal of duplicate values, 5343 mRNAs were selected for further analysis. To identify DEGs in LUAD, GSE130779 of eight LUAD patients was selected (Table 1). Using logFC ≥1 and P < 0.05 as the cut-off criteria, 107 upregulated DEGs were identified in lung cancer tissues and non-tumor lung tissues (Table S6). The intersection of DEGs was used to predict mRNAs, and 20 mRNAs were upregulated (Figure 2D).

Construction of the ceRNA Network

The findings were used to build a circRNA-miRNA-mRNA network using Cytoscape (Figure 3).

Figure 3 CeRNA networks of circRNA-miRNA-mRNA interactions in LUAD. Green indicate circRNA, red indicates miRNAs, and yellow indicates mRNAs.

MRNA Analysis and Establishment of a Sub-Network

To analyze the mRNAs selected by the ceRNA network, survival analysis was performed using UALCAN and GEPIA. The analysis identified 19 mRNAs (TMEM52, IL24, POF1B, KIF1A, NHS, HIST2H2BE, ABCC3, PYCR1, CD79A, IGF2BP3 ANKRD17, GTSE1, MKI67, CLSPN, PLAU, LUC7L, MAGIX, GPATCH4, ABAT) that were upregulated in LUAD based on TCGA samples. The expression of LBH was lower in tumor tissues than in normal tissues (Figure 4). We also compared the relative level of mRNAs in LUAD tissues, and we found that the relative level of PYCR1 was the highest (Figure 5). The expression of IGF2BP3 (P = 0.001), MKI67 (P = 0.0046), CD79A (P = 0.012) and ABAT (P = 0.035) had prognostic value in LUAD patients using UALCAN (Figure 6). We then assessed the correlation between the pathological stage of LUAD patients and the four mRNAs (IGF2BP3, MKI67, CD79A and ABAT), and we found a significant correlation between pathological stage and the expression of IGF2BP3 (P=0.00901), MKI67 (P=0.000305), CD79A (P=2.5e-5) and ABAT (P=0.000807) (Figure 7). The prognostic value of the four mRNAs in LUAD patients was evaluated by using GEPIA for further verification, we found that LUAD patients with low levels of IGF2BP3 (P = 0.0017) and MKI67 (P =0.027) were associated with longer overall survival, but high levels of CD79A (P = 0.024) and ABAT (P = 0.00074) were associated with longer overall survival (Figure 8). Finally, a potential circRNA-miRNA-mRNA sub-network was successively established in LUAD (Figure 9).

Figure 4 The expression of 20 DEGs in LUAD using UALCAN. “p<0.05” represents significant difference.

Figure 5 The relative level of 20 mRNAs in LUAD (GEPIA).

Figure 6 The expression level of the four DEGs in LUAD patient survival using UALCAN. “p<0.05” represents significant difference.

Figure 7 Correlation between the four mRNAs and the pathological stage of LUAD patients. (A) IGF2BP3, (B) MKI67, (C) CD79A, and (D) ABAT, p<0.05.

Figure 8 The prognostic value of 4 mRNAs in LUAD patients in the disease overall survival curve (GEPIA). (A) IGF2BP3, (B) MKI67, (C) CD79A, and (D) ABAT, p<0.05.

Figure 9 CeRNA sub-network of circRNA-miRNA-mRNA interactions in LUAD. Green indicate circRNA, red indicates miRNAs, and yellow indicates mRNAs.

qRT-PCR Analysis

Hsa_circ_0072088 was found upregulated between LUAD tissues with paired adjacent normal tissues by qRT-PCR (Figure 10).

Figure 10 Expression of hsa_circ_0072088 in tissues. *P<0.05.

Discussion

CircRNA is highly stable because of its covalent closed-loop structure. Many studies have suggested that circRNAs may be involved in tumorigenesis. The main way for circRNAs regulating tumors progression is to act as sponges for miRNAs and then further regulating the expression of downstream mRNAs. Although studies have reported the construction of circRNA-miRNA-mRNA regulatory networks in LUAD,^17,18 the current understanding of circRNA-related ceRNA networks in LUAD is limited, and further study is needed. In this study, we identified a potential circRNA-miRNA-mRNA ceRNA regulatory network in LUAD.

Differential expression analysis led to the selection of one circRNA for further investigation. Hsa_circ_0072088, which is associated with cancer cell proliferation and progression, regulates the expression of miR-532-3P and FOXO4 in colorectal cancer.¹⁹ Liu et al demonstrated that hsa_circ_0072088 contributes to lung cancer cell proliferation and invasion by sponging miR-4302 and promoting ZNF121-dependent MYC expression.²⁰

CircRNAs act as miRNA sponges and regulate the expression of miRNA target transcripts.^9,21 For instance, the proliferation and invasion of HCC are dependent on circMAST1 binding to miR-1299 and inhibiting CTNND1.²² Circ_100565 upregulates HMGA2 expression by sponging miR-506-3p to promote migration, invasion, and proliferation.²³ Hsa_circ_0043265 inhibits NSCLC progression by sponging miR-25-3p to promote FOXP2 expression.²⁴ Here, we used Circinteractome and the intersection of DEMs to identify hsa-miR-532-3p and hsa-miR-942-5p.

MiRNAs regulate gene function by negatively downregulating gene expression.^25,26 Therefore, we explored the molecular mechanisms underlying the roles of the hsa_circ_0072088-hsa-miR-532-3p and hsa_circ_0072088-hsa-miR-942-5p axis in LUAD. Firstly, target genes of hsa-miR-532-3p and hsa-miR-942-5p were predicted using the StarBase database. Then, the DEGs between LUAD tissues and normal tissues were screened. Twenty upregulated mRNAs (TMEM52, IL24, POF1B, KIF1A, NHS, LBH, HIST2H2BE, ABCC3, PYCR1, CD79A, IGF2BP3 ANKRD17, GTSE1, MKI67, CLSPN, PLAU, LUC7L, MAGIX, GPATCH4, ABAT) were selected according to the intersection of the DEG set and predicted gene set.

It is important to develop genomic detection methods for LUAD to predict the survival rate of patients. In this study, survival analysis was performed using UALCAN. The results identified 19 mRNAs (TMEM52, IL24, POF1B, KIF1A, NHS, HIST2H2BE, ABCC3, PYCR1, CD79A, IGF2BP3, ANKRD17, GTSE1, MKI67, CLSPN, PLAU, LUC7L, MAGIX, GPATCH4, ABAT) showing differential expression in LUAD. The differential expression levels of IGF2BP3, MKI67, CD79A, and ABAT indicated that they may have prognostic value for predicting survival in LUAD patients. Moreover, we found that IGF2BP3, MKI67, CD79A and ABAT increased with the tumors progressed. LUAD patients with low levels of IGF2BP3 and MKI67 were associated with longer overall survival, but high levels of CD79A and ABAT were associated with longer overall survival. A potential circRNA-miRNA-mRNA sub-network was established. According to the ceRNA mechanism, the relationship between circRNAs and mRNAs should be positive.²⁷ Hsa_circ_0072088 expression was positively correlated with the expression of IGF2BP3, MKI67, CD79A and ABAT in LUAD. In this study, the hsa_circ_0072088 was identified to be overexpressed through experimental verification in the tissue of LUAD. In fact, IGF2BP3 and mKI67 have been experimentally confirmed, but CD79A and ABAT still need further experimental verification. Guo et al demonstrated that IGF2BP3 is abnormally highly expressed in LUAD tissue, and can lead to worse overall survival.²⁸ Li et al found that MKI67 was significantly upregulated in human LUAD tissues was higher than that in the control tissues.²⁹

The present study had several limitations. First, the findings need to be confirmed in additional cohorts and further experimental verification is needed. Second, the sample size was too small, as there were few samples in each microarray dataset. Third, the DAVID (https://davidd.ncifcrf.gov/) was used to the function enrichment analysis with gene count, but there were not clustered.

Conclusions

The present study revealed a potential hsa_circ_0072088-mediated circRNA-miRNA-mRNA ceRNA network in LUAD, and intersection analysis, differential expression analysis, and survival analysis were performed. But the sample size is too small, hsa_circ_0072088, the miRNAs (hsa-miR-532-3p and hsa-miR-942-5p), and the mRNAs (IGF2BP3, MKI67, CD79A, and ABAT) need further broader investigation to explore its potential to serve as prognostic biomarkers in LUAD.