Abstract
Thirdhand smoke (THS) refers to residual tobacco smoking pollutants that can be adsorbed to indoor surfaces and dust and persist for years after active smoking. THS-related chemicals such as N-nitrosonornicotine (NNN) and 4-(methylnitrosamino)-1-(3-pyridyl)-1-butanone (NNK) are tobacco-specific lung carcinogens that involved in lung cancer development and progression. In this study, we computed the differentially expressed genes (DEGs) between THS and paired control samples. THS-related overexpressed genes (OEs) were overlapped with OEs of lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC). Survival analyses of these overlapped genes were performed using LUAD and LUSC data. 6 genes were selected for validation based on their expression levels and prognostic value. Hematological and neurological expressed 1 (HN1) was further selected due to its novelty in LUAD research. The potential roles of HN1 in LUAD were explored in several ways. In summary, HN1 is overexpressed in THS samples and is associated with the prognosis of patients with LUAD. It may promote cancer progression through several pathways and could serve as a potential therapeutic target especially for THS-related LUAD. In-depth mechanistic studies and clinical trials are warranted.
1. Introduction
Thirdhand smoke (THS) refers to residual tobacco smoking pollutants that can adsorbed to indoor surfaces, clothes, and dust and persist for years after active smoking [1]. These deposited chemicals such nicotine could reemit into the air and react with nitrous acid to form carcinogenic tobacco-specific nitrosamines (TSNAs) including N-nitrosonornicotine (NNN) and 4-(methylnitrosamino)-1-(3-pyridyl)-1-butanone (NNK), which are tobacco-specific lung carcinogens that are involved in the development of lung cancer [2]. The presence of THS in indoor environment is widespread, and it is a severely underestimated health hazard. Currently, there are no effective ways to eliminate THS.
A randomized clinical trial identified 389 genes that were differentially expressed in human respiratory epithelium in response to acute THS exposure or clean air for 3 hours. Gene ontology analysis indicated that these genes were enriched in cell stress and survival-related signaling pathways such as respiratory electron transport chain, DNA repair, and activation of cell viability [3]. It showed that human respiratory epithelium could respond rapidly to THS. Importantly, Hang et al. reported that THS could induce lung cancer development and increase lung cancer incidence in mice [4, 5]. However, the association between THS and lung cancer is still poorly understood due to limited numbers of studies to date.
In this study, we explored the potential carcinogenesis roles of THS in human lung cancer by investigating overlapped-overexpressed genes in THS and lung cancer. This research will provide a mechanistic link between THS and human lung cancer, which would pave the way for future further investigations.
2. Results
2.1. THS-Related Genes Are more Associated with Lung Adenocarcinoma (LUAD) Not Lung Squamous Carcinoma (LUSC)
Of 382 THS-related overexpressed genes (OEs), 12 genes are overexpressed both in LUAD and LUSC, 8 genes are overexpressed in LUAD, and 9 genes are overexpressed in LUSC. Of 7 THS-related downexpressed genes (DEs), one gene is downexpressed both in LUAD and LUSC while another gene is only downexpressed in LUSC (Figure 1, upper left). As shown in the heat map (Figure 1, upper right), 21 genes are associated with poor overall survival (OS) of LUAD. Even genes are overexpressed in LUSC and genes are exclusively overexpressed in LUSC based on filter of , and only KRT8 correlated with poor OS of LUSC while MFAP2 associated with better OS. Detailed survival analysis results of 31 common DEGs in THS and lung cancer are shown in Table 1. This is very interesting since THS-related DEGs are mainly related to prognosis of LUAD.
Gene expression levels are very important for their biological functions, the 21 OS-related genes were further filtered based on their expression in LUAD, and 6 genes with expression were obtained (middle-left). These 6 genes might play critical roles in THS-associated LUAD. Literature mining confirmed that all these 6 genes are associated with cancer. KRT8 and SNRBP are overexpressed and correlated with poor prognosis of NSCLC; ROMO1, GSTP1, and ALDOA are oncogenic in NSCLC. HN1 can promote cancer progression in other cancer types, but there is no report on its association with NSCLC.
2.2. Validation of the Prognostic Value of the 6 Genes in LUAD and Other Cancer Types
We validated the prognostic value of the 6 genes using a merged LUAD dataset from KmPlot. Survival analyses results confirmed that these 6 genes are all associated with poor prognosis of patients with LUAD (Figure 2, red line indicates high expression; black line indicates low expression). Specifically, HN1 (, logrank ); KRT8 (, logrank ); ROMO1 (, logrank ); GSTP1 (, logrank ); ALDOA (, logrank ); SNRPB (, logrank ).
We also explored the expression and prognostic value of these 6 genes in other types of cancer using GEPIA 2.0 webtool. As shown in Supplementary Figure 1, 6 genes were overexpressed in most cancer types, such as pancreatic cancer and colorectal cancer (heat map of log2 transformed fold changes, red box: high expression in cancer tissue against normal control, ; blue box: low expression in cancer tissue, ). Supplementary Figure 2 is the hazard ratio heat map of 6 genes in 29 cancer types. The upper panel is overall survival (OS), and the lower panel is recurrent free survival (RFS). Red box indicates association with poor prognosis while blue box indicates association with better prognosis (). As we can see in Supplementary Figure 2, these 6 genes are associated with poor prognosis of several cancer types. These results demonstrated that THS-related genes not only important in LUAD but may also play critical roles in many other cancer types.
2.3. Expression and Function Analyses of HN1 in LUAD
Literature mining showed that HN1 is the only one in the 6 genes that has not been reported its association with lung cancer. Thus, we selected HN1 to explore its potential function and oncogenic roles in LUAD. Further analyses show that HN1 expressions in different stages of LUAD were all significantly higher than in normal control (Supplementary Figure 3), and its high expression was associated with poor prognosis of patients with LUAD (Supplementary Figure 4). GSEA analyses were performed using DEGs between HN1 high and low groups. Result shows that HN1 expression is positively correlated with embryonic stem cell signatures, metastasis, invasiveness, and hypoxia pathways (Figure 3, , , , and , respectively). Pathway enrichment analysis using clusterProfiler shows that HN1-related DEGs were enriched in cell adhesion molecules (CAMs), Hedgehog signaling pathway, proteoglycans in cancer, cellular senescence, PI3K-Akt signaling pathway, etc. (Supplementary Table 1). Visualization of enrichment pathway network is presented in Figure 4.
Figure 5 shows the correlation between HN1 and RPPA protein abundances. Figure 5(a) is the box plot of spearman computed by using LUAD-HN1 expression data and RPPA data. HN1 is positively correlated with cell cycle and metastasis-related proteins such as CCNB1, CCNE1, FN1, and FOXM1, while it is negatively correlated with PI3K-AKT pathway-related proteins such as AKT pS473 and PRAS40 pT246. Figure 5(b) is the correlation map of HN1 and RPPA proteins (blue represents positive correlation while red represents negative correlation, correlations with value < 0.05 were presented in the map). These protein level analysis results are in accordance with the above pathway enrichment analyses.
(a)
(b)
2.4. Association Analyses of HN1 and Cancer Immunotherapy-Related Factors
The possible association between HN1 and cancer immunotherapy-related factors such as immune checkpoint blockade (ICB) therapy efficacy-related genes and abundance of cancer microenvironment cells were also investigated. Figure 6 shows the correlation between HN1 and ICB efficacy-related genes (Figure 6(a) is box plot of spearman value; Figure 6(b) is correlation map). Figure 6(a) shows that HN1 expression is positively correlated with LDHA, LDHB, TNFSF9, TNFRSF18, IFNG, etc., while negatively correlated with JAK1, JAK2, CD40LG, etc., indicating the possibility that HN1 might serve as a prognostic biomarker or potential therapeutic target for cancer immunotherapy. Figure 6(b) is the correlation map of HN1 and ICB efficacy-related molecules. It shows the correlation clusters of these molecules in LUAD.
(a)
(b)
The correlation between HN1 and cancer microenvironment cells is presented in Figure 7 (Figure 7(a) is box plot of spearman value; Figure 7(b) is correlation map). As we can see from Figure 7(a), HN1 is most correlated with Th1 cells, Th2 cells, CD8 naive T cells, etc., while negatively correlated with CD4 naive T cells, fibroblasts, mast cells, etc. Figure 7(b) is the correlation map of HN1 expression and abundance of stroma cells. It shows some intensity modules or clusters of different cell types. Figure 8 is the box plot of cancer infiltrating immune cells in HN1 high and low groups. The abundance of CD4 T cells, iTreg cells, B cells, etc. were lower in HN1 high samples, while Th1 cells, nTreg cells, etc. were higher in HN1 high samples. Correlation map and box plot of spearman are shown in Supplementary Figure 5. These results indicated that HN1 might play important roles in cancer progression by interacting with ICB-related factors and regulating biological behaviors of cancer microenvironment cells.
(a)
(b)
3. Discussion
Lung cancer is the leading cause of cancer death and the second most commonly diagnosed cancer, with an estimated 2.2 million new cancer cases and 1.8 million deaths worldwide in 2020 [6]. Big progresses have been achieved in clinical practices in recent years. KRAS and EGFR-targeted drugs were applied for lung cancer patients with specific mutations [7, 8]. Immune checkpoint blockade (ICB) therapies are also promising in the management of lung cancer [9]. However, therapeutic resistance emerges rapidly after treatment posing formidable obstacles to cancer therapeutics. Clarify the molecular mechanism lung cancer carcinogenesis and searching for novel prognostic biomarkers and therapeutic targets are of crucial importance.
THS is the residual of tobacco smoke that remains in the environments after active smoking. Several reports show that THS is a great public health hazard in indoor environment [10, 11]. For instance, THS could induce damage in human DNA and stimulate high levels of inflammatory cytokines and may involve in chronic obstructive pulmonary disease and asthma [5, 12]. Previously, Hang et al. reported that early exposure to THS was associated with increased lung cancer incidences in mouse model [4].
High-dimensional biological data resource is very important in screening new biomarkers and potential therapeutic targets [13]. In this study, the potential molecular mechanisms underlying the association between THS and LUAD were also explored using public available databases. Interestingly, we found THS-related DEGs were mainly associated the prognosis of patients with LUAD not LUSC. These results were in accordance with previous reports that THS was associated with increased lung adenocarcinoma incidences in mouse model.
Based on analysis results and literature mining, we choose HN1 for further investigation. HN1 was found to be overexpressed in LUAD and could serve as a prognostic biomarker for poor prognosis of patients with LUAD. Further analysis indicated that HN1 may promote LUAD progression by modulating tumor microenvironment and immune-related pathways. HN1, also named as JPT1, plays critical roles in the regulation of cell cycle and cell adhesion [14] and could negatively modulate AKT-mediated GSK3B signaling [15]. It has been reported that HN1 can promote tumorigenesis and metastasis in several cancer types such as breast cancer [16], prostate cancer [14], liver cancer [17], cervical cancer [18], and thyroid cancer [19]. However, this is the first report to show the association between HN1 and LUAD.
In summary, this is the first report using human data to show that THS-related transcriptional responses of respiratory epithelium are mainly associated with the development and progression of LUAD not LUSC. Exposure to THS is a significant health threat for nonsmokers, especially for children. The public should pay more attention to the potential risks of THS exposure, and potential buyers or renters of houses should be notified if there were THS risks in order to avoid unnecessary exposure to a potent lung carcinogen. Policy makers should also take THS into consideration when developing future environmental and health policies.
4. Materials and Methods
4.1. Ethics Statement
All the data were obtained from public datasets. The Research Ethics Committee of Zhejiang People’s Hospital waived the requirement for ethical approval.
4.2. Data Sources
Gene expression data and protein data of LUAD and LUSC were downloaded from The Cancer Genome Atlas (TCGA: http://cancergenome.nih.gov/). THS-related gene expression data were obtained from Gene Expression Omnibus (GEO, accession no. GSE129959) [3, 20]. This dataset contains 4 participants that receive the clean air exposure first and THS exposure second. Abundance data of stroma and immune cells were obtained from xCell [21] and ImmuCellAI [22] databases.
4.3. Bioinformatics and Statistical Analyses
Survival analyses of 6 genes were analyzed online and downloaded from KmPlot [23] using a merged LUAD dataset. Differential expression and survival analyses of LUAD and LUSC TCGA data was performed using GEPIA version 2.0 [24]. UALCAN webtool was used for HN1 expression analysis in different stages of LUAD and normal control [25]. R 3.4.4 (R Foundation for statistical computing [http://www.r-project.org/]) was used for computing differentially expressed genes. Venn diagram analysis was performed using an online tool (http://bioinformatics.psb.ugent.be/webtools/Venn/). Pathway enrichment visualization was performed using NetworkAnalyst version 3.0 [26]. Gene set enrichment analysis (GSEA) was performed using clusterProfiler version 4.0 [27] and GSEA v4.0.3 (https://www.gsea-msigdb.org/gsea/downloads.jsp). Correlation analysis and visualization were performed using “corrplot” package. Other statistical analyses were performed using GraphPad Prism 5.01 (GraphPad Software, Inc. (http://www.graphpad.com)). Adjust value was corrected for multiple comparisons using the Benjamini and Hochberg’s false discovery rate [28]. Statistical significance was defined as a value < 0.05.
Data Availability
Data used in this study are all from publicly available datasets, and all data sources are noted in the manuscript.
Conflicts of Interest
The authors declare no competing financial interests.
Authors’ Contributions
Jiamin Peng and Jia Shuopeng contribute equally to this work.
Supplementary Materials
Supplementary Figure 1: expression of HN1, KRT8, ROMO1, GSTP1, ALDOA, and SNRPB in other types of tumors in comparison to normal control. Red: high expression; blue: low expression. Supplementary Figure 2: prognostic value of HN1, KRT8, ROMO1, GSTP1, ALDOA, and SNRPB in other types of tumors. Red: good prognosis; blue: poor prognosis. Supplementary Figure 3: HN1 expression in different stages of LUAD is significantly higher in comparison to normal control. Supplementary Figure 4: high expression of HN1 is associated with poor prognosis of patients with LUAD. Supplementary Figure 5: correlation between HN1 and the abundance of immune cells from ImmuCellAI database. (A) Box plot of spearman value from correlation analysis between HN1 and abundance of immune cells. (B) Correlation heat map of HN1 and abundance of immune cells. Supplementary Table 1: GSEA pathway enrichment of HN1-correlated genes. (Supplementary Materials)