Common variation at 6q16 within HACE1 and LIN28B influences susceptibility to neuroblastoma.

Neuroblastoma is a cancer of the sympathetic nervous system that accounts for approximately 10% of all pediatric oncology deaths. Here, we report a genome-wide association study of 2,817 neuroblastoma cases and 7,473 controls. We identified two new associations at 6q16, the first within HACE1 (rs4336470; combined P=2.7×10(-11); odds ratio 1.26, 95% confidence interval (CI) 1.18-1.35) and the second within LIN28B (rs17065417; combined P=1.2×10(-8); odds ratio 1.38, 95% CI 1.23-1.54). Expression of LIN28B and let-7 miRNA correlated with rs17065417 genotype in neuroblastoma cell lines, and we observed significant growth inhibition upon depletion of LIN28B, specifically in neuroblastoma cells that were homozygous for the risk allele. Low HACE1 and high LIN28B expression in diagnostic primary neuroblastomas were associated with worse overall survival (P=0.008 and 0.014, respectively). Taken together, these data show that common variants in HACE1 and LIN28B influence neuroblastoma susceptibility and indicate that both genes likely have a role in disease progression.

Neuroblastoma is a malignancy derived from the developing sympathetic nervous system. The median age at diagnosis is 17 months, and the survival rate for the most aggressive subset remains approximately 50% despite intensive multi-modal cytotoxic therapy[1] and recent advances in immunotherapy[2]. The genetic etiology of familial neuroblastoma, which accounts for approximately 1% of cases, has recently come into focus[3-6]; however, the genetic and environmental factors that cause sporadic neuroblastoma remain largely unknown. We have recently reported common SNPs within or upstream of LINC00340 and FLJ44180, BARD1, LMO1, DUSP12, HSD17B12 and DDX4-IL31RA and a common CNV within NBPF23 as each being highly associated with neuroblastoma[7-11]. Collectively, however, these variants still account for only a small portion of neuroblastoma heritability, and it is likely that additional predisposition loci remain to be discovered.

To identify additional variants associated with neuroblastoma, we expanded our previous genome-wide association study (GWAS) discovery cohort and analyzed 2,101 neuroblastoma cases accrued through the North American-based Children’s Oncology Group (Supplementary Table 1) with 4,202 control subjects of European ancestry matched genetically and by genotyping array version to minimize genomic inflation (see Online Methods and Supplementary Figs. 1 and 2). All subjects were genotyped using the Illumina HumanHap550 or Quad610 Beadchip. We restricted our analysis to the SNPs present on both platforms that passed our quality control metrics; the genomic control inflation factor was 1.14 (Supplementary Fig. 3). Evaluation of the first three principal components in cases and controls confirmed that the slightly high inflation factor was not due to gross population stratification (Supplementary Fig. 4). Clusters of SNPs from six genomic loci reached genome-wide significance (P-values ranged from 7.8 × 10−16 to 4.8 × 10−8; Supplementary Fig. 1 and Supplementary Table 2), including three SNPs within LINC00340 and FLJ44180 at 6p22 (P-values ranged from 7.8 × 10−16 to 1.7 × 0−14), ten SNPs within or near BARD1 at 2q35 (P-values ranged from 4.1 × 10−14 to 3.7 × 10−8), two SNPs within LMO1 at 11p15 (P-values ranged from 1.2 × 10−13 to 3.8 × 10−10) and one SNP within HSD17B12 at 11p11 (P = 4.8 × 10−8), further confirming our previous reports7-10. In addition, we identified one SNP (rs4696715) at chromosome 4p16 and two SNPs (rs4336470 and rs9404576) at 6q16 that have not been reported previously (P-values ranged from 1.8 × 10−8 to 3.4 × 10−8; Table 1). Several genotyped SNPs in strong linkage disequilibrium (LD) with rs4696715 at 4p16 did not show evidence for association with neuroblastoma; therefore, this SNP was not considered further here.

Table 1

Significantly associated genotyped SNPs at the HACE1 and LIN28B loci at 6q16

		Discovery cohortaEuropean ancestry			Italian replicationa(TaqMan)			CHOP replicationaAfrican American			Combined
SNP	A1/A2(protective/risk allele)	Freq A1cases(n = 2,101 )	Freq A1controls(n = 4,202)	P b	Freq A1cases(n = 351 )	Freq A1controls(n = 780)	P b	Freq A1cases(n = 365)	Freq A1controls(n = 2,491 )	P c	MetaPd	ORe(95% CIf)
HACE1
rs4336470	T/C	0.30	0.35	1.8 × 10−8	0.30	0.34	0.060	0.62	0.70	1.4 × 10−3	2.7 × 10−11	1.26(1.18–1.35)
rs9404576	G/T	0.30	0.35	3.4 × 10−8	-	-	-	0.62	0.70	1.3 × 10−3	1.8 × 10−10	1.27(1.18-1.36)
rs4079063	G/A	0.43	0.47	4.0 × 10−5	-	-	-	0.70	0.78	1.3 × 10−4	1.3 × 10−7	1.20(1.12-1.29)
rs2499663	C/T	0.43	0.47	4.5 × 10−5	-	-	-	0.70	0.78	1.5 × 10−4	1.6 × 10−7	1.21(1.13-1.29)
rs2499667	G/A	0.43	0.47	2.6 × 10−5	-	-	-	0.71	0.78	2.6 × 10−4	1.2 × 10−7	1.21(1.13-1.29)
LIN28B
rs17065417	C/A	0.08	0.11	1.8 × 10−7	0.08	0.11	0.033	0.09	0.11	0.129	1.2 × 10−8	1.38(1.23-1.54)

No deviations from Hardy-Weinberg equilibrium were observed (P > 0.001) in all cohorts.

P-values were calculated by allelic test.

P-values were calculated by logistic regression with percent African admixture as covariate[12].

Meta-analysis P-value calculated using METAL[14].

OR, odds ratio of risk allele based on meta-analysis.

CI, confidence interval.

Closer examination of the 6q16 locus identified four additional SNPs that showed association (P < 1 × 10−4) with neuroblastoma (Table 1 and Supplementary Table 3). Three of these SNPs mapped to introns of the HECT domain and ankyrin repeat containing, E3 ubiquitin protein ligase 1 gene (HACE1) and exhibited a moderate degree of LD with rs4336470. The fourth SNP (rs17065417; P = 1.8 × 10−7) mapped to an intron of the lin-28 homolog B gene (LIN28B) and showed very little evidence for LD with rs4336470 in multiple HapMap populations (Supplementary Fig. 5). To ensure these results were not influenced by subtle substructure in the discovery phase, we included the first 20 principal components as covariates in a logistic regression analysis. This reduced the inflation factor to 1.04 but did not alter our conclusions regarding regions reaching genome-wide significance, and the P-values for SNPs at the newly identified 6q16 locus were essentially unchanged (Supplementary Table 4). Finally, we conditioned the analysis of the 6q16 SNPs on rs4336470 to investigate if there might be more than one independent association signal at 6q16. As expected, SNPs in modest LD with rs4336470 were no longer statistically significant after conditioning and clearly represent one signal (Supplementary Table 5). In contrast, while the signal at rs17065417 was attenuated, it remained significant (P = 2.5 × 10−4; Supplementary Table 5), suggesting that rs4336470 and rs17065417 may contribute independently to neuroblastoma risk. No association was observed between rs4336470 or rs17065417 genotypes and clinical or biological covariates (Supplementary Tables 6 and 7), and we observed only weak evidence for epistasis at the other significant loci (Supplementary Tables 8 and 9), suggesting that each may act independently to confer risk.

We next sought to replicate the rs4336470 and rs17065417 associations in an Italian cohort of 351 cases and 780 controls using PCR-based genotyping. Both SNPs showed evidence for association in the same direction seen in the discovery effort (Table 1 and Supplementary Table 3). To assess whether these variants influence susceptibility in other ethnic groups, and to seek additional replication, we analyzed a third independent case series comprised of 365 African American neuroblastoma cases and 2,491 genetically matched controls, all genotyped on the Illumina HumanHap550 or Quad610 BeadChip. After accounting for African admixture in logistic regression analysis, the genomic inflation factor was 1.02, as described previously[12] (Online Methods and Supplementary Fig. 6). Consistent with the lower incidence rate of neuroblastoma in African Americans[13], the rs4336470 protective allele (T) is actually the major allele in African Americans, whereas it is the minor allele in individuals of European ancestry. Despite these observed differences in allele frequencies, all SNPs mapped to HACE1 replicated robustly in the African American cohort (P-values ranged from 1.3 × 10−4 to 1.4 × 10−3; Table 1 and Supplementary Table 3), and all associations were in the same direction. Allele frequencies for the LIN28B SNP rs17065417 were comparable across ethnic populations and showed a trend toward association in African Americans in the same direction as the European American discovery and Italian replication cohorts (Table 1 and Supplementary Table 3). Combined analysis using the weighted inverse-variance method in METAL[14] demonstrated that all six genotyped SNPs had P values beyond or approaching the conservative Bonferroni adjusted genome-wide significance threshold (P-values ranged from 2.7 × 10−11 to 1.6 × 10−7; Table 1).

To identify variants at the 6q16 locus not assayed directly on the Illumina SNP arrays, we performed genotype imputation in our discovery cohort using data from the 1000 Genomes project. This analysis identified ten additional genome-wide significant SNPs (P-values ranged from 4.4 × 10−9 to 5.7 × 10−8; Fig. 1). Five of these imputed SNPs were in strong LD with rs4336470 (r2 > 0.8 in 1000 Genomes EUR population); three mapped to introns of HACE1 and the other two were located just downstream of HACE1 (Fig. 1a). The remaining five imputed SNPs reaching genome-wide significance were in strong LD with rs17065417 (r2 > 0.8 in 1000 Genomes EUR population) and were located within introns of LIN28B (Fig. 1b). In total, we identified 46 imputed SNPs that showed strong evidence for association (P < 1.0 × 10−6, Supplementary Table 10). To further evaluate whether two independent association signals exist at 6q16, we conditioned the regional association analysis on either rs4336470 or rs17065417 and found that neither could fully account for the observed associations (Supplementary Fig. 7). The signal at 6q16 was only abolished across the entire region after conditioning on both rs4336470 and rs17065417 (Supplementary Fig. 8). While we cannot rule out the possibility that both SNPs may be tagging the same underlying risk variant, these data are consistent with the presence of two independent association signals at 6q16, one implicating HACE1 and the other LIN28B.

Figure 1

Regional association plots at the HACE1 and LIN28B loci. Regional association plots including both genotyped and imputed SNPs for the HACE1 and LIN28B loci generated by LocusZoom[36]. Plotted are the significance of association (-log10 transformed P values) and the recombination rate. SNPs are color-coded based on pair-wise linkage disequilibrium (r2) with the most significant genotyped SNP in the EUR (European) 1000 Genomes Interim Phase I release genotypes. The most significant genotyped SNP and associated P-value are labeled, and the SNP is shown in purple. (a) HACE1 locus (b) LIN28B locus.

The HACE1 gene encodes an E3 ubiquitin-protein ligase that was first identified in a sporadic Wilms tumor that harbored a t(6:15)(q21;q21) translocation[15]. The rearrangement was associated with decreased HACE1 expression, and further study showed that HACE1 is silenced in the majority of Wilms tumors via hypermethylation of two CpG islands upstream of the transcriptional start site[15]. Similar epigenetic silencing has been reported in advanced colorectal cancer[16] and gastric carcinoma[17]. Indeed HACE1 is down-regulated in multiple human tumors and maps to a region of common deletion or LOH, consistent with a tumor suppressor function. Hace1 null mice form spontaneous tumors in a wide array of tissues and are susceptible to additional cancer triggers, both genetic and environmental[18]. Through its E3 ubiquitin ligase function, HACE1 has been shown to suppress cell growth and anchorage independence of human tumor cells, including the neuroblastoma cell line IMR32 (ref. 18). Studies suggest that HACE1 inhibits cell cycle progression during stress via regulation of cyclin D1 degradation[18] and that HACE1 also regulates retinoic acid receptor (RAR) activity[19]. In addition to multiple somatic alterations involving the HACE1 locus, a constitutional t(5;6)(q21;q21) translocation that disrupts HACE1 was recently identified, making it a putative Wilms tumor susceptibility gene[20].

The LIN28B gene encodes a developmentally regulated RNA binding protein and is a key repressor of the let-7 family of microRNAs[21]. Both LIN28A and LIN28B function as oncogenes that promote cellular transformation when ectopically over-expressed[22-24], and high levels of expression of LIN28A or LIN28B have been observed in several human cancers and correlate with low let-7 levels[22,24]. Post-transcriptional regulation of let-7 by LIN28A is required for normal development and contributes to the pluripotent state by preventing let-7–mediated differentiation of embryonic stem cells[21,25,26]. Over-expressing LIN28 or inhibiting let-7 with antisense RNA promotes reprogramming of human and mouse fibroblast to pluripotent stem cells[26,27]. In a panel of 60 pediatric cancer cell lines[28], LIN28B was consistently expressed at high levels in neuroblastoma (Supplementary Fig. 9). Recently, LIN28B and let-7 were also identified as key regulators of glucose homeostatis[29] and hematopoiesis[30]. To date, GWAS studies have identified variants in LIN28B associated with human height[31] and the age of onset of puberty[32] and menarche[33]. In addition, a recent candidate gene study identified a putative association with epithelial ovarian cancer[34].

To investigate the functional relevance of neuroblastoma-associated SNPs within HACE1 and LIN28B, we first analyzed a set of 12 neuroblastoma cell lines with matched genome-wide SNP genotyping and mRNA expression data. No correlation between rs4336470 and HACE1 mRNA expression was observed (P = 0.30), but LIN28B expression was significantly higher in cell lines homozygous for the rs17065417 risk allele (A) compared to heterozygous cell lines (P = 0.02; Fig. 2a). No cell lines tested were homozygous for the rs17065417 protective allele (C). LIN28A was not expressed in neuroblastoma cell lines, consistent with recent reports that LIN28A and LIN28B may be mutually exclusive in terms of their expression in cancer cells[21]. In addition to a correlation between rs17065417 genotype and LIN28B expression, we observed a strong positive correlation with MYCN expression (P = 0.0003; Fig. 2b). We next confirmed LIN28B protein levels by western blot in four cell lines at the extremes of LIN28B mRNA expression (Fig. 2c). Cell lines expressing high LIN28B showed lower let-7 expression across the entire miRNA family (Fig. 2d). Transient knockdown of LIN28B resulted in significant growth inhibition in neuroblastoma cells homozygous for the rs17065417 risk allele and with high LIN28B expression (Fig. 3a-d), but not in the heterozygous cell line SKNAS with low LIN28B expression (Fig. 3e,f). Taken together, these data are consistent with the hypothesis that LIN28B promotes neuroblastoma tumorigenesis in part through repression of let-7 miRNAs, and that the risk alleles are associated with growth advantage through increased LIN28B expression; however, additional studies with larger sample sizes are needed to confirm this.

Figure 2

LIN28B risk alleles correlate with increased LIN28B expression and decreased let-7 miRNA expression. (a) LIN28B mRNA expression is significantly higher in neuroblastoma cell lines homozygous for the rs17065417 risk allele (AA) compared to neuroblastoma cell lines heterozygous for the risk allele (CA). Given the MAF of rs17065417, we did not identify any cell lines homozygous for the protective allele (CC). (b) Neuroblastoma cell lines homozygous for the rs17065417 risk allele show high expression of MYCN. (c) Western blot confirms increased protein expression of LIN28B in neuroblastoma cell lines homozygous for the rs17065417 risk allele. (d) Neuroblastoma cell lines homozygous for the rs17065417 risk allele and with high LIN28B expression show decreased or absent expression of the let-7 miRNAs.

Figure 3

Transient knockdown of LIN28B influences neuroblastoma cell growth in an expression-specific manner. (a-d) In cells homozygous for neuroblastoma risk alleles and with higher LIN28B expression levels, LIN28B knockdown leads to significant growth inhibition. (e) In cells heterozygous for the risk allele (carrying one protective allele) and with low LIN28B expression, LIN28B knockdown does not affect cell growth. (f) LIN28B knockdown as measured by quantitative reverse transcription PCR and western blot for experiments a–e.

To examine the relevance of HACE1 and LIN28B in tumor samples, we assayed expression of both genes by real-time quantitative RT-PCR in a representative set of 87 primary tumors obtained at time of diagnosis (Supplementary Table 1). HACE1 expression was significantly lower (P = 0.002; Fig. 4a) and LIN28B expression was significantly higher (P = 0.032; Fig. 4b) in the high-risk tumors. Accordingly, low HACE1 expression was associated with a worse overall survival (P = 0.008; Fig. 4c), as was high LIN28B expression (P = 0.015; Fig. 4d). Analysis of four independent mRNA expression array datasets comprising 517 neuroblastoma tumors provided robust replication of these observations (Fig. 4e-h and Supplementary Figs. 10-12). These data support the hypothesis that HACE1 may function as a tumor suppressor, and LIN28B as an oncogene, in advanced neuroblastomas.

Figure 4

HACE1 and LIN28B expression are associated with advanced disease and poor outcome in neuroblastoma. (a) HACE1 mRNA expression is significantly lower in high-risk neuroblastoma tumors. Bar graph is shown for Children’s Oncology Group (COG) risk groups. (error bars: s.e.m.). (b) LIN28B mRNA expression is significantly increased in high-risk neuroblastoma tumors. Bar graph is shown for Children’s Oncology Group (COG) risk groups. (error bars: s.e.m.). (c) Decreased HACE1 expression in primary tumors obtained at diagnosis is associated with worse overall survival. Kaplan-Meier analysis is shown, with patients grouped by tertiles of HACE1 expression. Log rank P-value is shown. (d) Increased LIN28B expression in primary tumors obtained at diagnosis is associated with worse overall survival. Kaplan-Meier analysis is shown, with patients grouped by tertiles of LIN28B expression. Log rank P-value is shown. (e) Replication of decreased HACE1 in advanced stage neuroblastoma using published Affymetrix U133 plus v2 array data (R2 bioinformatics tool). Bar graph is shown for INSS stage 1-4 (error bars: s.e.m.). (f) Replication of increased LIN28B in advanced stage neuroblastoma (R2 bioinformatics tool). Bar graph is shown for INSS stage 1-4 (error bars: s.e.m.). (g) Decreased expression of HACE1 is associated with worse outcome in dataset from e. Kaplan-Meier is analysis shown, with patients grouped by tertiles of HACE1 expression. Log rank P-value is shown. (h) Increased expression of LIN28B is associated with worse outcome in dataset from f. Kaplan-Meier analysis is shown, with patients grouped by tertiles of LIN28B expression. Log rank P-value is shown. ****P < 0.0001, ***P < 0.001, **P < 0.01, *P < 0.05.

In conclusion, here we have identified common variants within HACE1 and LIN28B that are associated with neuroblastoma. Similar to other genes identified in our ongoing GWAS efforts[8,9,35], it is likely that the germline susceptibility variants identified here play an important role not only in tumor initiation but also in disease progression via cis-effects on major cancer genes. Further study of HACE1 and LIN28B in neuroblastoma is warranted and may lead to new insights into the genetic and epigenetic mechanisms underlying an aggressive clinical phenotype.