Forty paired LUAD and adjacent nontumorous tissues were obtained from patients undergoing surgery at the First Affiliated Hospital of Henan University and Puyang Hospital of traditional Chinese medicine, China, between 2018 and 2020. This study was approved by the Ethics Committee of Medical School of Henan University, China. All methods in this study were performed following the approved guidelines. Written informed consent was obtained from each patient before sample collection.

All NSCLC cell lines (HCI-H1299 and A549), immortalized lung epithelial cell line BEAS-2B and HEK293T were obtained from the Cell Bank of Type Culture Collection of the Chinese Academy of Sciences (Shanghai, China). The cells were cultured in DMEM medium (Corning, USA) with 10% fetal bovine serum (Pan biotechnology, Germany) at 37°C and 5% CO₂ (all cells were cultured under the same conditions) and were not passaged more than 25 times after thawing. Cells growing at an exponential rate were used for experiments. Cells were periodically evaluated to confirm Mycoplasma-negative status, and cell lines were authenticated by growth characteristics, examination of morphology and short tandem repeat analysis.

CCT6A in H1299 and A549 cells were stably knockdown by transduction with pre-made lentiviral short hairpin RNA (shRNA) (TranSheepBio, Shanghai, China). The shRNA vectors included TRCN0000062514 named as sh1CCT6A (target sequence: 5’- CGTGTCATTAGAGTATGAGAA-3’) and TRCN0000062515 named as sh2CCT6A (target sequence: 5’-CCAGAACATCTCTTCGTACTA-3’). The sequences of scramble control shRNA (shNC) were 5’- TTCTCCGAACGTGTCACGTAAT-3’. According to the manufacturer’s instruction, the recombinant lentiviral vector, envelope plasmid (pMD2.G) and packaging vector (psPAX2) were co-introduced into HEK293T cells through Lipofectamine 2000. After 48 h transduction, the culture supernatant was collected and filtered with 0.45μm filter, and then was used to infect H1299 and A549. The infected cells were screened by puromycin at 9.18 μmol/L. The efficiency of gene silencing was detected by Western blot analysis.

Cell proliferation and migration were assessed by colony formation and wound-healing assay. Details of the relevant contents have been described previously (13)

Tissues or cells were lysed on ice in RIPA lysis buffer containing 500 mM NaCl, 50 mM Tris pH 8.0, 1mM EDTA, 1% NP-40 and 1×cocktail of protease inhibitors (Roche, Lewes, UK). Protein concentrations were quantified according to the manufacturer’s instructions (Pierce, Rockford, IL). Protein lysates were separated by SDS-PAGE and transferred to polyvinylidene fluoride membranes. Details of the relevant content have been described previously (13). Information of antibodies used was provided in Supplementary Table 1 .

Raw sequencing data of LUAD-related RNAs expression and complete clinical data of the corresponding patients were downloaded from the TCGA database (https://portal.gdc.cancer.gov/). After homogenizing the TCGA raw data using the trimmed mean of M-values (TMM) method, the expression level of RNA was converted to a log2 value. Then, the R package of edgeR was used to screen for differentially expressed lncRNAs, miRNAs and mRNAs in LUAD and adjacent nontumorous tissues. Statistical significance was defined as P<0.05 and absolute Log Fold Change ≥1. Heatmaps and volcano plots of differentially expressed RNAs (DERNAs) were drawn using ggplots and heatmap software packages.

The miRcod database (http://www.mircode.org/) was used to predict the interaction between differentially expressed lncRNA and miRNA (14). The mRNA targeted by differentially expressed miRNA were searched from the miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/), miRDB (http://www.mirdb.org/), and TargetScan (http://www.targetscan.org/) databases using the Perl program (version:5.26.1) (15–17). A lncRNA-miRNA-mRNA ceRNA network based on the “ceRNA hypothesis” was established and visualized by Cytoscape software (http://cytoscape.org/).

The survival-related node genes in the ceRNA network were extracted by univariate Cox regression analysis with the threshold of Hazard Ratio ≠1 and P<0.05. To minimize overfitting, the R package “glmnet” was used to further extract prognostic genes for following multivariate Cox regression by a least absolute shrinkage and selection operator (LASSO) regression analysis. Then, the R package “caret” was used to randomly divide TCGA_LUAD patients into training and testing cohorts (18). A prognostic signature based on the ceRNA network was constructed in the training cohort by multivariate Cox regression (19). The prognostic gene signatures were shown as risk score=sum [gene expression×coefficient]. The risk score for each patient was calculated. We divided the LUAD patients into the high-risk and low-risk groups with the median risk scores as our cut-off. To verify the prognostic value of the above ceRNA-related genes signature, survival and ROC curve analyses were performed using the testing cohort and entire cohort as the validation set. Principal components analysis (PCA) was used to explore the distribution patterns of the different risk groups. A prognostic nomogram including risk scores and clinical features for predicting the likelihood of 1-, 3-, and 5-year OS was developed by R “rms” package. The calibration curves and C-index were used to evaluate the predictive accuracy of the nomogram. Finally, the mRNA expression levels of the prognostic signature genes were validated using Oncomine (http://www.oncomine.org/), TCGA and GSE32863 databases. The protein expression of levels was further verified by western blot analysis using 40 pairs of LUAD and adjacent nontumorous tissues.

To investigate the potential biological pathways and processes of the ceRNA network-related genes signature, we conducted the KEGG pathway analysis through Gene Set Enrichment Analysis (GESA) for the training cohort, with FDR<0.05 as a threshold for the significant pathways.

The CIBERSORT algorithm was utilized to estimate the fraction of 22 immune cell types in the LUAD samples from gene expression data in the training cohort. Samples with a CIBERPORT output value of P<0.05 were considered to meet the conditions for further analysis. The difference of immune cells in the proportion between the high- and low-risk groups was analyzed by Wilcoxon rank sum test. Additionally, it is worth noting that immune checkpoints are biomarkers for selecting LUAD patients for immunotherapy. Therefore, in this study, we analyzed the correlation between the ceRNA network-associated genes signature and key immune checkpoints (PD-1, PD-L1, CTLA-4, LAG3, TIM-3, TIGIT, CD80, CD276, TNFSF4 and VTCN1) using the R package “lmma”.

Chemotherapy is one of the effective methods for treating advanced patients with LUAD. The clinical response of each LUAD patient in high- and low-risk groups to chemotherapy was estimated according to the Genomics of Drug Sensitivity in Cancer (GDSC) data. Nine commonly used chemotherapy drugs, cisplatin, docetaxel, doxorubicin, erlotinib, etoposide, gemcitabine, paclitaxel, vinorelbine and cytarabine, were selected for the chemotherapeutic response prediction through the ridge regression using the “pRRophetic” R package. The half-maximal inhibitory concentration (IC50) predicted of each TCGA_LUAD patient was used to assess differential chemotherapeutic response (20).

SPSS21.0 (SPSS Inc., Chicago, IL, USA) and R software (version 3.6.0) were used for all statistical analyses. The difference between the two groups was assessed using the Student’s t-test. The continuous data are expressed as the mean ± standard deviation (SD). All statistical tests were two-tailed, and statistical significance was set at P<0.05.

We analyzed the differential expression of lncRNA, miRNA, and mRNA in LUAD by R software package using TCGA database containing 497 LUAD and 54 paracancerous samples, with P<0.05 and absolute Log2 Fold Change of 1 as the threshold. A total of 214 lncRNA (51 down- and 163 upregulated), 198 miRNA (87 down- and 111 upregulated) and 2989 mRNA (1344 down- and 1645 upregulated) were identified between LUAD and adjacent nontumorous samples ( Supplementary Table 2 ). The differentially expressed lncRNA, miRNA and mRNA distribution were visualized by volcano plots and heatmaps ( Supplementary Figure 1 ).

To explore the potential regulatory mechanism of ceRNA in LUAD, we tried to establish the ceRNA network for LUAD based on the ceRNA hypothesis (19). Using miRcode, miRTarBase, and miRDB databases, a total of interactions of 119 miRAN-mRNA pairs and 15 lncRNA-miRNA pairs were identified ( Supplementary Table 3 ). Finally, a LUAD-specific lncRNA-miRNA-mRNA ceRNA network, consisting of 7 lncRNAs, 15 miRNAs, and 95 mRNAs, including 117 nodes and 134 edges, was constructed and visualized ( Supplementary Figure 2 ).

Because this ceRNA network constructed above is composed of many genes and their interactions, it isn’t easy to clarify its diagnostic and prognostic significance. Therefore, a univariate Cox regression analysis (Hazard Ratio ≠1, p<0.05) was performed to screen the node genes related to the overall survival (OS) of LUAD patients in the ceRNA network using the TCGA_LUAD database. The results showed that 24 node genes were significantly associated with the OS ( Figure 1A and Supplementary Table 4 ). Subsequently, to reduce the complexity of the risk model, Lasso regression analysis was used to remove genes with relatively lower correlation, and 9 of the 24 prognostic genes were screened out ( Figures 1B, C ). Then, we randomly divided 462 TCGA_LUAD patients with survival data into training and testing cohorts using the R package “caret”. A prognostic signature model was developed in the training cohort based on multivariate Cox regression analysis ( Supplementary Table 5 ). Then, four candidate signature genes were identified, namely, GPI, IL22RA1, CCT6A, and SPOCK1. Figures 1D , E presented the forest plots and heatmap of the four prognostic genes. Furthermore, a 4-gene signature-based risk score formula was constructed as Risk score=0.314∗GPI+0.127∗ IL22RA1+0.330∗CCT6A+0.104∗SPOCK1. According to the median risk scores, 232 LUAD patients in the training cohort were divided into high-risk and low-risk groups. As depicted in Figures 2A, B , the increase of risk score was related to the poor OS of patients with LUAD. Kaplan-Meier curve showed that OS decreased in patients in the high-risk group ( Figure 2C ). Time-dependent ROC analysis observed that compared with each of the above genes, the four-gene prognostic signature had larger AUC values ( Figure 2D and Supplementary Figure 3 ). Moreover, principal component analysis (PCA) clearly identified a significantly different distribution between the two risk groups ( Figure 2E ).

Finally, to verify the predictive value of the four-gene signature, we used the testing cohort (n=230) and the entire cohort (n=462) as the validation set to evaluate the findings from the training cohort. Similar to the results of the training cohort, the KM curves of the two validation sets demonstrated that patients in the low-risk group exhibited better OS ( Figures 2F, I ). The AUC of the four-gene signature was 0.633, 0.635, 0.665, 0.673, 0.672 and 0.668 at the 1-year, 3-year, and 5-year timepoints in the testing cohort and the entire cohort, respectively ( Figures 2G, J ). In addition, PCA also showed similar results as the training cohort ( Figures 2H, K ).

We performed external validation of four prognostic signature genes. The Oncomine database analysis observed that compared to normal lung tissues, the expression levels of GPI, IL22RA1, CCT6A, and SPOCK1 mRNA in LUAD tissues were significantly higher ( Figure 3A ). Moreover, we further confirmed mRNA expression of these four genes in LUAD from TCGA’ paired samples and GSE32863 databases ( Figures 3B, C ). Additionally, we investigated the protein expression levels of these four genes in 40 pairs of LUAD tissues and adjacent nontumorous tissues by Western blot analysis. The results showed that GPI, IL22RA1, CCT6A, and SPOCK1 expression in the tumor tissues (T) was markedly higher than in the control tissues (N) (n=40, P<0.05; Figures 3D, E and Supplementary Figure 4 ).

To further evaluate the prediction performance of the ceRNA network-related genes signature, we selected four other published gene signatures obtained from Li’s (21), Mo’s (22), Sun’s (23) and Zhang’s (24) for comparison. The risk score of each patient was calculated according to the corresponding genes in these four models using the same method (multivariable Cox regression analysis) in the training cohort and then evaluated the time-dependent ROC. Figure 2D and Figures 4A –D revealed that the AUC of the ceRNA network-related gene signature for 5-year OS was 0.678, which was significantly larger than that of Li’s (0.616), Mo’s (0.576), Sun’s (0.603) and Zhang’s (0.622) gene signatures. The C-index of all prognostic signatures calculated by the restricted mean survival (RMS) package showed that our model had the highest C-index with 0.668 ( Figure 4E ). Moreover, the RMS time curve of all five prognostic models also demonstrated that our 4-gene signature performed best at a time period greater than 8 years ( Figure 4F ). These results suggested that the ceRNA network-related genes signature might provide better prognosis evaluation performance for LUAD.

Univariate Cox regression analysis showed that tumor stage (HR=1.706, 95% CI=1.353-2.150, P<0.001), recurrence (HR=3.133, 95% CI=1.908-5.143, P<0.001), and the risk score (HR=2.144, 95% CI=1.563-2.941, P<0.001) were closely correlated with OS in training cohort ( Figure 5A ). Multivariate Cox regression analysis further confirmed the above results ( Figure 5B ). Therefore, these three factors were combined to construct a compound nomogram for predicting the OS of patients with LUAD at 1-, 3- and 5-year ( Figure 5C ). The calibration plot of the nomogram for predicting 3-year OS of LUAD patients in the training cohort showed great consistency between actual observation and nomogram prediction ( Figure 5D ), and the nomogram model’s C-index for prediction of OS was 0.778 (95% CI=0.728-0.829, P=1.97e-27). Moreover, the AUC value of the nomogram for predicting 3-year OS was larger than that of the stage, recurrence and risk score, suggesting that using this nomogram to predict OS might bring more net benefit ( Figure 5E ). Finally, we further assessed the predictive value of the four-gene prognostic nomogram using the testing cohort and the entire cohort. The calibration plots ( Figures 5F, H ) and the time-dependent ROC curves of risk score ( Figures 5G, I ) were consistent with the results derived from the training set.

To further explore the potential biological pathways and processes related to the four-gene signature, we conducted gene set enrichment on samples within the training cohort by GSEA. We found that critical pathways associated with tumorigenesis, including cell cycle, DNA replication, P53 signaling pathway, proteasome, and spliceosome, were significantly enriched in the high-risk group ( Figure 6A and Supplementary Table 6 ). Additionally, we investigated whether the risk model was related to the tumor immune microenvironment. We plotted the heatmap of 22 tumor-immune cell types, showing the distribution of these immune cells ( Figure 6B ). Then, Wilcoxon rank-sum test was used and revealed that high-risk LUAD patients had significantly higher proportions of T cells CD4 memory activated, NK cells resting, macrophages M0 and macrophages M1 and lower proportions of T cells CD4 memory resting, monocytes, and Mast cells resting ( Figure 6C ). Furthermore, we further analyzed the correlation between the risk groups and the expression of immune checkpoint molecules, including programmed death-1 (PD-1), programmed death-ligand 1 (PD-L1), cytotoxic T-lymphocyte-associated protein 4 (CTLA-4), lymphocyte activation gene-3 (LAG3), T-cell immunoglobulin and mucin-domain containing-3 (TIM-3), T cell immunoreceptor with Ig and ITIM domains (TIGIT), CD80, CD276, tumor necrosis factor superfamily member 4 (TNFSF4), and V-set domain-containing T-cell activation inhibitor-1 (VTCN1). The results showed that compared with the low-risk group, the high-risk group had markedly higher expression levels of TNFSF4, CD274, PD-L1, and LAG3 ( Figure 6D ). Therefore, the heterogeneity of immune cell infiltrations and immune checkpoint molecules expression observed in these results may provide potential prognostic indicators and targets for immunotherapy in patients with LUAD.

Besides immune checkpoint blockades therapy, chemotherapy is still an effective treatment for advanced LUAD patients. Thus, we attempted to investigate the response to the common chemotherapeutic drugs in low- and high-risk patients with LUAD in the entire cohort. The IC50 values of the high- and low-risk groups were calculated based on the GDSC data. The results demonstrated that high-risk LUAD patients showed increased sensitivity to cisplatin, docetaxel, doxorubicin, erlotinib, etoposide, gemcitabine, paclitaxel and Vinorelbine, while there was no significant difference in IC50 value of cytarabine between the high-risk group and low-risk group, which indicated that the four-gene risk model might act as a potential predictor for chemosensitivity ( Figure 7 ).

Finally, to obtain more mechanistic insights into the impact of GPI, IL22RA1, CCT6A and SPOCK1 on the progression of LUAD, Gene Set Enrichment Analysis (GESA) was performed using TACG_LUAD database to compare the high expression and low expression of GPI, IL22RA1, CCT6A and SPOCK1, respectively. We observed that many important regulatory genes involved in the PI3K-AKT-mTOR signaling pathway, which plays an important role in regulating various oncogenic processes (25), were markedly enriched in cells with high GPI, IL22RA1, CCT6A and SPOCK1 expression ( Figure 8A ). In view of the highest normalized enrichment score (NES) of CCT6A (NES) (NES=2.21, P= 6.19e-4), we mainly chose CCT6A to explore the relevant mechanism. We first examined the expression level of CCT6A in BEAS-2B, A549 and H1299 cell lines by Western blot ( Figure 8B ). Due to the high expression level of CCT6A in two non-small cell lung cancer (NSCLC) cell lines, A549 and H1299, we used short hairpin RNA (shRNA) targeting CCT6A to silence CCT6A. The results showed that CCT6A knockdown substantially reduced proliferation and migration of A549 and H1299 ( Figures 8C–F ). Moreover, Western blot analysis revealed that CCT6A silencing significantly inhibited epithelial-mesenchymal transition (EMT) of A549 and H1299, evidenced by a marked reduction in the protein levels of N-cadherin and marked increase of E-cadherin. Notably, we also observed that CCT6A knockdown markedly decreased the protein levels of p-PI3K and p-AKT ( Figure 8G ), suggesting that CCT6A may affect the EMT of LUAD cells by activating PI3K/AKT pathway and then affect the malignant of LUAD and prognosis of LUAD patients. In the future, we will continue to focus on the specific mechanism of GPI, IL22RA1 and SPOCK1 in LUAD.