HSD17B4, ACAA1, and PXMP4 in Peroxisome Pathway Are Down-Regulated and Have Clinical Significance in Non-small Cell Lung Cancer

To explore the potential functions and clinical significances of peroxisomes during lung cancer development and progression, we investigated the expressional profiles of peroxisome pathway genes and their correlations with clinical features in non-small cell lung cancer (NSCLC). The RNA-seq data of NSCLC including lung squamous carcinoma (LUSC) and lung adenocarcinoma (LUAD) patients with their clinical information were downloaded from The Cancer Genome Atlas (TCGA). Gene expression comparisons between tumor and normal samples were performed with edgeR package in R software and the results of the 83 peroxisome pathway genes were extracted. Through Venn diagram analysis, 38 common differentially expressed peroxisome pathway genes (C-DEPGs) in NSCLC were identified. Principal components analysis (PCA) was performed and the 38 C-DEPGs could discriminate NSCLC tumors from the non-tumor controls well. Through Kaplan-Meier survival and Cox regression analyses, 11 of the C-DEPGs were shown to have prognostic effects on NSCLC overall survival (OS) and were considered as key C-DEPGs (K-DEPGs). Through Oncomine, Human Protein Atlas (HPA) and the Clinical Proteomic Tumor Analysis Consortium (CPTAC), three K-DEPGs (HSD17B4, ACAA1, and PXMP4) were confirmed to be down-regulated in NSCLC at both mRNA and protein level. Their dy-regulation mechanisms were revealed through their correlations with their copy number variations and methylation status. Their potential functions in NSCLC were explored through their NSCLC-specific co-expression network analysis, their correlations with immune infiltrations, immunomodulator gene expressions, MKI67 expression and their associations with anti-cancer drug sensitivity. Our findings suggested that HSD17B4, ACAA1, and PXMP4 might be new markers for NSCLC diagnosis and prognosis and might provide new clues for NSCLC treatment.


INTRODUCTION
Lung cancer is the most frequently diagnosed cancer and the leading cause of cancer death worldwide, with about 2.1 million new lung cancer cases and 1.8 million lung cancer deaths every year (Bray et al., 2018). Almost 80-85% of the lung cancer cases were non-small cell lung cancer (NSCLC) which included two major histological types: lung squamous carcinoma (LUSC) and lung adenocarcinoma (LUAD) (Travis et al., 2015). For the NSCLC patients at an early stage (stage I-II), surgery is the recommended treatment and the 5-year survival is about 53-92% (Vansteenkiste et al., 2014;Goldstraw et al., 2016). However, for the patients at a locally advanced or metastatic stage (stage III or IV), the operation is not amendable and the 5-year-survival is only 4-20% (Siegel et al., 2017;Osmani et al., 2018). To improve the lung cancer outcomes, it is crucial to find new markers for its early diagnosis and prognostic predication.
Peroxisomes are ubiquitous cellular organelles which can be found in nearly all eukaryotic cells. They were first described by Johannes Rhodin in 1954 and termed as "microbodies." As there were a lot of hydrogen peroxide metabolizing enzymes in the microbodies, they were then called "peroxisomes" functionally (De Duve and Baudhuin, 1966). In fact, besides involving in the synthesis and turnover of reactive oxygen species (ROS) (Bonekamp et al., 2009;Antonenkov et al., 2010), the peroxisomes also play important roles in fatty acid Aand/or β-oxidation (Poirier et al., 2006;Wanders and Waterham, 2006), the catabolism of purines (Sprecher et al., 1995) and the biosynthesis of glycerolipids (Wanders et al., 2010) which are indispensable to human health and development. Metabolic reprogramming is one of the core traits of the cancer cells, due to the multiple alterations during the multi-step process of tumorigenesis (Soga, 2013;El Hassouni et al., 2019). In recent years, with the deepening understanding of the important functions of the peroxisomes, more and more attention were paid to the peroxisomes research (Tanner et al., 2013;El Hassouni et al., 2019) and they were demonstrated to be implicated in innate immunity (Ferreira et al., 2019), signal transduction (Schrul and Schliebs, 2018), aging (Deori et al., 2018), and cancer (Dahabieh et al., 2018).
In fact, dy-regulations of peroxisomal enzymes and/or their effects were shown in numerous tumor types including prostate cancer (Ko et al., 2018), colorectal carcinoma (Lakis et al., 2010), liver cancer (Lu et al., 2015;Chen et al., 2018), oral squamous cell carcinoma (Lai et al., 2018), pancreatic cancer (Deplanque et al., 2015), breast cancer (Keller et al., 1993), and lymphoma (Zheng et al., 2019). However, there was no systemic study of peroxisomes in lung cancer. Here, we identified 38 common differentially expressed genes in peroxisome pathway in NSCLC from The Cancer Genome Atlas (TCGA) project and evaluated their prognostic effects in LUSC and LUAD, respectively. Interestingly, 11 of the genes were shown to have prognostic effects but only in LUSC or LUAD individually. However, three of them were confirmed to be consistently under-expressed in LUSC and LUAD at both mRNA and protein level. We then constructed the NSCLCspecific co-expression network of the three genes, explored their potential functions, evaluated their correlations with the immune infiltrations, immunomodulator gene expressions, proliferation marker (MKI67) expression and the sensitivities of anti-cancer drugs in NSCLC. These results might provide new clues for the values of peroxisomes in NSCLC diagnosis, prognosis and treatment.

Available Data From TCGA Database
The RNA-seq data of 501 TCGA-LUSC patients (including 501 primary tumor and 49 normal tissue samples) and 513 TCGA-LUAD patients (including 513 primary tumor and 59 normal tissue samples) with their clinical information were downloaded from Genomic Data Commons (GDC) data portal 1 . The clinical features of the patients were listed in The information of the peroxisome pathway genes (n = 83) was investigated with KEGG database 2 . EdgeR package in R software (R3.5.2) was used for expressional comparisons of the genes between tumor and normal tissues in TCGA-LUSC and TCGA-LUAD datasets and the expressional differences of the peroxisome pathway genes were extracted. The genes with false discovery rate (FDR) < 0.01 were considered to be statistically significant differentially expressed peroxisome pathway genes (DEPGs). The intersection of the two sets of the DEPGs in LUSC and LUAD were considered common DEPGs (C-DEPGs) which were consistently up-or downregulated in the two subtypes. Principal components analysis (PCA), which was applied widely for effective dimension reduction and exploratory visualization, was confirmed to be useful to correct the possibility of false association and show the difference between case and control clearly (Price et al., 2006;Zhang and Castello, 2017). In this study, through GEPIA 3 , PCA was performed to evaluate the discriminating power of the C-EDPGs in differentiating NSCLC from non-tumor lung tissues. In GEPIA, the Genotype-Tissue expression (GTEx) normal data was used to solve the imbalance between the tumor and normal data which can cause inefficiency in various differential analyses and the TCGA and GTEx gene expression data were all Trans Per Million (TPM) normalized from the raw RNA-Seq data by the UCSC Xena project based on a uniform pipeline (Tang et al., 2017).
To evaluate the prognostic effects of the C-DEPGs on overall survival (OS) of the NSCLC patients, with SPSS 18.0, Kaplan-Meier survival analysis with log rank test was performed in LUSC and LUAD, respectively, and the hazard ratios (HRs) were obtained from univariate Cox proportional hazard models. For the above analyses, the median expression of each gene was set as the threshold and the patients were divided into low expression and high expression groups. The genes with significant prognostic effects (p < 0.05) were considered as the key C-DEPGs (K-DEPGs).

Validations of the Expressional Differences of the K-DEPGs in NSCLC
At mRNA level, the K-DEPGs were compared between tumor and normal lung tissues in other LUSC and LUAD datasets via Oncomine database. For the comparisons, the filters were used as follows: analysis type: lung adenocarcinoma vs. normal analysis, squamous cell lung carcinoma vs. normal analysis; data type: mRNA; p-value: 0.05. At protein level, the immunohistochemical (IHC) data of the K-DEPGs in lung cancer and normal lung tissues were investigated from the Human Protein Atlas 2 https://www.kegg.jp/ 3 http://gepia.cancer-pku.cn/ (HPA) database and the results were validated by the lung cancer proteomic data from Clinical Proteomic Tumor Analysis Consortium (CPTAC) 4 which included 111 unique LUAD tumor samples and 102 unique normal samples analyzed by global proteomic mass spectrometry using the 10-plexed isobaric tandem mass tags (TMT-10) following the CPTAC reproducible workflow protocol. The clinical features of the lung cancer dataset from CPTAC were shown in Supplementary Table S1. The expressional differences of the proteins between early stage (stage I) and late stage (stage II/stage III/stage IV) and between well/moderate differentiation (grade 1/grade 2) and poorly differentiation (grade 3) were also evaluated to investigated their associations with NSCLC progression. The non-parametric, two-independent-samples Wilcoxon test was used to evaluate the protein differences and p < 0.05 was considered significant.

Influences of Copy Number Variations and Methylation Values on the Expressions of the Confirmed Genes in NSCLC
To further uncover the potential mechanisms of the dyregulation of the confirmed genes, their correlations with copy number variations (CNVs) in TCGA-LUSC and TCGA-LUAD were analyzed through cBioPortal 5 , a publicly accessible resource providing visualization and analysis tools for more than 5,000 tumor samples from 232 cancer studies in the TCGA pipeline. The correlations between the methylation status of the genes and their expressions were investigated via MEXPRESS 6 , an online tool for visualization of DNA methylation and expression data from TCGA. Spearman's correlation and Pearson's correlation were used for the analyses and the absolute value of correlation coefficient >0.1 with p < 10 −5 was considered significant.

NSCLC Specific Co-expression Network Construction of the Confirmed K-DEPGs
Co-expressed genes were shown to be involved in similar biological processes and functionally related (Stuart et al., 2003). Co-expression network analysis was confirmed to be effective to explore new functions for specific genes (Zhong et al., 2017;Dugo et al., 2018). Here, The NSCLC specific co-expression network of the three confirmed K-DEPGs including HSD17B4, ACAA1, and PXMP4 were constructed through a database for tissue and cancer specific biological networks (TCSBN) 7 . The co-expressed genes in the network were then analyzed through Metascape 8 to explore the potential functions of the three confirmed genes during NSCLC development.

Further Insight of Potential Roles of the Confirmed K-DEPGs in NSCLC Immunoregulation and Tumor Proliferation
The abundances of six immune infiltrates (B cells, CD4 + T cells, CD8 + T cells, Neutrophils, Macrophages and Dendritic cells) in LUSC and LUAD were investigated Via TIMER 9 , a comprehensive resource for systematical analysis of immune infiltrates across diverse cancer types. The immunomodulator genes (n = 91) including 24 immunoinhibitors, 46 immunostimulators, and 21 major histocompatibility complex (MHC) related genes were downloaded from The Cancer Immunome Atlas (TCIA) 10 . To investigate their associations with immune processes, the confirmed C-DEPGs (HSD17B4, ACAA1, and PXMP4) expressions were analyzed for their purity-corrected Spearman correlations with the immune infiltrations and the immunomodulator gene expressions in NSCLC through TIMER.
MKI67, also known as Ki-67 protein, is expressed in all the proliferating cells and widely used as a proliferation marker (Gerdes et al., 1983;Scholzen and Gerdes, 2000;Booth et al., 2014). MKI67 expression profiles in NSCLC were investigated through TIMER at mRNA level and through HPA 9 https://cistrome.shinyapps.io/timer/ 10 https://tcia.at/ and CPTAC at protein level. To evaluate the potential roles of HSD17B4, ACAA1, and PXMP4 in NSCLC proliferation, their correlations with MKI67 expression were also evaluated both at mRNA level and protein level through TIMER and CPTAC, respectively. Wilcoxon test and Spearman correlation analysis was used for expressional comparisons and correlation evaluation, respectively. For all the analyses, p < 0.01 was considered statistically significant.
Further Insights of the Potential Effects of the Confirmed K-DEPGs on Anti-cancer Drug Sensitivity/Resistance in NSCLC Cell Lines The relative expressions of HSD17B4, ACAA1, and HSD17B4 in NSCLC cell lines and the pharmacologic profiles for 24 anticancer drugs across 504 cell lines were downloaded from Cancer Cell Line Encyclopedia (CCLE) database 11 . The half maximal inhibitory concentrations (IC50) of 24 anti-cancer drugs in 89 NSCLC cell lines were extracted. Spearman's correlation analysis was also applied to investigate the associations between the three gene expressions and the IC50 values of the drugs in the cell lines. In addition, for each drug, the NSCLC cell lines were divided into sensitive group and resistant group with the drug IC50 of 8 µM as the threshold according to a previous study (Xiang et al., 2019) and the expressional differences of the three genes between 11 https://portals.broadinstitute.org/ccle/  Frontiers in Genetics | www.frontiersin.org the two groups were evaluated with Mann-Whitney U tests. The analyses were performed through SPSS 18.0 and p < 0.05 was considered significant.

Identification of C-DEPGs Between Tumor and Normal Samples in NSCLC
With the filter FDR < 0.01, there were 57 and 53 DEPGs in LUSC and LUAD, respectively (Figures 1A,B). Among them, 38 genes including 16 down-regulated ( Figure 1C) and 22 up-regulated genes ( Figure 1D) were consistently differentially expressed in LUSC and LUAD and considered as C-DEPGs. Upon the PCA analysis through GEPIA (Figure 1E), the 38 C-DEPGs could discriminate the NSCLC tumors and non-tumor tissues effectively, indicating their diagnostic power in differentiating NSCLC and normal controls.

Prognostic Effects of C-DEPGs
Through Kaplan-Meier survival analysis (Supplementary Table  S2 Table S1) and meta-analysis could not be performed. (E) PXMP4 was significantly down-regulated in LUSC. (F-J) significant overall down-regulation of HSD17B4, ACAA1, CAT, HMGCLL1, and PXMP4 in LUAD. NSCLC, non-small cell lung cancer; LUSC, lung squamous carcinoma; LUAD, lung adenocarcinoma; C-DEPGs, common differentially expressed peroxisome pathway genes. For all the analyses, p < 0.05 was considered significant. was shown to have prognostic effects on LUSC OS while 10 other genes (Figures 2B-K) were indicated to have significant prognostic effects on LUAD OS while no significance was shown for the other C-DEPGs (Supplementary Figures S1, S2). The 11 genes with significant prognostic effects were considered as K-DEPGs and their HRs were shown in Figure 2L. Notably, the down-regulated (ACAA1, CAT, HMGCCL1, and PXMP4) and up-regulated genes (NUCT19, PEX26, PEX5L, PEX6, SLC25A17, and SOD1) were shown to have favorable and unfavorable effects on LUAD OS, respectively. However, for HSD17B4, which was down-regulated in NSCLC, its unfavorable effects on LUSC OS were shown. It was indicated that although the genes were differentially expressed in both LUSC and LUAD, their prognostic effects were different in the two subtypes.

Validation of the 11 K-DEPGs in Other NSCLC Datasets
With the filters, 12 NSCLC datasets in Oncomine database were selected for expressional comparisons of the K-DEPGs between the NSCLC tumors and normal lung tissues. There were six and ten comparisons for LUSC and LUAD datasets, respectively. Except HMGCLL1 and NUDT19 in LUSC (Hou lung), there were two and above comparisons of the K-DEPGs each in LUSC or LUAD datasets (Supplementary Table S3) and meta analyses were applied to obtain the overall differences of the gene expressions. As shown in Figure 3, HSD17B4, ACAA1, CAT, HMGCLL1, and PXMP4 were all confirmed to be downregulated in LUSC (Figures 3A-E) and LUAD (Figures 3F-J), respectively. Among the six up-regulated K-DEPGs, only NUDT19 were confirmed to be significantly increased in both LUSC ( Figure 4A, Hou lung: fold change = 1.535, p < 0.01) and LUAD (Figure 4G) while no significance of PEX26 (Figures 4B,H, p > 0.05) and PEX5L (Figures 4D,J, p > 0.05) was shown in LUSC and LUAD. PEX6 (Figures 4C,I) and SOD1 (Figures 4F,L) were only confirmed to be up-regulated in LUAD while SLC25A17 (Figures 4E,K) was confirmed to be significantly increased in LUSC.
At protein level, for the six genes which were confirmed to be under-expressed (HSD17B4, ACAA1, CAT, HMGCLL1, and PXMP4) or over-expressed (NUDT19) both in LUSC and LUAD in above analyses, there was no IHC data available for HMGCLL1 and NUDT19 in HPA database and their expressional differences could not be obtained. For CAT, with negative staining both in lung cancer tumor cells and normal pneumocytes, no significant expressional difference was shown. However, for the other three genes, the IHC data indicated their down-regulation in NSCLC ( Figure 5). HSD17B4, ACAA1, and PXMP4 were all strong  expressed in lung (pneumocytes) while there was a negative to moderate expression of HSD17B4 and ACAA1 in the cytoplasm or membrane of LUSC and LUAD tumor cells. For PXMP4, although not all, most (9/12) of the NSCLC tumors also presented its negative to moderate expression and its decrease in the tumors was also obvious. In addition, as shown in Figure 6, the under-expression of HSD17B4 (Figure 6A), ACAA1 (Figure 6B), and PXMP4 ( Figure 6C) was also shown in the CPTAC lung cancer dataset, consistent with their expression profiles in above analyses. Furthermore, through their comparisons between different grades (Figures 6D-F), a lower trend of ACAA1 expression (p = 0.084, 0.05 < p < 0.1) and a significant lower expression of PXMP4 (p = 0.02, p < 0.05) was shown in the poor differentiation tumors than the ones with well-moderate differentiation while no significant difference of HSD17B4 was shown. However, none of the three was shown to be statistically differentially expressed between the tumors between early and late stages (Figures 6G-I, p > 0.05).

The Influences of Copy Number Variations and Methylation Values on the Expressions of HSD17B4, ACAA1, and PXMP4
As shown in Figure 7, through cBioPortal, HSD17B4 and ACAA1 expressions were shown to be positively correlated with  their copy number values in both LUSC and LUAD. However, positive correlation between PXMP4 expression and its copy number value was observed in LUSC ( Figure 7E) while not in LUAD ( Figure 7F). Table S4), interestingly, positive correlations were shown between HSD17B4 expression and its DNA methylation status of six CpG sites (cg23314948, cg06903010, cg24537512, cg16261704, cg01229506, and cg13432928) both in LUSC and LUAD. In contrast, significant negative correlations were shown between ACAA1 expression and its methylation value of CpG site cg10548708 in LUSC (Pearson r = −0.234, p = 4.412E-06) and LUAD (Pearson r = −0.411, p = 1.002E-18). For PXMP4, eight CpG sites (cg20588982, cg06231372, cg12297619, cg27361727, cg25092328, cg27194921, cg18669346, and cg24270031) in its promoter region were shown to have significant negative correlations between their methylation values and PXMP4 expression in the LUSC (Pearson r ranging from −0.374 to 0.457, p < 10 −5 ) and LUAD (Pearson r ranging from −0.509 to 0.572, p < 10 −5 ).

Through MEXPRESS analyses (Supplementary
NSCLC Specific Network of HSD17B4, ACAA1, and PXMP4 and Their Potential Functions in Immunoregulation and Tumor Growth Regulation in NSCLC Through TCSBN and Cytoscape 3.6, the NSCLC co-expression network of HSD17B4, ACAA1, and PXMP4 was constructed ( Figure 8A). There were 70 nodes and 1431 edges in the network and they represented the genes and their coexpression scores, respectively. The genes were applied to enrichment analysis and they were significant enriched in 11 terms ( Figure 8B). Interestingly, except peroxisome, several immune-related processes including phagocytosis, regulation of adaptive immune response, regulation of lymphocyte mediated immunity, and ECM affiliated were enriched, consistent with the associations of peroxisomes with immune response and inflammation reported in previous studies (Vijayan et al., 2017;Di Cara et al., 2019), indicating the implication of HSD17B4, ACAA1, and PXMP4 in these processes.
In addition, via TIMER, the purity-corrected Spearman's correlations between the expressions of HSD17B4, ACAA1, and PXMP4 and the immune infiltrations (Table 2) in NSCLC were obvious. Notably, HSD17B4 and ACAA1 were shown to be negatively correlated with LUSC purity while no significant correlation was found between the three gene expressions and LUAD purity ( Table 2). However, both in LUSC and LUAD, HSD17B4 was found to have positive correlations with infiltrations of B cell, CD8 + T cell, Macrophage and dendritic cell. In contrast, ACAA1 was shown to be negatively correlated with CD8 + T cell and neutrophil infiltrations only in LUAD but not in LUSC. For PXMP4, its positive correlations with macrophage in LUSC while its significant correlations with B cell and neutrophil in LUAD were indicated. Consistent with their correlations with the immune infiltrations, the three gene expressions were shown to be significantly correlated with most (in LUSC: 77/91, in LUAD: 77/91, in both LUSC and LUAD: 67/91) of the immunomodulator gene expressions with 116 and 145 significant positive/negative correlations in LUSC ( Figure 9A) and LUAD (Figure 9B), respectively. Although there were more positive correlations than the negative ones, the three gene expressions were positively or negatively correlated with both immunostimulatory and the immunoinhibitory gene expressions, indicating their potential roles in the immune inhibition and/or immune stimulation of NSCLC. Interestingly, most of the significant correlations (LUSC: 36/38, LUAD: 44/47) between the three gene expressions and MHC-related gene expressions were positive NSCLC, non-small cell lung cancer; LUSC, lung squamous carcinoma; LUAD, lung adenocarcinoma. *p < 0.05; **p < 0.01. Spearman's purity-corrected correlation analysis was used for and p < 0.01 was considered significant. Partial cor., partial correlation.
ones. Considering the crucial roles of MHC molecules in presenting antigens to lymphocytes and that (Lenz, 2018), the positive roles of HSD17B4, ACAA1, and PXMP4 in regulation of antigen presentation were also indicated. In addition, among the top 10 correlations (with the lowest p-values) ( Supplementary Table S5), all the 10 correlation in LUSC and 7 of 10 correlations in LUAD were associated with HSD17B4 and all the correlations were positive ones, indicating although all the three genes were associated immune response, the immunoregulatory potential of HSD17B4 in NSCLC was highlighted. Through TIMER, higher MKI67 expression was shown to be in LUSC and LUAD tumors than their normal controls (Supplementary Figures S3A,B, p < 0.01). And MKI67 expression was negatively correlated with expressions of HSD17B4, ACAA1, and PXMP4 in both LUSC and LUAD (Supplementary Figures S3C,D, p < 0.01). At protein level, positive expression of MKI67 in NSCLC tumors while negative expressed in normal lung tissues was shown in HPA (Supplementary Figures S3E-G) and MKI67 overexpression in NSCLC tumors was also confirmed in CPTAC lung cancer dataset (Supplementary Figure S3H, p < 0.01). Although no significant correlation between HSD17B4 and MIKI67 expression (Supplementary Figure S3I, p > 0.05) was shown FIGURE 9 | Significant correlations between HSD17B4, ACAA1, and PXMP4 expressions and immunomodulator gene expressions. (A) There were 114 significant positive or negative correlations between HSD17B4, ACAA1, and PXMP4 expressions and immunomodulator gene expressions in LUSC. (B) There were 145 significant positive or negative correlations between HSD17B4, ACAA1, and PXMP4 expressions and immunomodulator gene expressions in LUAD. The nodes and the edges represented the genes and the correlations, respectively. Purity-corrected Spearman correlation was used and p < 0.01 was considered significant.
FIGURE 10 | Gene expressional differences of ACAA1 and PXMP4 between NSCLC cell lines with different anti-cancer drug sensitivities. (A) For anti-cancer drugs AZD0530, AZD6244, and ZD6474, ACAA1 was higher expressed in drug-sensitive NSCLC cell lines than the drug-resistant ones. (B) For Erlotinib and Lapatinib, there was a higher trend (0.05 < p < 0.1) of ACAA1 expression in the sensitive cell lines than the resistant ones. (C) PXMP4 was higher expressed in the Topotecan-sensitive NSCLC cell lines than the resistant ones. NSCLC, non-small cell lung cancer. Two independent samples Mann-Whitney U test was used for the comparisons and p < 0.05 was considered significant.
in CPTAC lung cancer, the other two proteins (ACAA1 and PXMP4) were confirmed to be negatively correlated with MKI67 expression (Supplementary Figures S3J,K, p < 0.01), indicating their potential roles in regulation of NSCLC growth. Table S6), ACAA1 expression was shown to be negatively correlated with the IC50 values of five anti-cancer drugs including AZD0530 (r = −0.261, p = 0.014), AZD6244 (r = −0.281, p = 0.008), Erlotinib (r = −0.211, p = 0.047), Lapatinib (r = −0.234, p = 0.028), and ZD-6474 (r = −0.285, p = 0.007) while no significant correlation was shown for HSD17B4 or PXMP4 expression with the IC50 values of all the 24 drugs (p > 0.05). The expressional differences of the three genes between the sensitive cell lines and the resistant ones were shown in Table 3, consistent with the negative correlations above, for anti-cancer drugs AZD0530, AZD6244, and ZD-6474, there was a significant higher expression of ACAA1 in the sensitive cell lines than in the resistant ones ( Figure 10A). For Erlotinib and Lapatinib, a higher trend of ACAA1 expression in the sensitive cell lines were shown (Figure 10B, 0.1 > p > 0.05). Interestingly, although no significant correlation was shown between PXMP4 expression and Topotecan IC50 values, there was a significant lower expression of PXMP4 in the Topotecan-sensitive NSCLC cell lines than the resistant ones ( Figure 10C), indicating that the expressional difference of PXMP4 might be due to the differences between the two groups but not due to the drug itself.

DISCUSSION
Reprogramming of metabolic pathways are implicated in the process of tumor development to facilize its unregulated growth and metastatic dissemination (Hanahan and Weinberg, 2011). Recently, the involvement of peroxisome pathway in cancer was demonstrated (Dahabieh et al., 2018). However, although dysfunctions of peroxisomes were shown in many tumors (De Craemer, 1995;Valenca et al., 2015), their roles in lung cancer were rarely explored. Here, we focused on the peroxisomes in NSCLC and identified 38 C-DEPGs in LUSC and LUAD. Among the 38 genes, HSD17B4, ACAA1, and PXMP4 were highlighted for their down-regulation both at mRNA level and protein level, their prognostic effects, their correlations with immune infiltrations, immunomodulator gene expressions, the proliferation marker (MKI67) expression and/or anti-cancer drug sensitivities in NSCLC. HSD17B4 encodes 17β-hydroxysteroid dehydrogenase type 4 (HSD17B4), a 80 kDa multifunctional enzyme localized in peroxisomes, also known as D-specific bifunctional protein (DBP) and multifunctional protein 2 (MFP-2)(de Launoit and Adamski, 1999). In normal cells, HSD17B4 plays important roles in the peroxisomal β-oxidation of long-and branchedchain fatty acids (Baes et al., 2000;Breitling et al., 2001), bile acid biosynthesis (Ferdinandusse et al., 2005), and sterol metabolism (Poutanen et al., 1995). Over the last decade, HSD17B4 has been reported to be involved in the tumorigenesis and progression of many tumors. However, the specific roles of HSD17B4 seems to be inconsistent in different contexts. For example, its overexpression was associated with the tumor cell proliferation in liver cancer (Lu et al., 2015, prostate cancer and colon cancer (Frasinyuk et al., 2017) while exerted tumor suppressive functions in adrenocortical carcinoma (Ding et al., 2017). Here, HSD17B4 was shown to be under-expressed in NSCLC, indicating its diagnostic potential. Consistent with its positive correlation with inflammation in liver cancer tissues (Pan et al., 2018), here, its significant correlations with the immune cell infiltrations and immunomodulator gene expressions were obvious in NSCLC, indicating the important roles of HSD17B4 in immune response regulation. Considering the important roles of immune dysfunction during lung cancer occurrence and progression (Datar and Schalper, 2016;Mascaux et al., 2019), the potential effects of HSD17B4 in immunoregulation might provide new clues for NSCLC immunotherapy. However, its unfavorable prognostic effects were shown in LUSC while not in LUAD, indicating the heterogeneity of its roles in different lung cancer subtypes. Further study is needed for investigation of its specific activities in NSCLC.
ACAA1 gene encodes acetyl-CoA acyltransferase 1 (ACAA1), another important enzyme during peroxisomal β-oxidation of fatty acids (Antonenkov et al., 1997). In previous studies, it was demonstrated to be under-expressed in liver cancer (Yan et al., 2017) and kidney renal clear cell carcinoma (Zhang et al., 2019). Here, similar to HSD17B4, ACAA1 was also found to be down-regulated and negatively correlated with MKI67 expression in NSCLC, indicating its anti-tumor potential. Notably, considering the importance of HSD17B4 and ACAA1 in peroxisomal β-oxidation of very long fatty acids, the dysfunction of the process in NSCLC was indicated, consistent with the anti-tumor function of β-oxidation in lung cancer in a previous study (Srivastava et al., 2014). In addition, as one of the innate immunity genes, the single nucleotide polymorphisms (SNPs) of ACAA1 have been reported to be associated with pathogenesis of childhood asthma (Sharma et al., 2012) and the protective effects of its exposure to endotoxin (Sordillo et al., 2011). In this study, interestingly, besides its negative correlations with CD8 + T cell and neutrophil infiltrations in LUAD, its significant correlations with immunomodulator gene expressions was presented in both LUSC and LUAD, indicating its associations with innate and adaptive immune response. Furthermore, lower ACAA1 expression was shown to be associated the resistance of NSCLC cell lines to five different anti-cancer drugs including AZD0530, AZD6244, ZD-6474, Erlotinib and Lapatinib, providing new clues for NSCLC chemotherapy.
PXMP4, also termed PMP24, a 24-kDa peroxisomal integral membrane protein, was first isolated from rat liver peroxisome membranes (Reguenga et al., 1999). Although it was known to be able to bind chaperone/membrane transporter PEX19 (Sacksteder et al., 2000), its function in peroxisomes was unclear. However, its dysfunction in prostate cancer was reported and its under-expression was associated with the DNA hypermethylation of its CpG island (Wu and Ho, 2004;Zhang et al., 2010). Here, we also found the down-regulation of PXMP4 and the negative correlations between its expression and the methylation values of its CpG sites in LUSC and LUAD. However, PXMP4 expression was shown to be positively correlated with its copy numbers in LUSC while not in LUAD, indicating the difference between the differences of PXMP4 regulation between the two subtypes. Furthermore, its negative correlation with MKI67 expression and LUAD differentiation indicated the tumor-suppressor activities. In addition, besides their prognostic effects in LUSC or LUAD, since there was no expressional difference of HSD17B4, ACAA1, and PXMP4 between NSCLC tumors of early and late stages, their prognostic and diagnostic potential was shown.
In summary, dy-regulation of peroxisome pathway was common in NSCLC and more than 30% of peroxisome pathway genes were consistently down-or up-regulated in LUSC and LUAD, indicating their crucial roles in NSCLC. To our knowledge, this is the first study for the systemic investigation of peroxisomes in lung cancer. The C-DEPGs might provide new clues for the study of peroxisomes in NSCLC. HSD17B4, ACAA1, and PXMP4 might be new markers for NSCLC diagnosis and prognosis and/or new therapeutic targets for NSCLC treatment. In addition, HSD17B4 and ACAA1 down-regulation highlighted the potential effects of dysfunction of peroxisomal β-oxidation of fatty acids in NSCLC. However, considering the differences in their prognostic effects, regulations, correlations with tumor purity and immune infiltrations effects between the two subtypes, further study is needed to investigate the specific functions of HSD17B4, ACAA1, and PXMP4 in LUSC and LUAD, respectively.

AUTHOR CONTRIBUTIONS
XZ and LD conceived and designed the study. XZ, HY, and FG collected and analyzed the data. XZ and JZ interpreted the data. XZ and HY drafted the manuscript. DL reviewed and revised the manuscript. All authors read and approved the final manuscript.