All genomic and clinical data were obtained from the Gene Expression Omnibus (GEO) database under accession numbers GSE83666, GSE75309 and GSE103352. To evaluate the impact of ceRNA networks expression on the survival of NSCLC patients, another lncRNA, mRNA and miRNA expression data set of NSCLC samples was downloaded from The Cancer Genome Atlas (TCGA) database.

Data on microarray analyses of lncRNA, miRNA and mRNA expression in NSCLC cells from patients, which were resistant or sensitive to EGFR-TKIs were obtained from the GEO database. All expression data were standardized and analyzed using limma (Ritchie et al., 2015). We identified lncRNAs, miRNAs and mRNAs differentially expressed between sensitive or resistant cancer, based on the criteria of -log10 (fold change) < 0.01 and p < 0.05.

Differentially expressed lncRNAs, miRNAs and mRNAs were analyzed based on Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways using clusterProfiler. Enrichment was defined as p < 0.05 (Yu et al., 2012).

Interactions between differentially expressed miRNAs and mRNAs in NSCLC resistant to EGFR-TKIs were identified by screening the RNAInter database (Lin et al., 2020). And, 104,857,6 interactions pairs are extracted from the database. A similar procedure was followed to identify interactions between differentially expressed lncRNAs and miRNAs.

The lncRNA-miRNA pairs and miRNA-mRNA pairs, identified as described above, were merged together into a ceRNA network based on lncRNAs, miRNAs and mRNAs related to NSCLC resistance to EGFR-TKIs. Pearson correlation analysis was carried out among the species in the network, and interactions associated with a correlation coefficient R > 0.75 and p < 0.01 were considered reliable. The potential network of biological processes enriched in NSCLC resistant to EGFR-TKIs was built using Pathway Enrichment Analysis ClueGO (version 2.3.2) (Bindea et al., 2009).

To identify associations between the expression of nodes in the ceRNA network and overall survival of patients, Kaplan‐Meier survival analysis was performed using the ‘survfit’ function in the ‘survival’ package in R (https://CRAN.Rproject.org/package=survival) (Xu et al., 2014). We also screened ceRNA regulatory sub-networks of interacting lncRNAs, miRNAs and mRNAs in the TCGA data in order to identify those significantly associated with survival (p < 0.05).

The NSCLC cell lines PC9, sensitive to Gefitinib, and PC9-GR, resistant to Gefitinib, were obtained from Heilongjiang Cancer Institute (Harbin, China). As described, an EGFR-TKI resistant cell line (PC9-GR) were established by treating EGFR-TKI-sensitive PC9 cell with continuous gefitinib culture (Huang et al., 2020b). These cells were cultured in RPMI 1640 medium (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% fetal bovine serum (Gibco; Thermo Fisher Scientific, Inc.) and 1% penicillin-streptomycin at 37°C in a humidified atmosphere of 5% CO₂.

This quantitative real-time PCR analysis was conducted as described using the following primers, all synthesized by Sangon (Shanghai, China) (Wang et al., 2019a): ABTB1, 5′-GGC​GGG​ATT​ACT​ATG​ACG​AC-3′ (forward) and 5′-GCC​TGA​GAA​CCA​CGA​CAC​TC-3′ (reverse); NPTN, 5′-GCT​CCT​AAA​GCA​AAC​GCC​ACC​A-3′ (forward) and 5′-TCT​TGC​GCC​ATA​TCC​AGT​CTG​G-3′ (reverse); SERPINE1, 5′-CTC​ATC​AGC​CAC​TGG​AAA​GGC​A-3′ (forward) and 5′-GAC​TCG​TGA​AGT​CAG​CCT​GAA​AC-3′ (reverse); and hsa-miR-150-5p-R, 5′-CCG​TCT​CCC​AAC​CCT​TGT​AC-3′ (forward) and 5′-CAG​TGC​AGG​GTC​CGA​GGT-3′ (reverse).

Cell viability assay was performed as previously demonstrated (Wang et al., 2019a). The cells were seeded in 96-well plates with 2*10³ cells per well with or without tiplaxtinin (40 μm) (HY-15253, MedChem Express) incubated for 24 h before seeding on the 96-well culture plate. Subsequently, gefitinib was added to the wells incubating for 48 h at the indicated concentrations. After gefitinib treatment, the viability of cells was measured.

All statistical analyses were carried out using SPSS 24.0 software (IBM Corp, Armonk, NY, United States ) and GraphPad Prism software 8.0 (GraphPad Software, San Diego, CA, United States). Data were expressed as mean ± standard deviation. Differences between samples resistant or sensitive to EGFR-TKIs were assessed for significance using Student’s t test. The ability of the core ceRNA networks to predict resistance to EGFR-TKI was assessed based on the area under curves (AUCs) of the receiver operating characteristic (ROC) curves. The AUC is obtained by summing the areas of the parts under the ROC curve. Survival was analyzed using the Kaplan-Meier approach, and survival was compared between groups using the log-rank test. Differences were considered significant if they were associated with a two-tailed p < 0.05.

Expression profiles for NSCLC cells sensitive or resistant to EGFR-TKI were acquired from the GEO databases, then standardized and analyzed using limma. The flow chart summarizing the framework is exhibited in Figure 1. The GSE83666 and GSE75309 databases were combined and divided into resistant group and non-resistant group. 1755 up-regulated mRNAs, 1,436 down-regulated mRNAs, up-regulated 4,598 lncRNAs and 4,391 down-regulated lncRNAs were screened with p value <0.05. GSE103352 data set was divided into resistant and non-resistant groups. 505 up-regulated miRNAs and 578 down-regulated miRNAs were screened with p value <0.05. DEGs were depicted in a Manhattan chart (Figure 2A) and a heatmap (Figures 2B–D), revealing 6,858 genes that were upregulated in resistant cancer, comprising 1,755 mRNAs, 4,598 lncRNAs, and 505 miRNAs; as well as 6,405 downregulated genes, comprising 1,436 mRNAs, 4,391 lncRNAs, and 578 miRNAs.

To gain insights into the function of DEmRNAs identified in the GSE83666 and GSE75309 datasets, we analyzed their enrichment in GO processes and KEGG pathways. A total of 968 biological processes and 31 KEGG pathways were identified with the filter criteria of adjust p-value< 0.05, of which 20 biological processes and 15 KEGG pathways have previously been linked to lung cancer resistance to EGFR-TKIs (Figures 3A,B). In addition, Gene Set Enrichment Analysis (GSEA) software was applied to identify GO-BP and KEGG pathway of the candidate hallmarks (Subramanian et al., 2005). GSEA identified seven biological processes and nine KEGG pathways significantly involved in resistance to EGFR-TKIs, including cellular amino acid metabolic processes, p53 signaling and proteasome pathways (Figures 4A,B). A potential network of biological processes enriched in NSCLC resistant to EGFR-TKIs was built using Pathway Enrichment Analysis ClueGO (Supplementary Figure S1).

We constructed a ceRNA regulatory network linking lncRNAs, miRNAs and mRNAs by identifying miRNA–mRNA interactions (Figure 5A) and lncRNA–miRNA interactions (Supplementary Figure S2). The correlation function in Hmisc package was used for correlation analysis, and 329 ceRNA sub-networks were obtained. Pearson correlation analysis of lncRNAs and mRNAs in the network were screened for those satisfying the criteria of R > 0.75 and p < 0.01, leading to 31 ceRNA regulatory sub-networks (Figure 5B). When enriched KEGG pathways were added to the network, p53 signaling and peroxisome pathways emerged as potentially linked to resistance (Figure 5C). In this way, our study provides valuable leads for future experiments to explore the molecular mechanisms of NSCLC resistance to EGFR-TKIs.

To get the core ceRNA sub-networks, we screened these pathways for those satisfying the criteria R > 0.85 and p < 0.01, identifying 12 core ceRNA regulatory sub-networks (Figure 6A; Table 1). These core ceRNA regulatory sub-networks (triple 1–12) were expressed at significantly different levels between sensitive or resistant samples (Figure 6B). The ability of these regulatory sub-networks to predict resistance was assessed based on AUCs. Triple 3, 7, 9, 10 and 12 showed better predictive performance than other ceRNA regulatory sub-networks (Supplementary Figure S3).

Numerous researches have reported that ceRNA may be involved in EGFR-TKI resistance, and may influence the progression of EGFR-mutant NSCLC. To examine whether ceRNAs are related to prognosis of NSCLC patients carrying EGFR mutations, Kaplan-Meier analysis and Cox regression modeling were performed based on TCGA data on progression-free interval (PFI) and overall survival (OS). Triple 9 and 12 were associated with PFI and OS (Figures 7, 8).

In addition, the Kaplan-Meier plotter online tool (http://kmplot.com/analysis/) was used to verify the effect of SERPINE1 on the survival of lung cancer. As shown in Supplementary Figure S4, the lung cancer patients with elevated SERPINE1 expression levels had poor OS than patients with corresponding low expression levels (p < 0.001). Consistently, as demonstrated in Supplementary Figure S5, high SERPINE1 expression predicts a poor prognosis in NSCLC using the web-based tools in The Human Protein Atlas (https://www.proteinatlas.org) based on the TCGA database.

To verify the reliability of the bioinformatics results, we investigated the expression of the core ceRNA triple 9 and 12 using qRT-PCR in NSCLC cell lines. Expression of the two lncRNAs ABTB1 and NPTN as well as the hsa-miR-150–5p was lower in the PC9-GR cell line than in the PC9 line (Figures 9A–D). Conversely, expression of the mRNA SERPINE1 was higher in the PC9-GR cell line than in the PC9 line. These results were consistent with the bioinformatics analysis. As shown in Figure 9E, inhibition of SERPINE1 can affect the sensitivity of EGFR-TKI drug-resistant cells to gefitinib at different concentrations.