Bioinformatics Identification of ZNFs/LINC00520/miR-181d/BCL2 Axis as a Novel Network in Cisplatin-Resistant Lung Adenocarcinoma Cells ()
1. Introduction
Lung cancer ranks high in the incidence and mortality of malignant tumors worldwide [1] [2]. Approximately, 85% of lung cancer is Non-Small Cell Lung Cancer (NSCLC), and in this type, Lung Adenocarcinoma (LUAD) has the worst prognosis [3]. Currently, chemotherapy is still the major treatment method, but drug resistance often leads to treatment failure or tumor recurrence. The potential mechanisms of chemoresistance are as follows: 1) the overexpression of efflux transporters such as ATP-Binding cassette transporters may lead to increased drug efflux and/or decreased drug influx [4]; 2) tumor cells may have defects in the apoptotic machinery by overexpressing anti-apoptotic proteins or under-expressing pro-apoptotic proteins [5] [6]; 3) there may be dysregulation of DNA damage repair pathways [7], including base excision repair, nucleotide excision repair, single-strand break repair, and DNA mismatch repair [8]; 4) there may be activation of detoxifying systems, including glutathione and glutathione S-transferase P1 [9]; 5) there may be an epithelial-to-mesenchymal transition, which is conducive to chemoresistance [4]; 6) cancer stem cells, as a subgroup of tumor-initiating cells, can self-renew and may cause resistance to conventional chemotherapy [10]. However, the mechanisms of chemoresistance in lung cancer have not been fully elucidated and further studies are needed.
Recently, substantial evidence has suggested that autophagy, a conserved metabolic process, which promotes orderly degradation and recycling of the intracellular organelles and cytoplasmic proteins [11], plays a vital role in chemoresistance. For instance, galectin-1 could lead to chemoresistance against cisplatin (DDP) treatment by inducing autophagic flux in hepatocellular carcinoma cells, which could be reversed by an autophagy inhibitor or knockdown of Autophagy-Related Gene 5 (ATG5) [12]. Kim et al. [13] demonstrated that DDP treatment induced GDNF Family Receptor α 1 (GFRA1) expression in human osteosarcoma cells. On the other hand, the induction of GFRA1 reduced DDP-induced apoptotic cell death by regulating AMPK-dependent autophagy. Tan et al. [14] indicated that the Thioredoxin-Related Protein of 14 kDa (TRP14) induced autophagy and consequently DDP resistance in ovarian cancer cells via the AMPK/mTOR/p70S6K signaling pathway. Moreover, autophagy and DDP resistance were decreased or increased with knockdown or overexpression of TRP14, respectively, and this effect was reversed by treatment with Rapa or ATG5 knockdown. These studies strongly suggest that regulation of autophagy activity might be an efficient strategy for reversing chemoresistance.
A competitive endogenous RNA (ceRNA) hypothesis, which describes that transcripts could regulate each other at the post-transcription level by competing for shared microRNA (miRNA) response elements, was proposed in 2011 [15]. Recently, a growing body of evidence has indicated that long non-coding RNA (lncRNA)-miRNA-mRNA ceRNA regulatory networks are closely associated with the occurrence, development, invasion, and metastasis of NSCLC [16] [17] [18] [19], as well as chemoresistance. Xiao et al. [20] reported that lncRNA CASC2, as a ceRNA, bound to miR-18a and repressed its function to release target gene IRF-2 in NSCLC cells, thus strongly inhibiting the proliferation, migration, and invasion of DDP-resistant NSCLC cells (H226/DDP and A549/DDP) in vitro and tumor growth in vivo. Moreover, a handful of studies have described interactions among lncRNA, miRNA, and mRNA, which play crucial roles in autophagy-related chemoresistance in NSCLC. For example, Sun et al. [21] demonstrated that the knockdown of lncRNA-XIST restored the chemosensitivity of DDP-resistant A549 cells to DDP via regulating autophagy by the ceRNA pathway of lncRNA-XIST/miR-17/ATG7. Huang et al. [22] found that BLACAT1 was upregulated in DDP-resistant NSCLC cells, and it promoted autophagy and chemoresistance of NSCLC cells through the miR-17/ATG7 signaling pathway. Therefore, it is necessary to construct more holistic and systematic ceRNA networks that regulate chemoresistance depending on autophagy in NSCLC.
In the present study, we aim to identify Autophagy-Related Differentially Expressed mRNAs (ARDEmRNAs), differentially expressed miRNA (DEmiRNAs), and differentially expressed lncRNAs (DElncRNAs) between DDP-sensitive A549 cells and DDP-resistant A549/DDP cells by bioinformatics analysis. Furthermore, the ceRNA regulatory network related to autophagy was constructed, and functional annotation and network analysis were subsequently performed on the genes in the network. In addition, the subnetwork was extracted based on critical nodes. The relationships with prognoses were further analyzed utilizing the clinical data from The Cancer Genome Atlas (TCGA). The independent prognostic factors LINC00520, miR-181d, and BCL2 were screened out to construct the risk model. Meanwhile, Transcription Factors (TFs) located upstream of LINC00520 were predicted. This study will contribute to understanding the molecular mechanism of chemoresistance in LUAD and provide new therapeutic targets for the treatment of LUAD chemoresistance.
2. Materials and Methods
2.1. Data Acquisition
Three primitive microarray datasets including expression profiles of mRNAs, miRNAs, and lncRNAs in DDP-sensitive A549 cells and DDP-resistant A549/DDP cells, i.e. GSE108214 [23], GSE157692 [24], and GSE43493, were downloaded from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) [25]. The basic information for the three datasets was listed in Table 1.
2.2. Data Processing and Differential Expression Analysis
To obtain gene expression matrix data, the R program was applied to normalize the microarray data and annotate the probe. The limma R package was used to screen Differentially Expressed Genes (DEGs). Significantly DEmRNAs, DEmiRNAs and DElncRNAs were screened byP < 0.05 and |log2(fold-change)| > 1.
2.3. Screening of ARDEmRNAs
A total of 232 ATGs were extracted from Human Autophagy Database (HADb, http://www.autophagy.lu/index.html). Then, the intersection of 232 ATGs and significantly DEmRNAs in GSE108214 dataset was taken to obtain ARDEmRNAs using Venn software (v2.0; http://www.bioinformatics.com.cn/static/others/jvenn/index.html).
![]()
Table 1. Information on the three chips from gene expression omnibus database.
2.4. Construction of Autophagy-Related ceRNA Network
The interactions between miRNAs and 23 ARDEmRNA were predicted using two online analysis tools: miRDB (http://www.mirdb.org/miRDB) and TargetScan (http://www.targetscan.org), and the miRNAs that were predicted in two databases simultaneously were obtained. Additionally, the DElncRNAs-miRNAs interactions were predicted by miRcode (http://www.mircode.org). The intersection of these two miRNA sets and the DEmiRNAs was taken using Venn software. Based on the “ceRNA hypothesis”, a ceRNA regulation network must meet the following criteria: 1) The lncRNA was co-expressed positively with mRNA. 2) The lncRNA was co-expressed negatively with miRNA. 3) The mRNA was co-expressed negatively with miRNA. The interactive networks of DElncRNAs, DEmiRNAs, and ARDEmRNAs were thus established. The Cytoscape (Version 3.6.0) software was used to visualize the ceRNA network.
2.5. Functional Enrichment Analysis
Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the DAVID v6.8 online tool (https://david-d.ncifcrf.gov/). GO functional categories included Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). Statistically significant terms were determined using a two-sided P < 0.05.
2.6. Protein-Protein Interaction (PPI) Analysis
A PPI network of the ARDEmRNAs included in the ceRNA network was constructed using the STRING online database (https://string-db.org/), with the cut-off criterion of an interaction score > 0.4. Then, the TSV file of the PPI network was imported into Cytoscape v3.6.0. CytoHubba, a plugin in Cytoscape, was utilized to select hub genes of the PPI network. Subsequently, significant modules of the network were screened out using the Molecular Complex Detection (MCODE), another plugin in Cytoscape. P < 0.05 was considered statistically significant.
2.7. Extraction of Subnetwork
The ceRNA network topological properties were calculated using the NetworkAnalyzer Tool in Cytoscape v3.6.0. Based on degree, betweenness centrality, and radiality, the top 10% nodes were taken as critical nodes, which were used to further construct a subnetwork. Additionally, the TSV file of the whole ceRNA network was imported into Cytoscape v3.6.0. The hub genes were defined by the Maximal Clique Centrality (MCC) algorithm in CytoHubba.
2.8. Survival Analysis
Gene expression and clinical data of 586 LUAD patients were obtained from the TCGA dataset by using R program software. Patients were classified into groups with up-regulation and down-regulation based on the optimal cut-off value. Kaplan-Meier survival curves were plotted using SangerBox online tools (http://sangerbox.com/). Furthermore, univariate and multivariate Cox regression analyses were conducted using SPSS 19.0, and P < 0.05 indicated a significant difference. The significant prognosis-related genes were used to construct a prognostic risk score model for predicting Overall Survival (OS). The formula of the risk score model was as follows: risk score = expression level of Gene 1 × β1+ expression level of Gene 2 × β2+ ··· + expression level of Gene N × βN, where β is the regression coefficient calculated by multivariate Cox regression analysis [26]. According to the median risk score, patients were classified into high-risk and low-risk groups. To verify the predictive performance of the risk score in correctly classifying patients into deaths and lives, Receiver Operator Characteristic (ROC) curves were plotted. The univariate and multivariate Cox regression analyses were conducted again to assess whether the risk score is independent of the clinical features in LUAD patients.
2.9. Screening of TFs
Firstly, the promoter sequences of lncRNAs were obtained from the NCBI database (https://www.ncbi.nlm.nih.gov/gene/) and stored in FASTA format. The promoter region was defined as the 2000 bp upstream and 100 bp downstream of the transcriptional start site. Next, TFs that could bind to the promoter region were predicted by the UCSC Genome Browser (http://genome.ucsc.edu/index.html). The JASPAR2018 TFBS hg38 Track minimum score was set to 800. Finally, possible TF binding sites in the promoter region were predicted by the JASPAR database (https://jaspar.genereg.net/), and higher scores indicated higher reliability.
2.10. Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)
A549 were cultured in RPMI medium 1640 and A549/DDP cells were cultured in 10% FBS-RPMI medium 1640 and 2 µg/mL DDP (Sigma, United States) to maintain resistance. All cells were cultured in an incubator at 37˚C in an atmosphere containing 5% CO2. Total RNA was extracted from each cell line using TRIzol Reagent (Invitrogen, United States) according to the manufacturer’s instructions. The resulting RNA (2 µg) was used as a template for reverse-transcribing to first-strand cDNA synthesis of mRNA and lncRNA by using a PrimeScript™ RT reagent Kit (TaKaRa, China) in a total volume of 20 µL. RT-qPCR was performed on a MX3000P Quantitative PCR System (Agilent Technologies, Japan) using a TB Green® Premix Ex TaqTMII kit (TaKaRa, China) referring to the protocol. Glyceraldehyde-3-Phosphate Dehydrogenase (GAPDH) was used as the internal control, and mRNA and lncRNA values were normalized to that of GAPDH. The miRNA sequence-specific reverse transcription qRT-PCR for miR-181d and endogenous control U6 were performed according to Hairpin-itTM miRNAs qRT-PCR quantization kit and U6 snRNA real-time PCR normalization kit (GenePharma, Shanghai, China). The relative expression of RNA was calculated by the 2−∆∆Ct method. Each sample was repeated in triplicate. The primers were purchased from Shanghai Sangong Bioengineering Co., Ltd. A list of primers is shown in Table 2.
2.11. Statistical Analysis
The differential expression levels of genes were compared by using two independent samples Student t-test between the A549 cells and A549/DDP cells in the R program. The log-rank test was used to evaluate the difference in survival. The correlations between hub genes, seed genes, or clinical characteristics and prognoses were performed by univariate and multivariate Cox regression analyses. P < 0.05 was considered statistically significant.
3. Results
3.1. Identification of DEmRNAs, DEmiRNAs, and DElncRNAs
Three microarray datasets (GSE108214, GSE157692, and GSE43493) were
![]()
Table 2. Oligonucleotide primers used for RT-PCR validation of candidate noncoding RNAs and mRNAs in the present study.
downloaded from the GEO database. After data standardization, the gene expression levels were compared between A549 cells and A549/DDP cells by a limma R package. A total of 3179 DEmRNAs (1684 upregulated and 1495 downregulated) in GSE108214 dataset, 180 DEmiRNAs (69 upregulated and 111 downregulated) in GSE157692 dataset, and 160 DElncRNAs (75 upregulated and 85 downregulated) in GSE43493 dataset were identified using the standard of P < 0.05 and |log2(fold-change)| > 1. The volcano maps showed the overall distribution of DEmRNAs, DEmiRNAs, and DElncRNAs between A549 cells and A549/DDP cells (Figures 1(A)-(C)). The top 50 genes of DEmRNAs, DEmiRNAs and DElncRNAs were displayed in the heatmaps (Figures 1(D)-(F)).
3.2. Screening of 35 ARDEmRNAs
Since autophagy is involved in chemoresistance, the intersection between
![]()
Figure 1. Identification of DEmRNAs, DEmiRNAs, and DElncRNAs. Volcano plots of differentially expressed RNAs between A549 cells and A549/DDP cells. The red points in the plot represent up-regulated RNAs and the green points represent down-regulated RNAs with statistical significance. (A) mRNAs. (B) miRNAs. (C) lncRNAs. Heatmaps of top 50 differentially expressed RNAs between A549 cells and A549/DDP cells. The colors indicate the expression level: red indicates up-regulation, green means down-regulation and black indicates normal expression. Furthermore, each column represents a sample and each row represents a differentially expressed gene. (D) mRNAs. (E) miRNAs. (F) lncRNAs.
DEmRNAs in GSE108214 dataset and 232 ATGs from HADb was analyzed, and 35 ARDEmRNAs were identified (Figure 2).
3.3. Construction of Autophagy-Related ceRNA Network
Firstly, 174 miRNAs targeting 35 ARDEmRNAs were predicted by both miRDB and TargetScan. Meanwhile, 262 miRNAs targeting 160 DElncRNAs were predicted by miRcode. After analyzing the intersection of 174 miRNAs, 262 miRNAs, and the 180 DEmiRNAs in the GSE157692 dataset, the 16 common miRNAs were obtained, as well as 16 common miRNA-DEmRNA interactions, and 16 common miRNA-DElncRNA interactions. According to the regulatory relationships in the ceRNA hypothesis, 10 DEmRNAs (Table 3), 9 overlapped DEmiRNAs, and 14 DElncRNAs (Table 4) were included in the ceRNA network. The specific steps for constructing an autophagy-related ceRNA network were shown in the flow chart (Figure 2), and a visualization of the ceRNA network
![]()
Figure 2. Flow chart of the autophagy-related ceRNA network construction.
![]()
Table 3. Putative differentially expressed miRNAs that may target autophagy related differentially expressed mRNAs.
![]()
Table 4. Putative differentially expressed miRNAs that may be targets of differentially expressed lncRNAs.
was shown in (Figure 3(A)). Meanwhile, the expression levels of 14 DElncRNAs, 10 DEmRNAs, and 9 DEmiRNAs in A549 and A549/DDP cells were visualized with box plots using GraphPad Prism Software (Version 8.0) (Figures 3(B)-(D)).
3.4. Functional Enrichment Analysis of ceRNA
The functional enrichment of the 33 DEGs (10 ARDEmRNAs, 9 DEmiRNA, and 14 DElncRNA) in the ceRNA network was performed by DAVID and the top 10
![]()
Figure 3. Construction of autophagy-related ceRNA network. (A) Visualization of the ceRNA network containing 10 ARDEmRNAs, 9 DEmiRNAs, and 14 DElncRNAs. The circles represent ARDEmRNAs; diamonds represent DEmiRNAs; squares represent DElncRNAs; black represents upregulation; gray represents downregulation; and, edges indicate lncRNA-miRNA-mRNA interactions. (B) The expression levels of 9 DEmiRNAs in the ceRNA network. (C) The expression levels of 14 DElncRNAs in the ceRNA network. (D) The expression levels of 10 ARDEmRNAs in the ceRNA network. Gray represents A549 while black represents A549/DDP. *P < 0.05, compared to A549.
GO terms and KEGG pathways were displayed in (Figure 4). The GO terms of “Single-multicellular organism process”, “Multicellular organismal process” and “Negative regulation of biological process” were enriched in BP (Figure 4(A)). In CC (Figure 4(B)), “Cytoplasm”, “Intracellular membrane-bounded organelle” and “Intracellular organelle” terms were enriched. In MF (Figure 4(C)), the main enriched GO terms were “Identical protein binding”, “Molecular function regulator” and “Protein homodimerization activity”. KEGG analysis results revealed that the genes were predominantly enriched in “HIF-1 signaling pathway”, “Proteoglycans in cancer” and “PI3K-Akt signaling pathway” (Figure 4(D)).
3.5. PPI Analysis of 10 ARDEmRNAs
Based on the STRING database, we established the PPI network to reveal the functional interaction of 10 ARDEmRNAs. Ten nodes and 7 edgeswere selected in the PPI network, with PPI enrichment p-value < 3.79 × 10−4 (Figure 5(A)). According to the MCC algorithm, the top 5 hub genes were identified using the CytoHubba plugin of Cytoscape, which were the nodes with higher degrees in the network, including Vascular Endothelial Growth Factor Receptor (VEGFR), ribosomal protein S6 kinase B1 (RPS6KB1), FKBP prolyl isomerase 1A (FKBP1A), inositol 1,4,5-trisphosphate receptor type 1 (ITPR1), and Nicotinamide Phosphoribosyltransferase (NAMPT) (Figure 5(B)). MCODE analysis was further
![]()
Figure 4. Functional enrichment results of ceRNA network. (A) Top 10 significantly enriched GO terms in BP. (B) Top 10 significantly enriched GO terms in CC. (C) Top 10 significantly enriched GO terms in MF. (D) Top 10 significantly enriched pathways in the KEGG pathway.
![]()
Figure 5. PPI network analysis of 10 ARDEmRNAs. (A) The PPI regulation network was constructed by the STRING database. (B) The top 5 hub genes in the PPI network were identified using Cytoscape. (C) A module was extracted from the PPI network using MCODE.
performed to obtain a significant module, which consisted of 3 genes and 3 edges from the whole network (Figure 5(C)). Moreover, the top ranked hub gene VEGFR was also a seed gene in the module.
3.6. Extraction of Subnetwork
The top 10% - 20% of critical nodes ranked by topological properties play a key role in the whole ceRNA network. Once the critical nodes are disrupted, the network function is abnormal, thereby affecting tumorigenesis. Therefore, the topological properties of the ceRNA network were analyzed and 33 genes including 10 DEmRNAs, 9 DEmiRNAs, and 14 DElncRNAs were ranked according to betweenness centrality, degree, and radiality (Table 5). Based on the top 2 critical nodes of miR-125b and miR-181d, which were also the top 2 hub genes defined by the MCC algorithm in the CytoHubba (Figure 6(A)), the subnetwork was constructed (Figure 6(B)).
3.7. Survival Analysis of the Critical Genes in the ceRNA Subnetwork
Of the 33 critical genes, 4 genes including LINC00520, miR-125b, miR-181d, and BCL2 were found to be related to the prognosis of LUAD patients by Kaplan-Meier survival analysis and a log-rank test. Furthermore, upregulated LINC00520 and downregulated miR-125b, miR-181d, and BCL2 predicted poor prognosis (Figure 7(A)). Based on univariate (Figure 7(B)) and multivariate Cox regression analysis (Figure 7(C)), LINC00520, miR-181d and BCL2 were selected to construct a prognostic risk score model as follows: risk score = −0.653 × ExpressionLINC00520 + 0.475 × ExpressionmiR-181d + 0.605 × ExpressionBCL2. Then, with the median-risk score as the boundary, the LUAD patients were divided into the high-risk group (score higher than the median-risk score) and the low-risk group (score lower than the median-risk score). The risk score distribution plot, survival status scatter plot, and heatmap of risk gene expression in the two groups were presented in (Figure 7(D)). The risk curve and scatterplot indicated that samples with higher risk scores had a higher risk of
![]()
Figure 6. Construction of the ceRNA subnetwork. (A) The top 10 hub genes in the ceRNA network were identified using Cytoscape. (B) Visualization of the ceRNA subnetwork containing 3 ARDEmRNAs, 2 DEmiRNAs, and 2 DElncRNAs. The circles represent ARDEmRNAs; diamonds represent DEmiRNAs; squares represent DElncRNAs; black represents upregulation; gray represents downregulation; and, edges indicate lncRNA-miRNA-mRNA interactions.
![]()
Figure 7. Overall survival analysis of LUAD patients in TCGA cohort. (A) Kaplan-Meier survival curves of the critical genes in ceRNA subnetwork. (B) The forest plot of univariate Cox regression analysis. (C) The forest plot of multivariate Cox regression analysis. (D) Risk score analysis. Upper panel: risk score distribution of patients in the prognostic model. Middle panel: survival status scatter plots for patients in the prognostic model. Bottom panel: heatmap of LINC00520, miR-181d, and BCL2 expression in LUAD samples. (E) The 1-, 3-, and 5-year ROC curves of risk score. (F) Kaplan-Meier survival curves of the risk score model. (G) Univariate Cox regression analysis of the risk score and clinical characteristics. (H) Multivariate Cox regression analysis of the risk score and clinical characteristics.
![]()
Table 5. Topological properties of the 33 genes in the whole ceRNA network.
mortality. The heatmap revealed that LINC00520 might have a positive effect on LUAD, while BCL2 and miR-181d might have adverse effects. Furthermore, ROC curve analysis indicated that the average Area Under the Curve (AUC) at 1, 3, and 5 years were 0.63, 0.66, and 0.74, respectively (Figure 7(E)). Meanwhile, the Kaplan-Meier survival curve revealed that the high-risk group was closely associated with poor survival status (Figure 7(F)). To verify the independence of the risk score, univariate and multivariate Cox regression analyses were performed on age, gender, clinic stage, and TNM stage (Table 6), and risk score. The p-value of the risk score in univariate Cox regression analysis was less than 0.05, indicating that it was a well-predicting model (Figure 7(G)); however, multivariate Cox regression analysis revealed that the risk score was not independent prognostic factor (Figure 7(H)). Additionally, univariate Cox regression analysis demonstrated that the clinical variables, including clinic stage, T stage, N stage, and M stage, were significantly associated with patient survival, while multivariate Cox regression analysis indicated that only the M stage was an independent prognostic factor.
3.8. Prediction of TFs and Construction of TF Regulatory ceRNA Subnetwork
First, JASPAR was utilized to predict the possible TFs that could bind with the promoter of LINC00520 (Figure 8(A)). Five ZNF family members with more
![]()
Table 6. Clinical characteristics of lung adenocarcinoma patients in the Cancer Genome Atlas database.
![]()
Figure 8. Prediction of Transcription Factors (TFs) and construction of TF regulatory ceRNA subnetwork. (A) TFs for LINC00520. (B) The binding sites of 5 ZNF family members with more than 12 scores for LINC00520. (C) A ZNF family regulatory ceRNA subnetwork.
than 12 scores, including ZNF384, ZNF24, ZNF263, ZNF341, and ZNF148, were predicted as TFs of LINC00520 (Table 7), and the site motifs were shown in (Figure 8(B)). Furthermore, a ZNF family regulatory ceRNA subnetwork, ZNFs/LINC00520/miR-181d/BCL2 axis, was constructed (Figure 8(C)), which might provide new insight into the molecular mechanism of chemoresistance in LUAD.
3.9. Validation of Regulatory Axis Gene Expression
To verify the accuracy of the lncRNA-miRNA-mRNA regulatory network we constructed, the expression of these genes in cisplatin-resistant and parental A549 cells was tested in vitro. The RT-PCR results showed that both LINC00520 and BCL2 were significantly upregulated in A549/DDP cells relative to A549 cells (P < 0.05; Figure 9(A) and Figure 9(C)), while the expression of miR-181d
![]()
Table 7. Ten members of ZNF family with more than 10 score as transcription factor for LINC00520.
![]()
Figure 9. Validation of regulatory axis gene expression by qRT-PCR analysis (A) LINC00520. (B) miR-181d. (C) BCL2. (D) ZNF384. (E) ZNF24. (F) ZNF263. (G) ZNF341. (H) ZNF148. *P < 0.05, compared to A549.
was significantly downregulated in cisplatin resistance cells relative to the parental cells (P < 0.05; Figure 9(B)). At the same time, the expression of 5 ZNF family members, including ZNF384, ZNF24, ZNF263, ZNF341, and ZNF148 (P < 0.05; Figures 9(D)-(H)) were also upregulated in A549/DDP cells relative to A549 cells. These results were consistent with the conclusions of the bioinformatics analysis.
4. Discussion
Chemoresistance is a huge obstacle in the treatment of lung cancer. Accumulating evidence has suggested that autophagy is an important mechanism of chemoresistance [27] [28], which could be regulated by ceRNA networks. LncRNA MALAT1 might act as a ceRNA for miR-23b-3p and miR-30b, through which inhibitory effects on ATG12 and ATG5 were attenuated, leading to chemo-induced autophagy and chemoresistance in gastric cancer [29] [30]. Wang et al. [31] indicated that lncRNA H19 acted as a ceRNA of miR-194-5p, which is a suppressive miRNA of Sirtuin-1 (SIRT1), and enhanced 5-Fluorouracil (5-FU) chemoresistance via miR-194-5p/SIRT1/autophagy pathway in colorectal cancer cells. Additionally, Shi et al. [32] demonstrated that lncRNA HANR increased autophagy-related sorafenib resistance by inhibiting the miR-29b/ATG9A axis in hepatocellular carcinoma cells. Although some studies have suggested that ceRNA plays a role in autophagy-related chemoresistance in NSCLC [21] [22], the more potentially important RNA molecules and their regulatory networks have not yet been identified. Bioinformatics has a great advantage in cancer prediction. Cancer biomarkers can be screened by bioinformatics analysis and used for the accurate evaluation and management of the disease. In addition, specific bioinformatics tools have been used to reduce the difficulty of analysis, and speed up the study schedule and efficiency. In the present study, we identified 3179 DEmRNAs, 180 DEmiRNAs, and 160 DElncRNAs from three datasets of GSE108214, GSE157692, and GSE43493 by bioinformatic analysis, respectively. Given the crucial role of autophagy in chemoresistance, 35 autophagy-related DEmiRNAs were obtained by taking the intersection between ARGs and 3179 DEmRNAs. Then, the interactions among the DElncRNAs, DEmiRNAs, and ARDEmRNAs were predicted. In the end, a comprehensive autophagy-related ceRNA regulatory network inDDP resistance LUAD cell line was constructed, in which 14 DElncRNAs, 10 DEmRNAs, and 9 DEmiRNAs were involved. This might provide a direction for further research on chemoresistance in LUAD.
To explore the biological function of the DEGs in the ceRNA network, GO and KEGG enrichment analyses were performed. GO is an important initiative to unify the representation of gene and gene product attributes across all species. The enrichment analysis is to test whether a GO term is statistically enriched for 33 DEGs, but regrettably, no meaningful information about chemoresistance could be taken into consideration. It is worth noting that the KEGG pathway enrichment analysis indicated various combinations of these DEGs were involved in 5 pathways, including “HIF-1 signaling pathway”, “Proteoglycans in cancer”, “PI3K-Akt signaling pathway”, “MTOR signaling pathway”, and “Cholinergic synapse”. Moreover, most of them have been confirmed to be closely associated with chemoresistance. Hypoxia-Inducible Factors (HIFs) are a kind of transcription factor activated under a hypoxic environment to maintain cell activity [33]. The hypoxic and ischemic environment may trigger the HIF-1 signaling pathway, which can activate cell autophagy and further promote drug resistance to tumor [34]. Kitajim et al. [35] elucidated that HIF-1α knockdown could improve the 5-FU sensitivity of gastric cancer cells. Méndez et al. [36] also found that the activation of HIF signaling pathway could promote chemoresistance and was critical in maintaining the characteristics of cancer stem cells in pancreatic cancer. Therefore, the HIF-1 signaling pathway is closely associated with multidrug resistance. As is well known, PI3K/Akt signaling pathway and mTOR signaling pathway are classical autophagy pathways, which are closely related to chemoresistance [37] [38]. PI3K and Akt suppression can inhibit mTOR phosphorylation at the Ser2448 site, thereby enhancing the expression of autophagy-related proteins and inducing autophagy and chemoresistance [39]. Although there is no direct evidence to support the correlation between “Proteoglycans in cancer”, “Cholinergic synapse”, and autophagy or chemoresistance, a recent study implied that proteoglycans exerted diverse functions in the occurrence of cancer, mainly related to binding, hydration, transport, and resistance to drug [40], while autophagy has been reported to interact with synaptic plasticity, including structural changes in the synapse number, shape, size, and composition [41] [42].
Considering the importance of the key drivers in autophagy-related ceRNA network [43], a novel ceRNA subnetwork was extracted according to topological properties calculated by the NetworkAnalyzer Tool, including 2 lncRNAs (HOTTIP, LINC00520), 2 DEmiRNAs (miR-125b, miR-181d), and 3 DEmRNAs (BCL2, FKBP1A, and RPS6KB1). HOTTIP, miR-125b, and BCL2 have been studied for their roles in autophagy-related chemoresistance in cancer. Zhao et al. [44] provided results suggesting that HOTTIP was markedly upregulated in the cancer tissues of gastric cancer patients who were treated with gastrectomy and DDP chemotherapy, and HOTTIP expression up-regulation induced DDP resistance by regulating the miR-216a-5p/BCL2/Beclin1/autophagy pathway, which provided a novel strategy to overcome chemoresistance in gastric cancer. Several studies have revealed that miR-125b is also involved in autophagy and chemoresistance. For example, Yu et al. [45] reported that miR-125b functioned as an important downstream mediator upon the activation of the CXCL12/CXCR4 axis to confer 5-FU resistance in colorectal cancer probably through increasing autophagy both in vitro and in vivo. Wang et al. [46] found that miR-125b was downregulated and Forkhead box P3 (Foxp3) was upregulated in both thyroid cancer tissues and cell lines, and demonstrated a novel mechanism that miR-125b might promote autophagy and enhance the efficacy of DDP in thyroid cancer by negatively regulating Foxp3. Another study has also demonstrated that MAP kinase interacting serine/threonine kinase 2 enhanced chemoresistance of ovarian cancer by suppressing autophagy via miR-125b [47]. Moreover, the antiapoptotic protein BCL2 is frequently upregulated in acquired chemo-resistant cancer cells, which blocks drug-induced apoptosis and enhances autophagy [48]. The underlying mechanisms of LINC00520, miR-181d, FKBP1A, and RPS6KB1 in chemoresistance have been rarely reported. They might be novel and pivotal molecules driving the disease process and should be further studied.
We selected 7 genes from the novel ceRNA subnetwork and explored their role in the prognosis of LUAD patients using the Kaplan-Meier survival analysis. The results suggested that LINC00520, miR-181d, miR-125b, and BCL2 were identified to be associated with prognosis. Of them, high expression of LINC00520 and low expression of miR-125b, miR-181d, and BCL2 were related to poor OS in LUAD patients. Similar results have been reported by others. Xia et al. [49] confirmed that LINC00520 could predict poor prognosis and promote the progression of lung cancer by inhibiting miR-3175 expression. Luan et al. [50] reported that high expression of LICN00520 was a risk factor for the prognosis of melanoma patients and LINC00520 exerted its oncogenic role by competitively binding to miR-125b to promote Eukaryotic initiation factor 5A2 expression. Additionally, Wang et al. [51] demonstrated that LINC00520 expression level was elevated in NSCLC tissues and cells, and LINC00520 up-regulation was distinctly correlated with advanced tumor stage and shorter survival time in NSCLC patients. Zeybek et al. [52] showed that expression of miR-125b was significantly up-regulated in stage I LUAD tissues and the AUC was 0.875, which indicates that miR-125b might be a prognostic and diagnostic biomarker for LUAD. Zhao et al. [53] revealed the prognostic value of BCL2 in NSCLC patients, i.e. high expression of BCL2 protein in carcinoma cells did predict good disease-free survival and goodOS in NSCLC. Our results are consistent with these previous results, suggesting that the bioinformatics analysis of our study has high reliability and accuracy.
It is worth mentioning that the association between the expression level of miR-181d and prognosis with LUAD has not been reported until now. However, miR-181d, LINC00520, and BCL2 were identified as independent prognostic factors of LUAD using univariate and multivariate Cox regression analyses in this study. Based on this prognosis-associated LINC00520/miR-181d/BCL2 ceRNA, we established a risk score signature for survival prediction in LUAD patients, in which patients with high-risk scores exhibited a poor survival rate according to the result of the Kaplan-Meier survival curve and the risk score distribution plot. Unfortunately, subsequent univariate and multivariate Cox regression analyses found that this novel ceRNA signature could not independently predict OS in LUAD patients, and the most likely reason is the small sample size.
Emerging evidence has indicated that TF plays crucial role in regulating the aberrant expression of lncRNAs in human cancers [54]. In this study, JASPAR was employed to predict TFs for LINC00520, and 5 ZNF family members were identified as candidates to regulate LINC00520. The roles of ZNF family members in autophagy-related chemoresistance in cancer have been confirmed. Cui et al. [55] demonstrated that ZNF263 increased chemoresistance by activating Endoplasmic Reticulum Stress (ERS) related autophagy, functioning as a functional ERS-related tumor activator and a potential target for hepatocellular carcinoma therapy. In addition, Gao et al. [56] found that ZNF148 and TOP2A regulated each other through ceRNA regulatory mechanism in colorectal cancer, which has biological effects on cell proliferation. These results suggested that ZNF family members might affect autophagy and chemoresistance in tumors by regulating target genes through ceRNA regulatory mechanisms. Finally, we constructed a ZNF regulatory ceRNA subnetwork, which might provide novel therapeutic strategies for LUAD and new directions for further study of LUAD.
5. Conclusion
In this study, we successfully constructed an autophagy-related ceRNA regulatory network and explored the differential expression of lncRNAs, miRNAs, and autophagy-related mRNAs in DDP-resistant LUAD cells for the first time. In addition, the prognosis risk score model was established based on the independent prognostic factors (LINC00520/miR-181d/BCL2). Furthermore, the ZNF regulatory ceRNA subnetwork was constructed according to the regulatory roles of TF-lncRNA crosstalk. This study might contribute to our understanding of the molecular mechanism of LUAD chemoresistance. However, to shed lighter on the novel biomarkers and therapeutic targets of LUAD, further experimental studies are warranted.
Acknowledgements
This study was supported by the National Natural Science Foundation of China (No. 81703001), the Natural Science Foundation of Hebei Province (No. H2021406021), Hebei Province Medical Science Research Project (No. 20210247; 20221335), the government-funded training program for outstanding clinical medical talents, Chengde Medical University Scientific Research Major Projects (No. KY2020005), Hebei Province Key Research and Development Projects (No. 19277783D), and Hebei Province Pathogenic Biology Emphasis Subject Projects.
NOTES
*Ying Xu and Na Guo contributed equally to this work.
#Corresponding author.