To date, 23 PATs, also known as ZDHHC proteins, expressed in mammalian cells have been discovered. Given that the function of ZDHHC in LUAD tumorigenesis is variable, we reasoned that finding suitable ZDHHC is important for studying the modulatory mechanisms of palmitoylation. To identify potential ZDHHC targets in LUAD, we analyzed the RNA-seq data from TCGA database to compare differentially expressed genes between LUAD cases and normal lung tissues. As shown in Figure 1A, 12 ZDHHCs were significantly upregulated in LUAD, and 7 ZDHHCs were significantly downregulated (fold change > 1, P < 0.05). We further interrogated MS-based proteomic data from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) Confirmatory/Discovery cohorts to investigate ZDHHC expression at the protein level, and identified 9 upregulated and 2 downregulated (fold change > 1, P < 0.05) ZDHHCs (Figure 1B). Heat maps (Figure 1C) showed the hazard ratio (HR) of 23 ZDHHCs in LUAD based on overall survival (OS) and disease-free survival (DFS) results. We found that high expression of ZDHHC5 and ZDHHC15 was associated with an unfavorable prognosis and a favorable outcome, respectively. To our surprise, ZDHHC5 was the only ZDHHC with differential expression and prognostic value in LUAD (Figure 1D).

We compared the transcriptional level of ZDHHC5 in LUAD with that in normal samples using 3 independent studies from the Gene Expression Omnibus (GEO) database. In Figure 2A, the mRNA expression of ZDHHC5 was also significantly upregulated (fold change > 1, P < 0.05). The mass spectrometry results of ZDHHC5 from the CPTAC database (Figure 2B) revealed the over expression of ZDHHC5 protein in LUAD. There are 5 known immune subtypes within LUAD tumors (Thorsson et al., 2018), including C1 (wound healing), C2 (IFN-γ dominant), C3 (inflammatory), C4 (lymphocyte depleted), and C6 (TGF-β dominant); overall survival analysis stratified by immune subtype indicated that C3 had the most favorable prognosis, and C6 and C4 had the poorest outcome. This immunogenomics pipeline was used to characterize these heterogeneous tumors, and we used the resulting data to explore the distribution of ZDHHC5 in immune subtypes. ZDHHC5 was slightly downregulated in C3 and was overexpressed in C4 (Figure 2C), which is in accordance with the proposed malignant role of ZDHHC5 in LUAD. We also explored the Spearman correlation of ZDHHC5 with 3 types of immune regulators and tumor-infiltrating lymphocytes in LUAD, but there was no significant association (Supplementary Figure 1A). TCGA Consortium proposed a transcriptional data-based molecular classification for LUAD, which includes 3 molecular subgroups into the clinical classification of LUAD (Cancer Genome Atlas Research Network, 2014). We performed gene expression analysis of ZDHHC5 in the 3 molecular subtypes (Supplementary Figure 1B), but there was no significant difference (P > 0.05). As shown in Figure 2D, the mRNA expression level of ZDHHC5 was higher in TP53-mutant LUAD cases than in the non-mutant cases.

The 2015 World Health Organization (WHO) classification of lung tumors results from an integrated multi-disciplinary approach involving radiology, molecular biology, and medical oncology (Travis et al., 2015). We analyzed the RNA-seq data of ZDHHC5 in LUAD pathological types (Figure 2E). Immunohistochemical (IHC) results of LUAD and normal lung tissue from the Human Protein Atlas (HPA) database (Uhlén et al., 2015) were obtained to be determined (Figure 2F). Basic annotation parameters were scored for the HPA samples by pathologists. ZDHHC5 was stained in over 75% of tumor cells, including medium and low expression LUADs. In lung normal tissues, ZDHHC5 was highly stained in macrophages, with a high fraction (>75%) of stained cells. However, ZDHHC5 showed medium-staining intensity in only 25–75% alveolar cells, which was lower than in the LUAD. Different pathological types represent different prognoses: lepidic represent good, acinar, and papillary standing for intermediate, and micropapillary and solid representing worse prognosis. The clinicopathological diagnosis of LUAD reports the percentage of all types in increments of 5%. Thus, we evaluated the histological subtype of the IHC results and reported the percentage of all identifiable patterns in 5% increments. This classification indicated ZDHHC5 protein expression varied within different histological subtypes, which further supported the heterogeneity of the growth patterns among LUAD. However, the IHC data was not sufficient for statistical analysis.

To confirm the significance of ZDHHC5 in the prognosis of LUAD, we performed OS analysis of ZDHHC5 in 17 independent cohorts. The survival analysis results of ZDHHC5 in these 17 study cohorts was moderately heterogeneous (I² = 45%). As shown in Figure 2G, the HR of ZDHHC5 corrected by meta-analysis was 1.19 (P < 0.01), which indicated ZDHHC5 was associated with an unfavorable prognosis in LUAD. We also performed survival analysis of ZDHHC5 in 3 molecular subtypes of LUAD (Figure 2H). In the overall survival, high expression of ZDHHC5 only showed lower survival probability in the PP (Proximal-proliferative, formerly magnoid) subgroup (P < 0.05), whereas the PI (Proximal-inflammatory, formerly squamoid) and TRU (Terminal respiratory unit, formerly bronchioid) subgroup did not (P > 0.05). However, the Kaplan-Meier curves of disease-free survival indicated ZDHHC5 expression was not associated with disease-free survival in the 3 subtypes (P > 0.05).

We performed Pearson’s correlation analysis of ZDHHC5 in LUAD. There were 179 positive associations and 143 negative associations (FDR < 0.01; R > 0.3 for positive, R < −0.3 for negative). Heatmaps (Figure 3A) revealed the top 50 positively and negatively associated genes, respectively. ZDHHC5 maintains GSC self-renewal capacity, promotes GSC tumorigenicity, and is essential for preserving the migratory, invasive, and angiogenic potentials of GSCs (Chen et al., 2017). However, the stemness of ZDHHC5 in LUAD has not been reported yet. Previously, we found a set of 242 LUAD CSC genes based on the mRNAsi index (Zhang et al., 2020). We mapped the CSC genes with ZDHHC5 association results, and INCENP was the only intersecting gene (Figure 3B). The INCENP gene expression was up-regulated in LUAD cases and different among LUAD subtypes; the overall survival analysis results indicated an unfavorable prognosis of INCENP in LUAD (Supplementary Figure 2). It has been reported that ZDHHC5 depletion suppresses cell proliferation and colony formation of some NSCLC cell lines but not of human bronchial epithelial cell (HBEC3) (Tian et al., 2015). We analyzed the RNA-seq data of ZDHHC5 and INCENP in the LUAD cell line of A549 and HBEC3. As shown in Figure 3C, the mRNA levels of INCENP decreased in A549 compared with that in HBEC3, but ZDHHC5 did not. The scatter plots of the Pearson’s association analysis results revealed that the correlation between ZDHHC5 and INCENP in LUAD increased by clinical stages (Figure 3D). Next, we generated a set of INCENP correlated genes in LUAD, and used the Stouffer-based meta-analysis method to compare the associations results of ZDHHC5 and INCENP (Supplementary Table 1). Subsequently, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using the results from the meta-analysis using gene set enrichment analysis (GSEA). To reduce redundancy, enriched gene sets were post-processed by the method of Weighted set cover (Vasaikar et al., 2018). Among the clustered results (Figure 3E), the cell cycle pathway (hsa04110) had the highest normalized enrichment score (NES) among the GSEA results (Figure 3F).

As the Pearson’s correlation revealed only the potential physical association between gene expression levels rather than the direct analysis of the interplay between two proteins, we investigated the potential palmitoylation with the amino acid sequences of the two proteins. INCENP can bind to the oncogenes survivin and borealin through its N-terminal domain (INCENP_N) form a trimer to exert its cancer-promoting function (Klein et al., 2006; Vader et al., 2006). The conserved motif encompassing amino acids 32–44 has been proven to be essential for targeting INCENP to centromeres during mitosis, and the first 26 amino acids in the protein sequence of INCENP has also been reported to facilitate this functional targeting centromeres (Ainsztein et al., 1998). However, effects of the deletion of the INCENP protein structure have not been excluded. CSS-Palm 4.0 is a software for in silico predicting palmitoylation sites based on proteomic analysis, including a fourth-generation Group-based Prediction System (GPS) algorithm (Ren et al., 2008). The latest training data set contains 583 palmitoylation sites from 277 proteins, experimentally proved and reported in scientific literature. Particle Swarm Optimize (PSO) was also integrated to improve the convergence speed and training accuracy. The prediction performance and system robustness of CSS-Palm 4.0 was evaluated by the leave-one-out validation and 4-, 6-, 8-, 10-fold cross-validations. To further uncover the mechanism underlying the association between INCENP and ZDHHC5 expression, we analyzed the S-palmitoylation sites of INCENP using CSS-Palm 4.0. As shown in Figure 4A, INCENP has 7 cysteine sites predicted to be S-palmitoylated. Among them, Cys¹⁵ had the highest prediction score of 23, and the scores of other sites were less than 8. Actually, the first 26 amino acids are a non-conserved sequence. This peculiarity makes it difficult to predict its effect on protein structure. In this regard, we analyzed the Cys¹⁵ site, and found that it was a conserved amino acid site among species, as illustrated in Figure 4B. Taken together, the results indicated that Cys¹⁵ may be the key site for palmitoylation of INCENP, as well as for targeting centromeres.

To confirm the function of INCENP in CSCs, we extracted genes that were directly associated with INCENP in the LUAD stem cell gene set. Then we constructed a protein-protein interaction (PPI) interaction network using STRING, and conducted a Biological Process (Gene Oncology, GO) enrichment analysis (Figure 4C). The well-known function of INCENP is associated with centromeres, which was responsible for many kinetochore functions, including microtubule binding, motor activity, and cell cycle signaling (Ainsztein et al., 1998). The enrichment analysis results also indicated that the role of INCENP in CSCs was mainly attributable to the cell cycle.

Since the proteins of the ZDHHCs family all have the function of catalyzing palmitoylation, and ZDHHC5 is not the only ZDHHC that is highly expressed in LUAD, we asked why only ZDHHC5 was associated with poor prognosis in LUAD. It is widely established that palmitoylated proteins are (almost) always localized at the membrane (except ZDHHC23) and palmitoylation influences protein trafficking between the cytomembrane and the cytosol (or Golgi apparatus). However, over the years, several nuclear proteins have been identified as targets of palmitoylation. Further, INCENP is also a nuclear protein, which suggests a potential role for ZDHHC5 in the control of INCENP in the nucleus. In order to determine the subcellular localization of ZDHHCs, we predicted the ZDHHCs nuclear localization signals (NLSs) specific to the importin αβ pathway using cNLS Mapper (Figure 5A). Briefly, a GUS-GFP reporter protein fused to an NLS with a score of 8, 9, or 10 indicates exclusive localization to the nucleus; scores of 7 or 8 indicate partial localization to the nucleus; a score of 3, 4, or 5 indicates localization to both the nucleus and the cytoplasm; and scores of 1 or 2 localization to the cytoplasm. We identified 11 ZDHHCs predicted to be localized in the nuclear compartment and the cytoplasm (scores between 3 and 5). To verify these in silico results, we obtained immunofluorescence images of 20 ZDHHC proteins in human cell lines from the HPA database (Supplementary Table 2). We found that ZDHHC5, ZDHHC8, ZDHHC12, ZDHHC14, ZDHHC15, ZDHHC16, and ZDHHC23 localized in the nucleus (Supplementary Figure 3). Figure 5B showed that ZDHHC5 was mainly localized to the plasma membrane and nucleoplasm in A549, in A431 (an epidermoid carcinoma cell line), and in U-2 OS (an osteosarcoma cell line). There were 4 ZDHHCs predicted and verified to be located in nucleus, including ZDHHC5, ZDHHC14, ZDHHC16, and ZDHHC23.

Since ZDHHC5 is not the only ZDHHC protein with nuclear localization, we also performed correlation analysis to define the specificity of the effects of ZDHHC5 and INCENP (Figure 5C). The Pearson’s correlation analysis showed that ZDHHC5 was the only ZDHHC protein positively correlated with INCENP gene expression. Although ZDHHC16 also showed nuclear localization, its gene expression was negatively associated with INCENP. Accumulating evidence has indicated that the nuclear translocation of ZDHHC5 plays an important role in the regulation of many biological processes including S-palmitoylation, and its down-regulation is almost certainly involved in LUAD tumorigenesis.

The RNA-seq results and clinical data for LUAD patients (n = 515) and normal samples (n = 59) was obtained from the LUAD dataset of the TCGA database.¹ The GEO database² is an open-access platform, providing gene expression results performed by microarray, including clinical data. The CPTAC is a database, designed for large-scale proteome and genome analysis. Protein expression analysis in this study used mass-spectrometry-based proteomic data from the CPTAC Confirmatory/Discovery cohorts. The HPA (Uhlén et al., 2015) is an open access knowledge resource, with the aim to map all the human proteins in cells, tissues and organs. The HPA consists of six separate parts, and the Tissue Atlas (Uhlén et al., 2005), Pathology Atlas (Uhlen et al., 2017), Cell Type Atlas and Cell Atlas (Thul et al., 2017) were used in this study.

UALCAN³ is a web resource providing open access to TCGA and CPTAC data. We performed mRNA and protein gene expression analysis in UALCAN. RNA-seq results of mRNA expression was normalized as Transcript per million (TPM), and the expression data is log₂ (TPM + 1) transformed. As there was an observed high variability in within-profile expression median or standard deviation across samples, the log2 spectral count ratio values from CPTAC were first normalized to standard deviations from the median value within each sample profile, and then normalized across the sample. Z-values represent standard deviations from the median across samples for LUAD. We used the median to calculate the deviation between LUAD and the normal cases. Significant differences were estimated by Student’s t-test, and a P-value < 0.05 was considered to define differentially expressed genes (DEGs).

The Single Cell Type Atlas, stores expression data of protein-coding genes in single human cell types. Normalized RNA expression (NX) data of lung cell lines was analyzed by the Cell Type Atlas of the HPA and visualized in bar plot.

The survival map and Kaplan-Meier curve were generated from Gene Expression Profiling and Interactive Analyses vision 2 (GEPIA2).⁴ The OS or DFS analysis was based on gene expression and survival time; the median was selected as the threshold for splitting the high-expression and low-expression cohorts (Cut-off = 50%); the HR based on Cox proportional hazards Model were calculated; P-value < 0.05 was considered statistically significant using the log-rank test.

The Lung Cancer Explorer (LCE) (Cai et al., 2019) is an open access web resource, which collects data from 56 lung cancer studies that include over 6700 patients. In this study, LCE was employed to evaluate the prognostic value of ZDHHC5 and INCENP using meta-analysis. Forest plots generated by the LCE meta-analysis module indicate heterogeneity test using the I² statistic, which refers to the percentage of variation across datasets that is due to heterogeneity. All data was obtained from 19 independent LUAD studies from GEO and TCGA, listed as Tomida_2009 (GSE13213), Hou_2010 (GSE19188), Wikerson_2012 (GSE17710), Staaf_2012 (GSE29016), Kuner_2009 (GSE10245), Rousseaux_2013 (GSE30219), Okayama_2012 (GSE31210), Bild_2006 (GSE3141), Girard_N_c (GSE31548), Botling_2013 (GSE37745), Jones_2004 (GSE1037), Sato_2013 (GSE41271), Tang_2013 (GSE42127), Der_2014 (GSE50081), Schabath_2016 (GSE72094), Zhu_2010 (GSE14814), Takeuchi_2006 (GSE11969), Shedden_2008 (Consortium) (Shedden et al., 2008), and TCGA_LUAD_2016. In this study, I² = 45% (ZDHHC5) and 54% (INCENP), which indicates intermediate different heterogeneity of LUAD datasets.

The Tumor–Immune System Interactions Data Base (TISIDB) is an integrated repository portal for tumor and immune system interaction (Ru et al., 2019). Distribution of gene expression across immune subtypes (Thorsson et al., 2018) using TCGA data was also analyzed in TISIDB.

The LinkFinder module of LinkedOmics⁵ (Vasaikar et al., 2018) was used to search and visualize genes associated with ZDHHC5 and INCENP mRNA expression in LUAD cases from the TCGA database (515 patients). The results were analyzed using Pearson’s correlation analysis. R represents the Pearson correlation coefficient. Genes with the result of | R | > 0.3 and FDR < 0.05 were considered statistically significant.

The HPA contains images of histological sections from normal (Tissue Atlas) and cancer tissues (Human Pathology Atlas) obtained by IHC. Since specimens are derived from surgical material, normal was defined as non-neoplastic and morphologically normal. Tissue microarrays were used to show ZDHHC5 antibody (Atlas Antibodies Cat#HPA014670, RRID:AB_2257442) staining in samples from 3 individuals defined as normal lung tissue, and samples from 6 cancer patients corresponding to 12 LUAD samples. Each sample was represented by 1-mm tissue cores. Basic annotation parameters included an evaluation of (i) staining intensity (negative, weak, moderate, or strong), (ii) fraction of stained cells (<25, 25–75, or > 75%), and iii) subcellular localization (nuclear and/or cytoplasmic/membranous). The histological subtype was classified by a clinical pathologist in our group, and semiquantitative estimation of the different histological patterns present in 5% increments according to the 2015 WHO Classification of Lung Tumors (Travis et al., 2015). LUAD is classified according to the most representative pattern in the cross-sectional area of the tissue section (the so-called main pattern), and pathological types are introduced according to the main pattern (Kuhn et al., 2018).

The Cell Atlas of the HPA database displays high-resolution, multicolor images of proteins labeled by indirect immunofluorescence. The standard immunostaining protocol for immunofluorescence can be found in the open access repository for science methods at protocols. The cells were stained with the HPA antibody of ZDHHCs protein (green) and DAPI for the nucleus (blue).

The prediction of S-palmitoylation sites in INCENP was performed by GPS-Palm (CSS-Palm 4.0), a deep learning-based graphic presentation system for the prediction of S-palmitoylation sites in proteins. We used InterPro to predict domains and important sites of INCENP.

We predicted the ZDHHCs nuclear localization signals (NLSs) specific to the importin αβ pathway by cNLS Mapper, based on the amino acid sequence of each protein. cNLS Mapper extracted the putative NLS sequences with a score equal to or more than the selected cut-off score. Higher scores indicated stronger NLS activities.

The PPI network was retrieved from STRING v.11.0.⁶ The minimum required interaction score was set at medium confidence (0.400). Active interaction sources included text-mining, experiments, databases, co-expression, neighborhood, gene fusion, and co-occurrence. Only individual networks with 10 or more nodes were included for further analysis. Disconnected nodes were hidden in the network. The edges indicated both functional and physical protein associations. Each node represented an input protein, and the line thickness indicated the strength of data support.

We performed Gene Oncology (GO) functional enrichment with genes in each PPI network using STRING. FDR indicated P-values corrected for multiple testing within each category using the Benjamini-Hochberg procedure. FDR was log₁₀ transformed to describe the significance of the enrichment and visualized in a bar-plots according to the BP, biological process; MF, molecular function; and CC, cellular component.

The Link-Interpreter module of LinkedOmics was used to perform enrichment analysis of the meta-analysis results of ZDHHC5 and INCENP associated genes. In the Link-Interpreter module, the meta data was signed and ranked by the meta P-value, and Gene Set Enrichment Analysis (GSEA) (Subramanian et al., 2005) was used to generate analyses based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. The minimum number within per gene size was set as 5, and 1,000 simulations were performed.

The mRNA expression-based stemness index (mRNAsi index) (Malta et al., 2018) was used to assess the degree of oncogenic dedifferentiation, which was calculated using an innovative OCLR machine-learning algorithm (Sokolov et al., 2016). Higher mRNAsi scores were associated with malignant biological processes in CSCs and greater tumor dedifferentiation, according to the histopathological grades. In our previous study (Ruiz-Cordero and Devine, 2020), we calculated a set of LUAD stemness genes based on the mRNAsi index. Using WGCNA, stemness genes were clustered into modules. Key genes were screened from the blue module based on the gene significance for mRNA stemness indices (GS.mRNAsi), gene significance for the epigenetically regulated-mRNAsi (GS EREG-mRNAsi) and module membership (MM) scores. In this study, we used genes with the GS mRNAsi/GS EREG-mRNAsi score > 0.5 and P < 0.05 as stem genes for further analysis.