The RNA sequencing data of LUAD in TCGA database were downloaded from the UCSC Xena online server (https://xenabrowser.net/). The gene expression matrices of GSE30219 and GSE31210 were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Data were transformed as follows: i=2ⁱ−1 for subsequent analysis.

Endocytosis-associated genes were collected from the Gene Ontology (GO) (https://geneontology.org/), Kyoto Encyclopedia of Genes and Genomes (https://www.kegg.jp/) database and the Reactome Knowledgebase (www.reactome.org). The search term ‘endocytosis’ was used. In total, 24 gene sets were obtained, and were intersected for endocytosis-associated genes.

Based on the clinical information of the LUAD cohort from TCGA, the RNA sequencing data of 60 tumor and paired paracancerous tissues were used for differential expression analysis via the DESeq2 R package (version 1.30.0) (19). DEGs were identified using the criteria of absolute value of Log fold-change (FC) >1.2, adjusted P≤0.05 and visualized as volcano plots. The intersection between DEGs and endocytosis-associated genes was used as a candidate for univariant Cox regression analysis.

Cox regression and Kaplan-Meier survival analyses were performed against candidate genes using the Survival R package (version 3.5–5; http://github.com/therneau/survival) to measure the association between expression level and overall survival in TCGA-LUAD, GSE30219 and GSE31210 cohorts. Genes with significant P-values were selected for Lasso Cox regression analysis using the glmnet R package (version 4.1–3; http://glmnet.stanford.edu/), and the refined prognostic model was constructed. The risk scores of each patient were calculated using the following formula: RiskScore=Σexp_i × β_i, where exp_i was the expression of gene i and β_i referred to the coefficient calculated by univariate Cox regression analysis of gene i.

Metascape (https://metascape.org) is an online tool for gene function annotation analysis, and it was used to analyze the GO enrichment of the 18 genes in the novel signature. Results were visualized as network plots.

The tumor purities of TCGA-LUAD, GSE30219 and GSE31210 cohorts were calculated using the ESTIMATE R package (version 1.0.13; http://r-forge.r-project.org/projects/estimate/) (20). Briefly, gene symbols and their corresponding expression levels were used as input data. After transformation and calculation by the ‘filterCommonGenes’ and ‘estimateScore’ functions, purity data were obtained. The Pearson's correlation coefficient was used for investigating the correlation between the riskScore and tumor purity was visualized as scatter plots.

Immune cell compositions were acquired via the CIBERSORTx (https://cibersortx.stanford.edu/) online analytical tool. By performing a series of steps to input the signature matrix file and expression file with default settings, the results were downloaded after the analysis was finished. The risk score distributions of each patient in the high and low groups were visualized as boxplots.

The human lung epithelial cell line BEAS-2B, and the NSCLC cell lines A549, H1299 and H1975 were obtained from the American Type Culture Collection. The NSCLC cell line PC9 was purchased from the National Collection of Authenticated Cell Cultures. The BEAS-2B and A549 cell lines were cultured in Dulbecco's Modified Eagle Medium (Thermo Fisher Scientific, Inc.) supplemented with 10% fetal bovine serum (FBS) (Gibco; Thermo Fisher Scientific, Inc.). The H1299, H1975 and PC9 cell lines were cultured in Roswell Park Memorial Institute-1640 (Thermo Fisher Scientific, Inc.) supplemented with 10% FBS. All cells were cultured at 37°C and 5% CO₂.

Cells were trypsinized to a single-cell suspension and counted. A total of 3,000 cells per well were seeded into 96-well culture plates. CQ (cat. no. S6999; Selleck Chemicals) was dissolved in DMSO to a storage concentration of 10 mM and diluted in culture media to a working concentration of 20 µM. After either 24, 48, 72, 96 or 120 h of CQ treatment, Cell Counting Kit-8 (cat. no. CK04; Dojindo Laboratories, Inc.) reagent was added, followed by a 2 h-incubation. The OD450 absorbance was measured using a microplate reader.

Cultured cells were seeded into glass bottom culture dishes. After 48 h of treatment with 20 µM CQ, EdU (cat. no. C0078S; Beyotime Institute of Biotechnology.) incorporation and staining were performed according to the manufacturer's protocol. Briefly, after pretreatment with 10 µM EdU for 2 h, cells were fixed with 4% paraformaldehyde solution for 10 min at room temperature and permeated using PBS with 0.2% Triton X-100 for 10 min at room temperature. Subsequently, Azid Alexa Fluor 594 was used to label EdU, and Hoechst 33342 was employed to stain nuclei for 10 min at room temperature protected from light. Samples were imaged using a laser confocal microscope (Olympus Corporation). The proportion of EdU-positive cells was calculated by counting the EdU-positive cells and the total number of cells using ImageJ (version 1.53; National Institutes of Health).

After 48 h of 20 µM CQ treatment, total RNA was isolated from cells using TRIzol^® (Invitrogen; Thermo Fisher Scientific, Inc.) reagent according to the manufacturer's protocol, and reverse-transcribed to cDNA by using the RevertAid First Strand cDNA Synthesis Kit (cat. no. K1622; Thermo Fisher Scientific, Inc.) at 42°C for 60 min, 70°C for 5 min, and then stored at 4°C. The SYBR Green Mix (cat. no. Q711; Vazyme Biotech Co., Ltd.) was used in a thermocycler instrument (QuantStudio 3; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions as follows: Predenaturation at 95°C for 5 min, followed by 40 cycles of denaturation at 95°C for 15 sec, and annealing and extension at 60°C for 30 sec. GAPDH served as the internal control. Primer sequences are summarized in Table SI. The mRNA levels were quantified using the 2^−ΔΔCq method (21).

Cell apoptosis was stained using Annexin V-FITC apoptosis detection kit (cat. no. C1062M; Beyotime Institute of Biotechnology) according to the manufacturer's protocol. Briefly, after 48 h of 20 µM CQ treatment, total cells in the supernatant medium and on the plate, were harvested and centrifuged at 500 × g for 10 min at room temperature. Cell pellets were resuspended in PBS buffer and counted. A total of 1×10⁵ cells were centrifuged at 500 × g for 10 min at room temperature, and resuspended in binding buffer (cat. no. C1062M; Beyotime Institute of Biotechnology) containing Annexin V-FITC and propidium iodide and incubated for 15 min at room temperature in the dark. Subsequently, cell apoptosis was detected using an Agilent NovoCyte flow cytometer (Agilent Technologies, Inc.). The data were analyzed using FlowJo (version 10.6.2; FlowJo LLC).

All statistical analyses were performed using GraphPad Prism (version 8, Dotmatics) and R (version 4.0.2; http://www.r-project.org/about.html). A paired Student's t-test was used for the analysis of two paired groups, and an unpaired independent-sample Student's t-test was used for the comparison of two experimental groups. One-way analysis of variance was performed to compare ≥3 groups, and Dunnett's post hoc test was used for multiple comparisons. The Pearson's correlation coefficient was used to investigate the correlation between risk score and either stromal or immune score, or tumor purity. Statistical results are presented as mean ± standard deviation. Univariate and multivariate Cox regression analysis were used to analyze the effect of risk factors on survival. Kaplan-Meier plots were used to assess the correlation between risk score and survival. P<0.05 was considered to indicate a statistically significant difference.

The analysis flow chart of the present study is presented in Fig. 1. First, LUAD data were downloaded from TCGA database via the UCSC XENA online server, and the transcriptome sequencing data of paired cancerous and adjacent normal samples were selected for the identification of DEGs using the DESeq2 R package. A total of 2,324 upregulated and 1,409 downregulated DEGs were identified using the filtering criteria of absolute value of LogFC >1.2 and P-adjust ≤0.05 (Table SII), and were visualized as a volcano plot (Fig. 2A). To select the DEGs correlated with the function of endocytosis, the list of DEGs and the endocytosis-associated gene set were intersected (Table SIII), and 138 candidate genes were obtained (Fig. 2B). Subsequently, a univariate Cox regression analysis was performed to evaluate the association between the expression levels of these genes and the overall survival of patients. A total of 33 genes were indicated to be significantly correlated with survival times, of which the expression levels of five genes were associated with poor prognosis, and the others were associated with a favorable prognosis (Fig. 2C and Table SIV). The expression levels of these genes in tumor tissues (n=525) and normal tissues (n=60) in the LUAD cohort from TCGA containing a total of 585 samples were compared and visualized as boxplots (Fig. 2D). Compared with the normal tissues, the expression levels of RAB27B, ATP6V0A4, LOXL2, HTR2B, SYT2, F2RL1 and IGHM were significantly increased in tumor tissues, whereas the expression levels of the other genes, excluding EREG, were significantly decreased in the tumor tissues.

Based on the aforementioned identified candidate genes, Lasso Cox regression analysis was used to screen the optimal genes to construct the refined prognostic model. The coefficient results calculated by Lasso Cox regression analysis are presented in Fig. 3A. According to the calculation of risk factors using Cox regression, and Kaplan-Meier survival analyses, a novel gene signature of 18 genes, including CFTR, RAB27B, ADRB1, DPYSL2, ATP6V0A4, EREG, SFTPD, LOXL2, HTR2B, SYT2, ALOX15, F2RL1, IL7R, PCSK9, ADRB2, GATA2, IGHM and MRC1, was involved (Fig. 3B and Table SV). The GO analysis results revealed the enriched signaling pathways and molecular events, among which the functions associated with ‘endocytosis’ and ‘regulation of vesicle-mediated transport’ were the most significantly enriched (Fig. 3C and D).

The risk score of every patient in the LUAD cohort from TCGA was calculated based on the expression of 18 genes and the distributions in the clinical information were examined. As indicated in Fig. 3E-H, the risk score of the novel signature in different pathologic tumor (T)/nodal (N)/metastasis (M) grades and tumor stages in LUAD was significantly increased compared with that in normal tissues. As the grades or stages increased, the risk score also had an upward trend. Furthermore, univariate and multivariate Cox regression analyses were performed to evaluate the prognostic value of the novel gene signature. As presented in Table I, the pathologic M [M1 vs. M0; hazard ratio (HR), 2.112; 95% confidence interval (CI), 1.235–3.612], pathologic N (N1 vs. N0; HR, 2.392; 95% CI, 1.700–3.365; N2 vs. N0; HR, 3.046; 95% CI, 2.084–4.452), pathologic T (T2 vs. T1; HR, 1.491; 95% CI, 1.046–2.125; T3 vs. T1; HR, 2.971; 95% CI, 1.765–4.999; T4 vs. T1; HR, 3.004; 95% CI, 1.547–5.834; TX vs. T1; HR, 4.794; 95% CI, 1.153–19.939) and tumor stage (stage II vs. stage I; HR, 2.451; 95% CI, 1.710–3.513; stage III vs. stage I; HR, 3.492; 95% CI, 2.390–5.102; stage IV vs. stage I; HR, 3.813; 95% CI, 2.203–6.599) indicated statistical significances in the univariant Cox regression analysis, whereas these were not statistically significant in the multivariant Cox regression analysis. Both univariant and multivariant Cox regression analyses of the risk score of the novel gene signature exhibited poor survival for patient prognosis (univariant; HR, 1.476; 95% CI, 1.336–1.630; multivariant; HR, 1.418, 95% CI, 1.272–1.580), demonstrating that this novel gene signature could be an independent prognostic marker.

To further validate the efficacy of the risk score in the prognostic prediction of patients with LUAD, TCGA-LUAD cohort (training set) was divided into two groups according to the median value of the risk score. The patients in the high-risk score group had a shorter survival time and an increased ratio of ‘DEAD’ status compared with the patients in the low-risk score group (Fig. 4A). The same results were illustrated in the two independent validation sets, GSE30219 and GSE31210, obtained from the GEO database (Fig. 4B and C). Heatmaps of the 18 genes in the training set and validation sets are presented in Fig. 4D-F.

Furthermore, Kaplan-Meier plots were used to examine the ability of the novel signature to predict prognosis. As presented in Fig. 5A-C, an increased probability of a longer survival time was positively correlated with a low-risk score of the novel gene signature in the training set and the two validation sets.

Tumor cells interact with tumor-infiltrating immune cells such as macrophages to induce the immunoinhibitory phenotype (22,23). Considering that these genes were endocytosis-associated genes, the present study attempted to investigate the association between the novel signature and tumor-infiltrating immune cells of LUAD. The ESTIMATE algorithm was used to calculate tumor purity scores, and these were compared with the risk scores of the novel signature in TCGA-LUAD cohort. As presented in Fig. 6A-C, the risk score was negatively correlated with the stromal and immune scores, but positively correlated with tumor purity, indicating that patients with a low-risk score had a more complicated immune microenvironment. Immune cell compositions in different groups with high or low risk scores were further examined using the CIBERSORTx algorithm. B cells, plasma cells, CD4⁺ T cells, monocytes, dendritic and mast cells had an increased composition in the low-risk score group compared with the high-risk score group; however, the NK cells and M0 and M1 macrophages were increased in the high-risk score group compared with the low-risk score group (Fig. 6D). In addition, the correlation between the risk score and marker genes of different immune cells was calculated. As presented in Table II, the majority of genes were negatively associated with these genes. These findings indicated that the novel gene signature could hint at the complexity of the tumor microenvironment.

To further investigate the expression levels of endocytosis-associated genes in the novel signature, total RNA was isolated from the lung epithelial cell line BEAS-2B and the lung cancer cell lines A549, H1975, H1299 and PC9, and then reverse transcribed for qPCR assays. The result was visualized as a heatmap in Fig. 7A. Compared with the levels in BEAS-2B cells, the expression levels of RAB27B, ATP6V0A4, LOXL2, HTR2B, SYT2, F2RL1 and IGHM were increased in lung cancer cells, and the other genes had decreased expression levels in tumor cells, which was consistent with the expression patterns in TCGA-LUAD dataset (Fig. 2D).

CQ is a conventional drug for treating amebiasis. CQ has been reported to induce a blockade of endocytosis via pH changes in lysosomes in cells (24). Given the roles of endocytosis in LUAD progression and the value of the endocytosis-associated signature in the prognostic prediction of patients with LUAD, CQ was used to treat lung cancer cell lines as well as normal lung epithelial cell lines. CQ treatment decreased the expression levels of the genes that were positively correlated with poor prognosis, and increased the expression levels of genes that were negatively associated with short survival time in lung cancer cells such as A549, H1299, H1975 and PC9 (Fig. 7B). However, the CQ-induced expression level changes in these genes were not consistently observed in BEAS-2B (Fig. 7B), indicating that CQ had a preferential effect on tumor cells. In addition, CQ treatment decreased the viability of A549, H1975, H1299 and PC9 cells, but no such effect was observed in BEAS-2B cells (Fig. 8A). EdU staining indicated reduced proliferation of lung cancer cell lines in the CQ treatment group compared with that in the DMSO group (Fig. 8B). Additionally, the level of apoptosis after stimulation with CQ was examined in these cells. Flow cytometry results revealed that the percentage of apoptotic cells was increased, and that of live cells was decreased after CQ treatment (Fig. 9). Taken together, these results suggest that CQ treatment resulted in an inhibited proliferation and increased apoptosis of lung cancer cells, presumably due to the reversal of gene expression patterns in tumors, demonstrating that targeting endocytosis could benefit the clinical treatment of patients with lung cancer, including patients with a high-risk score.