Inherited variations in human pigmentation-related genes modulate cutaneous melanoma risk and clinicopathological features in Brazilian population

Ultraviolet light exposure and cutaneous pigmentation are important host risk factors for cutaneous melanoma (CM), and it is well known that inherited ability to produce melanin varies in humans. The study aimed to identify single-nucleotide variants (SNVs) on pigmentation-related genes with importance in risk and clinicopathological aspects of CM. The study was conducted in two stages. In stage 1, 103 CM patients and 103 controls were analyzed using Genome-Wide Human SNV Arrays in order to identify SNVs in pigmentation-related genes, and the most important SNVs were selected for data validation in stage 2 by real-time polymerase-chain reaction in 247 CM patients and 280 controls. ADCY3 c.675+9196T>G, CREB1 c.303+373G>A, and MITF c.938-325G>A were selected for data validation among 74 SNVs. Individuals with CREB1 GA or AA genotype and allele “A” were under 1.79 and 1.47-fold increased risks of CM than others, respectively. Excesses of CREB1 AA and MITF AA genotype were seen in patients with tumors at Clark levels III to V (27.8% versus 13.7%) and at III or IV stages (46.1% versus 24.9%) compared to others, respectively. When compared to others, patients with ADCY3 TT had 1.89 more chances of presenting CM progression, and those with MITF GA or AA had 2.20 more chances of evolving to death by CM. Our data provide, for the first time, preliminary evidence that inherited abnormalities in ADCY3, CREB1, and MITF pigmentation-related genes, not only can increase the risk to CM, but also influence CM patients’ clinicopathological features.

. (A) Pigmentation regulation by alpha-melanocyte stimulating hormone (MSHα) and G-proteins from melanocortin receptor 1 (MC1R): MSHα/MCR1 can trigger the activation of the adenylate cyclase (ADCY) and 3′-5′-cyclic adenosine monophosphate AMP (cAMP). The cAMP signals activate the protein kinase A (PKA) that phosphorylates and activates cAMP responsive element binding protein (CREB) transcription factor, which induces the expression of melanocyte inducing transcription factor (MITF) and induction of proliferation and differentiation of melanocytes. G-proteins: β, γ, and α; ATP: adenosine triphosphate. (B) Kaplan-Meier (K-M) curves for progression-free survival according to ADCY3 c.675+9196T>G genotypes, where patients with TT genotype presented lower survival than those with TG or GG genotype. (C) K-M curves for melanoma-specific survival according to MITF c.938-325G>A genotypes, where patients with GA or AA genotype presented lower survival than those with GG genotype.
Data and specimen collection. Clinical information of individuals (age at diagnosis, gender, skin color, skin phototype, sun exposure, type of sun exposure, and number of nevi) was obtained by specific questionnaires. Skin phototype was defined using reported criteria 23 . Individuals exposed to the sun for more than 2 h per day and for more ten years were considered positive for sun exposure 23 . Sun exposure was classified as intermittent in cases of recreational activities performed less than 50% of the week or holidays, or chronic, activities at home or work under sunlight exposure during more than 50% of the time 24 .
The diagnosis of CM was established by histopathological evaluation of tumor fragments embedded in paraffin and stained with hematoxylin and eosin. Pathological aspects of the tumor (tumor location, histological type, Breslow thickness, Clark level, and tumor stage) were obtained from medical records of patients 25 . Tumor stage was identified using the TNM classification of the American Joint Committee on Cancer, where T describes the size of tumor, N describes spread of tumor to nearby lymph nodes, and M describes distant metastasis 2 . Patients with desmoplastic, acro-lentiginous and amelanotic melanomas were excluded from the study.
Surgical excision (n = 217) was the primary treatment for patients with localized tumor 26 . Sentinel lymph node biopsy (n = 41) was recommended in patients with tumor measuring more than 1 mm (mm) and lymphadenectomy (n = 24) was performed in patients with clinically positive lymph nodes or lymph nodes with tumor infiltration on histopathological evaluation. Patients with operable single metastasis or relapse (n = 30) underwent surgical resection 27 . Those patients with inoperable relapse or multiple metastases (n = 30) received chemotherapy with dacarbazine 28 . Radiotherapy was also used in the local treatment of patients with surgical impossibility (n = 4), particularly in bleeding lesions, bone or brain metastases 29 .
Stage 1: screening of SNVs, candidate genes choice and SNVs selection. DNA from leukocytes of peripheral blood of CM patients and controls were genotyped for a total of 906,660 SNVs using the Affymetrix Genome-Wide Human SNV Arrays 6.0 (AFFYMETRIX, USA), according to the manufacturer's recommended protocols. The intensities resulting from the arrays scanning process were made available via CEL files, one per DNA sample with total quality control higher than 90% (AFFYMETRIX, USA). Tools from the Bioconductor (https ://www.bioco nduct or.org) were used to process the CEL files. The genotyping was performed applying the corrected robust linear mixture model (crlmm) algorithm 30 .
The genes previously reported as involved in the pigmentation pathway were selected for study. The pathway analysis was performed using the Database for Annotation, Visualization and Integrated Discovery (https ://david .ncifc rf.gov 31 and Kyoto Encyclopedia of Genes and Genomes pathway maps (https ://www.kegg.jp) 32 .
Each pigmentation related-gene was analyzed using the in silico method by the Human Splicing Finder algorithm (version 3.1) (https ://www.umd.be/HSF3/index .html 33 in order of identifying SNVs in splicing regulatory sequences. For analysis, wild-type (ancestral) allele was taken as reference. SNVs showing deviation from the HWE and those with the minor allele frequency less than 10% were excluded from the selection 34 . SNVs that potentially alter expression or function of the encoding proteins 13,14 were selected for further validation.
Stage 2: validation of selected SNVs in risk and characteristics of melanoma. DNA from leukocytes of peripheral blood of CM patients and controls was analyzed by real-time polymerase chain reaction with TaqMan SNV genotyping assays (APPLIED BIOSYSTEMS, USA) for ADCY3 (rs11900505, assay ID: C_7868411_20), CREB1 (rs10932201, assay ID: C_2859093_20) and MITF (rs7623610, assay ID: C_29012190_10) SNVs, following manufacturer instructions. Twenty percent of genotype determinations were carried out twice in independent experiments with 100% of concordance.

Statistical analysis.
Association between disease statuses, CM patients versus controls, and genotypes for study's stage 1 was performed using logistic regression model, and analyses were adjusted by age at diagnosis, skin color, and sun exposure. SNVs that presented raw p-values below the 0.001 thresholds were selected for further inspection. These analyses were implemented in R software (version 3.3.0) (https ://www.r-proje ct.org).
The Hardy-Weinberg equilibrium (HWE) was tested using chi-square (χ 2 ) statistics for the goodness-to-fit test, and logistic regression model served to obtain age, skin color, sun exposure, and number of nevi statusadjusted crude odds ratios (ORs) with 95% confidence intervals (CI) in comparisons evolving patients and controls for study's stage 2. To evaluate the robustness of risk estimates, the false discovery rate (FDR) was computed, which reflects the expected ratio of false-positive findings to the total number of significant findings; the differences revealed were considered statistically significant at FDR values < 0.05 25 . χ 2 and Fisher's exact tests were used to evaluate associations between clinicopathological features and genotypes of selected SNVs. Bonferroni method was used to adjust values of multiple comparisons in patients stratified by tumor aspects 36 . For ADCY3, CREB1 and MITF expression analysis, data sets were probed for normality using Shapiro-Wilk's test. Because data sets assume normal distribution, analysis of variance performed comparisons of groups 36 .
For survival analysis, the progression-free survival (PFS) was calculated from the date of surgery until the date of first recurrence, or the date of progression of disease, or the date of death by any cause, or the date of last follow-up. The melanoma-specific survival (MSS) was calculated from the date of diagnosis until the date of death by the disease or last follow-up. PFS and MSS were calculated using Kaplan-Meier estimates, and differences between survival curves were analyzed by log-rank test 25 . The impact of age at diagnosis, gender, tumor location, Breslow thickness 37 , Clark level, TNM stage and genotypes of each analyzed SNV in survival of patients were evaluated using univariate Cox proportional hazards ratio (HR) regression. At a second time, all variables with p < 0.20 were included in the multivariate Cox regression. The significant results of Cox analysis were internally validated using a bootstrap resampling study to investigate the stability of risk estimates (1,000 replications) 25 .
All tests were done using the SPSS 21.0 software (SPSS INCORPORATION, USA). Significance was two sided and achieved when p values were ≤ 0.05.

Study population. The clinicopathological features of patients and the clinical features of controls included
in stage 1 and stage 2 of the study are presented in Table 1. Controls were younger than patients, and CM patients presented more white skin color, referred more sun exposure and presented more nevi than controls, and all differences were corrected in comparisons involving patients and controls by appropriate statistical analysis. Similar clinicopathological features were observed in patients and controls analyzed in both stages of the study.
The most significant melanoma associated SNVs identified in stage 1 (p < 0.0001) are presented in Table S1 Supplement. Seventy-four SNVs in 28 human pigmentation-related genes were found to be involved with CM risk (Table S2 Supplement).
In accord with results of the in silico analysis, the variant allele "G" of ADCY3 c.675+9196T>G may abolish a potential branch point site and an exonic splicing enhancer (ESE), and this variation may create a site of ligation for SRp55 and 9G8 splicing proteins. Besides, new sites for an exon-identity element (EIE) and an intron-identity element (IIE) may be created. The variant allele "A" of CREB1 c.303+373G>A may abolish a splice donor site (5′ end of the intron), an exonic splicing silencer (ESS), an IIE site, and a binding site for the hnRNP A1, and this variation may create a new branch point, an EIE site, and a putative exonic splicing enhancer. The variant allele "A" of MITF c.938-325G>A may create a splice donor site, an ESS, and a binding site for the hnRNP A1, and this variation may break a potential branch point site, an EIE, and silencer motifs and an IIE site (Table S3  Supplement). ADCY3 c.675+9196T>G, CREB1 c.303+373G>A, and MITF c.938-325G>A SNVs were selected for analysis in stage 2 of the study due to their potential effects on encoding proteins 13 www.nature.com/scientificreports/ CREB1 GA or AA genotype and allele "A" were more common in patients than in controls; carriers of the above genotypes and allele were under 1.79 and 1.47-fold increased risks for CM than those with the GG genotype and allele "G", respectively ( Table 2). No associations between ADCY3, CREB1 and MITF SNVs combined genotypes were seen in CM patients and controls (Table S4 Supplement).
No associations of studied SNVs genotypes were seen in CM patients stratified by age, gender, and skin color (Table S5), phototype, sun exposure, type of sun exposure, and number of nevi (Table S6). However, CREB1 AA genotype was more common in patients with tumors located in limbs than in head or trunk (31.7% versus 15.9%, p = 0.009) and tumors with Clark levels III to V than in those with tumors of I or II Clark levels (27.8% versus 13.7%, p = 0.012), and MITF AA genotype was more common in patients with III or IV tumor stage than in those with tumors at 0 to II stages (46.1% versus 24.9%, p = 0.007). These results were significant even after Bonferroni correction (corrected p value: 0.0125) (Table 3).

Association of clinicopathological aspects and genotypes with patients' survival.
We obtained consisted survival data from 210 CM patients. The median follow-up time of patients enrolled in the survival analysis was 97 months (range 5-228 months). The patient's final status was established on January 2020, when 136 patients were alive (132 without disease, 4 with disease) and 74 patients had died ( (Table 4).

Discussion
In this study, we investigated and identified intronic SNVs ADCY3 c.675+9196T>G, CREB1 c.303+373G>A, and MITF c.938-325G>A in pigmentation-related genes in association with CM risk and clinicopathological features. After screening SNVs (stage 1), we found more than 6,000 SNVs associated with CM risk in introns of genes, according to previous studies [14][15][16][17][18] , and we selected three SNVs involved in the splicing regulatory sequences of pigmentation-related genes for data validation, due to their potential roles in determining abnormalities in production and/or function of the respective encoded proteins 13,14 .
In fact, previous GWAS have shown that the majority of disease-associated variants reside in the non-coding regions of the genome, suggesting that gene regulatory changes contribute to disease risk 18 . On the other hand, splicing comprises a two-step reaction of intron removal and exon ligation and is essential for gene expression: pre-mRNA splicing is catalyzed by the spliceosome, a large complex of ribonucleoproteins (RNPs), and this complex recognizes the target sequences and assembles on pre-mRNA 20 .
After SNVs validation (stage 2), we observed that CREB1 GA or AA genotype and allele "A" were more common in CM patients than in controls, and that individuals with referred genotypes and allele were under 1.79 and 1.47-fold increased risks of CM than others, respectively.
CREB1 was highly expressed in tumor cells, such as human gastric cell lines and knockdown of CREB1 inhibited human gastric cancer cells growth 38 . CREB1 has also been seen as an important gene in CM development 8 , and analysis of common network from cancer type-specific RNA-Seq co-expression data showed CREB1 as a melanoma-associated gene 39 . To the best of our knowledge, there are no previous studies focusing the roles of ADCY3 c.675+9196T>G, CREB1 c.303+373G>A, and MITF c.938-325G>A SNVs in risk of CM, and therefore the association of CREB1 GA or AA genotype and allele "A" with CM risk seen in the present study is a new finding. The search for potential splicing regulatory elements using in silico algorithm in this study indicated Table 3. ADCY3 c.675 + 9196 T > G, CREB1 c.303 + 373G > A, and MITF c.938-325A > G genotypes in 247 patients with cutaneous melanoma stratified by tumor features. Values are expressed as number and percentage. * The numbers of patients were not the same included in the study (n = 247) because no consistent information could be obtained from some individuals. P values < 0.05 are presented in bold letters. ** Significant even after Bonferroni correction for multiple comparisons (corrected p value = 0.0125). www.nature.com/scientificreports/ that gene variants induce the creation or abrogation of binding sites 33 . The allele "A" of CREB1 c.303+373G>A may alter binding sites for splicing elements, such as the hnRNP A1, and possibly increases CREB1 activity due to altered efficiency of splicing 19,20,33 . Since CREB1 is a transcription factor that stimulates the MITF activity, the increase of its activity may in turn increase MITF activity, having proliferation of abnormal melanocytes and increased risk for CM as consequence 8 . When genotypes were analyzed in patients stratified by clinicopathological aspects, we noted that CREB1 AA variant genotype was more common in patients with tumors located in limbs than in patients with tumors located in head or trunk and with tumors at Clark level III to IV than in patients with tumors at I or II level. In addition, an excess of MITF AA genotype was found in patients with tumors at stage III or IV than in those with stage I or II tumors.
It was already described that acquisition of metastatic phenotype in CM involved the gain in expression of CREB/activating transcription factor-1 (CREB/ATF-1) 40 and MITF amplification 41 . However, how far our knowledge reaches, this study is the first to describe the influence of CREB1 c.303+373G>A and MITF c.938-325G>A SNVs on clinicopathological features of CM. Indeed, CREB1 promotes tumorigenesis by increasing cell migration, proliferation, and invasiveness, through its effects on the MITF pathway 8 . The in silico analysis showed that the variant allele "A" of MITF c.938-325G>A may create a site of ligation for splicing factors, including the Table 4. Clinicopathological aspects and genotypes in survival of 210 cutaneous melanoma patients. HR, hazard ratio; CI, confidence interval; NC, characteristic not computed in multivariate analysis. * The total numbers of individuals differed from the total quoted because it was not possible to obtain consistent information about characteristics in some individuals. a P bootstrap = 0.01; b P bootstrap = 0.001. c P bootstrap = 0.001. d P bootstrap = 0.02; e P bootstrap = 0.002. f P bootstrap = 0.006. g P bootstrap < 0.0001 in multivariate analysis. Significant differences between groups are presented in bold letters. www.nature.com/scientificreports/ hnRNP A1, possibly determining increase in gene expression 33 . Thus, we postulate that CREB1 AA and MITF AA genotypes may increase abnormal melanocytes proliferation and consequently improve aggressiveness of CM. We also noted that ADCY3 c.675+9196 TT genotype was associated with shorter PFS while MITF GA or AA genotype was associated with shorter MSS in CM patients, when compared to the remaining genotypes.
Up-regulation of ADCY3 increased the tumorigenic potential of gastric cells 42 and predicted shorter overall survival in patients with pancreatic cancer 43 . Overexpression of ADCY2 was previously associated with aggressive behavior of CM 44 , and MITF amplification predicted worst survival of CM patients 41 . The in silico analysis showed that the ADCY3 c.675+9196T>G variant may alter sites of ligation for splicing factors, including the SRp55 and 9G8, with a possible increase in the efficiency of splicing and gene expression 19,20,33 . Since ADCYs participate in CREB activation, and CREB regulates the expression of MITF 8 , the increase in ADCY3 activity in CM patients with the TT genotype may favor proliferation of abnormal melanocytes leading to relapse or death by CM effects. Again, the possible increased activity of MITF in patients with GA or AA genotype may have contributed to this clinical unfavorable outcome.
It is also worth to comment that pigmentation-related genes have been seen as potential therapeutic targets. Previous studies showed that increased ADCY expression generated resistance to MAPK inhibitions and up regulates MITF in melanoma cells 8 , and the suppression of MITF expression by the CH6868398 agent caused melanoma cell growth inhibition 45 . Inhibition of p300 acetyltransferase transcriptional coactivator of MITF by p300/CBP complex had growth inhibitory effects in melanoma cells expressing MITF 46,47 , and Kazinol U reduced melanogenesis by inhibition of MITF in melanoma cells 48 . Since response to new agents depends on ADCY3 and MITF expressions, it is possible that patients with distinct genotypes of these genes present differentiated responses to therapies.
At this time, we draw attention to the fact that no differences in ADCY3, CREB1 and MITF expressions were identified in leukocytes of peripheral blood of individuals with the distinct genotypes of ADCY3, CREB1 and MITF SNVs. It is possible that the sample size evaluated was not enough to identify differences in gene expression among individuals or, alternatively, these variants may determine gene expression abnormalities only in tumor tissue or only protein functional changes.
In summary, we described for the first time the potential importance of ADCY3 c.675+9196T>G, CREB1 c.303+373G>A, and MITF c.938-325G>A SNVs in the pigmentation-related genes in CM risk and clinicopathological features in Brazilian individuals. We recognize that the present study has limitations: it was conducted on a relatively small number of individuals and only quantitative analysis of gene expression in normal leukocytes was performed. Thus, we believe that our results will require confirmation in a further larger epidemiological study in our population and others, and quantitative and functional analyses of ADCY3, CREB1 and MITF SNVs in melanoma cells. If these findings are confirmed, they might help to identify individuals with high risk for CM who deserves to receive additional recommendations for CM prevention and early tumor detection and/ or differentiated treatment, perhaps including the targeting lineage specific MC1R signalizing pathway agents.