Meta-analysis identifies four new loci associated with testicular germ cell tumor.

We conducted a meta-analysis to identify new susceptibility loci for testicular germ cell tumor (TGCT). In the discovery phase, we analyzed 931 affected individuals and 1,975 controls from 3 genome-wide association studies (GWAS). We conducted replication in 6 independent sample sets comprising 3,211 affected individuals and 7,591 controls. In the combined analysis, risk of TGCT was significantly associated with markers at four previously unreported loci: 4q22.2 in HPGDS (per-allele odds ratio (OR) = 1.19, 95% confidence interval (CI) = 1.12-1.26; P = 1.11 × 10(-8)), 7p22.3 in MAD1L1 (OR = 1.21, 95% CI = 1.14-1.29; P = 5.59 × 10(-9)), 16q22.3 in RFWD3 (OR = 1.26, 95% CI = 1.18-1.34; P = 5.15 × 10(-12)) and 17q22 (rs9905704: OR = 1.27, 95% CI = 1.18-1.33; P = 4.32 × 10(-13) and rs7221274: OR = 1.20, 95% CI = 1.12-1.28; P = 4.04 × 10(-9)), a locus that includes TEX14, RAD51C and PPM1E. These new TGCT susceptibility loci contain biologically plausible genes encoding proteins important for male germ cell development, chromosomal segregation and the DNA damage response.

In the United States, testicular germ cell tumors (TGCT) are the most common cancers in young men, with a peak incidence among those aged 25 to 34 years. The incidence of TGCT has more than doubled among white men in the United States over the past 30 years; similar increases in incidence rates have been observed in other populations of European ancestry[1-3]. Of note, the incidence of TGCT varies widely between populations and is much higher in individuals of European compared to African ancestry[2]. Established risk factors for TGCT include family history, cryptorchidism, adult height and prior TGCT history; several recent studies also have implicated marijuana use[4-7]. First degree relatives of affected men have been shown consistently to have an increased TGCT risk (5- to 19-fold for brothers and 2- to 4-fold for fathers)[8-11], the highest for any cancer. Further, the estimated heritability of TGCT is third among all cancers, with genetic effects estimated to account for 25% of TGCT susceptibility[12]. These observations, coupled with twin studies[13-15], support a strong genetic component contributing to TGCT susceptibility.

Despite the greatly increased relative risk of TGCT in family members, candidate gene and linkage approaches yielded little progress in identifying specific genetic risk factors. Initially, two independent genome wide association studies (GWAS) identified allele variation within KITLG on 12q22 as the strongest genetic risk factor for TGCT, with a per allele odds ratio (OR) greater than 3[16,17]. Variants on 5p15.33 (TERT- two independent loci), 5q31.3 (SPRY4), 6p21.3 (BAK1), 9p24.3 (DMRT1- two independent loci), and 12p13.1 (ATF7IP) also have been associated with TGCT risk[16-21]. The per allele ORs for the identified TGCT susceptibility alleles are in large part higher than those identified for other cancers, which may be due, in part, to the homogeneity of the disease, as all TGCT are thought to arise from the primordial germ cell[22,23]. Multiple additional loci are expected to contribute to susceptibility as has been shown for cancers of lower heritability[24]. Combining multiple GWAS represents a step to increase power to detect additional genetic risk factors failing to reach genome-wide significance in individual studies.

We performed a meta-analysis of the most promising 340 SNPs (after excluding previously reported loci) observed in the adjusted pooled analysis of the combined NCI scan (STEED, US Servicemen’s Testicular Tumor Environmental and Endocrine Determinants Study; and FTCS, NCI Familial Testicular Cancer Study) with the previously reported University of Pennsylvania (UPENN) TGCT scan (Online Methods). Allelic ORs for known loci are shown in Supplementary Table 1 for the combined NCI scan. Forty SNPs from nine loci had P values below 10−4, of which 12 localized to the MAD1L1 gene locus (7p22.2) (details of correlation between NCI and UPENN study for top 40 SNPs in Supplementary Table 2). The most significant SNP marker from each of nine loci, plus eight additional markers were selected for replication (n=17). An in silico analysis of these 17 SNPs was performed in the GWAS data from the University of Southern California (USC) and the UK Testicular Cancer Collaboration (UKTCC)[18], followed by genotyping in four additional TGCT case-control studies from: Fred Hutchinson Cancer Center (Adult Testicular Lifestyle and Blood Specimen [ATLAS] study), University of Pennsylvania (Testicular Cancer in Philadelphia Area Counties [TestPAC] study), Oslo University Hospital-Radium Hospital, Norway (OUHRH), and MD Anderson Cancer Center (MDA). Details of each study are included in the Supplementary Note. The combined analysis included 4,142 TGCT cases and 9,566 controls (Supplementary Table 3). In the combined meta-analysis, we observed four new loci significantly associated with TGCT (P value < 5×10−8) (Table 1; Supplementary Table 4).

Table 1

Meta-analysis identified novel loci associated with testicular germ cell tumor

SNP1	Nearby genes	Study2	Cases	Controls	EAF3	Allelic OR (95% CI)	P value	P for heterogeneity
rs17021463	HPGDS	NCI	582	1055	0.475	1.33 (1.14–1.55)	2.12 × 10−4
G|T	4q22.2	UPENN	349	919	0.433	1.28 (1.07–1.53)	6.67 × 10−3
		Discovery	931	1974		1.31 (1.16–1.47)	6.09 × 10−6	0.750
		UKTCC	979	4947	0.420	1.19 (1.08–1.31)	5.66 × 10−4
		USC	358	258	0.428	1.09 (0.82–1.46)	0.464
		OUHRH	798	377	0.399	1.16 (0.98–1.39)	0.092
		TestPAC	267	575	0.417	1.14 (0.92–1.41)	0.225
		ATLAS	297	664	0.420	1.11 (0.91–1.35)	0.292
		MDA	236	351	0.412	1.01 (0.79–1.27)	0.960
		Replication	2935	7172		1.15 (1.07–1.23)	7.01 × 10−5	0.862
		Combined	3866	9146		1.19 (1.12–1.26)	1.11 × 10−8	0.583
rs12699477	MAD1L1	NCI	582	1056	0.412	1.31 (1.13–1.52)	4.64 × 10−4
T|C	7p22.3	UPENN	349	919	0.407	1.37 (1.15–1.63)	4.25 × 10−4
		Discovery	931	1975		1.34 (1.19–1.50)	7.16 × 10−7	0.704
		UKTCC	979	4947	0.380	1.17 (1.06–1.30)	1.34 × 10−3
		USC	358	258	0.352	1.32 (0.99–1.74)	0.024
		TestPAC	266	573	0.375	1.09 (0.88–1.35)	0.440
		ATLAS	298	671	0.373	1.08 (0.89–1.32)	0.480
		Replication	1901	6449		1.16 (1.07–1.25)	2.41 × 10−4	0.645
		Combined	2832	8424		1.21 (1.14–1.29)	5.59 × 10−9	0.318
rs4888262	RFWD3	NCI	582	1052	0.552	1.36 (1.17–1.59)	9.12 × 10−5
T|C	16q22.3	UPENN	349	919	0.498	1.39 (1.16–1.67)	3.12 × 10−4
		Discovery	931	1971		1.37 (1.22–1.54)	1.39 × 10−7	0.858
		UKTCC	957	4940	0.458	1.22 (1.10–1.34)	1.02 × 10−4
		USC	358	258	0.440	1.41 (1.09–1.83)	5.14 × 10−3
		TestPAC	260	575	0.465	1.12 (0.91–1.38)	0.2876
		ATLAS	295	649	0.495	1.16 (0.95–1.41)	0.1328
		Replication	1870	6422		1.21 (1.12–1.31)	1.62 × 10−6	0.568
		Combined	2801	8393		1.26 (1.18–1.34)	5.15 × 10−12	0.397
rs9905704	TEX14	NCI	582	1054	0.719	1.37 (1.16–1.61)	1.88 × 10−4
G|T	17q22	UPENN	349	919	0.663	1.30 (1.08–1.56)	5.46 × 10−3
		Discovery	931	1973		1.33 (1.19–1.52)	3.44 × 10−6	0.674
		UKTCC	979	4947	0.680	1.28 (1.16–1.43)	3.65 × 10−6
		USC	358	258	0.649	1.35 (1.08–1.72)	0.017
		OUHRH	802	382	0.695	1.23 (1.02–1.49)	0.028
		TestPAC	259	575	0.669	1.14 (0.91–1.43)	0.253
		ATLAS	300	666	0.669	1.09 (0.88–1.33)	0.403
		MDA	234	351	0.672	1.15 (0.89–1.47)	0.273
		Replication	2932	7179		1.23 (1.15–1.32)	1.37 × 10−8	0.628
		Combined	3863	9152		1.27 (1.18–1.33)	4.32 × 10−13	0.668
rs7221274	PPM1E	NCI	582	1056	0.660	1.39 (1.19–1.61)	3.56 × 10−5
G|A	17q22	UPENN	349	919	0.598	1.19 (1.00–1.43)	0.053
		Discovery	931	1975		1.30 (1.16–1.47)	7.79 × 10−6	0.197
		UKTCC	979	4947	0.620	1.23 (1.12–1.37)	3.08 × 10−5
		USC	358	258	0.558	1.37 (0.87–2.13)	9.43 × 10−3
		OUHRH	802	380	0.630	1.22 (1.02–1.47)	0.029
		TestPAC	243	538	0.612	1.01 (0.81–1.27)	0.890
		ATLAS	301	671	0.604	1.00 (0.83–1.22)	0.989
		MDA	215	351	0.615	1.05 (0.83–1.33)	0.682
		Replication	2898	7145		1.16 (1.09–1.25)	3.27 × 10−5	0.228
		Combined	3829	9120		1.20 (1.12–1.28)	4.04 × 10−9	0.130

SNP genotype depicted as reference allele|effect allele

NCI depicts combined analysis results of the two GWAS scans STEED and FTCS performed at NCI

Discovery depicts initial meta-analysis of NCI and UPENN and Replication depicts meta-analysis of the rest

EAF depicts effect allele frequency in control

The most significant 4q22.2 SNP marker, rs17021463, is located within the intron of the hematopoietic prostaglandin D synthase gene, HPGDS (P = 1.11×10−8, OR 1.19, 95%CI 1.12–1.26) (Figure 1A, Table 1). In mice, hpgds is expressed in the early embryonic male gonad and appears to regulate nuclear localization of the sox9 protein[25]. Disruption of hpgds leads to modification of the phenotype of apcMin/+ mice[26]. Seventy-one surrogate markers highly correlate with HPGDS rs17021463 (r ≥ 0.8, 1000 Genomes CEU data, Supplementary Table 5). Notably, rs35744894 (r = 0.87) changes a DMRT2 binding motif (Supplementary Table 6); variation in DMRT1 has been associated with TGCT risk[19].

Figure 1

Recombination plot and linkage disequilibrium structure for the four new TGCT susceptibility regions at 4q22.2, 7p22.3, 16q22.3 and 17q22 (a–d)

Regional plots of association results, recombination hotspots and linkage disequilibrium for the (a) 4q22.2–22.3:94,904,738–95,514,609, (b) 7p22.3:1,651,900–2,479,029, (c) 16q23.1:74,179,928–74,812,676 and (d) 17q22–23.1:56,083,934–57,680,480 TGCT susceptibility loci. (a–d) Combined meta-analysis results are shown as red diamonds with rs numbers labeled, and the NCI scan in gray. For each plot, −log10P values (y axis, left) of the SNPs are shown according to their chromosomal positions (x axis). Linkage disequilibrium structure based on NCI controls (n=1,188) was visualized by snp.plotter software. The line graph shows likelihood ratio statistics (y axis, right) for recombination hotspot by SequenceLDhot software and five different colors represent 5 tests of 100 controls from NCI without resampling. Physical locations of each region are based on NCBI Build 37 of the human genome.

Fifty-three of 71 surrogate markers that were highly correlated with HPGDS rs17021463 (r ≥ 0.8, 1000 Genomes CEU data, Supplementary Table 5) across a 200kb window mapped within or near an adjacent gene, SMARCAD1 (SWI/SNF-related, matrix-associated actin-dependent regulator of chromatin, subfamily a, containing DEAD/H box 1). SMARCAD1 is a chromatin remodeler, which restores silenced heterochromatin domains in dividing cells and participates in DNA damage response[27,28]. Homozygous mutant mice display developmental defects, including impaired fertility[29]. Surrogate markers included one nonsynonymous substitution, rs7439869, at codon 301 (r = 0.93, 1000 Genomes CEU, T>C, Val>Ala). Although it is predicted to be tolerated by PolyPhen2[30], it changes an OCT4 (POUF5F1) and SOX4 binding motif (Supplementary Table 6). OCT4 is a transcription factor, which regulates pluripotency in a number of cell types, including primordial germ cells, and is expressed in TGCT[31-36].

We identified a locus on 7p22.3, harboring mitotic arrest deficient-like 1 (MAD1L1) gene, which encodes MAD1. The most significant SNP without study heterogeneity, rs12699477, localized in intron 17 (P = 5.59 × 10−9, OR 1.21, 95%CI (1.14–1.29)) (Figure 1B). Of note, the risk allele (C) at rs12699477 is more prevalent in populations of European (29%) than those of African ancestry (8%) in 1000 Genomes[37]. MAD1 is a spindle assembly checkpoint protein that delays the onset of anaphase in the mitotic cell cycle until all sister chromatids achieve proper alignment and microtubule attachment, thereby preventing aneuploidy and maintaining genomic stability[38].

Among the 35 SNPs that are highly correlated with MAD1L1 rs12699477 (r ≥ 0.7, 1000 Genomes CEU data, Supplementary Table 5), rs1801368 is a missense mutation at codon 558 (G>A, Arg>His) that resides in the MAD1L1 second leucine zipper domain. Arg558His has been reported to be associated with lung cancer risk[39] and may lead to reduced binding of MAD2 to MAD1, resulting in decreased proficiency in enforcing mitotic arrest[40]. We observed additional statistically significant associations with TGCT for neighboring SNPs in the MAD1L1 region, including rs10275045 (P=3.78×10−10, OR 1.20, 95%CI (1.13–1.27)) and rs3778991 (P=6.73×10−10, OR 1.21, 95%CI (1.14–1.28)). However, both displayed significant study heterogeneity (Supplementary Table 4). The r between our strongest signal at rs12699477 and these markers is 0.66 and 0.50, respectively, in the STEED controls. A conditional analysis resulted in a marked attenuation of the signal, supporting a single TGCT susceptibility locus across MAD1L1 on 7q22.3 (Supplementary Table 7).

We observed a significant TGCT association with rs4888262 on 16q22.3 (P = 5.15×10−12, OR 1.26, 95%CI 1.18–1.34), which is a synonymous SNP in codon 404 (G>A, Thr) of the ring finger WD domain 3 (RFWD3) (Table 1, Figure 1C). RFWD3 is an E3 ubiquitin ligase that positively regulates p53 stability by forming a RFWD3-MDM2-p53 complex, thereby protecting p53 from degradation by MDM2 polyubiquitination[41,42]. Within the LD interval are SNPs that map to two additional genes, the golgi glycoprotein 1 (GLG1) and mixed lineage kinase domain-like (MLKL); the latter of which has been recently identified as a key of mediator of TNF-induced necrosis, downstream of receptor interacting protein kinase 3 (RIP3)[43,44] (Figure 1C). We note that rs3851729, which is highly correlated with rs4888262 (r = 0.77, 1000 Genomes CEU), maps to a highly conserved sequence in the 3′ UTR of GLG1; similarly, rs4072222 (r = 0.87, 1000 Genomes CEU) maps to an intron of MLKL (Supplementary Table 5). Both susceptibility variants are cis-eQTLs that influence MLKL and RFWD3 expression in monocytes[45].

We identified two highly correlated SNPs (r = 0.74 in the STEED controls) on 17q22, rs9905704 (P = 4.32×10−13, OR 1.27, 95%CI 1.18–1.33) and rs7221274 (P = 4.04×10−9, OR 1.20, 95%CI 1.12–1.28) (Table 1, Figure 1D). In a conditional analysis, the signal at one SNP was markedly attenuated by the other, indicating a single 17q22 TGCT susceptibility locus (Supplementary Table 7). Within this LD block are at least six plausible candidate genes: RAD51C (RAD51 homolog C [S. cerevisiae]), TEX14 (testis expressed 14), PPM1E (protein phosphatase, Mg2+/Mn2+ dependent, 1E), SEPT4 (septin 4), TRIM37 (tripartite motif containing 37), and SKA2 (spindle and kinetochore associated complex subunit 2) (Figure 1D). Proteins encoded by these candidate genes, except for SKA2, have been implicated as having roles in spermatogenesis[46-51]. RAD51C is a DNA repair gene, in which rare mutations confer susceptibly to ovarian cancer[52,53]. Of male rad51c mice, approximately one-third were found to be infertile due to impaired spermatogenesis[49]. TEX14 is an essential component of germ cell intercellular bridges, evolutionarily conserved structures from invertebrates to humans that allows clonal development of daughter cells in syncytium; targeted disruption of Tex14 results in male sterility in mice[48]. TEX14 also has been implicated as an important component of kinetochores (KTs) and interacts with MAD1[54]. PPM1E encodes a phosphatase that dephosphorylates to switch-off CaMK4 (calcium/calmodulin-dependent protein kinase IV), deficiency of which causes infertility in mice[50,55]. TRIM37 encodes a RING-B-box-coiled-coil protein; rare mutations in this gene cause the autosomal recessive disease mulibrey nanism (MUL; MIM 253250)[56], in which adult males have testicular failure[57]. Three SNPs - rs8077332, rs11652713, and rs9898048 - map within TRIM37 and are in perfect LD with rs7221274 (r = 1, 1000 Genomes CEU, Supplementary Table 5); all are cis-eQTL affecting RAD51C expression in monocytes[45]. Thus, fine mapping and functional studies will be required to elucidate the biological basis of the association signal in this interval on 17q22.

In our meta-analysis of GWAS studies, we have identified four new TGCT susceptibility loci at 4q22, 7q22, 16q22.3, and 17q22. In total, 10 loci now have been conclusively associated with TGCT susceptibility. The four newly identified susceptibility alleles account for 2% of the risk to the brothers and 3% of risk to the sons of TGCT patients, increasing the cumulative total of 12 susceptibility alleles (two susceptibility alleles from TERT-CLPTM1L [5p15] and two from DMRT1 locus [9p24]) to 14% and 21% of the risk to brothers and sons, respectively. Based on the high heritability of TGCT, more than one hundred additional loci are expected to be discovered[24]. Notably, the allelic ORs associated with these novel loci are in the range of 1.2 to 1.3, continuing the trend of identifying loci with higher odds ratios for TGCT than for other cancer types[23].

Interestingly, each locus harbors biologically plausible candidate genes implicating several pathways – most strikingly, spermatogenesis and male germ cell development (HPGDS, SMARCAD1, SEPT4, TEX14, RAD51C, PPM1E, TRIM37), chromosomal segregation (MAD1L1, TEX14, SKA2), and DNA damage response (SMARCAD1, RFWD3, RAD51C). None of the four newly identified loci have been previously implicated in GWAS of other cancers, further supporting that there are distinct pathways and regions implicated in TGCT susceptibility; however rare mutations in RAD51C have been implicated in ovarian cancer susceptibility[53]. TGCT susceptibility is particularly unique in that many of the associated genes affect male germ cell development and differentiation, thus emphasizing the potential detrimental effect that inherited variation in this developmental process can have on the tumorigenic potential of the primordial germ cell.

This study is the first to implicate variation within genes involved in chromosomal segregation as associated with cancer susceptibility. TGCT karyotypes are unique among cancers, in that nearly all carry the same chromosomal aberration, a gain of 12p, most often in the form of an isochromosome, which is considered essential for tumor development[58-60]. Variation in these genes may lead to chromosomal instability and facilitate the development of aneuploidy. Numerous potential regulatory SNPs were identified, suggesting that newly identified associations might be mediated by plausible genes within each locus, which warrant further fine-mapping and functional studies to elucidate the biological bases of the TGCT susceptibility regions. Studies of the genetic basis of TGCT continue to provide novel insights into this unique disease with high heritability.

ONLINE METHODS

Studies

Detailed characteristics of the study populations are described in both the Supplementary Note and Supplementary Table 3. Subjects used in the current study are all of European descent and data from each study were collected and analyzed in accordance with local ethical permissions and informed consent. Three studies (STEED, FTCS, and UPENN) were included in the discovery meta-analysis, and six studies contributed to replication by de novo genotyping (TestPAC, ATLAS, OUHRH, and MDA) or in silico look-up in existing data (UKTCC and USC).

Genotyping and quality control

Genotype quality control metrics for the reported GWAS scans (UPENN and UKTCC) were previously described[18,19]. Genotype quality control metrics for STEED, FTCS, and USC are described in Supplementary Note[61]. OUHRH and MDA studies were genotyped using the 5′ exonuclease assay (TaqMan™) and the ABI prism 7900HT sequence detection system, all according to the manufacturer’s instructions, across several genotyping centers. Primers and probes were supplied directly by Applied Biosystems as Assays-By-Design™. Technical validation was performed in the HapMap samples (n=270) with greater than 99% genotype concordance. TestPAC and ATLAS studies conducted genotyping using the iPLEX mass array platform (Sequenom, Inc.) following manufacturer’s protocol. Assays at all genotyping centers included at least four negative controls and 2–5% duplicates on each plate. Standard quality control protocol was implemented; SNP call rate > 95%, no deviation from Hardy-Weinberg equilibrium in controls at P<0.00001, <2% discordance between genotypes in duplicate had to be fulfilled and cluster plots for SNPs that were close to failing any of the QC criteria were re-examined centrally.

Statistical analysis

Two genome-wide scans from the National Cancer Institute (STEED and FTCS) were analyzed as a combined dataset using a logistic regression model for trend effect adjusted for age, study, and additionally for one eigenvector (only one with p < 0.05) to account for population stratification in this European population. From the top 500 SNPs by trend P values from the NCI scan excluding previously reported ones, 340 SNPs were selected based on the availability of surrogates (r > 0.6) in the previous TGCT GWAS scan from the University of Pennsylvania. Since SNP content differs between the Illumina and Affymetrix platforms, the best correlated surrogate per each marker was paired to perform a discovery meta-analysis (111 SNPs, direct match; 229 SNPs, surrogate match). From the discovery meta-analysis, 17 of 40 SNPs with P values < 10−4 were selected for follow up in the remaining studies. In silico follow-up was done in the USC and UKTCC scans, whereas additional genotyping was done in TestPAC, ATLAS, OUHRH and MDA studies (Supplementary Table 3). Not all markers were available for replication efforts from all sites (see Supplementary Table 4).

The meta-analysis was conducted using the suite of tools in GLU (Genotyping Library and Utilities) software, combining study-specific odds ratio (OR) estimates using a fixed effects model, which used the inverse-variance method to estimate the combined OR and its 95% confidence intervals (CIs). To assess existence of heterogeneity among studies, Cochran’s Q statistic was used to calculate P for heterogeneity.

Recombination hotspots were identified in the vicinity of the novel TGCT associated loci using SequenceLDhot[62], a program that uses the approximate marginal likelihood method[63] and calculates likelihood ratio statistics at a set of possible hotspots. We tested five unique sets of 100 control samples drawn from STEED. PHASE v2.1 program was used to calculate background recombination rates[64,65] and LD heatmap was visualized in r using snp.plotter program[66].

The relative risk attributable to a set of SNPs (λ) was estimated using the following formula[67]

where q is the minor allele frequency of SNPi and p = 1 − q. SNP specific risks for rare homozygotes, heterozygotes, and common homozygotes are denoted by r0, r1, and r2, respectively. The NCI controls (n=1,140) were used to estimate minor allele frequencies and odds ratio estimates from SNP association analyses were used to estimate relative risks. This formula assumes the effects of all SNPs in the set are multiplicative. The proportion of familial risk attributable to a set of SNPs was calculated as , where λ0 is the familial relative risk estimated from TGCT epidemiological studies (λ0 =4 for affected father, λ0 =8 for affected brother)[68].

Genomic annotation

Genomic annotation on high LD surrogates (r ≥ 0.8, 1000 Genomes CEU) of 5 SNPs (rs17021463, rs12699477, rs4888262, rs9905704, and rs7221274) from the four TGCT susceptibility loci identified in the current study was conducted using ENCODE tools – HaploReg[69] and RegulomeDB[70] (Supplementary Table 5). rs12699477 did not have surrogates with r ≥ 0.8 threshold, thus we lowered the threshold to 0.7 for surrogates, and then conducted annotation. All surrogates were queried in RegulomeDB browser to cross examine predicted regulatory DNA elements such as regions of DNase hypersensitivity, binding sites of transcription factors, and promoter regions that have been biochemically characterized to regulation transcription. Summaries of each SNP analysis by RegulomeDB browser expressed in scores are added to Supplementary Table 5. To predict potential regulatory SNPs, we assessed SNPs that meet one of the following criteria - 1) conserved (GERP and/or Siphy); 2) present in a promoter or DNase hypersensitivity region; or 3) predicted to have a cis eQTL or having a RegulomeDB score of ≤ 3. Twenty-nine SNPs that passed one of these criteria also changed a motif, and are annotated further with the motif of interest and their log-odds (LOD) motif score for the specific SNP of interest in Supplementary Table 6. Two SNPs in 3′-UTR regions were evaluated using SNP Function Prediction for changes in miRNA binding sites and are included in Supplementary Table 6.