Association of genetic variants of TMEM135 and PEX5 in the peroxisome pathway with cutaneous melanoma-specific survival.

BACKGROUND: Peroxisomes are ubiquitous and dynamic organelles that are involved in the metabolism of reactive oxygen species (ROS) and lipids. However, whether genetic variants in the peroxisome pathway genes are associated with survival in patients with melanoma has not been established. Therefore, our aim was to identify additional genetic variants in the peroxisome pathway that may provide new prognostic biomarkers for cutaneous melanoma (CM).
METHODS: We assessed the associations between 8,397 common single-nucleotide polymorphisms (SNPs) in 88 peroxisome pathway genes and CM disease-specific survival (CMSS) in a two-stage analysis. For the discovery, we extracted the data from a published genome-wide association study from The University of Texas MD Anderson Cancer Center (MDACC). We then replicated the results in another dataset from the Nurse Health Study (NHS)/Health Professionals Follow-up Study (HPFS).
RESULTS: Overall, 95 (11.1%) patients in the MDACC dataset and 48 (11.7%) patients in the NHS/HPFS dataset died of CM. We found 27 significant SNPs in the peroxisome pathway genes to be associated with CMSS in both datasets after multiple comparison correction using the Bayesian false-discovery probability method. In stepwise Cox proportional hazards regression analysis, with adjustment for other covariates and previously published SNPs in the MDACC dataset, we identified 2 independent SNPs (TMEM135 rs567403 C>G and PEX5 rs7969508 A>G) that predicted CMSS (P=0.003 and 0.031, respectively, in an additive genetic model). The expression quantitative trait loci analysis further revealed that the TMEM135 rs567403 GG and PEX5 rs7969508 GG genotypes were associated with increased and decreased levels of mRNA expression of their genes, respectively.
CONCLUSIONS: Once our findings are replicated by other investigators, these genetic variants may serve as novel biomarkers for the prediction of survival in patients with CM. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

Skin cancers, including cutaneous melanoma (CM), are among the most common cancers worldwide. In the United States, their incidence continues to rise dramatically (1), having doubled in the last 4 decades. In 2019, approximately 96,480 patients were diagnosed with CM in the United States, accounting for 5.5% of all new cancer cases. As the most lethal skin cancer, deaths from CM account for 1.2% of cancer-related mortality (2). Although Breslow thickness, ulceration, and the mitotic rate are considered to be important prognostic indicators, many primary CM cases unexpectedly metastasize after complete surgical resection (3); consequently, CM has the highest mortality rate among all skin cancers. Therefore, it is of great urgency to identify and better understand additional prognostic biomarkers that may provide CM patients with more appropriate and individualized treatment options.

CM is a multifactorial disease originating from melanocytes. Host factors, such as fair skin pigmentation, and environmental factors, such as ultraviolet radiation, can contribute to the development of CM (4). Studies have also shown that genetic susceptibility is involved in both the development and progression of CM (5,6). In recent genome-wide association studies (GWASs), genetic variants have been found to be associated with both risk and survival of CM (7). However, few of the identified genetic variants, particularly single-nucleotide polymorphisms (SNPs), were biologically functional, nor did these studies provide support for the biological plausibility of their findings. To identify truly functional variants, more sophisticated and combined analyses in the post-GWAS era are needed, including the combination of hypothesis-based pathway analysis, meta-analysis, and functional analysis.

Peroxisomes are ubiquitous and dynamic organelles that are involved in the metabolism of reactive oxygen species (ROS) and lipids (8,9). They also act as intracellular signaling platforms in inflammatory, redox, lipid, and innate immune signaling (10-13). Emerging studies suggest that peroxisomes play roles in aging, nonalcoholic fatty liver disease, and neurodegenerative disorders (14,15). Other studies have suggested that peroxisomes play an important role in the progression of several cancers, and abnormal peroxisome metabolism is one of the established characteristics of cancer cells (16-19). However, few studies have investigated the role of peroxisomes in disease progression of CM, and little is known about the effect of genetic variation in the peroxisome pathway genes on the clinical outcomes of CM patients.

Given that some genetic variants are known to be associated with CM prognosis, identifying additional genetic variants in some specific signaling pathways may provide novel prognostic biomarkers of CM. Considering the importance of the peroxisome pathway in cancers, it is likely that genetic variants in peroxisome pathway genes are associated with survival in CM patients. Therefore, we tested this hypothesis by using genotyping datasets from 2 previously published GWASs.

We present the following article in accordance with the TRIPOD reporting checklist (available at http://dx.doi.org/10.21037/atm-20-2117).

Methods

Study populations

In the present study, we used the data of 858 patients from a CM GWAS conducted at The University of Texas MD Anderson Cancer Center (MDACC) (20) as a discovery dataset. Another GWAS of 409 CM patients from 2 Harvard University cohort studies, the Nurses Health Study (NHS) and the Health Professionals Follow-up Study (HPFS) (21), was used as a replication dataset. All the CM patients were non-Hispanic Whites. Informed consent was obtained from all patients. The access to the datasets used in the present study was approved by the Institutional Review Boards of the MDACC (No. LAB03-0048), Brigham and Women’s Hospital (No. 2007-P-000616), and Harvard T.H. Chan School of Public Health (No. 2007-P-000616), and conducted in accordance with the Declaration of Helsinki (as revised in 2013).

The discovery dataset was requested from the Database of Genotypes and Phenotypes (dbGaP: http://www.ncbi.nlm.nih.gov/gap), with an accession number of phs000187.v1.p1 (6). Genotyping was performed with the Illumina HumanOmni-Quad_v1_0_B array (Illumina Inc, San Diego, CA) and subsequent genome-wide imputation (imputation quality r2≥0.8) was performed with the MaCH software using sequencing data from the 1000 Genomes Project, phase I v2 CEU data (March 2010 release) (22). The replication dataset was generated with the Illumina HumanHap610 array, and genome-wide imputation (imputation quality r2≥0.8) was also performed with the MaCH software using sequencing data from the 1000 Genomes Project CEU data (March 2012 release) (23). Compared with the NHS/HPFS dataset, patients in the MDACC dataset had relatively complete clinical information, including age, sex, Breslow thickness, metastasis, mitotic rate, ulceration, and survival outcome, in addition to individual genotyping data, while patients in the NHS/HPFS replication dataset had information only for age, sex, and survival outcome, in addition to genotyping data.

Gene and SNP extraction

We selected 90 peroxisome pathway genes, according to the Kyoto Encyclopedia of Genes and Genomes (KEGG), BioCarta, Reactome, and Gene Ontology (GO) databases in the Molecular Signatures Database (http://software.broadinstitute.org/gsea/msigdb/index.jsp). Because females carry 2 copies of the X chromosome (males are hemizygous), and no standard statistics have been established for sex-specific analysis, genes on the X chromosome were excluded from further analyses. Finally, after the exclusion of 2 genes (ABCD1 and ACSL4) on the X chromosome (Table S1), 88 genes located on the autosomes remained as candidate genes. SNPs within these 88 genes including 2 kb upstream and downstream were extracted from the MDACC GWAS dataset. The quality control criteria for genotyped and imputed SNPs were: (I) a minor allelic frequency (MAF) ≥0.05; (II) a genotyping success rate ≥95%; and (III) Hardy-Weinberg equilibrium (HWE) P≥1×10−5.

Statistical analysis

CM-specific survival (CMSS) was defined as the time from diagnosis with CM to CM-related death or the date of last follow-up. Multivariate Cox proportional hazards regression analysis was performed to assess associations between SNPs and CMSS with an additive model using the GenABEL package of R software. The Bayesian false-discovery probability (BFDP) approach was used, instead of the false discovery rate, for multiple testing correction with a prior probability setting of 0.1 and an upper detected hazards ratio of 3 (24). Only SNPs with a BFDP value <0.8 in both the discovery and replication datasets were chosen for further independent analysis with adjustment for clinical covariates and previously published 40 SNPs from the same MDACC GWAS dataset.

A meta-analysis was further performed using the results of the SNPs from both the discovery and replication datasets using PLINK 1.90. Cochran’s Q statistic and the I2 index were used to assess inter-study heterogeneity. A Cochran’s Q test P value >0.1 and I2<50% indicated no substantial heterogeneity between 2 studies; in such cases, a fixed-effects model was employed. Otherwise, a random-effects model was used. Potential functions of the identified SNPs and those with high linkage disequilibrium (LD) (r2≥0.8) were predicted by RegulomeDB (http://www.regulomedb.org/) and HaploReg4.1 (25). Additionally, Haploview v4.2 (26) was applied to construct a Manhattan plot, and LocusZoom (27) was used to produce regional association plots. A Kaplan-Meier curve was used to estimate survival function, and log-rank tests were performed to compare the effects of different genotypes on CMSS. To estimate the prediction accuracy of the model with demographic variables and the identified SNPs, time-dependent receiver operating characteristic (ROC) analysis was carried out with the “timeROC” package in R (28). Other statistical analyses were performed using SAS software version 9.4 (SAS Institute, Cary, NC, USA), if not specified otherwise.

Results

Basic characteristics of the study populations

The basic characteristics of the 858 CM patients from the MDACC dataset and the 409 CM patients from the NHS/HPFS dataset have been described elsewhere (6). Briefly, the MDACC CM patients were enrolled in a hospital-based case-control study at a tertiary care center; they tended to be younger with more late-stage disease than patients undergoing follow-up in the NHS/HPFS cohort studies. In the MDACC study, patient age at diagnosis ranged between 17 and 94 years, with a mean age of 52.4±14.4 years, and 57.8% (496/858) of the patients were males. The median follow-up time was 81.1 months, and at the last follow-up, 133 patients had died of various causes, including 95 (11.1%) who had died from CM. In the NHS/HPFS studies, patient age at diagnosis ranged between 34 and 87 years, with a mean age of 61.1±10.8 years, and 33.7% (138/409) of the patients were males. The median follow-up time was 179.0 months, which was considerably longer than that of the MDACC study. However, the mortality rate of CM in the NHS/HPFS studies was 11.7% (48/409), which was similar to that in the MDACC patients.

Multivariate analyses of the associations between SNPs in the peroxisome pathway genes and CMSS

As summarized in the study flowchart shown in , we extracted 1,215 genotyped and 7,182 imputed SNPs in the peroxisome pathway genes from the MDACC dataset. Associations between these SNPs and CMSS were evaluated through multivariate Cox proportional hazards regression analysis with adjustment for age, sex, regional/distant metastasis, Breslow thickness, ulceration, and the mitotic rate. In the single-locus analysis of these 8,397 SNPs, 332 were found to be significantly associated with CMSS at P<0.05 in an additive genetic model, and after multiple testing correction, 227 SNPs, with BFDP <0.8, remained statistically noteworthy; these 227 SNPs were then replicated in the NHS/HPFS dataset (Figure S1A). After multivariate Cox regression analysis with adjustment for age and sex in the NHS/HPFS dataset (Figure S1B and Table S2), 27 SNPs in 4 genes, with P<0.05 and BFDP <0.8, remained significantly associated with CMSS. Among the 27 SNPs, there were 16 in TMEM135 (transmembrane protein 135), 7 in ACOX2 (acyl-CoA oxidase 2), 3 in PEX5 (peroxisomal biogenesis factor 5), and 1 in RAB8B (member RAS oncogene family). In the subsequent meta-analysis, these 27 SNPs remained significantly associated with CMSS, with no heterogeneity observed across the 2 datasets ().

Figure 1

Study workflow for SNPs in the peroxisome pathway. AUC, area under the receiver operating characteristic curve; BFDP, Bayesian false-discovery probability; CMSS, cutaneous melanoma-specific survival; GWAS, genome wide association study; HWE, Hardy-Weinberg equilibrium; HPFS, Health Professionals Follow-up Study; MAF, minor allele frequency; MDACC, The University of Texas MD Anderson Cancer Center; NHS, the Nurse Health Study; PEX5, peroxisomal biogenesis factor 5; ROC, receiver operating characteristic; SNP, single nucleotide polymorphism; TMEM135, transmembrane protein 135.

Table 1

Identification of 27 significant and replicated survival-associated SNPs in the peroxisome pathway genes

SNP	Allelea	Gene	Discovery-MDACC (n=858)					Validation-NHS/HPFS (n=409)					Combined-Meta-analysis (n=1,267)
			EAF	HR (95% CI)	Pb	BFDP		EAF	HR (95% CI)	Pc	BFDP		Phet	I2	HR (95% CI)	Pd
rs117285370*	C>T	RAB8B	0.05	2.40 (1.37–4.22)	0.002	0.337		0.06	1.95 (1.03–3.69)	0.039	0.788		0.632	0	2.19 (1.44–3.34)	2.63×10−4
rs7969508*	A>G	PEX5	0.30	1.43 (1.04–1.97)	0.029	0.780		0.24	1.86 (1.23–2.81)	0.003	0.362		0.325	0	1.58 (1.22–2.03)	4.14×10−4
rs7969635	A>G	PEX5	0.30	1.43 (1.04–1.97)	0.029	0.780		0.24	1.86 (1.23–2.81)	0.003	0.362		0.325	0	1.58 (1.22–2.03)	4.14×10−4
rs7969751	A>T	PEX5	0.30	1.43 (1.04–1.97)	0.029	0.780		0.24	1.86 (1.23–2.81)	0.003	0.362		0.325	0	1.58 (1.22–2.03)	4.14×10−4
rs6807633*	C>T	ACOX2	0.09	1.65 (1.04–2.60)	0.032	0.763		0.08	2.10 (1.17–3.76)	0.013	0.628		0.524	0	1.81 (1.26–2.59)	0.001
rs6790046	G>A	ACOX2	0.09	1.65 (1.04–2.60)	0.032	0.763		0.09	1.82 (1.17–2.85)	0.008	0.597		0.631	0	1.77 (1.23–2.53)	0.002
rs9790147	C>G	ACOX2	0.09	1.65 (1.04–2.60)	0.032	0.763		0.08	2.10 (1.17–3.76)	0.013	0.628		0.524	0	1.81 (1.26–2.59)	0.001
rs9789958	G>A	ACOX2	0.09	1.65 (1.04–2.60)	0.032	0.763		0.09	1.89 (1.04–3.40)	0.035	0.767		0.721	0	1.74 (1.21–2.49)	0.001
rs9647382	A>G	ACOX2	0.09	1.65 (1.04–2.60)	0.032	0.763		0.08	2.10 (1.17–3.76)	0.013	0.628		0.524	0	1.81 (1.26–2.59)	0.001
rs76039946	G>A	ACOX2	0.09	1.64 (1.04–2.59)	0.033	0.775		0.08	1.97 (1.08–3.59)	0.028	0.740		0.634	0	1.75 (1.22–2.52)	0.002
rs113351358	C>T	ACOX2	0.09	1.64 (1.04–2.59)	0.033	0.775		0.08	1.97 (1.08–3.59)	0.028	0.740		0.634	0	1.75 (1.22–2.52)	0.002
rs567403*	C>G	TMEM135	0.10	1.61 (1.03–2.51)	0.036	0.782		0.07	1.97 (1.03–3.74)	0.040	0.782		0.613	0	1.72 (1.19–2.48)	0.004
rs542279	C>T	TMEM135	0.10	1.60 (1.03–2.50)	0.038	0.793		0.07	2.57 (1.38–4.78)	0.003	0.386		0.224	32.46	1.88 (1.31–2.70)	6.24×10−4
rs631930	G>A	TMEM135	0.10	1.60 (1.03–2.50)	0.038	0.793		0.07	2.57 (1.38–4.78)	0.003	0.386		0.224	32.46	1.88 (1.31–2.70)	6.24×10−4
rs520718	C>T	TMEM135	0.10	1.60 (1.03–2.50)	0.038	0.793		0.07	2.57 (1.38–4.78)	0.003	0.386		0.224	32.46	1.88 (1.31–2.70)	6.24×10−4
rs35707386	G>T	TMEM135	0.10	1.60 (1.03–2.50)	0.038	0.793		0.06	2.48 (1.31–4.69)	0.005	0.492		0.269	18.24	1.85 (1.28–2.66)	9.71×10−4
rs556724	G>A	TMEM135	0.10	1.60 (1.03–2.50)	0.038	0.793		0.07	2.57 (1.38–4.78)	0.003	0.386		0.224	32.46	1.88 (1.31–2.70)	6.24×10−4
rs605658	T>C	TMEM135	0.10	1.60 (1.03–2.50)	0.038	0.793		0.07	2.57 (1.38–4.78)	0.003	0.386		0.224	32.46	1.88 (1.31–2.70)	6.24×10−4
rs618159	A>G	TMEM135	0.10	1.60 (1.03–2.50)	0.038	0.793		0.07	2.57 (1.38–4.78)	0.003	0.386		0.224	32.46	1.88 (1.31–2.70)	6.24×10−4
rs2847352	T>A	TMEM135	0.10	1.60 (1.03–2.50)	0.038	0.793		0.07	2.57 (1.38–4.78)	0.003	0.386		0.224	32.46	1.88 (1.31–2.70)	6.24×10−4
rs498059	A>G	TMEM135	0.10	1.60 (1.03–2.50)	0.038	0.793		0.07	2.57 (1.38–4.78)	0.003	0.386		0.224	32.46	1.88 (1.31–2.70)	6.24×10−4
rs578116	G>A	TMEM135	0.10	1.60 (1.03–2.50)	0.038	0.793		0.07	2.57 (1.38–4.78)	0.003	0.386		0.224	32.46	1.88 (1.31–2.70)	6.24×10−4
rs482944	G>A	TMEM135	0.10	1.60 (1.03–2.50)	0.038	0.793		0.07	2.57 (1.38–4.78)	0.003	0.386		0.224	32.46	1.88 (1.31–2.70)	6.24×10−4
rs34782550	T>C	TMEM135	0.10	1.60 (1.03–2.50)	0.038	0.793		0.07	2.57 (1.38–4.78)	0.003	0.386		0.224	32.46	1.88 (1.31–2.70)	6.24×10−4
rs34806361	A>G	TMEM135	0.10	1.60 (1.03–2.50)	0.038	0.793		0.07	2.57 (1.38–4.78)	0.003	0.386		0.224	32.46	1.88 (1.31–2.70)	6.24×10−4
rs500725	A>G	TMEM135	0.10	1.60 (1.03–2.50)	0.038	0.793		0.07	2.57 (1.38–4.78)	0.003	0.386		0.224	32.46	1.88 (1.31–2.70)	6.24×10−4
rs493708	C>T	TMEM135	0.10	1.60 (1.03–2.50)	0.038	0.793		0.07	2.57 (1.38–4.78)	0.003	0.386		0.224	32.46	1.88 (1.31–2.70)	6.24×10−4

a, reference allele/effect allele; b, adjusted for age, sex, Breslow thickness, distant/regional metastasis, ulceation and mitotic rate in the additive model; c, adjusted for age and sex in the additive model; d, meta-analysis in the fix-effects model. SNPs with * remained statistically significantly in independent analysis. SNP, single-nucleotide polymorphism; MDACC, MD Anderson Cancer Center; NHS, the Nurse Health Study; HPFS, the Health Professionals Follow-up Study; EAF, effect allele frequency; HR, hazards ratio; CI, confidence interval; BFDP, Bayesian false-discovery probability; Phet, P value for heterogeneity by Cochrane’s Q test; ACOX2, acyl-CoA oxidase 2; PEX5, Peroxisomal Biogenesis Factor 5; TMEM135, transmembrane protein 135; RAB8B, member RAS oncogene family.

Genetic variants in the peroxisome pathway genes as independent predictors of survival

To further identify the independent predictors for CMSS survival, initial stepwise multivariate Cox regression analyses were performed to assess the effects of the 27 replicated SNPs on the survival of CM patients. These analyses were performed in the MDACC dataset only, because the NHS/HPFS dataset did not have the same detailed genotyping data and clinical covariates. Four SNPs in 4 genes remained statistically significantly associated with CMSS in the presence of covariates (i.e., age, sex, Breslow tumor thickness, regional/distant metastasis, tumor cell mitotic rate, and ulceration of tumors). These SNPs were ACOX2 rs6807633 (imputed), TMEM135 rs567403 (imputed), PEX5 rs7969508 (genotyped), and RAB8B rs117285370 (imputed). Finally, following adjustment for the 40 additional previously published survival-associated SNPs with the same MDACC GWAS dataset, we found that 2 SNPs (TMEM135 rs567403 C>G and PEX5 rs7969508 A>G) remained significant (P=0.003 and 0.031, respectively) for further analysis ().

Table 2

Predictors of CMSS obtained from the stepwise Cox proportional hazards regression analysis with adjustment for the published 40 survival-associated SNPs in the same MDACC dataset

Parametera	Categoryb	Frequency	HR (95% CI)a	Pa	HR (95% CI)c	Pc
Age	≤50/>50	371/487	1.02 (1.00–1.04)	0.018	1.05 (1.03–1.07)	<0.0001
Sex	Female/male	362/496	1.48 (0.93–2.36)	0.100	1.15 (0.68–1.93)	0.610
Regional/distant metastasis	No/yes	709/149	4.04 (2.61–6.26)	<0.0001	13.71 (7.38–25.47)	<0.0001
Breslow thickness(mm)	≤1/>1	347/511	1.18 (1.12–1.25)	<0.0001	1.26 (1.17–1.36)	<0.0001
Ulceration	No/yes	681/155	2.96 (1.91–4.61)	<0.0001	4.53 (2.64–7.76)	<0.0001
Mitotic rate (mm2)	≤1/>1	275/583	2.46 (1.21–4.98)	0.013	2.14 (0.95–4.79)	0.065
TMEM135 rs567403 C>G	CC/CG/GG	693/154/11	1.77 (1.13–2.77)	0.013	2.31 (1.34–3.97)	0.003
PEX5 rs7969508 A>G	AA/AG/GG	406/388/64	1.45 (1.05–2.00)	0.023	1.57 (1.04–2.37)	0.031

a, Stepwise Cox analysis included age, sex, regional/distant metastasis, Breslow thickness, ulceration, mitotic rate and 27 replicated SNPs; b, the “category/” was used as the reference; c, adjustment for 40 previously published survival-associated SNPs in the same MDACC dataset. CMSS, cutaneous melanoma-specific survival; MDACC, The University of Texas MD Anderson Cancer Center; HR, hazards ratio; CI, confidence interval; PEX5, Peroxisomal Biogenesis Factor 5; TMEM135, transmembrane protein 135.

As shown in , in the MDACC dataset, the effects of the TMEM135 rs567403 G and PEX5 rs7969508 G alleles on the survival of CM patients were statistically significant (trend test: P=0.036 and 0.030, respectively), and similar results were observed in the NHS/HPFS studies (trend test: P=0.040 and 0.004, respectively) as well as in their combined dataset (trend test: P=0.023 and 0.001, respectively). We then depicted these associations by using Kaplan-Meier survival curves (). Regional association plots were also generated to display the LD between the identified SNPs and those in TMEM135 and PEX5, including the 50 kb regions flanking these 2 genes (Figure S2A,B).

Table 3

Associations between two independent SNPs in the peroxisome pathway genes and CMSS of the patients in the MDACC dataset, the NHS/HPFS dataset and the combined dataset

Genotype	MDACC (n=858)					NHS/HPFS (n=409)					MDACC + NHS/HPFS (n=1,267)
	All	Death (%)	HR (95% CI)a	Pa		All	Death (%)	HR (95% CI)b	Pb		All	Death (%)	HR (95% CI)c	Pc
TMEM135 rs567403 C>G
CC	693	72 (10.39)	1.00			357	38 (10.64)	1.00			1,050	110 (10.48)	1.00
CG	154	22 (14.29)	1.71 (1.04–2.78)	0.035		50	9 (18.00)	1.71 (0.83–3.53)	0.149		204	31 (15.20)	1.55 (1.04–2.31)	0.031
GG	11	1 (9.09)	1.71 (0.23–12.68)	0.602		2	1 (50.00)	9.66 (1.28–72.77)	0.028		13	2 (15.38)	1.88 (0.46–7.69)	0.379
Trend test				0.036					0.040					0.023
CG+GG	165	23 (13.94)	1.71 (1.05–2.78)	0.032		52	10 (19.23)	1.86 (0.93–3.73)	0.081		217	33 (15.21)	1.56 (1.06–2.32)	0.024
PEX5 rs7969508 A>G
AA	406	38 (9.36)	1.00			236	20 (8.47)	1.00			642	58 (9.03)	1.00
AG	388	48 (12.37)	1.63 (1.05–2.55)	0.031		148	21 (14.19)	1.74 (0.94–3.21)	0.079		536	69 (12.87)	1.51 (1.07–2.15)	0.020
GG	64	9 (14.06)	1.73 (0.79–3.78)	0.170		25	7 (28.00)	3.59 (1.52–8.52)	0.004		89	16 (17.98)	2.25 (1.29–3.92)	0.004
Trend test				0.030					0.004					0.001
AG+GG	452	57 (12.61)	1.65 (1.07–2.53)	0.023		173	28 (16.18)	2.00 (1.12–3.55)	0.019		625	85 (13.60)	1.61 (1.16–2.25)	0.005
Number of combined risk genotypesd
0	329	30 (9.12)	1.00			205	16 (7.80)	1.00			534	46 (8.61)	1.00
1	441	50 (11.34)	1.45 (0.90–2.32)	0.126		183	26 (14.21)	1.90 (1.02–3.55)	0.043		624	76 (12.18)	1.48 (1.02–2.13)	0.037
2	88	15 (17.05)	2.96 (1.57–5.57)	0.008		21	6 (28.57)	3.92 (1.52–10.12)	0.005		109	21 (19.27)	2.72 (1.62–4.56)	0.0002
Trend test				0.002					0.003					0.0003
0–1	770	80 (10.39)	1.00			388	42 (10.82)	1.00			1,158	122 (10.54)	1.00
2	88	15 (17.05)	2.39 (1.37–4.19)	0.002		21	6 (28.57)	2.76 (1.16–6.55)	0.022		109	21 (19.27)	2.17 (1.36–3.45)	0.001

a, adjusted for age, sex, Breslow thickness, distant/regional metastasis, ulceration and mitotic rate in the MDACC dataset; b, adjusted for age and sex in the NHS/HPFS dataset; c, adjusted for age and sex in the MDACC and NHS/HPFS combined dataset; d, risk genotypes included TMEM135 rs567403 CG+GG and PEX5 rs7969508 AG+GG. SNP, single-nucleotide polymorphism; CMSS, cutaneous melanoma-specific survival; MDACC, The University of Texas MD Anderson Cancer Center; NHS, the Nurse Health Study; HPFS, the Health Professionals Follow-up Study; HR, hazards ratio; CI, confidence interval; PEX5, Peroxisomal Biogenesis Factor 5; TMEM135, transmembrane protein 135.

Figure 2

Kaplan-Meier survival curves of cutaneous melanoma-specific survival (CMSS): TMEM35 rs567403 in dominant model in MDACC dataset (A), NHS/HPFS dataset (B), and the MDACC and NHS/HPFS combined dataset, (C) and PEX5 rs7969508 in dominant model in the MDACC dataset (D), NHS/HPFS dataset (E), and the MDACC and NHS/HPFS combined dataset (F). Kaplan-Meier survival curves of the combined risk genotypes in CMSS: 0-1 risk genotypes group and 2 risk genotype group in MDACC (G), NHS/HPFS (H), and MDACC and NHS/HPFS combined datasets (I). 1Univariate analysis; 2Multivariate analysis. SNP, single nucleotide polymorphism; CM, cutaneous melanoma; PEX5, peroxisomal biogenesis factor 5; TMEM135, transmembrane protein 135; MDACC, The University of Texas MD Anderson Cancer Center; NHS, the Nurse Health Study; HPFS, the Health Professionals Follow-up Study.

Combined genotype analyses of the 2 independent SNPs

To assess the joint effect of the 2 independent SNPs on CMSS, we further combined the risk genotypes of TMEM135 rs567403 CG+GG and PEX5 rs7969508 AG+GG into a genetic score as the number of risk genotypes (NRG), and categorized all the patients into 3 groups with 0, 1, and 2 risk genotypes. The trend test showed that a higher NRG was associated with a progressively worse survival in the MDACC dataset in a dose-response manner (P=0.002), as well as in the NHS/HPFS dataset (P=0.003) and the combined dataset (P=0.0003), with adjustments for the available covariables (). Next, we further dichotomized the NRG score into 2 groups: low-risk (0-1 risk genotype) and high-risk (2 risk genotypes). We found that patients in the high-risk group had a higher CM-death risk in the MDACC dataset (HR =2.39, 95% CI: 1.37–4.19, P=0.002) and the NHS/HPFS dataset (HR =1.97, 95% CI: 1.16–6.55, P=0.022), as well as in the combined dataset (HR =2.17, 95% CI: 1.36–3.45, P=0.001), compared to those in the low-risk group (). These associations were also depicted by Kaplan-Meier curves ().

Stratified analyses for the effect of combined risk genotypes on CMSS

In the previous analysis we found that patients in the high NRG score group showed a substantially greater risk of CM-related death in the presence of the available clinical variables. The risk was more evident in the subgroup aged ≤50, male patients, patients with ulceration or late-stage disease, and those with a tumor cell mitotic rate >1/mm2 in the MDACC dataset, as well as in the subgroup aged >50 and female patients in the NHS/HPFS dataset. However, no interactions were found (Table S3).

ROC curve and area under the curve (AUC) for prediction of CMSS

We further assessed the predictive effect of the 2 independent SNPs on CMSS using the time-dependent AUC of the ROC curve. The ROC of the 5-year CMSS showed that the inclusion of the 2 identified SNPs only improved the predictive performance of the models including clinical variables without and with risk genotypes by 1.6% in the MDACC dataset (AUC =65.88% vs. 67.48%, P=0.336) (Figure S3A,B). However, we observed that the model with age, sex, and genetic variants had a better performance than those only including demographic variables in the NHS/HPFS dataset (54.05% vs. 70.88%, P=0.003) (Figure S3C,D) and in the combined dataset (AUC =63.6% vs. 67.61%, P=0.073) (Figure S3E,F).

Functional predictions of the 2 independent SNPs

Functional prediction by RegulomeDB showed that TMEM135 rs567403 C>G and PEX5 rs7969508 A>G had RegulomeDB scores of 5 and 1f, respectively, and also indicated that these SNPs might be located at DNase I regulating sites or transcription factor binding sites. We further searched for SNPs with high LD (r2≥0.8) with these 2 independent SNPs and made functional predictions by using HaploReg. We found that TMEM135 rs567403 C>G might disrupt or change the motifs of Dobox4 and Irf, whereas PEX5 rs7969508 A>G, located in DNase I hypersensitive sites, may be a marker of promoter histone and enhancer histone in embryonic stem cell-derived tissue (ESDR), and may have a high linear correlation with mRNA expression of its corresponding gene PEX5 (Table S4). We further assessed the potential functions of these 2 independent SNPs by using data from the ENCODE Project. TMEM135 rs567403 C>G has no significant finding (Figure S4A), while PEX5 rs7969508 A>G is located in the H3K4Me1 region, which may have potential enhancer activity (Figure S4B).

Two independent SNPs and mRNA expression levels of the associated genes

To provide molecular support for the observed survival associations with, and predictions by, the genotypes, we further performed expression quantitative trait locus (eQTL) analysis to explore correlations between genotypes of the 2 independent SNPs and their associated mRNA expression levels in the publicly available RNA sequencing data of lymphoblastoid cell lines generated from 373 European descendants in the 1000 Genomes Project. The eQTL analysis revealed that rs567403 C>G was significantly correlated with elevated mRNA expression levels of TMEM135 in an additive model (P=0.003, ), while rs7969508 A>G demonstrated a significant correlation with decreased mRNA expression levels of PEX5 in all genetic models (P=1.38×10–10, P=3.91×10–10, and P=0.0018 for the additive, dominant, and recessive models, respectively; ). We also performed an eQTL analysis using data from the Genotype-Tissue Expression (GTEx) Project V7.p2 (http://www.gtexportal.org/home). The results showed that rs7969508 A>G was correlated with significantly decreased PEX5 mRNA expression in normal tissues from unexposed suprapubic skin (P=2.6×10–18), sun-exposed lower leg skin (P=8.6×10–23), and whole blood cells (P=4.1×10–14) (). These results were consistent with the findings of the 1000 Genomes Project. However, there were no significant correlations between rs567403 genotypes and mRNA expression levels of TMEM135 in normal skin tissues from unexposed suprapubic skin (P=0.260), sun-exposed lower leg skin (P=0.480), or whole blood cells (P=0.730) (Figure S5) from the GTEx project.

Figure 3

The expression quantitative trait loci (eQTL) analysis for genotypes of PEX5 rs7969508 and TMEM135 rs567403. (A) Correlation between TMEM135 mRNA expression levels and rs567403 in 373 Europeans from the 1000 Genomes Project in an additive model. Correlation between PEX5 mRNA expression and rs7969508 genotypes in 373 Europeans from the 1000 Genomes Project in additive (B), dominant (C), and recessive (D) models. Correlation between PEX5 mRNA expression and rs7969508 in unexposed skin (E), sun-exposed skin (lower leg) (F), and whole blood cells (G) in GTEx. GTEx, Genotype-Tissue Expression Project; PEX5, peroxisomal biogenesis factor 5; TMEM135, transmembrane protein 135.

Discussion

In the present study, we assessed the associations between SNPs in the peroxisome pathway-related genes and CM survival. We found that 2 independent and potentially functional SNPs (TMEM135 rs567403 C>G and PEX5 rs7969508 A>G) were independently or jointly associated with survival in CM patients. Patients with a higher NRG score of these 2 genetic variants had a progressively worse survival. Importantly, the effect was consistent across most subgroups. Moreover, the TMEM135 rs567403 risk G allele was found to be significantly correlated with increased mRNA expression levels of TMEM135, whereas the PEX5 rs7969508 risk G allele was found to be significantly correlated with decreased mRNA expression levels of PEX5. These findings suggest that SNPs in peroxisomal pathway genes may play biological roles in the prognosis of CM, suggesting some biological plausibility of the survival-associated SNPs.

Peroxisome homeostasis is achieved by mediating the responses of peroxisome biogenesis and degradation to environmental conditions. By maintaining the integrity and number of organelles under different environmental stresses, pexophagy, as the process of selective autophagy, is essential for maintaining the stability of cellular homeostasis, which is one of the mechanisms underlying peroxisome degradation (29). Peroxisomes are important sites of ROS production and degradation, and thus, play crucial roles in the development and progression of cancer (30). High levels of intracellular ROS can trigger cancer development by promoting pro-oncogenic mutations, and excessive ROS may contribute to progression of cancer cells (31). These may be the mechanisms through which peroxisomes affect cancer development and prognosis. However, no published studies have investigated the molecular mechanisms of the peroxisomes involved in CM survival.

TMEM135, located on chromosome 11q14.2, encodes transmembrane protein 135, which is one of peroxisomal proteins belonging to the Tim17 protein family. To date, the role of TMEM135 in the biological processes of energy expenditure and the related fatty acid metabolism have remained unclear (32,33). Few studies have investigated the roles of TMEM135 in cancer. One study that used sequencing analysis of independent cohorts for breast cancer found that TMEM135 was a potential novel gene (34), and another study reported that rs11235127 near TMEM135 had a statistically significant association with breast cancer risk (35). In prostate cancer patients, a novel fusion transcript in TMEM135-CCDC67 was found to be associated with prostate cancer metastasis, recurrence, and prostate cancer-specific death after operation (36). A genome-wide linkage analysis of Spanish melanoma-prone families found a locus at 11q14.1-q14.3 with significant linkage evidence from multiple pedigrees, and the subregion contained TMEM135 (37). In the present study, we found that the rs567403 G allele might up-regulate the expression of TMEM135, suggesting that TMEM135 may play an oncogenic role in CM biology. However, additional experimental investigations are required to determine how TMEM135 rs567403 C>G influences CMSS.

PEX5, located on chromosome 12p13.31, encodes the type-1 peroxisomal targeting signal (PTS1) receptor necessary for peroxisome biogenesis (38). Few studies have investigated the role of PEX5 in cancer. One study suggested that under starvation conditions, PEX5 depletion prevented the activation of autophagy (39), an important process in intracellular degradation that results in programmed cell death or the promotion of cell survival (40). Therefore, it is likely that the expression levels of PEX5 may alter the processes of autophagy or cell death, which then affects the development and prognosis of cancer. Studies have also reported that cell proliferation was inhibited after PEX5 depletion in human hepatocellular carcinoma cells (16), and PEX5 expression was significantly increased in colon carcinoma tissue compared with normal colon mucosa from 5 patients, as determined by northern blotting analysis (41). These findings suggest that PEX5 may play an oncogenic role. However, in the present study we found that the PEX5 rs7969508 G allele down-regulated the mRNA expression of PEX5, resulting in poor survival in CM patients. This suggests that PEX5 may play a suppressor role in CM biology, which appears to disagree with reports about this gene in the literature, although we did not have additional experimental data to explain this finding.

The limitations of the present study are mainly related to the 2 GWAS datasets that were available to us. The sample sizes of the 2 GWAS datasets were insufficient, and some valuable clinical information, such as treatment and response rate, were not available for adjustment in the analysis; thus, the associations between these SNPs and CM survival should be further investigated in larger studies. Furthermore, the prognosis prediction model was built in non-Hispanic white populations of the United States, which may not be generalizable to other populations. Finally, additional functional experiments for these 2 independent SNPs should be performed to identify the molecular mechanisms underlying the observed survival association in CM patients.

Conclusions

In summary, we analyzed SNPs in the peroxisome pathway genes for their associations with CM survival using data from 2 published GWAS datasets, both of which had strict quality-control standards. The 2 novel SNPs we identified were predicted to have potential biological functions. Once our findings are replicated by other investigators, these genetic variants may serve as new biomarkers for predicting survival in CM patients.

The article’s supplementary files as