A putative competing endogenous RNA network in cisplatin-resistant lung adenocarcinoma cells identifying potentially rewarding research targets.

Lung adenocarcinoma (LUAD) is the most common type of non-small cell lung cancer and has a poor 5 year survival rate (<10%). Cisplatin is one of the most effective chemotherapeutic treatments for LUAD, even though it is of limited overall utility due to acquired drug resistance. To identify possible genetic targets for the mitigation of cisplatin resistance, gene expression data from cisplatin-resistant cell lines were integrated with patient information. Expression data for cisplatin-resistant and cisplatin-sensitive A549 cell lines were obtained from the Gene Expression Omnibus database, while LUAD patient data was obtained from The Cancer Genome Atlas (TCGA) database. Differentially expressed mRNAs (DEmRNAs), microRNAs (DEmiRNAs) and long non-coding RNAs (DElncRNAs) were identified between the cisplatin-sensitive and cisplatin-resistant cells. Using the TCGA patient data, 33 DEmRNAs associated with survival were identified. A total of 74 DElncRNAs co-expressed with the survival-associated DEmRNAs, and 11 DEmiRNAs that regulated the survival-associated DEmRNAs, were also identified. A competing endogenous RNA (ceRNA) network was constructed based on the aforementioned results, which included 17 survival-associated DEmRNAs, 9 DEmiRNAs and 16 DElncRNAs. This network revealed 8 ceRNA pathway axes possibly associated with cisplatin resistance in A549 cells. Specifically, the network suggested that the lncRNAs HOXD-AS2, LINC01123 and FIRRE may act as ceRNAs to increase cisplatin resistance in human LUAD cells. Therefore, it was speculated that these lncRNAs represent potentially rewarding research targets. Copyright: © Li et al.

Introduction

Lung cancer is a common malignancy and is the leading cause of cancer-associated mortality worldwide (1,2). Non-small cell lung cancer (NSCLC) accounts for ~85% of all lung cancers (3). Lung adenocarcinoma (LUAD) is the most common pathologic subtype of NSCLC in non-smoking males, and in all females (both smokers and non-smokers) (4,5). Although numerous resources have been directed towards the development of novel LUAD treatments, the prognosis of patients with advanced LUAD remains unsatisfactory, with a 5 year survival rate <10% in 2018 (6). LUAD is relatively sensitive to primary chemotherapy, but tumors rapidly acquire chemoresistance, leading to death for most patients (7,8).

Cisplatin is one of the most effective chemotherapeutic drugs and is used to treat various tumors, including testicular cancer, ovarian cancer, cervix carcinoma, breast cancer, prostate carcinoma, bladder cancer, lung cancer, melanoma and head-and-neck cancer (9,10). Cisplatin has a broad-spectrum anticancer activity, but its use is limited due to it causing severe side effects and due to a number of tumors acquiring cisplatin resistance (9). Although the side effects caused by cisplatin have been mildly alleviated by newly-developed antagonists (11), cisplatin resistance, which commonly originates from multiple cellular self-defense adaptations, often results in disease recurrence (12). Thus, the development of cisplatin resistance remains a substantial challenge for chemotherapeutics.

A major impediment to a comprehensive understanding of the molecular mechanisms underlying cisplatin-induced drug resistance is that most currently available results were generated using isolated cell lines. These studies can be misleading when extended to in vivo experiments and clinical trials (12,13). However, the integration of cell line data with clinical information, especially overall survival (OS) time, may improve this issue. For example, Zhao et al (14) used The Cancer Genome Atlas (TCGA) database to demonstrate that patients expressing high levels of the long non-coding RNA (lncRNA) HOMEOBOX A11 antisense RNA (HOXA11-AS) have shorter survival rates compared to the low expression level group; mechanistic experiments subsequently showed that the microRNA (miRNA/miR) targeted by HOXA11-AS affects cisplatin resistance in LUAD cells. The aforementioned study thus provides a framework for the identification of additional miRNAs associated with cisplatin resistance in LUAD cells.

In the present study, the framework of Zhao et al (14) was used to identify miRNA targets that may be useful for the mitigation of cisplatin resistance. The present study aimed to: i) Identify differentially expressed (DE) mRNAs (DEmRNAs), DEmiRNAs and DElncRNAs between two LUAD cell lines, namely A549 (cisplatin-sensitive) and A549-DDP (cisplatin-resistant), using data from the Gene Expression Omnibus (GEO) database (15); ii) quantify the expression levels of these DEmRNAs in samples of patients with LUAD using data downloaded from the TCGA database; iii) construct a competing endogenous RNA (ceRNA) network based on the aforementioned data; and iv) assess the associations between the elements of the ceRNA network and patient OS time to identify potential research targets.

Materials and methods

A549/A549-DDP data retrieval

Two miRNA and mRNA expression datasets were downloaded from the GEO database (16): GSE43249 (17), which was derived from the GPL14613 (miRNA-2) Affymetrix Multispecies miRNA-2 Array, and GSE43493 (18), which was derived from the GPL15314 Arraystar Human LncRNA microarray V2.0 (Agilent_033010 Probe Name version). Each dataset contained six samples, three that were cisplatin-sensitive and three that were cisplatin-resistant.

A549/A549-DDP data pre-processing

The raw microarray data were read using the package affy v1.52.0 (19) in R v3.4.3 (http://www.bioconductor.org/packages/release/bioc/html/affy.html), and was standardized using the robust multi-array average (20,21) method, with background adjustment, quantile normalization and summarization on a log2 scale. Using the platform annotation file, the probe was annotated and the unmatched probe was removed. To map different probes to the same mRNA or miRNA data, the mean value of each different probe was used as the final expression, and the genes were divided into mRNAs and lncRNAs following the guidelines of the HUGO Gene Nomenclature Committee (22).

Identification of DEmRNAs, DEmiRNAs and DElncRNAs

The DEmRNAs, DElncRNAs and DEmiRNAs were identified in the GEO datasets using the R package limma v3.34.9 (23). The classical Bayesian test was used to calculate P-values. mRNAs, lncRNAs and miRNAs were considered significantly differentially expressed if |log2 (fold change)|≥1 and P<0.05. To visualize the DEmRNAs, DElncRNAs and DEmiRNAs, heat maps and volcano maps were generated using the R packages ggplot2 (24) and heatmap2 (25), respectively.

TCGA patient data retrieval

RNA sequence data and clinical information (specifically, cisplatin treatment status and OS time) for 576 patients with LUAD were retrieved from the TCGA database (https://www.cancer.gov/tcga; accessed on August 29, 2017). The use of TCGA data in the present study is in accordance with TCGA publication guidelines (https://cancergenome.nih.gov/publications/publicationguidelines). Since the patient data used originated from the TCGA database, no further ethical approval was required.

Identification of DEmRNAs associated with patient survival

The expression levels of each of the identified DEmRNAs were quantified in each patient with LUAD. For each DEmRNA, patients were divided into a low- and a high-expression group based on mean gene expression. Kaplan-Meier survival curves were generated, and the DEmRNAs that were significantly associated with OS were identified using a log-rank test. P<0.05 was considered to indicate a statistically significant difference.

Functions and interactions of the survival-associated DEmRNAs

The Database for Annotation, Visualization and Integrated Discovery v.6.8 (26) was used to identify the Gene Ontology (GO) (27) terms and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (28) pathways significantly enriched in the survival-associated DEmRNAs (i.e. those with P<0.05). The STRING database v10.5 (29) was used to predict protein-protein interactions (PPIs) of the survival-associated DEmRNAs. PPI scores ≤0.15 were considered of low confidence. Cytoscape v3.6.1 (30) was used to visualize the PPI network and to calculate node degrees.

Co-expression of DElncRNAs and survival-associated DEmRNAs

The Pearson correlation coefficient between each DElncRNA and each survival-associated DEmRNA was calculated. P-values were adjusted using the false discovery rate (FDR) to control for the effects of multiple comparisons. DElncRNAs and survival-associated DEmRNAs were considered to be co-expressed when |r|>0.95 and P<0.05 (FDR-adjusted). The functions of the co-expressed DElncRNAs were predicted based on the lncRNA-mRNA regulatory network; the R package clusterProfiler (31) was used to identify the pathways significantly enriched in the target genes of the co-expressed DElncRNAs. Pathways with Benjamini-Hochberg-adjusted P-values <0.05 were considered significantly enriched.

DEmiRNA regulatory networks and KEGG pathway enrichment

The target gene prediction module of miRWalk v2.0 (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/miRretsys-self.html) (32) was used to identify possible target genes of the DEmiRNAs in eight databases miRWalk (http://mirwalk.umm.uni-heidelberg.de), Microt4 (http://mirtarbase.mbc.nctu.edu.tw/php/index.php), MiRanda (http://www.microrna.org/microrna/home.do), miRDB (http://www.mirdb.org/miRDB/policy.html), miRMap (https://mirmap.ezlab.org), PITA (http://genie.weizmann.ac.il/pubs/mir07/mir07_dyn_data.html), RNA22 (https://cm.jefferson.edu/rna22) and Targetscan (http://www.targetscan.org/vert_71/). To increase the reliability of the search results, only genes identified in ≥5 databases were used to construct the miRNA control network, which was visualized with Cytoscape v3.6.1 (30). The KEGG pathway enrichment of the predicted DEmiRNA target genes was investigated using clusterProfiler (31).

Construction of a ceRNA regulatory network

lncRNAs associated with the DEmiRNAs were identified using the prediction module of DIANA-LncBase v2 (33); only lncRNAs with scores >0.75 were included. Subsequently, a ceRNA network based on several data sources was constructed: The lncRNA-miRNA regulatory network; the miRNA-target mRNA regulatory network; and the DElncRNAs that were positively co-expressed with survival-associated DEmRNAs.

Association between OS time and the expression levels of selected lncRNA targets

Preliminary results demonstrated that the lncRNAs HOXD-AS2, LNC01123 and FIRRE appeared in one or more ceRNA axes. Therefore, the expression levels of these lncRNAs were quantified, as well as those of the co-expressed lncRNAs and survival-associated DEmRNAs, in non-LUAD tumors using the Gene Expression Profiling Interactive Analysis (GEPIA) server (34); GEPIA analyses RNA expression in 9,736 tumors and 8,587 normal samples from the TCGA and the Genotype-Tissue Expression projects.

Statistical analysis

The classical Bayesian test was used to test differentially expressed mRNAs, lncRNAs and miRNAs. DEmRNAs that were significantly associated with OS time were identified using the log-rank test. Fisher's exact test was applied for the GO enrichment of DEmRNAs associated with OS time. All comparisons were between cisplatin-resistant A549-DDP cells and cisplatin-sensitive A549 cells. P<0.05 was considered to indicate a statistically significant difference, unless otherwise specified. The statistical analysis was performed with R v.3.4.3 (35).

Results

DEmRNAs, DEmiRNAs and DEIncRNAs in the A549 and A549-DDP cell lines

A total of 842 mRNAs were identified to be differentially expressed between the A549 and A549-DDP cell lines. Among these DEmRNA, 245 (29.10%) were upregulated in the A549-DDP cell line compared with the A549 cell line, while 597 (70.90%) were downregulated (Fig. 1). In addition, 90 DElncRNAs and 18 DEmiRNAs were identified. Among these DEmiRNAs and DElncRNAs, 37 DElncRNAs (41.11%) and 8 DEmiRNAs (44.44%) were upregulated in the A549-DDP cell line compared with the A549 cell line, while 53 lncRNAs (58.89%) and 10 miRNAs (55.56%) were downregulated (Fig. 1; Table SI).

Figure 1.

Heat maps and volcano plots showing differentially expressed long non-coding RNAs, microRNAs and mRNAs. (A and B) DELncRNAs. (C and D) DEmiRNAs. (E and F) DEmRNAs. Blue columns, cisplatin-sensitive group; red columns, cisplatin-resistant group; orange dots, upregulated genes; green dots, downregulated genes; blue dots, unchanged genes.

Survival-associated DEmRNAs

In the TCGA patient dataset, 86 patients treated with cisplatin were identified. These patients expressed 786 of the identified DEmRNAs. Among these, 33 DEmRNAs were significantly associated with OS time (Table I). Five upregulated DEmRNAs were associated with low OS time: Rh family B glycoprotein (RHBG), phosphatase orphan 2 (PHOSPHO2), activity regulated cytoskeleton associated protein (ARC), thioredoxin (TXN) and kinesin family member 26A (KIF26A; Fig. 2A-E; Table SII). Four other upregulated DEmRNAs were associated with high OS time: Zinc finger protein 417, neural cell adhesion molecule 1 (NCAM1), mediator complex subunit 12 (MED12) and ADP ribosylation factor 4 (ARF4; Fig. 2F-G; Table SII). These nine DEmRNAs comparing with the other 24 DEmRNAs, were more related to the prognosis of patients.

Table I.

Survival-associated mRNAs differentially expressed between cisplatin-resistant and cisplatin-sensitive cell lines.

Symbol	Log2 FC	Low median	High median	Description
MT1A	−2.63		2.55	Metallothionein 1A
VGF	−2.50	4.15	2.55	VGF nerve growth factor inducible
SARM1	−2.22	2.86	4.15	Sterile alpha and TIR motif containing 1
DPP4	−2.21	2.86	4.15	Dipeptidyl-peptidase 4
SIRT4	−1.85	2.87		Sirtuin 4
SH2B2	−1.65		2.86	SH2B adaptor protein 2
PER1	−1.65		2.86	PER1
FKBP1B	−1.60	4.15	2.60	FK506 binding protein 1B
FAM117A	−1.55	2.86		Family with sequence similarity 117 member A
DIRAS3	−1.50		2.60	DIRAS family GTPase 3
STAC3	−1.47		2.60	SH3 and cysteine rich domain 3
MAGEH1	−1.42	2.60		MAGE family member H1
RAB9B	−1.40	2.87		RAB9B, member RAS oncogene family
SLC17A9	−1.28	2.86		Solute carrier family 17 member 9
ADRA1D	−1.21	4.15	2.55	Adrenoceptor alpha 1D
ELOVL2	−1.19		2.87	ELOVL fatty acid elongase 2
DBP	−1.19	2.87		D-box binding PAR bZIP transcription factor
NR1D1	−1.10	2.87		Nuclear receptor subfamily 1 group D member 1
HSPA2	−1.07	4.15	2.55	Heat shock protein family A (Hsp70) member 2
GJA1	−1.04	4.15	2.21	Gap junction protein alpha 1
CEACAM6	−1.03	4.15	2.55	Carcinoembryonic antigen related cell adhesion molecule 6
ID4	−1.01	4.15	2.60	Inhibitor of DNA binding 4, HLH protein
NDST3	1.01	2.53	2.53	N-deacetylase and N-sulfotransferase 3
ZNF417	1.06	2.60		Zinc finger protein 417
PHOSPHO2	1.09		2.86	Phosphatase, orphan 2
ARC	1.11		2.55	Activity regulated cytoskeleton associated protein
TXN	1.13	4.15	2.60	Thioredoxin
ARF4	1.15	2.86		ADP ribosylation factor 4
RHBG	1.27	4.15	2.86	Rh family B glycoprotein (gene/pseudogene)
EDEM1	1.28	2.86		ER degradation enhancing alpha-mannosidase like protein 1
KIF26A	1.48		2.55	Kinesin family member 26A
NCAM1	2.20	2.86		Neural cell adhesion molecule 1
MED12	2.46	2.53	4.15	Mediator complex subunit 12

Genes highlighted in bold are the five highly-expressed DEmRNAs associated with low OS time. The column ‘Symbol’ contains the gene name/ID; the column ‘LOG FC’ contains the log2 FC of up/downregulated genes; the column ‘Low median’ contains the median OS time of the low expression group; the column ‘High median’ contains the median OS time of the high expression group; and the column ‘Description’ contains the full gene name. OS time, overall survival; FC, fold change. If genes are without low median the overall survival time of their low expression group is out of the end of follow-up time point (e.g. five years). If genes are without high median the overall survival time of their high expression group is out of the end of follow-up time point.

Figure 2.

Kaplan-Meier survival curves for the six most upregulated mRNAs in the cisplatin-resistant A549-DDP cell line compared with the cisplatin-sensitive A549 cell line. (A) RHBG. (B) PHOSPHO2. (C) ARC. (D) TXN. (E) KIF26A. (F) ZNF417. (G) NCAM1. (H) MED12. (I) ARF4. Red lines represent the overall survival of the patients expressing low levels of each mRNA; blue lines represent the overall survival of the patients expressing high levels of each mRNA. The cut-off point for high/low levels was the mean of gene expression. The shaded areas are the 95% CI of the corresponding groups. RHBG, Rh family B glycoprotein; PHOSPHO2, phosphatase, orphan 2; ARC, activity regulated cytoskeleton-associated protein; TXN, thioredoxin; KIF26A, kinesin family member 26A; ZNF417, zinc finger protein 417; NCAM1, neural cell adhesion molecule 1; MED12, mediator complex subunit 12; ARF4, ADP ribosylation factor 4.

Functional enrichment and PPIs of the survival-associated DEmRNAs

The GO terms most over-represented in the DEmRNAs annotations were ‘extracellular region’, ‘blood microparticle’, ‘extracellular space’ and ‘linoleic acid metabolic process’ (Fig. 3A; Table SIII). No KEGG pathways enriched in the DEmRNAs were identified (data not shown). The PPI network of the survival-associated DEmRNAs (Fig. 3B) contained 19 nodes and 26 interaction pairs, including 6 upregulated and 13 downregulated DEmRNAs.

Figure 3.

Survival-associated mRNAs differentially expressed between the cisplatin-resistant and cisplatin-sensitive A549 cell lines. (A) GO enrichment. Fisher's exact test was used and P≤0.01 was considered to indicate a statistically significant difference. (B) Protein-protein interactions. Red nodes represent upregulated genes; green nodes represent downregulated genes. *P≤0.01 cisplatin-resistant vs. cisplatin-sensitive A549 cell lines. GO, Gene Ontology.

A total of 168 positively co-expressed pairs of DElncRNAs and survival-associated DEmRNAs were identified (74 DElncRNAs and 32 DEmRNAs). According to the DElncRNA-DEmRNA network, the target genes of the co-expressed DElncRNAs were over-represented in three KEGG pathways: ‘Protein processing in endoplasmic reticulum’, ‘mineral absorption’ and ‘circadian rhythm’ (Fig. 4; Table SIV).

Figure 4.

KEGG pathway enrichment of the co-expressed lncRNAs differentially expressed between the cisplatin-resistant and cisplatin-sensitive A549 cells lines. Dot colors represent the significance of the enrichment of each miRNA in each pathway; dot size represents the enriched to un-enriched gene ratio. KEGG, Kyoto Encyclopedia of Genes and Genomes; miRNA/miR, microRNA.

DEmiRNA target gene prediction and functional enrichment analysis

Using miRWalk v2.0 (32), 11 DEmiRNAs targets, 17 survival-associated DEmRNA targets and 52 DEmiRNA/survival-associated DEmRNA pairs were identified in the DEmiRNA regulatory network (Fig. 5; Table SV). Several KEGG pathways, including ‘cell adhesion molecules (CAMs)’ and ‘circadian rhythm’, were enriched in the target genes (Fig. 6; Table SVI).

Figure 5.

miRNA-mRNA regulatory network. Circles represent genes and triangles represent miRNAs. Red shapes represent upregulated molecules and green shapes represent downregulated molecules. miRNA/miR, microRNA.

Figure 6.

KEGG pathway enrichment of the miRNAs differentially expressed between the cisplatin-resistant and cisplatin-sensitive A549 cells lines. Dot colors represent the significance of the enrichment of each lncRNA in each pathway; dot size represents the enriched to un-enriched gene ratio. KEGG, Kyoto Encyclopedia of Genes and Genomes; lncRNA, long non-coding RNA.

ceRNA regulatory network

Using Cytoscape (30), the DElncRNA-DEmiRNA regulatory network was combined with the DEmiRNA-DEmRNA network to obtain a DElncRNA-DEmiRNA-DEmRNA ceRNA network (Fig. 7A; Table SVII). The ceRNA network included 9 DEmiRNAs, 16 DElncRNAs, 17 target DEmRNAs and 87 pairs with a regulatory association. The present putative ceRNA network (Fig. 7B) included eight axes: HOXD-AS2/hsa-miR-152-3p/MED12, HOXD-AS2/hsa-miR-152-3P/NCAM1, LINC01123/hsa-miR-152-3p/MED12, LINC01123/hsa-miR-152-3P/NCAM1, LINC01123/hsa-miR-152-3P/ARF4, LINC01123/hsa-miR-762/NCAM1, LINC01123/hsa-miR-762/RHBG and FIRRE/hsa-miR-1231/ARF4. HOXD-AS2/has-miR-152-3p/ARF4 was not included in the axes list as HOXD-AS2/hsa-miR-152-3p/ARF4 didn't form a triangle and there was no line to connect HOXD-AS2 and ARF4 (Fig. 7B). The putative ceRNA axes included three lncRNAs (HOXD-AS2, LINC01123 and FIRRE), three miRNAs (hsa-miR-152-3p, hsa-miR-762 and hsa-miR-1231) and four genes (MED12, RHBG, NCAM1 and ARF4).

Figure 7.

Competing endogenous RNA (ceRNA) regulatory network. (A) Whole ceRNA network. (B) Subset of the ceRNA network showing the potential miRNA-lncRNA-mRNA axes. Circles represent genes, triangles represent miRNAs and diamonds represent lncRNAs. Red shapes represent upregulated molecules and green shapes represent downregulated molecules. Solid lines represent regulatory relationships, where the gene at one end is regulated by the gene/miRNA/lncRNA at the opposite end. Dashed lines connect co-expressed molecules. Arrows indicate the miRNA-mRNA regulatory network. T-shapes show the lncRNA-miRNA regulatory network. ceRNA, competing endogenous RNA; miRNA, microRNA; lncRNA, long non-coding RNA.

Association between OS time and the expression of selected lncRNA targets

The GEPIA analysis identified RHBG and NCAM1 as co-expressed genes of HOXD-AS2 and LNC01123, respectively. LNC1123 was positively co-expressed with RHBG in LUAD (R=0.18; Fig. 8A), lung squamous cell carcinoma (LUSC; R=0.14; Fig. 8B) and testicular germ cell tumors (TGCT; R=0.63; Fig. 8C). LNC1123 was positively co-expressed with NCAM1 in mesothelioma (MESO; R=0.27; Fig. 8D), pheochromocytoma and paraganglioma (PCPG; R=0.34; Fig. 8E), and TGCT (R=0.41; Fig. 8F). In addition, LNC1123 upregulation was associated with shorter patient OS time in head and neck squamous cell carcinoma (HNSCC; Fig. 9A), and in cervical squamous cell carcinoma and endocervical adenocarcinoma (P=0.057; Fig. 9B).

Figure 8.

Co-expression associations in various types of cancer. Co-expression of LNC01123 and RHBG in (A) lung adenocarcinoma, (B) lung squamous cell carcinoma and (C) TGCT. Co-expression of LNC01123 and NCAM1 in (D) mesothelioma, (E) pheochromocytoma and paraganglioma and (F) TGCT. TPM, transcripts per million; RHBG, Rh family B glycoprotein; NCAM1, neural cell adhesion molecule 1; TGCT, testicular germ cell tumors.

Figure 9.

Patient OS time with respect to lncRNA expression in various types of cancer. Survival curve based on the expression of LNC01123 in (A) head and neck squamous cell carcinoma and (B) cervical squamous cell carcinoma and endocervical adenocarcinoma. Survival curve based on the expression of HOXD-AS2 in (C) colon adenocarcinoma, (D) brain lower-grade glioma, (E) lung squamous cell carcinoma and (F) uveal melanoma. Red lines represent the OS time of the patients expressing high levels of each lncRNA; blue lines represent the OS time of the patients expressing low levels of each lncRNA. The criteria for being considered as high/low levels was the mean of expression of lncRNA. The dotted lines represent the 95% CIs of the corresponding groups. lncRNA, long non-coding RNA; OS time, overall survival; HR, hazard ratio (P-value of HR is presented in the right upper corner of the figure).

HOXD-AS2 upregulation was associated with shorter patient OS time in colon adenocarcinoma (COAD; Fig. 9C), brain lower-grade glioma (LGG; Fig. 9D), LUSC (Fig. 9E) and uveal melanoma (Fig. 9F). HOXD-AS2 was positively co-expressed with NCAM1 in COAD (R=0.15; Fig. 10A), LUAD (R=0.17; Fig. 10B), LUSC (R=0.09; Fig. 10C), TGCT (R=0.6; Fig. 10D), uterine corpus endometrial carcinoma (R=0.25; Fig. 10E) and uterine carcinosarcoma (R=0.37; Fig. 10F). No positive associations were identified between HOXD-AS2 and RHBG (data not shown).

Figure 10.

Co-expression of the long non-coding RNAs HOXD-AS2 and NCAM1 in (A) colon adenocarcinoma, (B) lung adenocarcinoma, (C) lung squamous cell carcinoma, (D) testicular germ cell tumors, (E) uterine corpus endometrial carcinoma and (F) uterine carcinosarcoma. NCAM1, neural cell adhesion molecule 1; TPM, transcripts per million.

FIRRE was positively co-expressed with NCAM1 in glioblastoma multiforme (R=0.47; Fig. 11A), liver hepatocellular carcinoma (LIHC; R=0.12; Fig. 11B), PCPG (R=0.29; Fig. 11C) and TGCT (R=0.49; Fig. 11D). FIRRE was positively co-expressed with RHBG in prostate adenocarcinoma (R=0.24; Fig. 11E) and TGCT (R=0.44; Fig. 11F). In addition, FIRRE upregulation was associated with shorter patient OS time in kidney renal clear cell carcinoma (Fig. 12A), kidney renal papillary cell carcinoma (Fig. 12B), LGG (Fig. 12C), LIHC (Fig. 12D), MESO (Fig. 12E) and pancreatic adenocarcinoma (Fig. 12F).

Figure 11.

Co-expression relationships in various types of cancer. Co-expression of FIRRE and NCAM1 in (A) glioblastoma multiforme, (B) liver hepatocellular carcinoma, (C) pheochromocytoma and paraganglioma, and (D) TGCT. Co-expression of FIRRE and RHBG in (E) prostate adenocarcinoma and (F) TGCT. TPM, transcripts per million; NCAM1, neural cell adhesion molecule 1; RHBG, Rh family B glycoprotein; TGCT, testicular germ cell tumor.

Figure 12.

Patient OS time with respect to FIRRE expression in various cancers. (A) kidney renal clear cell carcinoma, (B) kidney renal papillary cell carcinoma, (C) brain lower grade glioma, (D) liver hepatocellular carcinoma, (E) mesothelioma and (F) pancreatic adenocarcinoma. Red lines represent the OS time of the patients with high levels of FIRRE; blue lines represent the OS time of the patients with low levels of FIRRE. The criteria for being considered as high/low levels was the mean expression of FIRRE. The dotted lines represent the 95% CIs of the corresponding groups. OS time, overall survival; HR, hazard ratio (P-value of HR is presented on the right upper corner of the figure).

Discussion

The present analysis of gene expression patterns (based on the GEO datasets) identified 33 genes differentially expressed in cisplatin-resistant A549-DDP cells compared with in cisplatin-sensitive A549 cells. Among these, nine were upregulated in the cisplatin-resistant cells and 24 were downregulated. By cross-referencing these results with patient data from the TCGA dataset, five of these upregulated genes (PHOSPHO2, ARC, TXN, RHBG and KIF26A) were identified to be associated with poor OS time outcomes. These five genes may be useful potential targets for the reversal of cisplatin resistance in LUAD.

RHBG was identified as being of particular interest, as this gene also appeared in one of the axes of the putative ceRNA network generated in the present study. RHBG is a non-erythroid membrane glycoprotein of the Rh antigen family, and the mechanisms regulating RHBG expression remain poorly studied. Consistent with the results of the present study, RHBG has been demonstrated to be expressed in LIHC and COAD cell lines (36). Additionally, RHBG has been implicated in the growth of brain tumors in mice (37). RHBG deserves further study both as a possible maker of poor LUAD outcomes and as a potential target for cisplatin-resistance reversal therapy.

Furthermore, NCAM1, MED12 and ARF4 appeared in one or more ceRNA network axes, but increased expression levels of these genes were associated with improved OS time. NCAM1 encodes a cellular adhesion protein and is a well-known potential target of antibody-based cancer immunotherapies (38). In addition, NCAM1 has been identified as an immunohistochemical marker for lung neuroendocrine tumors (39), and it was recently proposed that the NCAM1-180 splice variant might be a useful marker for NSCLC (40). Furthermore, NCAM1 may be a useful biomarker and therapeutic target for acute myeloid leukemia (41), the follicular variant of papillary thyroid carcinoma (42) and breast cancer (43). Although NCAM1 has been associated with cisplatin resistance in ovarian cancer (44,45), the in vitro expression of NCAM improved the response of multiple myeloma cells to Bortezomib (Btz) treatment (46). Consistent with this, the present study revealed that NCAM upregulation was associated with improved patient OS time outcomes.

MED12 is a component of the CDK8 subcomplex. MED12 mutations are associated with tumorigenesis (47). Indeed, somatic mutations in MED12 exon 2 have been observed in uterine leiomyosarcoma, colorectal cancer (CRC) (47), uterine leiomyoma, breast fibroadenoma, phyllodes tumors and prostate cancer (48). Additionally, inhibition of MED12 expression has been associated with resistance to cisplatin and other chemotherapy drugs (49,50). This is consistent with the results of the present study, in which patients with high levels of MED12 had improved OS time.

ARF4 is a small guanine-binding protein that serves a role in vesicular trafficking (51). Although the results of the present study suggested that ARF4 upregulation was associated with improved patient outcomes, it has been previously reported that high expression levels of ARF4 in patients with breast cancer are significantly associated with increased risk of distant metastasis and shorter OS time. Conversely, ARF4 silencing reduces the colonization of the lung by metastatic breast cancer cells in vivo (51). These contradictory results suggest that the role of ARF4 in LUAD deserves further investigation.

The present putative ceRNA network included three miRNAs (hsa-miR-152-3p, hsa-miR-762 and hsa-miR-1231) across the eight axes. In several types of cancer (including prostate, ovarian and breast), miR-152 expression has been shown to reduce tumor cell viability and proliferation (52–54). In addition, the suppression of miR-152 biogenesis increases cisplatin resistance in epithelial ovarian cancer (55). However, the overexpression of miR-152 increases cisplatin resistance and proliferation of nasopharyngeal carcinoma cells (56), while the overexpression of miR-762 stimulates the development of various tumors, including ovarian (57) and breast cancer (58). Conversely, the expression of miR-762 (in combination with other miRNAs) leads to the apoptosis of breast cancer cells (59). In the present study, miR-152 and miR-762 were downregulated in the cisplatin-resistant LUAD cells. Overall, these results suggested that the behavior of these miRNAs may vary in different types of cancer.

By contrast, miR-1231 expression consistently negatively regulates the progression of various types of cancer, including glioma (60,61), pancreatic cancer (62) and papillary thyroid cancer (63). Additionally, miR-1231 has been identified as an independent prognostic factor; low expression of miR-1231 is associated with worse patient outcomes compared with high expression of miR-1231 in glioma and pancreatic cancer (60–62). Consistent with these results, the present study revealed that miR-1231 upregulation was associated with improved patient OS time and cisplatin sensitivity.

The present putative ceRNA axes included three lncRNAs (HOXD-AS2, LINC01123 and FIRRE). Each of these three lncRNAs has been shown to be upregulated in one or more types of cancer, and each one is commonly associated with poor patient prognosis. For example, LINC01123 is upregulated in intrahepatic cholangiocarcinomas (64) and is associated with poor prognosis in prostate cancer (65). Similarly, HOXD-AS2 is upregulated in glioma cells and is associated with poor prognosis (66). Consistent with these previous studies, the present study revealed that LINC01123 and HOXD-AS2 were upregulated in numerous types of cancer and were associated with reduced patient OS time. Importantly, the HOXD-AS2/hsa-miR-152-3P/NCAM1 and LINC01123/hsa-miR-762/RHBG axes in the present putative ceRNA network were supported by the co-expression results, which showed that LINC01123 was co-expressed with RHBG and that HOXD-AS2 was co-expressed with NCAM1.

FIRRE upregulation is associated with poor OS time in diffuse large B-cell lymphoma, CRC and HNSCC (67). However, FIRRE upregulation is also associated with significantly improved OS time in CRC (68). The present study revealed that FIRRE was upregulated in numerous types of cancer, possibly indicating that this lncRNA behaves differently under different circumstances.

In combination with the aforementioned results, the analyses of the current study revealed that the mRNAs, miRNAs and lncRNAs that form the potential axes in the present putative ceRNA network serve various important roles in cancer pathogenesis and progression. Importantly, a number of these molecules may serve different roles in different types of cancer. Thus, the results of the present study suggest these molecules as important targets for future studies focused on cancer diagnosis, prognosis and therapy. NCAM1 and miR-152 are particularly intriguing targets with respect to cisplatin resistance, as NCAM1 increases Btz sensitivity and miR-152 reduces cisplatin-induced effects (46). However, further investigations are necessary to determine the ceRNA mechanisms underlying cisplatin resistance in LUAD.

In addition, the present study presents some limitations, such as that the TCGA dataset included relatively few patients that met the set criteria and that the available clinical survival data was restricted to OS time. Future studies should recruit patients with lung cancer for cisplatin chemotherapy, collect lung lesion biopsy samples from patients with disease progression after three cycles of chemotherapy and then quantify the expression levels of the candidate lncRNAs (HOXD-AS2 and LINC01123), miRNAs (hsa-miR-152-3p and hsa-miR-762) and mRNAs (NCAM1, MED12 and ARF4) identified herein in the biopsy samples. Additionally, patients should be followed-up at 3 months, 6 months, 1 year and 3 years after chemotherapy to determine survival rates.

Despite these limitations, the results of the present study suggested that the integration of cell line experimental data with clinical information may be a valuable method to identify key cancer genes and potentially useful research targets.