Genetic Analysis in African American Children Supports Ancestry-Specific Neuroblastoma Susceptibility

Neuroblastoma is an important pediatric cancer characterized by a high degree of clinical heterogeneity. Although a substantial proportion of patients shows a favorable outcome and may even experience spontaneous tumor regression, patients with neuroblastoma classified as high-risk show an aggressive clinical course with widespread metastases to bone and bone marrow at diagnosis (1). This group of children still has survival rates of less than 35% despite very aggressive therapy (2).

The age-adjusted incidence rate of neuroblastoma in North Americans of European descent [European Americans (EA)] is approximately 11.5 per million. Neuroblastoma is reported to be rarer among African American (AA) children, with an age-adjusted incidence of 8.5 per million (3), and odds ratio for having a child with neuroblastoma by parental origin of 0.74 [95% confidence interval (CI), 0.56–0.96] for AA parents compared with EA parents (4). According to the last available data with both age and race/ethnicity information in the Surveillance, Epidemiology, and End Results database (https://seer.cancer.gov/), for ages less than 1 year, the incidence of neuroblastoma is 58 per million in EAs (95% CI, 52–64) versus 30 per million in AAs (95% CI, 21–42); however, incidence rates in the two groups are similar for older ages (20 vs. 21 for ages 1–4; 4 vs. 3 for ages 5–9; 1 vs. 2 for ages 10–14). Older age at diagnosis is an adverse prognostic factor in neuroblastoma, and AA children with neuroblastoma are more likely to have the high-risk form of the disease compared with EA patients (57% vs. 44%), with associated lower overall survival probability (5-year event-free survival, 56% in AAs and 67% in EAs; ref. 5).

Although outcome disparities might also reflect differences in access to care or delays in medical assistance or adherence to treatment (6), numerous studies have shown that in neuroblastoma, tumors with favorable biology are unlikely to progress to unfavorable biology tumors (7–9).

Global African genomic ancestry estimated from genome-wide genetic data was found to be significantly associated to high-risk neuroblastoma (10, 11). Limited data exist on the incidence of neuroblastoma in Africa, but reports suggest it is similar to that of AAs or lower (3). More recent analyses support a lower frequency overall but higher proportion of metastatic disease at presentation, characteristic of the high-risk form, in sub-Saharan Africa compared with Northern Africa, where epidemiologic data on neuroblastoma are similar to those from Europe and North America (12).

Genome-wide association studies (GWAS) have identified several loci containing common variants associated to risk of developing neuroblastoma (13–20), but attempts at replicating these findings in a cohort of AA neuroblastoma cases and controls have produced mixed results (21). In fact, only limited information exists on genetic predisposition to neuroblastoma in AAs (21), and this information has mostly come from attempts at replicating susceptibility loci identified in the EA GWAS.

We present here a detailed genetic analysis of neuroblastoma in the largest cohort of AA patients available so far, aimed at the characterization of genetic risk factors in an underserved population which has been so far mostly excluded from genetic studies, and is disproportionately affected by the more severe form of this important pediatric cancer.

We analyzed an AA cohort consisting of 629 neuroblastoma cases and 2,990 controls, all genotyped on Illumina SNP arrays (HumanHap550, Human610, and OmniExpress). DNA for all neuroblastoma cases was obtained through the Children's Oncology Group. Samples were annotated with clinical and genomic information, including age at diagnosis, event-free and overall survival, disease stage, MYCN oncogene copy number, DNA index, and risk group classification (low, intermediate, high, or unknown). Controls were children recruited from the Philadelphia region through the CHOP Health Care Network with no serious underlying medical disorder including cancer. For both AA cases and controls, ancestry was initially determined based on self-report or parental report, and then confirmed by principal component analysis (PCA) using reference data from the HapMap3 founder populations (22). Two EA cohorts were also used for polygenic score analyses. Specifically, we used results described by McDaniel and colleagues (20) in a GWAS performed in a discovery cohort (EA1) consisting of 2,101 cases and 4,202 controls assayed with Illumina HumanHap550 and Human610 SNP arrays, and a validation cohort (EA2) comprised of 1,279 cases and 1,443 controls genotyped with Illumina OmniExpress SNP arrays (Table 1; Supplementary Table S1). Using PLINK (23), we performed standard GWAS quality control to remove SNPs with low call rates (<97%), significant deviation from Hardy–Weinberg equilibrium in controls (P < 10⁻⁵), and significant difference in missing rates between cases and controls (P < 10⁻⁵), as well as individuals with low call rates (<97%), inconsistencies between reported sex and sex chromosome genotypes, and duplicates or relatives up to third degree determined from estimated proportions of identity-by-descent sharing (24). PCA was performed on a subset of independent SNPs (r² < 0.2) after merging our data with the HapMap3 data to confirm the genetic ancestry of all samples and identify potential outliers or batch effects. Proportion of genome-wide African and European ancestries were estimated by a supervised procedure with two founder populations and the 1000 Genome CEU and YRI populations as surrogates for the two ancestral groups using ADMIXTURE (25). Individuals with <10% African ancestry were removed from the AA cohort, and individuals with < 90% European ancestry were removed from the EA2 cohort.

We used the Michigan Imputation Server (26) to perform imputation on the AA and EA2 cohorts, using reference panels from Consortium on Asthma among African-ancestry Populations in the Americas (www.caapa-project.org) and 1000G Phase III v5 (https://mathgen.stats.ox.ac.uk/impute/1000GP_Phase3.html) respectively. Eagle (27) was used to perform phasing. Genotypes from Illumina OmniExpress arrays were imputed separately from the other Illumina arrays. Following imputation, variants with minor allele frequency <1% and imputation accuracy (R²) <0.3 were removed and the remaining variants were tested for association with neuroblastoma by logistic regression analysis using SNPTEST version 2.5.4-beta3 (28), including the genome-wide proportion of African ancestry as covariate from the admixture analysis with HapMap founder populations CEU and YRI. Association analysis was performed in AA first including all cases and then selecting the subset of high-risk cases only. After association analysis, SNPs with an info score <0.8 were removed, leaving 6,749,711 SNPs (8,524,982 on OmniExpress arrays) in AA and 7,710,745 in EA2.

The summary statistics from OmniExpress and from the other illumina arrays in AA were then merged through inverse variance weighted meta-analysis with the METAL software (version released on 2011–03–25; ref. 29). Only SNPs with the same direction of effect were considered after meta-analysis, leaving 3,073,914 and 3,064,236 variants in all cases and in the subset of high-risk cases only analysis respectively. For new loci, standard GWAS thresholds of P value < 5 × 10⁻⁸ for genome-wide and P value < 10⁻⁵ for suggestive significance were applied. Regional association plot at candidate associated SNPs were obtained using LocusZoom (30).

We used RegulomeDB (31) to annotate candidate neuroblastoma-associated SNPs and SNPs in high linkage disequilibrium (LD) with them (r² > 0.8) with known or predicted regulatory elements. RegulomeDB assigns each SNP a ranking score from 1 (high) to 7 (low) for evidence of regulatory function, as well as a probability score ranging from 0 to 1, with 1 being most likely to be a regulatory variant. We also investigated transcriptional activity through the H3K27Ac data generated from BE2C neuroblastoma cell lines (32).

Cox regression analysis was performed with the CRAN survival package (33) on the candidate associated lead SNPs in previously unknown loci to see whether the presence of the risk allele in the genotype had an effect on the duration of event-free survival or overall survival (survival information was available for 509 patients).

After genotype imputation, the most likely genotype at each SNP with posterior probability >0.9 assigned by the Michigan Imputation Server to each subject was also used in local admixture estimation. Haplotypes were phased using the Eagle software (27) and Rfmix v2.03 (34) was used to assign local ancestry to contiguous windows of SNPs using the 1000 Genome CEU and YRI populations as surrogates for the two ancestral populations. Rfmix was used separately on vcf files from samples genotyped with Illumina OmniExpress arrays and from the other Illumina SNP arrays. The most likely ancestry at each window predicted by Rfmix was used to identify genomic regions with significant difference in proportion of African versus European ancestry in cases versus controls and high-risk cases versus intermediate- and low-risk cases, using the two-tailed two-proportions z-test. In order to account for multiple testing, the standard threshold of statistical significance (0.05) was divided by the number of genomic windows (242), defined by the expected number of ancestry switches on each chromosome (α/n ∼ 2 × 10⁻⁴).

We used PRSice v 2.2.11.b (35) and results from the EA1 neuroblastoma discovery GWAS (20) to perform a polygenic risk score (PRS) analysis in the AA and EA2 datasets, including the genome-wide proportion of African ancestry as covariate from the supervised admixture analysis with HapMap3 founder populations CEU and YRI. Specifically, we calculated a best fit PRS by comparing results obtained including increasing number of SNPs at variable thresholds of significance (from P value = 5 × 10⁻⁸ to P value = 1), as well as a PRS based on all confirmed neuroblastoma loci (36, 37) regardless of any otherwise applied inclusion filter (such as minor allele frequency and info score) – namely 1q23: DUSP12, 2q34: SPAG16, 2q35: BARD1, 3p21: KIF15, 3q25: MLF1, RSRC1, 4p16: CPZ, 5q11: DDX4, IL31RA, 6p22: CASC15, NBAT-1, 6q16-q21: HACE1, 6q16-q21: LIN28B, 8p21: NFEL, 11p15: LMO1, 11p11: HDS17B12, 11q22: MMP20, 12p13: CDKN1B, 17p13: TP53. Only independent SNPs (r² < 0.1) were included in the PRS to avoid overrepresentation of significant loci due to multiple markers in linkage disequilibrium. Association of the PRS to neuroblastoma in the EA2 and AA cohorts was tested by logistic regression, and proportion of disease risk variance explained by the PRS was estimated by the regression coefficient R². Significance of the best fit PRS was calculated by permutation to avoid over-fitting.

The data used in this study are available through dbGaP accession number phs000124.v3.p1.

Results of the GWAS in AA are reported in Fig. 1A. Two SNPs reached the genome-wide significance threshold (P value = 5 × 10⁻⁸) in a single locus, and a total of 93 SNPs in 28 genomic regions reached the suggestive significance threshold (P value = 10⁻⁵; Supplementary Table S2). The top signal (rs114655292 and rs116331909, P value = 4.4 × 10⁻⁸) was located in 10q11.21 in the TMEM72-AS1 gene (Supplementary Fig. S1), followed by a signal in 6p22 (rs9466535, P value = 2.6 × 10⁻⁷; Supplementary Fig. S2), downstream of the known neuroblastoma-associated gene CASC15, and one in 2q35 (rs17489231, P value = 5.2 × 10⁻⁷) within the known neuroblastoma-associated gene BARD1 (Supplementary Fig. S3). In contrast to what was observed in the EA GWAS (20), where the association signal at the BARD1 locus extended to a large region upstream of BARD1, the association signal in AA was mostly included within BARD1, due to the smaller LD blocks characteristic of populations of African descent. Other known neuroblastoma loci (36, 37) that were nominally significant (P value < 0.05) in the AA cohort and showed the same direction of effect observed in the EA1 GWAS included RSRC1, LIN28B, and LMO1 (Supplementary Table S3; ref. 16). A separate analysis of the subset of 254 high-risk cases (Fig. 1B; Supplementary Table S4) yielded two genome-wide significant signals on chromosome 2, at 2q35 (Supplementary Fig. S4) within BARD1 (rs7587476, P value = 7.9 × 10⁻¹⁰), and 2p22 (Supplementary Fig. S5) within CRIM1 (rs17018667, P value = 6.2 × 10⁻⁹), in spite of the smaller case sample size. Using RegulomeDB (31), we annotated with regulatory elements the top SNPs from each of the loci identified in the AA GWAS, including SNPs in high LD with them (r² > 0.8). Besides the already known variants within the BARD1 locus, no other candidate neuroblastoma-associated SNP showed strong evidence of known or predicted regulatory functions (RegulomeDB ranking scores ≥ 4; probability scores ≤ 0.61). We also considered a recent dataset generated using H3K27Ac chromatin immunoprecipitation sequencing in BE2C cell lines (32) to further characterize the transcriptional role of the previously identified variants, but also in this case, we only found evidence of regulatory potential within the BARD1 locus. We finally performed Cox regression analysis on candidate neuroblastoma-associated lead SNPs in previously unknown loci (rs114655292 and rs116331909 in TMEM72-AS1 and rs17018667 in CRIM1) to see whether the presence of the risk allele in the genotype had an effect on event-free survival time or overall survival time, but no significant differences were found.

Global genomic African ancestry proportion in AAs was slightly higher in cases than in controls (0.77 and 0.76 respectively) and in high-risk than intermediate- and low-risk cases (0.78 and 0.74 respectively), but these differences were not statistically significant (P value > 0.05). Estimates of local admixture across the autosomes revealed only one location with significantly higher proportion of European ancestry in cases versus controls after correction for multiple tests (P value < 2 × 10⁻⁴). This region (Supplementary Fig. S6) on 4q31.22 at 4:148206926–148510474 (GRCh37/hg19 coordinates) encompasses the EDNRA gene. Comparing local ancestry in high versus intermediate- and low-risk cases, we detected one region with excess of African ancestry in high-risk cases at 17p13.1 included in the ALOX12 gene, upstream of the known neuroblastoma-associated gene TP53 (P value < 2 × 10⁻⁴ at 17:6837953–6840635 in GRCh37/hg19 coordinates, Supplementary Fig. S7).

We used summary statistics from our published GWAS in EA cases and controls (20), to derive a best fit PRS by comparing results obtained including increasing number of SNPs at variable thresholds of significance (from P value = 5 × 10⁻⁸ to P value = 1; Table 2). The best fit PRS in EA2 cases and controls (P value = 1.8 × 10⁻⁷³, empirical P value < 10⁻⁵) included all independent SNPs and explained 18.52% of the disease risk variance, and in general the significance of the PRS increased as the P value threshold became more inclusive. In contrast, the best fit PRS in AA (P value = 3.2 × 10⁻¹¹; empirical P value < 10⁻⁵) included only 22 independent SNPs (Supplementary Table S5) with association P value < 2.75 × 10⁻⁶ in the EA1 GWAS, and explained only 2.03% of neuroblastoma risk variance in AAs. A PRS based on all confirmed neuroblastoma loci explained 6% of neuroblastoma risk variance in EA2 (P value = 3.4 × 10⁻²⁷), but only 1.2% in AAs (P value = 3.3 × 10⁻⁷).

Over the past decades, GWAS have started to unravel the genetic basis of susceptibility to complex diseases. Thousands of risk loci have been identified and confirmed in independent case and control cohorts, thanks to the availability of increasing number of genotyped samples. However, GWAS, in general, have largely been performed in populations of European descent, with 78% of individuals included in the NHGRI-EBI GWAS Catalog having European ancestry, followed by 9% East Asian, and only 2.4% African (38). In spite of this, individuals of African descent have contributed to 7% of all positive, significant association findings, and therefore significantly more than in proportion to their representation in GWAS (38). Obtaining robust and replicable results in GWAS of minorities is particularly challenging for a rare disease such as neuroblastoma, which is even rarer in AA than EA children. Nonetheless, in our association analysis on a relatively small sized AA cohort of neuroblastoma patients, we were able to replicate one known neuroblastoma susceptibility locus on 2q35 (BARD1) which reached genome-wide significance in the subset of high-risk cases. At this locus, rs7587476 was the lead SNP in the AA high-risk GWAS and was first identified as a neuroblastoma susceptibility SNP by Capasso and colleagues (ref. 14; P value < 10⁻⁷) in a high-risk EA GWAS and later replicated in AAs by Latorre and colleagues (ref. 21; P value = 2.4 × 10⁻⁶). Two previously unknown loci also showed genome-wide significant association: on 10q11 in the GWAS that included all cases (in a region harboring lincRNA TMEM72-AS1); and on 2p22 in the GWAS that included the subset of high-risk cases only, within the CRIM1 gene, which encodes a transmembrane protein that plays a role in tissue development. Interestingly, the top SNPs at these two previously unknown neuroblastoma candidate loci (rs114655292 and rs116331909 in TMEM72-AS1 and rs17018667 in CRIM1) are not polymorphic in the 1000 Genome European population, suggesting that, if confirmed, these associations may be linked to neuroblastoma susceptibility loci of African origin: indeed, the EA1 and EA2 GWAS show no association signals in these regions.

Admixture mapping can be used to map susceptibility loci for disorders with different prevalence among founder populations that became recently admixed (39). Given that neuroblastoma is more frequent in European compared with African populations, we estimated the local proportions of African and European admixture along the autosomal genome in the AA cohort, and looked for regions that showed an increase in European ancestry among cases compared with controls and thus may harbor neuroblastoma risk loci. Only one genomic region at 4q31.22 (4:148206926–148510474), containing the EDNRA gene, reached statistical significance after correction for the number of genomic intervals tested. EDNRA encodes the receptor for endothelin-1, that was originally identified as a vasoconstrictor peptide and has later also been shown to be involved in the development of cranial/cardiac neural crest-derived tissues and organs (40) and to regulate multiple aspects of neural cell crest development through evolutionary conserved pathways (41). However, we did not find significant association of any SNP from this region to neuroblastoma in the EA1 GWAS.

A previous case-only study in a subset of the same patients found that African genetic ancestry was significantly associated with high-risk neuroblastoma phenotype (21). We thus hypothesized that genomic regions with an increase in proportion of African ancestry in high-risk compared with non–high-risk cases may harbor loci predisposing to high-risk neuroblastoma. From this analysis, we found only one statistically significant region at 17:6837953–6840635, harboring the ALOX12 gene, upstream of the known neuroblastoma-associated gene TP53. We did not find suggestive association (P value < 10⁻⁵) to neuroblastoma of any SNP from this region in either AA or EA1 GWAS.

At a global level, PRS analysis indicated limited shared genetics underlying neuroblastoma risk in European and AA patients: significantly associated SNPs have sometimes different allelic frequencies in EA and AA, and they belong to blocks with variable LD structure. Several studies have pointed out that polygenic scores have limited ability to accurately predict disease risk when applied to populations of different genetic ancestry, thus further supporting the need for GWAS to be performed in populations of all races and ethnicities (42, 43).

Together these results strongly support the hypothesis that additional common variants, besides those identified in GWAS performed in children of European ancestry, contribute to risk of neuroblastoma in AA children. Additional efforts to recruit and study these patients are warrant and will allow identification of new candidate loci that open the way to future studies to confirm and characterize their functional role in neuroblastoma development and progression, with the final goal of improving care and treatment of all neuroblastoma patients.

No disclosures were reported.

A. Testori: Software, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. Z. Vaksman: Data curation. S.J. Diskin: Resources. H. Hakonarson: Resources, funding acquisition, writing–review and editing. M. Capasso: Writing-review and editing. A. Iolascon: Funding acquisition. J.M. Maris: Funding acquisition, writing–review and editing. M. Devoto: Conceptualization, formal analysis, supervision, investigation, methodology, writing–original draft, project administration, writing–review and editing.

This study was supported by grants from the NIH R01 CA124709 (to J.M. Maris), R01 CA180692 (to J.M. Maris), and R35 CA220500 (to J.M. Maris), and from Associazione Italiana per la Ricerca sul Cancro (grant number 20757 to A. Iolascon).

This study was also supported by the Giulio D'Angio Endowed Chair (to J.M. Maris) and by the Children's Hospital of Philadelphia Endowed Chair in Genetics Research (to H. Hakonarson).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.