https://orcid.org/0000-0001-9294-6546
Abstract
The 10q24.33 locus is known to be associated with susceptibility to cutaneous malignant melanoma (CMM), but the mechanisms underlying this association have been not extensively investigated.
We carried out an integrative genomic analysis of 10q24.33 using epigenomic annotations and in vitro reporter gene assays to identify regulatory variants. We found two putative functional single nucleotide polymorphisms (SNPs) in an enhancer and in the promoter of OBFC1, respectively, in neural crest and CMM cells, one, rs2995264, altering enhancer activity. The minor allele G of rs2995264 correlated with lower OBFC1 expression in 470 CMM tumors and was confirmed to increase the CMM risk in a cohort of 484 CMM cases and 1801 controls of Italian origin. Hi-C and chromosome conformation capture (3C) experiments showed the interaction between the enhancer-SNP region and the promoter of OBFC1 and an isogenic model characterized by CRISPR-Cas9 deletion of the enhancer-SNP region confirmed the potential regulatory effect of rs2995264 on OBFC1 transcription. Moreover, the presence of G-rs2995264 risk allele reduced the binding affinity of the transcription factor MEOX2. Biologic investigations showed significant cell viability upon depletion of OBFC1, specifically in CMM cells that were homozygous for the protective allele. Clinically, high levels of OBFC1 expression associated with histologically favorable CMM tumors. Finally, preliminary results suggested the potential effect of decreased OBFC1 expression on telomerase activity in tumorigenic conditions.
Our results support the hypothesis that reduced expression of OBFC1 gene through functional heritable DNA variation can contribute to malignant transformation of normal melanocytes.
Introduction
Cutaneous malignant melanoma (CMM) is a cancer of transformed neural crest (NC) derived melanocytes, pigment-producing cells, where both genetic and environmental factors are involved. The risk of CMM is largely modified by factors influencing individual sensitivity to UV radiation and sunlight exposure; sunburns during childhood in particular are a major behavioral risk factor (1). The heritability, that is, the contribution of genetic factors to CMM risk, has been estimated to be ~40–50% (1). Besides deleterious pathogenic variants in well-established melanoma susceptibility genes such as CDKN2A and CDK4, or in other melanoma risk genes of more recent identification, such as BAP1, POT1, TERT, ACD, TERF2IP, MITF, ATM (1,2,3,4) which confer a high/moderate CMM risk; common variants with low effect size are likely to be also involved in melanoma susceptibility.
Genome Wide Association Studies (GWAS) approach have led to insights into the architecture of disease susceptibility through the identification of novel disease-causing genes and mechanisms improving our knowledge of the complex disease etiology (5). Over the past decade, more than 430 cancer associated common variants at 262 distinct genomic regions have been successfully identified by GWAS (6). Most of these alterations resides in non-coding portions of the genome and may have regulatory consequences on cancer susceptibility (7). However, the functional role of the identified risk loci in cancer pathogenesis remains poorly investigated.
Large melanoma GWASs have identified several loci associated with CMM risk in the general population: PARP1, SLC45A2, TYR, MC1R, ASIP, CDKN2A-MTAP, CYP1B1, PLA2G6, TERT, ATM, ARNT-SETDB1, CDKAL1, OCA2, CCND1, AGR3, CASP8, FTO, CDK10, TMEM38B, OBFC1 and MX2 (8–10). Despite GWAS facilitating the initial identification of a risk locus, this approach presents some limitations mainly due to the difficulty to discern the causal variants due to linkage disequilibrium (LD). Moreover, common genetic loci that usually are likely hidden among signals discarded by the multiple testing correction represents a restrictive step of GWAS analytical process (5).
Post-GWAS strategies are trying to overcome these limitations by leveraging different approaches such as imputation analysis, next-generation sequencing (NGS), fine mapping and cross-phenotype meta-analysis (CPMA) (11). Leveraging the concept of pleiotropy, a CPMA of CMM and nevus GWAS demonstrated that several risk loci might act through nevus development, in line with clinical evidence (12). Moreover, a recent GWAS meta-analysis identified additional risk loci for CMM including some associated with telomere length or located near prominent telomere maintenance genes, including POT1, TERC, RTEL1, MPHOSPH6 and OBFC1 (10). Recently, we performed a cross-disease meta-analysis of neuroblastoma and CMM GWAS, in which we found several neuroblastoma-CMM cross-associated loci. Among these, we further confirmed the association of the previously identified (8) 10q24.33 CMM risk locus (index SNP rs11591710) (13). However, at this risk locus, the identification of the causal genetic variants (i.e. those that actually contribute to the development of CMM) and the detection of the genes whose function is influenced by the same causal variants remain to be established.
Here, we performed an integrative genomic analysis of the 10q24.33 CMM risk locus that led us to identify a gene-regulatory variant (rs2995264) with enhancer features within an intronic region of the OBFC1 gene. In silico and in vitro studies demonstrated that the G risk allele of rs2995264 SNP correlated with a decreased expression of OBFC1 gene suggesting its role as tumor-suppressor in CMM, which has been experimentally confirmed in melanoma cell lines. Despite the need for follow-up functional studies, a preliminary model of pathways potentially important for the CMM development is emerging through this approach.
Results
Integrative genomic analysis and functional annotation of candidate SNPs
We previously found the risk locus 10q24.33, with the rs11591710 index SNP in an intronic region of OBFC1, reaching genome-wide significance of association in a cross-disease meta-analysis of neuroblastoma and CMM (13). Particularly, the minor allele of rs11591710 was associated with increased risk of developing CMM (13). Therefore, we decided to further functionally analyze the 10q24.33 locus. To highlight potentially functional variants, we annotated 32 SNPs in LD (0.5 < r2 < 1) with the lead SNP rs11591710 with the melanoma super-enhancer (SE), enhancer and promoter histone marks obtained by an analysis of H3K27ac ChIP-Seq data (Supplementary Material, Supplementary Information) derived from seven human CMM cell lines and two hNC cell lines (GSE90683) deposited in GEO database. To prioritize causal functional variants in CMM, we first selected those SNPs overlapping at least one histone marker in either CMM or hNCC. The analysis pipeline allowed us to classify acetylation peaks as SE, enhancers or promoters (Supplementary Material, Table S1; see Supplementary Material, Supplementary Information). The three top SNPs with highest number of histone markers annotations (Supplementary Material, Table S1) rs34685262 (r2 = 0.77 with the index SNP), rs2995264 (r2 = 0.68 with the index SNP) and rs35176048 (r2 = 0.68 with the index SNP) fell in the same H3K27ac peak, located in the intron 3 of OBFC1, called as SE or enhancer in nine or eight cell lines whereas rs4387287 (r2 = 0.65 respect to the index SNP) was annotated in the promoter region of the same gene in eight cell lines (Fig. 1A and Supplementary Material, Table S1). Among these three top SNPs, the rs35176048 SNP was located in an H3K27ac but not in an H3K4me1 peak in melanoma and other ENCODE cell lines (Supplementary Material, Table S1). Beyond potentially cis-regulatory variants, we also noted that two potentially functional variants, rs10786775 and rs2487999, are missense variants, however these were predicted to be benign (Supplementary Material, Table S1). Based on these observations rs35176048, rs10786775, rs2487999 SNPs were not considered for further analysis. Instead, the SNPs rs34685262, rs2995264 and rs4387287, located in a well-defined melanoma regulatory DNA region (Fig. 1A), were tested for the induction of enhancer activity through luciferase reporter assay in HEK293T and A375 melanoma cells in order to validate their regulatory properties (Fig. 1B). Only the construct containing rs2995264-G risk allele induced a significant decrease of enhancer activity compared to the construct containing rs2995264-A reference allele (Fig. 1B).
H3K27ac activity at 10q24.33, allele-specific functional SNP, and SNP genomic interactions. (A) From the top to the bottom of the figure, it is showed the list of SNPs in LD with the lead SNP rs11591710 (zoom-in showing 83 154 bp), the H3K27ac data of NCC cell lines (GSE90683) and of CMM cell lines (GSE75352 and GSE82332). Functionally relevant SNPs that were further validated are highlighted by colored rectangles. rs2995264 (red rectangle) showed significant results. (B) Luciferase report gene assay carried out in HEK293T and A375 cells shows that only for rs2995264 the GG altered genotype correlates with a lower luciferase activity. Luciferase activity of the constructs was relative to PGL3 empty vector. The datapoints represent the mean of technical duplicates of each of the three independent experiments. P-values obtained by two-tailed T-test. (C) RNAseq profiling of primary tumors demonstrates that lower OBFC1 expression correlates with GG-rs2995264 risk genotypes. The dark line in the middle of the boxes is the median. The bottom of the box indicates the 25th percentile. The top of the box represents the 75th percentile. The T-bars that extend from the boxes indicate the minimum and maximum values (P-value = 0.002; two-tailed T-test). (D) Genomic interactions within the interval chr10:105111147-106 194 301 (genomic coordinates, hg19) obtained by HiC of the COLO-829 cell line are showed on the right panel. Black bordered triangles represent TADs. Bolded triangles highlight interactions between the genomic bin containing the rs2995264 (chr10:105660000-105 670 000) and the regions up-stream (distance = 2058 bp) or down-stream (distance = 7327 bp) of OBFC1 (see Supplementary Material, Table S3). Significantly interacting genes are underlined in red.
In GWAS, including 12 874 melanoma cases and 23 203 controls (6), the allele G of potential causal SNP rs2995264 is associated with CMM risk (P = 8.5 × 10−7, OR = 1.14) (Supplementary Material, Table S2). We sought to replicate this genetic association (using the SNP rs2488001 in high LD, r2 = 0.98, with rs2995264) in an independent cohort of 484 CMM cases and 1801 controls of Italian origin performing PCR-based genotyping. These data confirmed the association between the minor allele G to increased risk of CMM onset in the Italian population (P = 0.01, OR = 1.30) (Table 1). Moreover, to investigate whether more than one association signal may exist at 10q24.33, we conditioned our analysis of locus 10q24.33 on the potential causal SNP rs2995264 using summary statistics of CMM GWAS (8). We confirmed that no evidence for a separate association signal was observed at 10q24.33 locus (Supplementary Material, Fig. S1).
The genetic association of rs2995264 SNP in Italian melanoma cases and controls
| Healthy controls | Melanoma patients | P | OR (CI:95%) | |||
|---|---|---|---|---|---|---|
| n = 1801 | n = 484 | |||||
| ^rs2995264 | n | % | n | % | ||
| Genotype | ||||||
| AA | 1394 | 77.4 | 355 | 73.3 | - | - |
| AG | 385 | 21.4 | 113 | 23.4 | 0.25 | 1.15 (0.90–1.46) |
| GG | 22 | 1.2 | 16 | 3.3 | 0.001 | 2.85 (1.48–5.49) |
| Allele | ||||||
| A | 3173 | 88.1 | 823 | 85 | - | |
| G | 429 | 11.9 | 145 | 15 | 0.011 | 1.30 (1.06–1.59) |
| Dominant | ||||||
| AA/AG | 1779 | 98.8 | 468 | 96.7 | - | - |
| GG | 22 | 1.2 | 16 | 3.3 | 0.0015 | 2.76 (1.38–3.80) |
| Recessive | ||||||
| AA | 1394 | 77.4 | 355 | 73.4 | - | - |
| AG/GG | 407 | 22.6 | 129 | 26.6 | 0.062 | 1.24 (0.98–1.56) |
| Healthy controls | Melanoma patients | P | OR (CI:95%) | |||
|---|---|---|---|---|---|---|
| n = 1801 | n = 484 | |||||
| ^rs2995264 | n | % | n | % | ||
| Genotype | ||||||
| AA | 1394 | 77.4 | 355 | 73.3 | - | - |
| AG | 385 | 21.4 | 113 | 23.4 | 0.25 | 1.15 (0.90–1.46) |
| GG | 22 | 1.2 | 16 | 3.3 | 0.001 | 2.85 (1.48–5.49) |
| Allele | ||||||
| A | 3173 | 88.1 | 823 | 85 | - | |
| G | 429 | 11.9 | 145 | 15 | 0.011 | 1.30 (1.06–1.59) |
| Dominant | ||||||
| AA/AG | 1779 | 98.8 | 468 | 96.7 | - | - |
| GG | 22 | 1.2 | 16 | 3.3 | 0.0015 | 2.76 (1.38–3.80) |
| Recessive | ||||||
| AA | 1394 | 77.4 | 355 | 73.4 | - | - |
| AG/GG | 407 | 22.6 | 129 | 26.6 | 0.062 | 1.24 (0.98–1.56) |
^The rs2488001 SNP, in completed LD (r2 = 0.98) with our causal SNP rs2995264, was typed since the SNP TaqMan assay was not available.
OR: odd ratio; CI: confidence interval.
The significant results are reported in bold.
The genetic association of rs2995264 SNP in Italian melanoma cases and controls
| Healthy controls | Melanoma patients | P | OR (CI:95%) | |||
|---|---|---|---|---|---|---|
| n = 1801 | n = 484 | |||||
| ^rs2995264 | n | % | n | % | ||
| Genotype | ||||||
| AA | 1394 | 77.4 | 355 | 73.3 | - | - |
| AG | 385 | 21.4 | 113 | 23.4 | 0.25 | 1.15 (0.90–1.46) |
| GG | 22 | 1.2 | 16 | 3.3 | 0.001 | 2.85 (1.48–5.49) |
| Allele | ||||||
| A | 3173 | 88.1 | 823 | 85 | - | |
| G | 429 | 11.9 | 145 | 15 | 0.011 | 1.30 (1.06–1.59) |
| Dominant | ||||||
| AA/AG | 1779 | 98.8 | 468 | 96.7 | - | - |
| GG | 22 | 1.2 | 16 | 3.3 | 0.0015 | 2.76 (1.38–3.80) |
| Recessive | ||||||
| AA | 1394 | 77.4 | 355 | 73.4 | - | - |
| AG/GG | 407 | 22.6 | 129 | 26.6 | 0.062 | 1.24 (0.98–1.56) |
| Healthy controls | Melanoma patients | P | OR (CI:95%) | |||
|---|---|---|---|---|---|---|
| n = 1801 | n = 484 | |||||
| ^rs2995264 | n | % | n | % | ||
| Genotype | ||||||
| AA | 1394 | 77.4 | 355 | 73.3 | - | - |
| AG | 385 | 21.4 | 113 | 23.4 | 0.25 | 1.15 (0.90–1.46) |
| GG | 22 | 1.2 | 16 | 3.3 | 0.001 | 2.85 (1.48–5.49) |
| Allele | ||||||
| A | 3173 | 88.1 | 823 | 85 | - | |
| G | 429 | 11.9 | 145 | 15 | 0.011 | 1.30 (1.06–1.59) |
| Dominant | ||||||
| AA/AG | 1779 | 98.8 | 468 | 96.7 | - | - |
| GG | 22 | 1.2 | 16 | 3.3 | 0.0015 | 2.76 (1.38–3.80) |
| Recessive | ||||||
| AA | 1394 | 77.4 | 355 | 73.4 | - | - |
| AG/GG | 407 | 22.6 | 129 | 26.6 | 0.062 | 1.24 (0.98–1.56) |
^The rs2488001 SNP, in completed LD (r2 = 0.98) with our causal SNP rs2995264, was typed since the SNP TaqMan assay was not available.
OR: odd ratio; CI: confidence interval.
The significant results are reported in bold.
Based on the above-reported results, we decided to further functionally characterize the SNP rs2995264, which was the most significant functional SNP at 10q24.33 and located in an enhancer region of CMM and hNCC cell lines.
Evaluation of rs2995264 allele-specific enhancer activity toward OBFC1 promoter
To evaluate functional role of rs2995264, we verified if the candidate variant affected gene expression, by performing cis-expression Quantitative Trait Loci (eQTL) analysis (Supplementary Material, Supplementary Information). The analysis of gene expression variation using 470 RNAseq and SNP array data of CMM tumors from TCGA project demonstrated that the SNP rs2995264 altered expression of OBFC1 gene. Particularly, the presence of the G risk allele significantly correlated with decreased OBFC1 mRNA expression (P = 0.002, Fig. 1C). These results were further confirmed in skin tissue (P = 0.03) and cultured fibroblasts (P = 0.00005) (GTEx portal data, Supplementary Material, Fig. 2) and by luciferase reporter assay (Fig. 1B). Since rs2995264 affected OBFC1 expression and based on the previous evidences, we hypothesized that this SNP is located in an enhancer that physically interacts with the OBFC1 promoter. To demonstrate this assumption, we interrogated the public database Enhancer Atlas (31) that provides a list of long-range chromatin interaction partners for the queried locus obtained from 105 different human cell/tissue types. According to the analysis of Enhancer Atlas database, in the human foreskin tissue, the OBFC1 promoter interacts with the enhancer where the SNP rs2995264 is located (Supplementary Material, Fig. S3). Subsequently, we confirmed this interaction using HiC sequencing data obtained on the COLO-829 CMM cell line (Fig. 1D). Our analysis showed that the genomic bin (10 Kb) containing rs2995264 significantly interacted with a total of 11 genes. In particular, the SNP rs2995264 strongly interacted with both the up-stream (distance = 2058 bp; FDR = 1.76 × 10−38) and the down-stream (distance = 7327 bp; FDR = 1.26 × 10−16) regions of OBFC1 (Fig. 1D and Supplementary Material, Table S3).
To further validate physical interactions between polymorphic enhancer containing rs2995264 and OBFC1 promoter, we performed Chromosome Conformation Capture (3C) analysis in A375 and CJM melanoma cells. Sequences that are held nearby by genic regulation factors in chromosomal structure, but which might be far distant from one another on the linear chromosome, can be ligated and subsequently detected by PCR. A schematic representation of our 3C experiment is given in Fig. 2A. In addition to the enhancer/SNP and the promoter region of OBFC1 gene, we examined a genomic region physically opposite to the enhancer containing rs2995264 and lacking typical characteristics of a regulatory element as a negative control (mock region). An artificial template consisting of promoter region linked to mock region and obtained by overlap extension PCR was used as a positive control to validate technical set up (Supplementary Material, Fig. S4). Specific products were amplified in both cell lines with primers targeting the restriction fragments of enhancer/SNP and OBFC1 promoter, in samples that had been cross-linked, but not in samples that were not cross-linked (Fig. 2B). Ultimately, the results of 3C in the analyzed CMM cell lines, confirmed the interaction between the OBFC1 promoter and the regulatory element associated to the genetic variant rs2995264.
The enhancer containing rs2995264 interacts with the OBFC1 promoter in CMM cells. (A) Schematic representation of Chromosome Conformation Capture (3C) experiment displays approximate positions of analyzed regions, direction of transcription of OBFC1 and EcoRI cutting sites within the area. Primers P1 and P2 were designed to amplify a novel ligation product formed between the restriction fragments that encode the promoter region and enhancer DNA, respectively, whereas P2 and P3 amplify a ligation product for an interaction with the promoter region that is not expected. (B) The interaction between the enhancer/SNP region and the promoter region (P1-P2) and between the promoter region (P2) and a distal element (P3) in A375 and CJM cells was assessed. The interaction frequency corresponds with the intensity of amplified PCR products analyzed gels are shown in Supplementary Material, Fig. S4 and Supplementary Material, Supplementary Information. Data are shown as mean ± SD. (C) Design of CRISPR-mediated enhancer deletion in HEK293T cells showing the sites of targeted deletion in the intron 3 of OBFC1. The target regions are indicated by the dashed lines (deletion A of 720 bp and deletion B of 756 bp) flanked by the pairs of single guides RNA (sgRNA2 yellow, sgRNA3 green and sgRNA6 orange). Agarose gel image with validation PCR results of heterozygous and homozygous deletions are shown in Supplementary Material, Fig. S5. (D) mRNA and protein extracts were collected from H293T selected isogenic lines. (qRT)-PCR and western blot analysis were performed to verify the effect of homozygous enhancer deletion on OBFC1 mRNA and protein levels. Data shown are the mean ± standard deviation from three independent (qRT)-PCR experiments. (*P-value <0.01; two-tailed T-test). (E) Allele-specific binding proteins were identified by mass spectrometry using CMM cell nuclear extract and biotinylated double-stranded oligonucleotides. The dimethyl-labeling ratios of proteins bound to A protective allele (orange) or G risk allele (blu) probes are plotted on the x and y axes. (F) MEOX2 preferentially bound to A protective allele of rs2995264 both in basal conditions (−) and after MEOX2-overexpression (+), as determined by ChiP-assay. Data shown are the mean ± SD from two independent (qRT)-PCR experiments, each done in triplicate; the enrichment scores are relative to Rabbit IgG, used as control isotype.
To further illustrate the importance of this regulatory element in inducing expression of OBFC1 gene, we deleted the enhancer region containing rs2995264 by CRISPR/Cas9 system (Fig. 2C) in the HEK293T cell line that is frequently used for genome editing due to its high efficiency. The regulatory element was targeted by two single guide RNA (sgRNAs) pairs that efficiently deleted the region overlapping the enhancer (deletion A of 720 bp and deletion B of 756 bp) (Supplementary Material, Fig. S5). We confirmed that the homozygous deletion of the enhancer region decreased OBFC1 expression levels compared to wild type HEK293T cells (Fig. 2D). Collectively, these results provide strong evidence for the role of the enhancer containing rs2995264 in regulating OBFC1 gene expression.
To identify proteins that bind the SNP rs2995264 in an allele-preferential manner, we used affinity purification mass spectrometry (AP-MS) (27,32): DNA pull down using nuclear A375 (A/A genotype for rs2995264) and UACC1816 (G/G) extracts identified rs2995264-A and rs2995264-G preferential interactors, and MEOX2 resulted as the most significant interactor (Fig. 2E). In view of this, ChiP experiments were performed to determine whether the sequence containing rs2995264 actually binds MEOX2 nuclear protein in allele-preferential manner in A375 (A/A) and UACC1816 (G/G) cell lines in basal conditions and after 48 h of transfection with MEOX2 expressing plasmid (Supplementary Material, Fig. S6). Figure 2F indicates a MEOX2 binding enrichment in the presence of rs2995264-A allele. Accordingly, MEOX2 transcription factor recognizes and binds specifically to the protective rs2995264-A allele, with greater affinity, compared to the rs2995264-G risk allele. Altogether, these data suggest that the rs2995264 polymorphism alters the binding of MEOX2 transcription factor, possibly leading to alteration of the OBFC1 transcriptional machinery.
OBFC1 has tumor-suppressor effect in CMM
To unravel the potential OBFC1 contribution to CMM development, we tested the gene expression in three independent mRNA expression array datasets (Supplementary Material, Supplementary Information). Expression profile of OBFC1 was significantly lower in melanoma when compared to histologically benign tumors (nevi) (Fig. 3A) and in metastatic melanoma when compared to primary melanomas (Fig. 3B and C). We found no significant correlation between OBFC1 expression and patient survival (Supplementary Material, Fig. S7). Together, these data provide evidence that OBFC1 might play a biological role in CMM initiation rather than progression. So, we planned to test whether decreased OBFC1 levels could lead to cellular transformation required in tumor onset. We selected A375 and CJM (A/A rs2995264 genotype) cell lines with high OBFC1 expression and UACC1816 (G/G rs2995264 genotype) with low OBFC1 expression, validated by (qRT)-PCR and western blot analysis (Supplementary Material, Fig. S8A and B). We thus examined the consequence of OBFC1 knocked down by using short interfering RNA (siRNA) against OBFC1. Compared with scrambled siRNA used as a control, three-pooled siRNA (siRNA_A, siRNA_B, siRNA_C) specific for OBFC1, significantly reduced OBFC1 mRNA and protein levels at 48 h post-transfection (Fig. 3D–F). Knockdown of OBFC1 markedly increased cell viability in A375 and CJM cell lines, carrying the AA-rs2995264 protective genotype and high OBFC1 expression, compared to the control cells (siScrambled) (A375: T24 P ≤ 1×10−3; T48 P ≤ 2.1×10−6; T72 P ≤ 1.2×10−7; CJM T24 P ≤ 6×10−3; T48 P ≤ 7×10−3; T72 P ≤ 5.5×10−7), whereas in UACC1816, with GG-rs2995264 risk genotype and low basal OBFC1 expression, we observed non-substantial differences in cell viability after OBFC1 transient silencing (Fig. 3G–I). Importantly, these findings indicate that OBFC1 has potential tumor-suppressor effect in CMM and its reduced expression due to disease-predisposing alleles may contribute to CMM progression by promoting tumor cell survival.
Low OBFC1 expression is associated with unfavorable histology in CMM and transient knockdown of OBFC1 influences CMM cell viability in a genotype-specific manner. (A–C) The dark line in the middle of the boxplots shows the median value of mRNA expression of OBFC1 in GSE3189 (P-value = 0.025), GSE112509 (P-value = 0.0007) and TCGA datasets (P-value = 0.08). The bottom of the box indicates the 25th percentile. The top of the box represents the 75th percentile. The T-bars that extend from the boxes indicate the minimum and maximum values. P-values obtained by two-tailed T-test. (D–F) OBFC1 siRNA knockdown as measured by (qRT)-PCR and western blot 48 h post transfections for experiments. The datapoints represent the mean of technical duplicates of each of the three independent experiments. *P-value ≤0.01; P-value obtained by two-tailed T-test. (G–H) In cells homozygous for rs2995264 protective A allele and with higher OBFC1 expression levels (Supplementary Material, Fig. S8), OBFC1 transient knockdown leads to significant increase in cell viability. (I) In cells homozygous for rs2995264 risk G allele and with lower OBFC1 expression, OBFC1 transient knockdown does not affect cell viability. (G–I) Cell viability was determined using MTT assay at 0, 24, 48 and 72 h post siRNA transfection (X-axis). The amount of MTT formazan is directly proportional to the number of living cells (Y-axis). Data shown are the mean ± standard deviation from two independent MTT experiments, each done in six-duplicate; P-value obtained by two-tailed T-test.
OBFC1 functions in telomeres maintenance by regulating telomerase activity
The 10q24.33 locus has been associated with telomere length and CMM (8,33). As OBFC1 belongs to CST complex which turns off telomerase activity by inhibiting its binding to telomeric DNA (34), we hypothesized that decreased OBFC1 expression could predispose to CMM risk allowing telomerase activity. To test our hypothesis, we measured telomerase activity after OBFC1 transient silencing in melanoma cell lines (Supplementary Material, Fig. S9A and B) with a PCR-based assay that permitted quantitation of telomerase enzymatic activity. Compared to control cells (siRNA Scrambled), OBFC1 knockdown markedly increased telomerase activity in melanoma cell lines, suggesting a potential loss of CST complex capability to suppress telomerase access to lengthening telomeres (Fig. 4A). Literature data showed that endogenous telomerase action at telomeres is restricted to the cell cycle S phase (35–37). Consistent with these findings, we expected that increased telomerase activity after OBFC1 depletion could coincide with a higher percentage of cells in S phase. Here, we examined the proportions of A375 and CJM, cells with high OBFC1 expression, at each stage of the cell cycle by flow cytometry after siRNA OBFC1 treatment (Supplementary Material, Fig. S9C–F) and we found a direct proportionality between levels of telomerase activity and the percentage of cells in S phase (Fig. 4B and C). It has been known that the cyclin gene with the highest transcription rate during S phase is CCNA2 (38). So, we decided to reinforce the previous evidence of S phase cells accumulation with evaluation of CCNA2 mRNA levels after OBFC1 depletion. The mRNA of CCNA2 was higher in OBFC1 silenced CMM cells, as expected (Fig. 4D).
OBFC1 transient knockdown enhances telomerase activity. (A) OBFC1 transient knockdown in CMM cells increases significantly telomerase activity. Data shown represent three independent experiments and the datapoints represent the mean among technical duplicates of each experiment; *P-value <0.03, obtained by two-tailed T-test. (B–C) Percentage of cell cycle distribution of CMM cells after transient OBFC1 knockdown. Data shown represent three independent experiments and the datapoints represent the mean among technical duplicates of each experiment; P-value obtained by two-tailed T-test (**P-value <0.05). (D) Evaluation of CyclinA2 mRNA expression levels, marker of S phase, in CMM cells analyzed by flow cytometry for cell cycle distribution. Data represent three independent experiments and the datapoints represent the mean among technical duplicates of each experiment; P value obtained by two-tailed T-test (***P-value <0.01).
Discussion
The OBFC1 locus (10q24.33) has been previously identified as a CMM and neuroblastoma susceptibility locus (8,13,39). However, at this locus, most of the functional variant(s) responsible for biological mechanisms accounting for the risk and genes involved in CMM pathogenesis have yet to be elucidated. To functional characterize the CMM risk variant(s) at 10q24.33 and to determine the genes affected by the same variants, we have carried out an integrative genomic analysis of 10q24.33 locus that reached significance in our neuroblastoma-CMM GWAS meta-analysis performed previously (13).
We developed a specific strategy based on epigenomic annotations of a large number of CMM and NCC cell lines to identify regulatory variants at 10q24.33 locus, which could affect transcriptional machinery. We found three putative functional SNPs (rs2995264, rs34685262 and rs4387287) enriched in active enhancers of NCC and CMM cells, but only rs2995264 induced enhancer activity. We also confirmed that the minor allele G of rs2995264 is associated with CMM risk in an independent cohort of Italian origin. HiC and 3C experiments confirmed the interaction between the SNP rs2995264 and the promoter region of OBFC1. Moreover, an isogenic model characterized by CRISPR-Cas9 deletion of the enhancer region containing the SNP rs2995264 confirmed the enhancer regulatory potential on OBFC1 transcription.
In line with the established rs2995264 functional properties, we have demonstrated that G-rs2995264 risk allele decreased the binding affinity of MEOX2, a homeobox transcriptional factor that seems to mediate carcinogenesis by altering the normal mechanisms of angiogenesis and cell proliferation (40–42). These results provide evidence that a functional DNA variant in the enhancer region of OBFC1 influences CMM susceptibility.
Our data are strengthened by NHGRI-EBI GWAS catalog which reports the association of four SNPs (rs2487999, rs4387287, rs9420907 and rs9419958) located in OBFC1 genomic region with telomere length (43), suggesting an involvement of OBFC1 gene in tumorigenic transformation.
Furthermore, we demonstrated that the G risk allele correlated with low OBFC1 expression in CMM tumors, and the OBFC1 expression is reduced in primary melanoma and metastatic tumors, thus indicating a potential tumor-suppressor effect of OBFC1 in melanoma. Our findings are in accord with the work of Phelan et al. showing that the minor allele of an epithelial ovarian cancer predisposing SNP, in complete LD with rs2995264, correlated with low OBFC1 expression in ovarian cancer tissues (44).
OBFC1 (OB Fold-containing Protein 1), a human homolog of yeast STN1, is a subunit of an alpha accessory factor that stimulates the activity of DNA polymerase α primase, the enzyme that initiates DNA replication (45). OBFC1 is also known to be a key component of telomere-associated CST complex that binds telomeric single-stranded DNA in vitro and localizes at telomeres in vivo (46). The contribution of OBFC1 to cancer susceptibility firstly emerged from the large genetic association analyses of patients with different histotypes of epithelial ovarian cancer that identified 10q24.33 as a risk locus associated with borderline serous epithelial ovarian cancer (44). Our work confirms the involvement of OBFC1 gene in cancer predisposition demonstrating that transient knockdown of OBFC1 resulted in significant increase of cells viability in CMM cells with AA-rs2995264 protective genotype and basal OBFC1 expression, but had little effect on cells with GG-rs2995264 risk genotype and lower basal OBFC1 expression. These findings suggest that that GG genotype leads to a decreased OBFC1 expression level that contribute to promote cancer initiation.
It is known that in physiological conditions the human CST (CTC1, OBFC1 and TEN1) complex inhibits telomerase activity (47), as demonstrated by the excessive telomerase activity resulted after CST depletion (34). Additional experimental evidence showed that ectopic overexpression of an OBFC1 truncation mutation also led to telomere length increase over time (46). Here we present preliminary data supporting a role of OBFC1 in telomerase homeostasis, with an observed enhancement of telomerase activity in CMM cells after OBFC1 depletion, and provide evidence for a putative molecular mechanism that confers genetic susceptibility to melanoma. One potential limitation of this study is the use of A375 cell line that may present different transcriptomic and epigenomic profiles to other cell lines derived from cutaneous melanoma. Therefore, in vitro validation experiments of our data using additional CMM primary cells more closely similar to CMM tumors are needed to be performed in future studies. Moreover, additional research efforts are required in order to elucidate the molecular mechanisms of OBFC1 in promoting telomere maintenance required to CMM malignant transformation.
In conclusion, our results support the hypothesis that decreased expression of OBFC1 gene through functional heritable DNA variation can contribute to malignant transformation of normal cutaneous cells. Moreover, we provide preliminary data suggesting the potential effect of OBFC1 on the telomerase activity in tumorigenic conditions. This study has demonstrated that post-GWAS strategies are a useful step for the identification of causal functional variants at previously identified cancer risk loci and for the elucidation of the key roles of genes involved in tumor biology.
Materials and Methods
Identification of causal variant at 10q24.33
We first selected the variants in LD with the rs11591710 lead SNP (0.5 < r2 ≤ 1) (total including the lead SNP, n = 33) in European population using LDlink (analysistools.cancer.gov/LDlink). Then, to identify potential functional variants, we annotated these SNPs with multiple sources of in silico functional annotation from public databases, as detailed in Supplementary Material, Supplementary Information. We obtained the genome binding/occupancy profiling by high-throughput sequencing of the epigenetic marker H3K27ac in seven human derived-melanoma (GSE75352) and two human neural crest cell (hNCC) lines (GSE90683) through the National Center for Biotechnology Information.
CMM replication in an Italian cohort
The genomic DNA of CMM patients and healthy controls was extracted from peripheral blood using a Maxwell® RSC Blood DNA Kit (Promega, Madison, WI, USA), and DNA concentration and purity were evaluated using a NanoDrop™ 8000 Spectrophotometer. The rs2488001 SNP, in near perfect LD (r2 = 0.98) with rs2995264 SNP was typed by TaqMan® SNP Genotyping Assay (Applied Biosystems by Thermo Fisher Scientific, Waltham, MA, USA) in an Italian cohort of 484 CMM cases and 1801 controls. To monitor quality control, three DNA samples per genotype were genotyped by Sanger sequencing (3730 DNA analyzer, Applied Biosystems) and included in each 384-well reaction plate; genotype concordance was 100%. This study was approved by the Ethics Committee of the Medical University of Naples (N. 76/13).
Statistical analysis
A comparison of the genotypic and allelic frequencies between the groups was performed using the chi square test. Statistical significance was established at P < 0.05. Hardy–Weinberg equilibrium was evaluated using the goodness-of-fit chi squared test in control and case subjects (P > 0.05). The COnditional and JOint (COJO) function from the GCTA software package (14) was used to run conditional analyses at the rs2995264 locus, and identify independent association signals based on melanoma GWAS meta-analysis summary statistics including 12 874 cases and 23 203 controls (8). We used COJO to perform secondary association scans conditioning on rs2995264 with 1000 Genomes EUR reference panel for LD estimation. The association plot was generated using LocusZoom (15). All LD calculations (r2 and D′) were performed using the LDlink suite (https://ldlink.nci.nih.gov/?tab=home) and data from the 1000 Genomes Project European ancestry populations.
eQTL analysis for the SNP rs2995264
We tested the equality of variance by Levene’s Test between the two groups with opposite genotypes. Since the results of test was not statistically significant, we assumed that we had homogenous variance between groups. The independent T-test was thus performed to test statistical differences between the means of two groups.
Hi-C data analysis
As detailed in Supplementary Material, Supplementary Information, the sequencing was performed on an Illumina® HiSeq platform. Paired-end reads with length of 150 bp were mapped to the reference genome (build hg19/GRCH37) with Bowtie2 (16). The alignment BAM file was then filtered to remove duplicates, re-ligation or self-circularization artifacts that can be introduced during Hi-C library preparation. Then we used HiCExplorer tool v3.5.1 (16) to (i) build the interaction matrix at a resolution of 10 Kb (bin size = 10Kb); (ii) normalize the observed interaction matrix; (iii) determine Topologically Associating Domains (TADs, self-interacting genome regions) and their boundaries; and (iv) plot the results. Subsequently, we extended our region of interest containing LD SNPSs of 1 Mb up- and down-stream and calculated the statistical significance of the interactions between bins with the FitHiC v2.0.7 program (17). P-values were corrected for multiple tests by Benjamini-Hotchberg method (false discovery rate, FDR) and the cutoff was set at 1%. Finally, we annotated those bins with ANNOVAR (18) in order to map genomic bins to gene coordinates.
In vitro functional study
A detailed description of the Luciferase reporter assays and the experiments performed to evaluate the OBFC1 effect on CMM cell lines phenotype is reported in Supplementary Material, Supplementary Information.
Cell culture
The human A375 cell line and the human CJM cell line are CMM cell lines derived from skin epithelia and were donated from Professor Nick Hayward (QIMR Berghofer Medical Research Institute, Australia). The human UACC1816, a CMM cell line isolated from skin, was donated from Dr Kevin Brown (Translational Genomics Research Institute TGen, AZ, USA). The human HEK293T epithelial cell line was obtained from the American Type Culture Collection (CRL-3216). HEK293T and A375 cell lines were grown in Dulbecco’s Modified Eagle Medium (DMEM; Sigma); CJM cell line was grown in RPMI-1640 Media (Sigma); UACC1816 cell line was grown in RPMI-1640 Media (Sigma) supplemented with 25 mM HEPES (Sigma). The mediums were supplemented with 10% heat-inactivated FBS (Sigma), 1 mmol/L L-glutamine, penicillin (100 U/ml), and streptomycin (100 mg/ml; Invitrogen). The cells were cultured at 37°C, 5% CO2 in a humidified atmosphere. The cell lines used for all the experiments were re-authenticated and tested as mycoplasma-free. Early-passage cells were used and cumulative culture length was <3 months after resuscitation. Total cellular RNA extraction, reverse transcription and Quantitative Real Time (qRT)-PCR were performed as previously described (19).
Chromosome conformation capture (3C)
3C procedure is described more in details in Supplementary Material, Supplementary Information. CMM cells were used, A375 and CJM respectively. In order to obtain a negative control for the 3C analysis, in addition to the non–cross-linked sample, we selected a fragment chr10:105693245–105 701 218 that did not show any characteristic that can be associated with a regulatory region and for which is not expected an interaction with the promoter region. To guarantee the correct setting of the PCR experimental conditions it has been also necessary to produce, through overlapping PCR, a PCR product consisting of the restriction fragments corresponding to those of the intronic region and promoter. PCR products were resolved on 2% agarose gels. In order to normalize 3C-PCR signals, we used a loading control (internal primers located in the GAPDH gene (20)). The amount of DNA input was first titrated, and bands analyzed semiquantitatively using ImageJ software; the background was subtracted, and data normalized to an internal region unaffected by the restriction digest (LC region) (21,22). Two biological replicates were prepared and analyzed in three technical repeats.
CRISPR-based enhancer deletion
We used the CRISPR/Cas9 system to generate HEK293T isogenic cell line with a deletion of the enhancer region (hg19/chr10:105668100-105 669 000), as confirmed by the peak in H3K27ac Chip-Seq of HEK293T cell line (Encode project track ENCR000FCH). To delete this enhancer element, following Bauer et al. and Ran, FA et al. guidelines (23,24), two pSpCas9 (BB)-2A-GFP (PX458) vectors expressing Cas9 (Addgene) and the desired sgRNAs, designed using CRISPOR (25) (http://crispor.tefor.net), CHOPCHOP (26) (https://chopchop.cbu.uib.no) and CRISPR (24) (http://tools.genome-engineering.org) tools, were cotransfected into HEK293T cell line using Transfectin Lipid Reagent (Bio-Rad). sgRNA2 (5′- ACAGGCCTGCGGTGAGTCAG-3′) and sgRNA3 (5′- CGGGATGAGTCAGTGCGAGC-3′) were paired with sgRNA6 (5′- CAGCTATGGGCAGTACACTG-3′) to make clones with different genomic deletions (720 bp and 756 bp, respectively). After 48 h of transfection, GFP-positive single cells were FACS-sorted by size into 96-well plates. To identify and distinguish both mono-allelic and bi-allelic deletions, a PCR using two primer pairs flanking the sgRNA cleavage sites (Primer-OUT-Forward: 5′-TGCGAGGTCATTCTGGTCTTG-3′; Primer-OUT-Reverse: 5′-AACTTTGTGACCAAGAGCGT-3′) was performed with KAPA HiFi HotStart PCR Kit (Roche), following manufacturer’s instructions. Other two primers falling into the deleted sequence (Primer-IN-Forward: 5′- GTGAGTCAGGGGAAGCAGAA-3′; Primer-OUT-Reverse: 5′- TCCAGCTATGGGCAGTACAC-3′) were used to confirm deletion occurred (Supplementary Material, Fig. S5). These two sets of primer were used to screen single cell-deriving clones to evaluate which ones was edited correctly.
Affinity purification mass spectrometry
Quantitative AP-MS following SNP DNA pulldown was performed on the basis of procedures described by Choi et al. (27–29). For DNA pulldown, 500 pmol of annealed, forward-strand 5′-biotinylated oligonucleotide probe was coupled to Streptavidin Sepharose beads (GE Healthcare). rs2995264-A and 2 995 264-G allele probe sequences are: 5′-TGTACTTTCTGTTTCAAAAGA-3′; 5′-TGTACTTTCTATTTCAAAAGA-3′.
Chromatin immunoprecipitation (ChiP)
A375 and UACC1816, both in basal condition and after transfection with MEOX2 origene expressing plasmid pCMV6-ENTRY MEOX2 were used, and a detailed description of the procedure is shown in Supplementary Material, Supplementary Information. In order to validate MEOX2-ChIP reaction, MEOX2 binding sites in individual genes were identified using JASPAR Web Tools (http://jaspar.genereg.net/). The promoter of p21 gene known for MEOX2 binding (30) was used as positive control and analyzed using (qRT)-PCR. NFT3 genomic region (chr12:5542580–5 542 721) with the highest H3K27me3 peak and lowest H3K27ac peak in skin keratinocytes Chip-seq experiments, available in Encode (ENCSR621FNM; ENCSR736ZEG; ENCSR709ABP; ENCSR793NQA), was used as negative control (Neg Ctrl).
Acknowledgements
The authors would like to thank Professor Nick Hayward (QIMR Berghofer Medical Research Institute, Australia) for the donation of A375 and CJM cell lines and Dr Kevin Brown (Translational Genomics Research Institute TGen, Arizona USA) for the donation of UACC1816 cell line.
Conflict of Interest statement. None declared.
Funding
Associazione Italiana per la Ricerca sul Cancro (Grant no. 19255 to M.C. and Grant no 20757 to A.I.); Fondazione Italiana per la Lotta al Neuroblastoma (to M.C.); Associazione Oncologia Pediatrica e Neuroblastoma (to M.C.); Ministero della Salute (Grant no. RF-2016-02362288 to P.G.); Regione Campania ‘SATIN’ Grant 2018-2020 (to M.C.). The M. Vermeulen lab is part of the Oncode Institute, which is partly funded by the Dutch Cancer Society (K.W.F.).
