Association between putative functional variants in the PSMB9 gene and risk of melanoma--re-analysis of published melanoma genome-wide association studies.

Although the number of nevi and sun exposure (sunburn sensitivity) as well as other pigmentation-related characteristics are the known risk factors for melanoma ( Titus-Ernstoff et al., 2005 ), numerous studies have suggested that genetic variants in genes involved in several biologic pathways also have an impact on melanoma susceptibility ( Curtin et al., 2005 ; Figl et al., 2010 ; Zhang et al., 2010 ). Recently, genome-wide association studies (GWASs) have identified several new promising causal genetic variants or loci for melanoma risk and verified a few genes that were identified in earlier studies of high-risk melanoma kindreds ( Barrett et al., 2011 ; Bishop et al., 2009 ; Brown et al., 2008 ; Falchi et al., 2009 ; Macgregor et al., 2011 ). However, these GWAS-level significant SNPs explain only a small proportion of heritable melanoma risk. To identify additional novel but low-penetrance causal variants, the pathway-based analysis of GWAS data has been widely applied. Such analyses appear to be effective in detecting SNPs that confer relatively small but measurable and potentially biologically significant risk ( Donnelly, 2008 ; Elbers et al., 2009 ; Hong et al., 2009 ; Luo et al., 2010 ), thereby helping us to better understand the mechanisms underlying melanoma etiology. The cell cycle consists of a series of events that occur during cell division and duplication. It can be briefly divided in two periods: interphase and mitosis (including G1, S, G2 and M phases, linked by cell-cycle checkpoints, G1/S transition, and M/G1 transition). The dysfunction of cell cycle often leads to disordered cell growth and ultimately carcinogenesis ( Malumbres and Barbacid, 2007 ). The mutation or abnormal expression of cell-cycle related genes have been reported in studies of many cancers, such as lung cancer, breast cancer and melanoma ( Malumbres and Barbacid, 2007 ). The M/G1 transition is one of the important cell-cycle processes, which includes exiting from mitosis and onset of the G1 phase. The genes function of this process have been involved in the development of melanoma and other cancers ( Anantha and Borowiec, 2009 ; Clarke et al., 2009 ; Du et al., 2004 ; Khan et al., 2008 ). Currently, genetic variants in these genes have also been studied in several cancers ( Cunningham et al., 2009 ; Frank et al., 2008 ; Ma et al., 2011 ; Xiong et al., 2009 ). However, most of these early published candidate-gene studies have investigated only a few genes or variants in this cell cycle process. Recent data from melanoma GWASs provide a unique opportunity for us to elucidate the impact of other unknown low-penetrance variants of the genes involved in the cell-cycle control pathway on melanoma development. In the present study, we first performed a pathway-based analysis using our published GWAS dataset to evaluate the association between genetic variants in 76 M/G1 transition genes and melanoma risk. Then, we selected significant SNPs with putative functions for validation in other two GWAS datasets. Additional validations of the most promising functional variants were further performed through bioinformatics and laboratory approaches to provide biological supports for our findings.

The distribution of the major melanoma risk factors between cases and controls is presented in Supplementary Table S1 , which shows that cases were more likely to have light skin / hair color, more moles, more frequent severe sunburns and dysplastic nevi than were controls ( P < 0.001 for all significant variables). In the present study, a total of 1,149 SNPs in 76 M/G1 transition-related genes were extracted from our GWAS dataset ( Supplementary Table S2 ). The gene-based test had been performed with the VEGAS method ( Liu et al., 2010 ), which revealed seven genes with P-value < 0.05. A list of the SNPs with P value < 0.01 in the discovery set and their assigned genes are shown in Supplemental Table S3 , including 68 SNPs in eight genes. There were 34 SNPs with P value < 0.05 after corrections for multiple testing by Benjamin and Hochberg FDR method ( Benjamini and Hochberg, 1995 ). Most of the 68 SNPs (57/68 = 83.8%) were mapped within the PSMB9 gene region on chromosome 6, and the gene-based P value of PSMB9 was 0.003 according to the VEGAS method. Therefore, we focused on this region by selecting 18 SNPs with some putative functions for the in silico replication ( Supplementary Figure S1 ). Validation results are shown in Table 1 that used actual genotyping data for all of the 18 SNPs in the three datasets. Two significant SNPs in the discovery dataset were replicated in the GenoMEL (UK) dataset: rs1351383 in the first intron of PSMB9 ( P discovery = 0.005, P replication in UK = 0.013, P joint = 2.29×10 −4 ), rs2127675 in the 3’ flanking of PSMB9 ( P discovery = 0.001, P replication in UK = 0.004, P joint = 1.29×10 −5 ) but replication failed in the Australian GWAS dataset. In the meta-analysis of the three datasets, we found that the P fix value for rs1351383 in the fixed effect model was 0.052 and P fix value for rs2127675 was 0.006 ( Table 1 ). However, no significance remained in the random model ( P (R) = 0.255 and P (R) = 0.163, respectively), likely due to large heterogeneity after combining with the Australian dataset). The regional association plot for the PSMB9 region in the discovery set is presented in Figure 1 with additional 163 imputed SNPs. We then applied four genetic models to these two SNPs in our discovery GWAS dataset. It should be noted that this might overestimate the genetic effect when just using the discovery dataset due to the “Winner’s course” ( Zollner and Pritchard, 2007 ). For PSMB9 rs2127675, subjects carrying the AG or GG genotype had an increased risk of melanoma ( P = 8.00 × 10 −4 , OR = 1.37, 95% CI: 1.12–1.68; P = 1.93 × 10 −3 , OR = 1.65, 95% CI: 1.19–2.27, respectively), when compared with those with the AA genotype. The association was more significant under the dominant model ( P = 1.65 × 10 −4 , OR = 1.42, 95% CI: 1.17–1.73). When stratified by skin color, nevi, and moles status, significant associations were found mainly in subgroups of light skin color or with moles ( P = 0.005 and 5.37×10 −4 , respectively). Similar results were found for SNP rs1351383, which may be due to the fact that these two SNPs are in the same block with a strong LD (r 2 = 0.79) ( Table 2 ). In addition, we evaluated the mRNA expression of PSMB9 by the genotypes of rs2127675 and rs1351383 in 270 lymphoblastoid cell lines derived from the HapMap populations ( Figure 2 ). The risk genotypes of rs2127675 AG/GG were shown to be associated with higher expression levels of PSMB9 ( P trend = 0.024 in the CEU population and P trend = 0.004 in all unrelated populations, respectively) than the common AA genotype. Similar results were found for rs1352383: risk AC/CC genotypes were associated with higher expression levels, compared with the common genotype AA ( P trend = 0.049 in the CEU population and P trend = 0.002 in all unrelated populations, respectively). Although rs2071480 was not included in the GenoMEL (UK) GWAS dataset, this SNP was in high LD with rs1351383 and rs2127675 (r 2 = 0.99 and r 2 =0.79, respectively) and was at −79bp upstream of the transcription start site of PSMB9 . This SNP was also predicted to be located at the putative transcription factor binding sites by SNPinfo. Therefore, we further examined whether rs2071480 could change the binding affinity of transcriptional factors by the electrophoresis migration shift assay (EMSA). As shown in Figure 3A , the nuclear proteins prepared from A375 melanoma cells were able to bind to both oligo probes containing either rare or common or alleles of this SNP ( Figure 3A , lanes 2 and 6). Compared with the common G allele-specific shifted band (DNA-Protein complex) in lane 6, the relative intensity of shifted band for the variant T allele (lane 1) was slightly increased (the intensity ratios of shifted/upshifted bands are 0.16 and 0.30, respectively), indicating that the nuclear protein binding activity with T-allele oligo was stronger than that with the G allele oligo. However no shifted band was observed when 10× or 50× unlabeled probes with either the T or G allele were added to compete with the labeled probes, which indicated there was not a large difference between the binding activities of probes with different alleles ( Figure 3A , lanes 2, 3, 4, 5 and lanes 7, 8, 9, 10). To further test for the effect of rs2071480 on the promoter activity of the PSMB9 gene, the promoter sequence containing either G or T allele was inserted into the pGL3 vector ( Figure 3B and 3C ). We selected two clones that were identical to each other except for the SNP site and compared their promoter activities among A375, Hela and HCT116 cancer cell lines. As shown in Figure 3D , the luciferase activities driven by the construct containing T allele increased 1.4–2.3 folds compared with those driven by G allele construct ( P < 0.01). The expression analysis also showed that GT/TT genotypes were associated with higher mRNA expression levels of PSMB9 in HapMap lymphoblastoid cell lines ( P trend = 0.047 for 84 CEU samples and P trend = 0.004 for 201 unrelated samples, respectively; Supplementary Figure S2 ). These results were consistent with that of the EMSA assay ( Figure 3.A ). To provide further evidence for the association between PSMB9 expression levels and melanoma risk, we mined the microarray data at NCBI's Gene Expression Omnibus (GEO: http://www.ncbi.nlm.nih.gov/geo/ ). There was one previous study that had detected the differences in global gene expression profiles among seven normal skin, 18 benign nevi and 45 primary melanoma tissues ( Talantov et al., 2005 ). By using that expression dataset (GEO accession number: GSE3189 ), we compared the difference in PSMB9 mRNA expression levels between melanoma and benign tissues. As shown in Supplementary Figure S3 , the PSMB9 expression levels in primary melanoma were significantly higher than those in normal skin ( P = 0.013) and benign nevi ( P = 0.005), supporting a role of PSMB9 over expression in skin carcinogenesis.

In this GWAS-based pathway study, we investigated the association between genetic variants in the M/G1 transition-related genes and melanoma risk, using three published GWAS datasets. In the discovery phase, we first found that multiple SNPs with P-value < 0.01 located in the PSMB9 region, and we then selected 18 putative functional SNPs in this region to perform validation in other two GWAS datasets. We found that two SNPs, rs1351383 and rs2127675, also showed significant associations with melanoma risk in the GenomMEL GWAS dataset, but the validation failed in the Australian GWAS dataset. Further mRNA expression analyses and functional assays provided evidence that these SNPs might be associated with melanoma risk by increasing PSMB9 mRNA expression. Stratified analysis indicated the risk associated with these SNPs was more evident in subgroups with moles and fair skin, suggesting these PSMB9 SNPs may play a role in susceptibility to melanoma associated with pigmentogenesis and mole development. To our knowledge, this report provides the first evidence from large GWAS datasets for associations between PSMB9 SNPs and melanoma risk in a non-Hispanic white population. PSMB9 , also known as LMP2 , is located at 6q21, the high-risk region of major histocompatibility complex. It encodes one subunit of immunoproteasome, which involves in the antigen processing and presentation. Down-regulation of antigen-processing molecules, found in a variety of cancers, may change the spectrum of peptides presented by MHC molecules and may be associated with immune escape of tumors ( Igney and Krammer, 2002 ). However, the association between PSMB9 expression levels and melanoma risk has been inconclusive. One previous study had reported down-expression of antigen processing molecules in malignant melanoma using the immunohistochemical method, but it did not find significant difference in PSMB9 expression levels between nevi and primary melanoma lesions ( Kageshita et al., 1999 ). Another study also reported that the down-regulation of LMP7 and TAP2 , but not PSMB9 , was correlated with levels of the HLA class I surface expression in melanoma cell lines ( Mendez et al., 2008 ). However, Dannull and colleagues found that down-regulation of immunoproteasome by siRNA transfection of LMP2 , LMP7 and MECL-1 could stimulate the enhanced antimelanoma CTL activity ( Dannull et al., 2007 ). By mining the GEO database, we found that mRNA expression levels of PSMB9 increased significantly in melanoma tissues, compared with that of the normal skin and benign nevi. This may provide some additional support for our finding of an association between genotypes and an elevated expression of PSMB9 , which would increase melanoma risk. Further functional validations are warranted. The roles of genetic variants of PSMB9 have been investigated in several non-melanoma cancers and autoimmune diseases. Up to now, most of these studies focused on limited number of exonic variants. One PSMB9 variant of Arg60His (rs17587), which may influence the gene’s functions ( Mishto et al., 2006 ), had been reported to be associated with hypertension in adolescents ( Honcharov et al., 2009 ), multiple sclerosis ( Mishto et al., 2010 ), acute anterior uveitis ( Maksymowych et al., 1997 ), Mycobacterium tuberculosis infection ( Lv et al., 2011 ), insulin-dependent diabetes mellitus ( Deng et al., 1995 ), and cervical carcinoma ( Deshpande et al., 2008 ). Another exonic variant of Ile32Val (rs241419) was found to be weakly associated with colorectal cancer risk in a UK study ( Webb et al., 2006 ), but it was not replicated in a study of German families ( Frank et al., 2008 ). We have investigated the association of these variants with melanoma risk in the present study but no significant association was observed (data was not shown). In addition, other cell cycle-related genes and SNPs reported before ( Choudhury et al., 2004 ; Darieva et al., 2010 ; Lan et al., 2006 ) were not replicated in the present study. It should be noted that we did not find association that could reach genome-wide significance level and none SNPs were replicated in the Australian GWAS dataset. Although there are several SNPs in PSMB9 that could pass the multiple comparison correction (FDR = 0.045), we also cannot exclude possibility that these SNPs are not be true causal SNPs of melanoma. Considering the confounding influence of sun exposure on melanoma risk, one possible reason for the replication failure might be due to the difference in sun-exposure between the three study populations living at different latitudes ( Chang et al., 2009 ). Further stratification analysis by sun-exposure may help to estimate genetic effects of these variants on melanoma risk in Australia population. Heterogeneity of age at onset of the disease may be another reason for the failure of replication. The median onset age for the patients included in MD Anderson study was 51 years ( Amos et al., 2011 ), while for the Australia study, there were nearly half of cases (1064 cases) with onset age less than 40 ( Macgregor et al., 2011 ). Previous studies have discussed the influence of birth cohort effect on melanoma incidence ( Jemal et al., 2001 ). In the present study, although we did not have enough information to compare the birth year of patients in both the discovery study and the Australia study, considering the different median age, we cannot exclude the possible influence of cohort effects on the validation results. We also compared the MAF of the two SNPs between the three studies. The MAF of rs1351383 in the controls of MD Anderson, GenoMEL and Australia melanoma studies was 0.38, 0.42 and 0.42 respectively, while the MAF of rs2127675 in the controls of the three studies was 0.32, 0.36 and 0.37 respectively. The increased allele frequency in the Australian population might have decreased the power of the replication studies, leading to the replication failure ( Greene et al., 2009 ). As an alternative to population replication, supportive evidence from functional assays may add additional biological plausibility to the weak associations and support identification of a true causal variant ( Khoury et al., 2009 ). Therefore, further studies on the mechanisms underlying the observed associations are warranted to assess whether these or other untyped variants within this region may directly contribute to melanoma risk. In conclusion, our findings suggest that functional genetic variants in PSMB9 may influence the development of melanoma by increasing gene expression. Further functional assays and replication in additional population are warranted to verify these results.

For the discovery phase, we used the melanoma GWAS dataset from the MD Anderson Cancer Center that has been recently described ( Amos et al., 2011 ). Briefly, the study participants consisted of 1804 non-Hispanic White patients and 1026 controls, as a part of an ongoing melanoma investigation, who were recruited at M.D. Anderson between March 1998 and August 2008. In addition, both GWAS data and risk-factor questionnaire data were available for 931 melanoma patients and 1,026 cancer-free controls (friends and relatives of other cancer patients visiting the clinics), who were genetically unrelated and frequency-matched on age and sex. The study protocol was approved by the Institutional Review Board at MD Anderson, and a written informed consent was obtained from all participants. One-time whole blood samples were used for DNA extraction by various methods (including Gentra, Qiagen, and phenol/chloroform). DNA samples were genotyped using the Illumina Omnil-Quad array and the genotypes were called using the BeadStudio algorithm at the John Hopkins University Center for Inherited Disease Research (CIDR). Standard quality control (QC) procedures were applied to both samples and SNPs which had been described in the published paper ( Amos et al., 2011 ). Briefly, SNPs were included, if they had a minor allele frequency (MAF) > 0.01, call rate ≥ 95%, and Hardy-Weinberg equilibrium in controls with P ≥ 1×10−5. We excluded the duplicated samples, related (IBD) samples or outliers identified by principle component analysis (PCA). One replication was performed in silico utilizing GWAS from the GenoMEL Consortium ( Barrett et al., 2011 ; Bishop et al., 2009 ). The GenoMEL GWAS utilizes samples collected from multiple centers across Europe and Israel in two phases. Phase 1 of the original GenoMEL GWAS consisted of samples collected from 8 centres across 6 different European countries. These were supplemented with controls from the Wellcome Trust Case-Control Consortium. Standard QC measures were applied to both samples and SNPs, giving a total of 1,353 cases and 3,571 controls. Phase 2 of the GenoMEL GWAS was collected across 10 centers (4 not in Phase 1) in 8 different European countries and Israel, supplemented again by samples from the Wellcome Trust Case Control Consortium. After QC, 1,450 cases and 4,047 controls remained. Most GenoMEL Phase 1 samples were genotyped on the llumina HumanHap300 BeadChip version 2 duo array (with 317k tagging SNPs), with the exception of the French cases, which were genotyped on the Illumina Humancnv370k array. The GenoMEL Phase 2 samples were genotyped on the Illumina 610k array. The Australian data used in Phase 1 were dropped, as these samples were included in the Australian GWAS as another replication set. SNP quality control were applied to each genotyping platform separately. SNPs were excluded, if they had a call rate < 97%, or Hardy-Weinberg equilibrium P < 1×10 −20 or recommendation for exclusion by Wellcome Trust Case-Control Consortium (WTCC). The imputed p-values from the Phase 1 data were available for replication. The second in silico replication included 2,168 melanoma cases selected from the Queensland, Australia study of Melanoma: Environment and Genetic Associations (Q-MEGA) and the Australian Melanoma Family Study (AMFS), for which 1242 patient samples were typed on the Omni1-Quad and 926 typed on the Hap610 arrays ( Macgregor et al., 2011 ). Three Australian Caucasian sample populations were used as controls (n=4,387), for which 431 were typed on the Omni1-Quad and 3956 were typed on the Hap610 arrays. Cases and controls were combined into a single dataset for quality control analysis (including principal component analysis for outlier removal) and imputation. SNPs were excluded for MAF < 0.01, call rate < 95%, or Hardy-Weinberg equilibrium in controls with P < 1×10 −6 . Imputation was based on the 1000 Genomes Project data, which helped recover the full sample size for SNPs that were only typed on a subset of the arrays ( Macgregor et al., 2011 ). In the re-analysis of M.D. Anderson Melanoma GWAS, genes involved in the M/G1 transition process of the cell cycle were selected based on the following criteria: genes that have been reported to be involved in the M/G1 transition process; genes that have been included in the M/G1 transition of cell cycle pathway in the Web of Amigo ( http://amigo.geneontology.org/ , as of 12/03/2010) ( Carbon et al., 2009 ); and genes that have been covered by the Illumina Omni1-Quad BeadChips (Illumina, San Diego, CA). As a result, there were a total of 1149 SNPs in 76 genes in the M/G1 transition of cell cycle pathway available from our GWAS database. The assignment of a SNP to a gene was defined by Illumina annotation file “Human 1M Quad_gene_annotation.txt”, which annotated all SNPs to their closest gene regardless how far a SNP is away from the gene. To search for putative functional SNPs for these with P < 0.01 obtained for the initial association analyses, we used SNPinfo ( http://snpinfo.niehs.nih.gov/snpfunc.htm ) ( Xu and Taylor, 2009 ) to identify any putative functional SNPs based on the HapMap phase II data. To obtain more confident prediction, we validated the functional findings of SNPinfo with other bioinformatics softwares. For example, for SNPs that were predicted to affect transcription factor binding sites in the promoter region, we further evaluated their effects on the transcription factor binding in TFSEARCH ( Heinemeyer et al., 1998 ) ( http://molsun1.cbrc.aist.go.jp/research/db/TFSEAR CH.html ). Only SNPs that were identified to be functional in two or more bioinformatics software were further validated by the laboratory functional assays. For replication, we selected the 18 most significantly associated SNPs ( P < 0.01) in the PSMB9 gene region with putative function and evaluated their associations with melanoma using the genotyping data of two additional GWASs from UK and Australia. Nuclear extracts from melanoma cell line A375 were prepared according to the method of Andrews and Faller ( Andrews and Faller, 1991 ). Complementary single-stranded oligonucleotides for rs2071480 of PMSB9 (5’-GCGCGCGGCGCTAACTTGTGTAGGGCAGATC −3’ for the T allele and 5’-GCGCGCGGCGCTAACGTGTGTAGGGCAGATC −3’ for the G allele) were biotin-labeled using the 3’-end biotin labeling kit (Thermo Scientific, Rockford, IL) and re-annealed before performing the DNA binding. The binding of DNA and protein was performed by using the LightShift Chemiluminescent EMSA kit (Thermo Scientific, Rockford, IL). The DNA-protein complexes were separated on 6% polyacrylamide gel, and the products were detected by Stabilized Streptavidin-Horseradish Peroxidase Conjugate (Thermo Scientific). The competition assays were performed 50-fold excess of unlabeled wild-type and mutation oligonucleotides, respectively. The intensity of shifted bands in scanned films was determined by using TotalLab™ program (Nonlinear Dynamics Ltd., UK). PCR fragments containing T allele of SNP rs2071480 were amplified from genomic DNA isolated from homozygous T carriers using the following primers: forward primer 5'-AA GCTAGC ATCTGAGAATCTCGGGAGCA −3', reverse primer 5'-TT AAGCTT GGTTTCCAACCTGGGACAG −3'. The PCR products were then cloned into the pGL3-Basic vector (Promega, Madison, WI, USA) between Nhe I and Hind III sites and verified by directly sequencing. The common allele G of rs2071480 was introduced into the recombination vector by QuickChange site-directed mutagenesis kit (Cat # 200518; Stratagene, La Jolla, CA) using the forward mutagenic primer 5’-TTTGCGCGCGGCGCTAAC G TGTGTAGGGCAGATCT-3’ and reverse mutagenic primer 5’-AGATCTGGCCTACACA C GTTAGCGCCGCCCGCAAA-3’ according to the manufacturer’s protocol. The clones containing the expected G allele of rs2071480 were verified by direct sequencing. Three cancer cell lines (A375 from melanoma, HeLa from cervical cancer and HCT116 from colon cancer) were placed on 24-well plates at 1.0×10 5 cells per well with DMEM or 1x RPMI 1640 culture medium containing 10% fetal bovine serum and allowed to grow for one day prior to transfection (50–70% confluence). Transfection experiments were performed using FuGENE HD (Invitrogen, Carlsbad, CA, USA). Each transfection was performed in triplicates. Dual Luciferase Kit (Promega) was used to detect the activity of firefly luciferase and Renilla luciferase. We also analyzed PSMB9 mRNA expression by genotypes of rs1351383, rs2127675, and rs2071480 based on the transcript expression profiling data of 270 lymphoblastoid cell lines from CEU and other HapMap samples (including 90 CEU samples, 90 CHB/JPT samples and 90 YRI samples) (NCBI GEO accession ID: GSE7792 ) ( Stranger et al., 2007 ). The expression data and genotyping data were available for rs1351383 in 81 CEU and 199 unrelated lymphoblastoid cell lines, for rs2127675 in 81 CEU and 198 unrelated lymphoblastoid cell lines, and for rs2071480 in 84 CEU and 201 unrelated lymphoblastoid cell lines. The distribution differences of demographic variables and known risk factors between melanoma patients and controls were assessed by the χ 2 test. The associations between alleles or genotypes of each tagSNP and melanoma risk were primarily evaluated using the allelic test in PLINK1.07 ( Purcell et al., 2007 ). Benjamini and Hochberg FDR method was used for the multiple testing corrections ( Benjamini and Hochberg, 1995 ). LD patterns among tagSNPs were evaluated by Haploview ( Barrett et al., 2005 ). We adjusted for the five largest principal components of genetic variation to control potential effects of population structure. For the luciferase assay, significant differences between groups were determined by Student’s t test. For the in silico meta-analysis, betas from each study were combined using the inverse variance method. T-test and Wilcoxon-Mann-Whitney test were applied to compare the difference in mRNA expression levels of PSMB9 between different genotypes or different tissues. Unless specified otherwise, all other statistical analyses were performed using SAS 9.1.3 (SAS Institute Inc., Cary, NC), and all statistical tests were two sided, with a P < 0.05 set as the level of statistical significance.