Transcriptomic landscape and cell type classification of NSCLC

To analyze the transcriptional characteristics of NSCLC, we performed high-precision scRNA-seq analysis with tissues from 19 primary lung cancer patients who underwent surgery, including 14 ADC patients, 3 SCC patients, 1 combined small cell lung cancer (C-SCLC) patient, and 1 patient with mixed adenocarcinoma and neuroendocrine carcinoma (MANEC) [40]. Cells from healthy normal tissues from a pulmonary chondroid hamartoma (PCH) patient were collected as a control (Fig. 1a; Additional file 1: Figs. S1, S2a). In total, we obtained the single-cell transcriptome of 9002 cells. After stringent filtering (Additional file 1: Fig. S2b), we retained 7364 (81.8%) high-quality individual cells for subsequent analyses (Additional file 1: Fig. S2a).

To classify major cell types, we performed a t-distributed stochastic neighbor embedding (t-SNE) analysis using SCENIC identified cell clusters [36] (Fig. 1b). Based on the expression patterns of known canonical cell type markers, we identified six major cell types, including epithelial cells (EPCAM), fibroblasts (THY1), B cells (CD79A), T cells (CD3D), myeloid cells (CD68), and mast cells (KIT) (Fig. 1c, d). We observed variations in the proportions of these six cell types in the 16 patients sequenced without cell preselection (Fig. 1e). To comprehensively analyze the molecular characteristics of cancer cells, we focused our study on epithelial cells. Finally, we totally obtained 3373 epithelial cells from 16 patients of normal and tumor tissues for subsequent analysis, and the cell distribution in each patient was shown in Additional file 1: Fig. S2c.

scRNA-seq uncovered mixed-lineage tumor cells in NSCLC

In this study, three different lung cancer subtypes were covered: ADC, SCC, and NET. To further distinguish these lung cancer lineages at the single-cell level, we classified the cancer subtypes based on the expression of clinically well-established markers (Fig. 2a). According to the gene expression patterns of cancer subtype-specific markers, we redefined cancer types for each tumor epithelial cell, namely, ADC, SCC_ADChigh, SCChigh_ADC, NET_ADC, and triple-positive (coexpressing ADC/SCC/NET markers) (Fig. 2a, the details are in the “Methods” section). Intriguingly, we found that 55–98% of cancer cells simultaneously expressed classical marker genes (NKX2-1, KRT7, and NAPSA for ADC; TP63, KRT5, and KRT6A for SCC; CHGB, SYP, and NCAM1 for NET) for two or even three different histologic subtypes of NSCLC in the same individual cell (defined as a mixed-lineage cell) in six patients (P5, P10, P11, P13, P19, and P22). If we extend the marker genes to include the non-classical ones such as MUC1 for ADC, SOX2 for SCC, and ASCL1 for NET, the mixed-lineage cancer cells were notably observed in all of the sixteen patients from whom we isolated epithelial cells in tumor tissues. The ratio of the redefined cancer subtype across every patient is demonstrated in Additional file 1: Fig. S3a. To the best of our knowledge, this is the first study to identify mixed-lineage cancer cells at the whole-transcriptome level.

Next, to investigate the molecular heterogeneity of mixed-lineage cancer cells, we performed principal component analysis (PCA) for 1400 mixed-lineage cancer cells (Additional file 1: Fig. S3b) and found that they were not separated by cancer subtypes; instead, different cancer subtypes were mixed together within the same patients, which indicating that heterogeneity across patient was greater than cancer types within the patient. Furthermore, for these three lineages of cancer cells, we selected the top 50 individual cells that showed the strongest lineage-specific marker gene expression signatures for each lineage and performed PCA (Additional file 1: Fig. S3c). The cells of these three different lineages were divided into three independent clusters accordingly, which indicated that these 150 selected cancer cells accurately represent the unique molecular features of the ADC, SCC, and NET lineages. Therefore, we performed differentially expressed gene (DEG) analysis for these cancer cells to identify new markers for each subtype (Additional file 1: Fig. S3d; Additional file 2: Table S1). We found that the cancer subtype markers (GRP, CHGB, NEUROD1, and CHGA for NET; KRT5, DSC3, KRT6B, TP63, and KRT6A for SCC; and NKX2-1, KRT7, NAPSA, and MUC1 for ADC) were specifically expressed in the corresponding cancer subtypes. Furthermore, genes such as TUBB3 and MEST for NET, TRIM29 and CSTA for SCC, and CEACAM6 for ADC can be used as candidate markers to distinguish these three subtypes.

To verify the identity of each individual cell, we entered all tumor epithelial cells into a ggtern plot with cells scored based on the expression levels of markers specific for the ADC, SCC, and NET lineages (Fig. 2b). Consistent with the above cell clustering and cancer type redefinition data (Fig. 2a), the ggtern plot also showed that some intermediate cells lay in the middle among the ADC, SCC, and NET lineages, indicating their mixed-lineage nature (Fig. 2b). Specifically, the majority of cells from patient P12 displayed the molecular features of all three lineages (Additional file 1: Fig. S3e). Next, the principal component scores of these three lineages for each patient were counted (Fig. 2c). Among the three patients diagnosed with SCCs, P12 showed a lower SCC score than P1 and P2. Meanwhile, ADC patients P5 and P10 displayed lower ADC scores than the other patients with ADCs. For patient P22 with MANEC, the ADC score was higher than the NET score, indicating that P22 had stronger ADC characteristics. All the results above showed that there were mixed-lineage cancer cells that coexpressed marker genes of different cancer lineages in NSCLC at the whole-transcriptome level. Moreover, we analyzed the expression of the lineage-specific marker genes from previously published datasets using the same method [20, 22]. 12.4% (237 out of 1910) and 11.8% (432 out of 3668) of the cells in the tumor tissue were identified as mixed-lineage tumor cells, respectively (Additional file 1: Fig. S3f). In the former dataset downloaded, the information about every single cell is from which individual patient was missing, and we cannot trace the origin of the mixed-lineage tumor cells to each individual patient. However, in the latter dataset, mixed-lineage cancer cells were identified in every one of the seven patients analyzed, with the ratio ranging from 7.3 to 24.8%. These results clearly confirmed the prevalent existence of mixed-lineage tumor cells in NSCLC patients.

NKX2-1 and TP63 are the best-known markers to discriminate between ADC and SCC. In our study, we found that four ADC patients (P5, P10, P11, and P13) had a large proportion of tumor cells coexpressing NKX2-1 and TP63 (Additional file 1: Fig. S4a). To further confirm this, we performed multiplex fluorescent immunohistochemistry (IHC) staining. For P5, the tumor areas on the section were strongly and diffusely positive for NKX2-1 and TP63, and 51% (35,654 out of 70,292 cells analyzed) of the individual cancer cells co-expressed these two markers simultaneously (Fig. 2d). In contrast, cells from adjacent normal tissues did not show double positivity for NKX2-1 and TP63. Therefore, multiplex fluorescent IHC staining verified our scRNA-seq results at the protein level. For patients P10 and P11, the staining results also showed the presence of NKX2-1 and TP63 double-positive cells in tumor tissues (Additional file 1: Fig. S4b). Surprisingly, we found that many cancer cells from P11 not only coexpressed NKX2-1 and TP63 but also showed high expression of CHGB and UCHL1, which are markers of neuroendocrine cells [41] (Additional file 1: Fig. S4a). Multiplex fluorescent IHC staining of TTF1 and UCHL1 on tumor tissues from P11 also verified the coexpression of these two markers in a high proportion of tumor cells (Additional file 1: Fig. S4c). Moreover, we confirmed that NKX2-1 and TP63 were coexpressed in the same tumor cells by RNA in situ hybridization (Fig. 2e; Additional file 3). In summary, we verified that mixed-lineage cancer cells are prevalent in many different patients with NSCLC.

Mixed-lineage and single-lineage tumor cells in the same patient originate from common tumor ancestor cells

Patients P19 and P22 possessed tumors with two combined components. We found that for patient P19, who had C-SCLC, most cells highly expressed not only the NET markers CHGA, CHGB, ASCL1, and NEUROD1 but also the ADC marker NAPSA. Only a small fraction of cancer cells solely expressed ADC markers. For patient P22 with MANEC, almost all cancer cells highly expressed NAPSA, and some of them coexpressed ASCL1. Therefore, unlike ADC patients P5, P10, P11, and P13, patients P19 and P22 mainly had NET and ADC mixed-lineage cancer cells (Fig. 3a). The mixed-lineage cancer cells showed mixed features of different subtypes of NSCLC, indicating that they were highly plastic.

Fig. 3

Mixed-lineage and single-lineage tumor cells in the same patient originate from common tumor ancestor cells. a Line plots showing the gene expression levels of clinical lineage-specific markers for single tumor cells from P19 and P22. b Heatmap showing the identified mitochondria mutations specific to tumor epithelial cells from P19. c PCA plot of single cells from P19, colored by redefined cell type, sample regions, and selected mitochondrial mutations. d Heatmap showing identified mitochondria mutations specific to tumor epithelial cells from P22. e PCA plot of single cells from P22, colored by cancer type, sample regions, and selected mitochondrial mutations. f PCA profiles of single cells from P19, colored by CNV clusters inferred by single-cell RNA-seq data. Only the cells with a corresponding CNV subtype with more than 10 cells were used for PCA analysis. g PCA profiles of single cells from P22, colored by CNV clusters inferred by single-cell RNA-seq data. Only the cells with a corresponding CNV subtype with more than 10 cells were used for PCA analysis. h Survival analysis of disease-free survival (DFS) and overall survival (OS) in ADC and SCC samples from TCGA. Mixed-lineage features were calculated based on the expression levels of the identified lineage marker genes for ADC, SCC, and NET. Samples with high scores and low scores represent high lineage-mixing features and low lineage-mixing features, respectively. i Mixed-lineage features evaluation for EGFR mutant samples and wild-type samples in ADC and SCC samples from TCGA, respectively. Mixed-lineage features were calculated based on the expression levels of the identified lineage marker genes for ADC, SCC, and NET. The higher the score indicates the higher lineage-mixed features

To determine the relationship of mixed-lineage cancer cells and single-lineage cancer cells in the same patient, we investigated the cancer phylogenetic structure based on our scRNA-seq data using mitochondrial mutation-based lineage tracking analysis and single-cell copy number variation (CNV) analysis following previously described methods [33, 34, 42]. We identified two tumor-specific mutations, 2645_G and 13226_G, in patient P19 (Fig. 3b). According to the cancer type identification above, tumor cells from P19 were mainly composed of the NET_ADC and triple-positive subtypes. Interestingly, the NET_ADC and triple-positive cancer subtypes shared these two mutations, suggesting that these two types of mixed-lineage cancer cells had common tumor ancestors in patient P19 (Fig. 3c). In addition, we identified more mitochondrial mutation sites specific to tumor cells for patient P22, such as 3558_T, 66_G, 3916_G, 5254_C, and 1985_T (Fig. 3d). These mutations were shared by all four cancer subtypes: ADC, NET_ADC, SCC_ADChigh, and triple-positive (Fig. 3d, e). Therefore, the mitochondrial mutation results indicated that the mixed-lineage cancer cells and ADC-based single-lineage cancer cells in the same patient had common tumor cell ancestors. Although we can not deduce the direction of the lineage changes, the most likely scenario is that in a NSCLC patient, a specific lineage of the normal epithelial cells was first transformed to single-lineage tumor cells during tumorigenesis, which further changed to mixed-lineage tumor cells. This does not exclude the possibility that at the late stage of tumorigenesis the mixed-lineage and single-lineage tumor cells can interchange easily with each other due to their plasticity and flexibility.

We next performed a single-cell CNV analysis based on the scRNA-seq data to further support the results of mitochondrial mutation analysis. We selected single-cell gene expression data from all normal epithelial cells as a control to calculate the CNVs of tumor epithelial cells. As demonstrated in Fig. 3f and Additional file 1: Fig. S5a, we found that two mixed subtypes, NET_ADC and triple-positive, in patient P19 had the same CNV patterns. To validate the accuracy of our single-cell CNVs, we assessed CNVs using the data generated from whole-genome sequencing (WGS) on bulk cells from the same patient and found that the results were consistent with these CNVs inferred by single-cell RNA-seq analysis (Additional file 1: Fig. S5b). For example, gain of chromosome 3 and chromosome 5 and loss of chromosome 4 in tumor cells were clearly evident in CNVs inferred by both scRNA-seq and bulk WGS. For patient P22, ADC cells and other mixed-lineage cancer cells also had the same CNV patterns (Fig. 3g and Additional file 1: Fig. S5c). As a result, we further confirmed that the mixed-lineage cancer cells and single-lineage cancer cells in the same patient had common tumor ancestors. Interestingly, we also found that a number of individual cells from adjacent normal tissues had copy number losses in chromosomes 4q and 8p (Additional file 1: Fig. S5a, c, d).

Intratumor heterogeneity contributes to clinical therapy failure and tumor progression [43]. Our scRNA-seq results highlighted the molecular heterogeneity and diversity in an individual tumor of NSCLC containing mixed-lineage cancer cells. These mixed cancer cells may have an aberrant cellular differentiation program or be associated with the transformation between the different subtypes. To analyze the relationship between mixed-lineage cancer cells and prognosis, we partitioned bulk RNA-seq samples of NSCLC from The Cancer Genome Atlas (TCGA) into two clusters based on the calculated lineage-mixing score. Patients with high lineage-mixing features (high score) were correlated with decreased survival, which indicates that a higher percentage of cells with mixed-lineage features in NSCLC predict poorer prognosis (Fig. 3h). Since data of TCGA were generated on bulk samples, it is possible that instead of expressing multi-cancer subtype markers in the same individual cell, marker genes of different cancer subtypes were sepatately expressed in different subpopulations of cancer cells in the tumor tissue. But the score can still reflect the general mixed trend of different cancer subtypes in the tumor tissues. To further investigate the potential connections between mixed-lineage cancer cells and EGFR mutation, we partitioned TCGA samples into two clusters based on EGFR mutation, EGFRWT and EGFRMut. We found that for ADC, the patients with EGFR mutation have lower mixed-lineage features than those with wild-type EGFR. In contrast, for SCC, the patients with EGFR mutation have higher mixed-lineage features than the patients with wild-type EGFR (Fig. 3i). Combined with the anticancer therapeutic effect of drugs against EGFR mutation, the ADC patients with mutated EGFR have better responses to epidermal growth factor receptor tyrosine kinase inhibitors (EGFR-TKIs) than the SCC patients. We speculate that the heterogeneity of mixed-lineage cancer cells may underlie variants in therapeutic responses to the EGFR-TKI target therapy between ADC and SCC patients.

Transcriptome dynamics analysis reveals gene regulation changes during tumorigenesis

To determine the transcriptional signatures of tumor cells, DEGs were identified between all normal epithelial cells and tumor epithelial cells. Eighty downregulated genes and 225 upregulated genes (log2fold change (tumor versus normal) > 2) were identified in tumor epithelial cells. Then, we constructed a pseudotime trajectory to uncover the transcriptome changes during tumorigenesis from the normal epithelium (left) to carcinoma (right) (Additional file 1: Fig. S6a). The genes whose expression levels changed along the trajectory were grouped into 4 distinct clusters based on the dynamic expression patterns (Additional file 2: Table S2). To better understand the biological significance of each cluster of genes, we performed a Gene Ontology (GO) analysis (Additional file 1: Fig. S6b). Genes in cluster 1 were dramatically downregulated in cancer epithelial cells. GO analysis showed that these genes were mainly involved in the regulation of defense (response to bacterium and antimicrobial humoral response), homeostasis, and activation of innate immunity by increasing cell chemotaxis and cytokine levels. Genes in clusters 2 and 3 were widely upregulated in tumor epithelial cells. Genes in cluster 2 were mainly enriched in the immune response to the virus and interferon signaling pathway, and genes in cluster 3 were enriched for the following GO terms: extracellular matrix organization, response to toxic substances, and P53 signaling pathway. Genes in cluster 4, which are mainly highly expressed in late-stage tumors, were strongly enriched in cell proliferation-related terms, including cell division, cell cycle, and mitosis, thus indicating why late-stage tumor cells are more likely to proliferate and metastasize.

AKR1B1 is necessary for tumor cell growth

To further investigate the molecular signatures that were involved in mixed-lineage features of tumor cells, we identified the DEGs among normal epithelial cells, ADC-based single-lineage tumor cells, and combined mixed-lineage tumor cells (Fig. 4a; Additional file 2: Table S3). We found that mesenchymal-related markers such as FN1, TGFBI, and COL1A1 were enriched in mixed-lineage tumor cells. It is known that tumor transformation may occur via EMT, during which process epithelial cells acquire mesenchymal-related features [44, 45]. In addition, epithelial cell differentiation regulation-related genes, such as AKR1B1, SPRR1B, and keratin genes KRT6A, KRT19, and KRT17 were also highly expressed in mixed-lineage tumor cells. Specially AKR1B1 displayed high expression in all four mixed-lineage tumor subtypes in tumor tissues (Fig. 4b). AKR1B1 is involved in the glucose transforming polyol pathway and has been reported to have the capacity to facilitate breast cancer tumorigenesis and metastasis via EMT process [46]. In addition, a study also showed that the expression of AKR1B1 was strongly correlated with EMT in lung cancer and colon cancer, during which process that tumor cells could obtain properties of cancer stem cells manifesting diverse plasticity [47]. Therefore, we speculated that AKR1B1 may be one of the master regulators of mixed-lineage tumor cells' plasticity. To further understand the underlying molecular mechanisms of the tumorigenesis of mixed-lineage tumor cells, we focused on functional analysis of AKR1B1. First, we performed siRNA knockdown experiments with two different siRNAs of AKR1B1 in the H2009 cell line (Additional file 1: Fig. S7a). Knockdown of AKR1B1 can significantly decrease the proliferation of H2009 cells compared with the non-targeting siRNA control (NC) (Fig. 4c). Moreover, knockdown of AKR1B1 strongly reduced the proportion of cells in the S phase and G2/M phase of the cell cycle (Fig. 4d).

Fig. 4

AKR1B1 is essential for the proliferation of lung tumor cells. a Heatmap showing differentially expressed genes among normal epithelial cells, ADC-based single-lineage tumor cells, and combined four mixed-lineage tumor cells. b Box plot in the left shows the single-cell gene expression level of AKR1B1 in normal epithelial cells and five cancer cell subtypes. Box plot in right shows the single-cell gene expression levels of AKR1B1 in epithelial cells from normal tissues, tumor tissues, and LyM tissues. LyM represents lymph node metastasis. c Proliferation analysis of H2009 cells after AKR1B1 was knockdown with two different siRNAs. Compared with non-targeting control (NC), siAKR1B1-1 and siAKR1B1-2 significantly reduced the proliferation of H2009 cells. *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001. p values were determined by t-test. d Cell cycle analysis after H2009 cells were transfected with two different siRNAs of AKR1B1 for 48 h. Compared with NC, both siAKR1B1-1 and siAKR1B1-2 significantly decreased the cell fractions of S and G2/M phases of the cell cycle. *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001. p values were determined by t-test. e Cell apoptosis analysis after H2009 cells were treated with DMSO or 100 μM epalrestat for 48 h. Epalrestat treatment significantly promoted the apoptosis of H2009 cells. ****p < 0.001. p values were determined by t-test. f Photograph of tumors treated with sterile water or epalrestat 36 days after injection. These 11 tumors in the control group were derived from 7 mice, and these tumors in the treatment group were derived from 6 mice. Scale bar, 10 mm. g Quantitation of tumor volumes. The tumor volumes were calculated by the following formula: volume = length × (width)2 × 0.5. The maximum and minimum detected in each time point were removed. ***p < 0.001; ****p < 0.0001. p values were determined by t-test. h Quantitation of tumor weight from tumors in f. Epalrestat treatment significantly inhibited tumor growth. ***p < 0.001. p values were determined by t-test

Next, we treated H2009 cells with 100 μM epalrestat, a specific inhibitor of AKR1B1 that has been approved for the treatment of diabetes complications [48]. As demonstrated in Additional file 1: Fig. S7b, we observed that the cell numbers in the epalrestat treatment group also decreased compared to those in the control group (treated with DMSO). To further exclude the potential side effects caused by the high concentration of epalrestat and to test the drug specificity, the colon cancer cell line SW480, which essentially does not express AKR1B1 according to the Cancer Cell Line Encyclopedia (CCLE) RNA-seq data (Additional file 1: Fig. S7c), was treated with the same concentration of epalrestat. After H2009 cells and SW480 cells were treated with 100 μM epalrestat in parallel for 48 h, the percentage of apoptotic cells was significantly higher in H2009 cells but not in SW480 cells, suggesting that epalrestat can specifically promote apoptosis by inhibiting AKR1B1 (Fig. 4e; Additional file 1: Fig. S7d). There is a minor possibility that the different tissue origins of H2009 cells and SW480 cells may also cause the differences in response to epalrestat between these two cancer cell lines. In addition, we performed bulk RNA-seq to analyze the gene expression in epalrestat treated cells and control (DMSO only) cells. We identified 1540 upregulated genes and 1737 downregulated genes (fold change (epalrestat versus DMSO) > 1.5, P value < 0.01) in epalrestat treated H2009 cells (Additional file 4: Table S4). As epalrestat is a non-competitive inhibitor of aldolase reductase, we found that when treating H2009 cells with epalrestat to inhibit the activity of aldolase reductase, the expression of AKR1B1 was slightly upregulated potentially though a negative feedback regulation mechanism. AKR1B1 could be involved in different metabolic and physiological processes and participate in a complex network of signaling pathways, such as inflammation, cell cycle, epithelial to mesenchymal transition, and mTOR pathway [49]. By analyzing the differentially expressed genes, we found that downregulated genes in epalrestat treated cells were mainly enriched in cell-cell adhesion, positive regulation of cell migration, negative regulation of cell differentiation, cell cycle, mitotic, prostaglandin biosynthetic process, and related metabolic processes, consistent with the proposed function of AKR1B1 (Additional file 1: Fig. S7e). Upregulated genes were mainly enriched in the cellular response to extracellular stimulus, response to nutrient levels, and cellular response to glucose starvation-related biological processes. In addition, positive regulation of the apoptotic process was also enriched in these upregulated genes, which was consistent with the observation that epalrestat treatment group samples had a higher proportion of apoptotic cells (Additional file 1: Fig. S7e). To further examine the effect of epalrestat on tumorigenesis in vivo, we injected H2009 cells into NOD-SCID mice. Then, the mice were intragastrically administered epalrestat daily. One month later, the tumors were collected for further analysis. Compared to control mice, epalrestat-treated mice exhibited significantly reduced tumor cell growth in vivo (Fig. 4f–h). In summary, these data demonstrated that AKR1B1 could play important roles in the tumorigenesis of lung cancer. However, as we did not detect the effectiveness of epalrestat on AKR1B1 target for tumor cells of mice, further evaluations on the potential off-target effectiveness of epalrestat in vivo should be carried out in the future.

To further understand the transcriptomic regulation of AKR1B1, we performed bulk RNA-seq to analyze the gene expression in H2009 cells with AKR1B1 knockdown. We identified 143 upregulated genes and 473 downregulated genes (fold change (siAKR1B1 versus NC) > 2, p value < 10−15) in H2009 cells with AKR1B1 knockdown (Additional file 1: Fig. S7f; Additional file 5: Table S5). Downregulated genes were mainly enriched in metabolism-related pathways, consistent with the known function of AKR1B1. In addition, cancer-related pathways, such as the p53 signaling pathway, the HIF-1 signaling pathway, focal adhesion, and cell cycle DNA replication, were also enriched in these downregulated genes, which is in line with the observation that knockdown samples had lower levels of proliferation characteristics (Additional file 1: Fig. S7f). Upregulated genes were mainly enriched in biosynthetic processes and negative regulators of cell proliferation (Additional file 1: Fig. S7g). To further explore the inhibition specificity of epalrestat toward AKR1B1, we compared these differentially expressed genes (DEGs) (fold change (epalrestat versus DMSO) > 1.5, p value < 0.01) of epalrestat-treated cells with those DEGs (fold change (siAKR1B1 versus NC) > 1.5, p value < 0.01) of siAKR1B1-treated cells. We identified 1179 upregulated genes and 1684 downregulated genes in siAKR1B1 treated cells. When we merged the DEGs of these two datasets, we found 84 overlapped upregulated genes and 170 overlapped downregulated genes (Additional file 6: Table S6). The overlapped downregulated genes were involved in many pathways, such as interferon alpha/beta signaling, response to decreased oxygen levels, positive regulation of cell migration, metabolism of carbohydrates, glucose metabolism, regulation of protein serine/threonine kinase activity, cell cycle, and cell population proliferation, which were correlated with the mentioned functions of AKR1B1 (Additional file 7: Table S7). The overlapped upregulated genes were mainly enriched in the metabolism and transport process, such as cellular amide metabolic process, tetrapyrrole metabolic process, long-chain fatty acid transport, and monocarboxylic acid transport. In addition, other GO terms such as cellular component morphogenesis, epithelial cell differentiation, and negative regulation of cell population proliferation were also enriched in these upregulated genes (Additional file 7: Table S7). Although analysis of common DEGs of epalrestat-treated cells and siAKR1B1-treated cells revealed that these overlapped genes enriched GO terms were consistent with the functions of AKR1B1, significant differences were also observed between these two datasets. Considering that our RNA sequencing results of siAKR1B1-treated cells were more reliable with 90% knockdown efficiency of AKR1B1, we speculated that epalrestat may have some off-target effects in the H2009 cell line. It has been reported that epalrestat also has a binding affinity toward AKR1B10, a member of aldo–keto reductase superfamily, and could suppress the enzymatic activity of AKR1B10. Therefore, we speculated that epalrestat may have some off-target effects in the H2009 cell line by targeting other members of aldo–keto reductase superfamily, such as AKR1B10, but the in-depth mechanism of AKR1B1 inhibition mediated by epalrestat in lung cancer need to be further investigated in the future. Collectively, our functional analysis of AKR1B1 showed that AKR1B1 was necessary for tumor growth and elucidated its transcriptomic regulation in NSCLC.