Exome-chip meta-analysis identifies novel loci associated with cardiac conduction, including ADAMTS6

RESEARCH Open Access Exome-chip meta-analysis identifies novel loci associated with cardiac conduction, including ADAMTS6 Bram P. Prins 1,2† , Timothy J. Mead 3† , Jennifer A. Brody 4 , Gardar Sveinbjornsson 5 , Ioanna Ntalla 6,7 , Nathan A. Bihlmeyer 8 , Marten van den Berg 9 , Jette Bork-Jensen 10 , Stefania Cappellani 11 , Stefan Van Duijvenboden 6,12 , Nikolai T. Klena 13 , George C. Gabriel 13 , Xiaoqin Liu 13 , Cagri Gulec 13 , Niels Grarup 10 , Jeffrey Haessler 14 , Leanne M. Hall 15,16 , Annamaria Iorio 17 , Aaron Isaacs 18,19 , Ruifang Li-Gao 20 , Honghuang Lin 21 , Ching-Ti Liu 22 , Leo-Pekka Lyytikäinen 23,24 , Jonathan Marten 25 , Hao Mei 26 , Martina Müller-Nurasyid 27,28,29 , Michele Orini 30,31 , Sandosh Padmanabhan 32 , Farid Radmanesh 33,34 , Julia Ramirez 6,7 , Antonietta Robino 11 , Molly Schwartz 13 , Jessica van Setten 35 , Albert V. Smith 36,37 , Niek Verweij 34,38,39 , Helen R. Warren 6,7 , Stefan Weiss 40,41 , Alvaro Alonso 42 , David O. Arnar 5,43 , Michiel L. Bots 44 , Rudolf A. de Boer 38 , Anna F. Dominiczak 45 , Mark Eijgelsheim 46 , Patrick T. Ellinor 47 , Xiuqing Guo 48,49 , Stephan B. Felix 41,50 , Tamara B. Harris 51 , Caroline Hayward 25 , Susan R. Heckbert 52 , Paul L. Huang 47 , J. W. Jukema 53,54,55 , Mika Kähönen 56,57 , Jan A. Kors 58 , Pier D. Lambiase 12,31 , Lenore J. Launer 51 , Man Li 59 , Allan Linneberg 60,61,62 , Christopher P. Nelson 15,16 , Oluf Pedersen 10 , Marco Perez 63 , Annette Peters 29,64,65 , Ozren Polasek 66 , Bruce M. Psaty 67,68 , Olli T. Raitakari 69,70 , Kenneth M. Rice 71 , Jerome I. Rotter 72 , Moritz F. Sinner 28,29 , Elsayed Z. Soliman 73 , Tim D. Spector 74 , Konstantin Strauch 27,75 , Unnur Thorsteinsdottir 5,76 , Andrew Tinker 6,7 , Stella Trompet 53,77 , André Uitterlinden 78 , Ilonca Vaartjes 44 , Peter van der Meer 38 , Uwe Völker 40,41 , Henry Völzke 41,79 , Melanie Waldenberger 29,64,80 , James G. Wilson 81 , Zhijun Xie 82 , Folkert W. Asselbergs 35,83,84,85 , Marcus Dörr 41,50 , Cornelia M. van Duijn 19 , Paolo Gasparini 86,87 , Daniel F. Gudbjartsson 5,88 , Vilmundur Gudnason 36,37 , Torben Hansen 10 , Stefan Kääb 28,29 , Jørgen K. Kanters 89 , Charles Kooperberg 14 , Terho Lehtimäki 23,24 , Henry J. Lin 48,90 , Steven A. Lubitz 49 , Dennis O. Mook-Kanamori 20,91 , Francesco J. Conti 92 , Christopher H. Newton-Cheh 34,93 , Jonathan Rosand 33,34 , Igor Rudan 94 , Nilesh J. Samani 15,16 , Gianfranco Sinagra 17 , Blair H. Smith 95 , Hilma Holm 5 , Bruno H. Stricker 96 , Sheila Ulivi 11 , Nona Sotoodehnia 97 , Suneel S. Apte 3 , Pim van der Harst 38,83,98 , Kari Stefansson 5,76 , Patricia B. Munroe 6,7 , Dan E. Arking 99 , Cecilia W. Lo 13 and Yalda Jamshidi 1,100* Abstract Background: Genome-wide association studies conducted on QRS duration, an electrocardiographic measurement associated with heart failure and sudden cardiac death, have led to novel biological insights into cardiac function. However, the variants identified fall predominantly in non-coding regions and their underlying mechanisms remain unclear. (Continued on next page) * Correspondence: [email protected] † Bram P. Prins and Timothy J. Mead contributed equally to this work. 1 Genetics Research Centre, Molecular and Clinical Sciences Institute, St George’s University of London, London SW17 0RE, UK 100 Genetics Research Centre, Molecular and Clinical Sciences Institute, St George’s University of London, London, UK Full list of author information is available at the end of the article © The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Prins et al. Genome Biology (2018) 19:87 https://doi.org/10.1186/s13059-018-1457-6 (Continued from previous page) Results: Here, we identify putative functional coding variation associated with changes in the QRS interval duration by combining Illumina HumanExome BeadChip genotype data from 77,898 participants of European ancestry and 7695 of African descent in our discovery cohort, followed by replication in 111,874 individuals of European ancestry from the UK Biobank and deCODE cohorts. We identify ten novel loci, seven within coding regions, including ADAMTS6, significantly associated with QRS duration in gene-based analyses. ADAMTS6 encodes a secreted metalloprotease of currently unknown function. In vitro validation analysis shows that the QRS-associated variants lead to impaired ADAMTS6 secretion and loss-of function analysis in mice demonstrates a previously unappreciated role for ADAMTS6 in connexin 43 gap junction expression, which is essential for myocardial conduction. Conclusions: Our approach identifies novel coding and non-coding variants underlying ventricular depolarization and provides a possible mechanism for the ADAMTS6-associated conduction changes. Keywords: Exome chip, Conduction, ADAMTS6, Meta-analysis, Background In the heart, the ventricular conduction system propagates the electrical impulses that coordinate ventricular chamber contraction. The QRS interval on an electrocardiogram (ECG) is used clinically to quantify duration of ventricular depolarization in the heart. Prolonged QRS duration is an independent predictor of mortality in both the general population [1–4] and in patients with cardiac disease [5–10]. QRS interval duration is a quantitative trait influenced by multiple genetic and environmental factors and is known to be influenced by both age and gender [11,12]. The heritability of QRS duration is estimated to be 35–55% from twin and family studies [13–16]. We previously performed a genome-wide association meta-analysis in 40,407 individuals and identified 22 genetic loci associated with QRS duration [17]. The QRS-associated loci highlighted novel biological processes such as kinase inhibitors, but also pointed to genes with established roles in ventricular conduction such as sodium channels, transcription factors, and calcium-handling proteins. However, the common risk variants identified in genome-wide association studies (GWAS) reside overwhelmingly in regulatory regions, making inference of the underlying causative genes difficult. Furthermore, as with most complex traits, the variants discovered to date explain only a small proportion of the total heritability (the “missing heritability”paradigm), suggesting additional variants are yet to be identified. In fact, the role of rare and low frequency variants, which cannot currently be detected using standard genome-wide single nucleotide polymorphism (SNP) chip arrays, have not been fully investigated. Here we used the Illumina HumanExome BeadChip to focus on rare (MAF < 1%), low frequency (MAF = 1–5%), and common (MAF ≥5%) putative functional coding variation associated with changes in ventricular depolarization. Results and Discussion We combined genotype data from 77,898 participants of European ancestry and 7695 of African descent participating in the Cohorts for Heart and Aging Research in Genomic Epidemiology (CHARGE) Exome-Chip EKG consortium (Additional file 1:TableS1).Atotalof 228,164 polymorphic markers on the exome-chip array passed quality control and were used as a basis for our analyses. Through single variant analysis in the combined European and African datasets, we identified 34 variants across 28 loci associated with QRS duration that passed the exome-chip-wide significance threshold (P<6.17× 10 −8 for single variants [Table 1, Additional file 2:Figure S1]). Eight of the identified loci were novel and five of these were driven by low frequency (MAF < 5%) and common (MAF ≥5%) non-synonymous coding variation. We confirmed 20 of the 29 previously identified QRS duration loci [14,17–19], the remaining loci were not covered by the Exome-Chip and/or did not pass quality control (QC) (Additional file 1: Table S2). As might be anticipated when combining two ancestries in association analyses, we detected heterogeneity of effects for one variant (Cochran’s heterogeneity P<1.47×10 −3 , a Bonferroni corrected P value of α=0.05/34 variants), Additional file 1: Table S2). We did not observe evidence for inflation of test statistics for any of the analyses (λ GC = 1.049, European and African ancestries, combined, Additional file 2: Figure S2, individual ancestry results, Additional file 2: Figures S3–S6). We next sought to replicate the 34 lead variants of our 28 loci in a replication meta-analysis of 111,874 individuals from the UK Biobank [20] and deCODE genetics [21] cohorts. In the replication meta-analysis, 30 lead variants for 25 loci replicated (P≤1.47 × 10 −3 = 0.05/34 variants), seven of which were novel, ten of which are known (Additional file 1: Table S2). The remaining four variants that did not replicate in UK Biobank encompass two previously established loci (one in locus SCN5A/SCN10A for which the other five variants replicated) and two novel Prins et al. Genome Biology (2018) 19:87 Page 2 of 17 loci (SENP2,IGF1R). This is likely due to differences in phenotype acquisition methods (UK Biobank having exercise ECGs measured), though effect size directions between discovery and replication remained consistent and Pvalues of non-replicating variants were all below nominal significance (P<0.05). Sex-specific associations with QRS duration Sex differences in QRS duration are well established (men have significantly longer QRS durations than women [22,23]), and might be attributable to differential effects of genetic variation in men and women. Therefore, we performed sex-stratified association analyses (Additional file 1: Table S3, Additional file 2: Figures S7 and S8). We included only those studies that had both male and female participants to mitigate potential bias due to contributions from single-sex cohorts. In total, up to 31,702 men and 39,907 women were included from both European and African ancestry studies. We found suggestive evidence for a sex-specific locus that was not identified in the combined analysis. The non-synonymous variant rs17265513 (p.Asn310Ser) in ZHX3 (zinc fingers and homeoboxes 3) showed a significant association only in men (P male =4.89×10 −8 ,β(SE) =−0.52(0.09)), whereas no effect was observed for women (P female =0.86,β(SE) = −0.01(0.08)); however, there was no Table 1 Lead SNPs for 28 loci identified for QRS duration in a combined European and African American ancestry meta-analysis Locus Band dbSNPID A1/A2 cMAF beta(se) Pn Nearest gene Annotation Novel loci 1 2q31.2 rs17362588 A/G 0.081 0.52 (0.08) 4.20 × 10 −11 85,593 CCDC141 Non-synonymous 2 3p22.2 rs116202356 A/G 0.015 −1.63 (0.17) 1.23 × 10 −20 85,593 DLEC1 Non-synonymous 3 3q27.2 rs6762208 A/C 0.357 −0.31 (0.05) 3.45 × 10 −12 85,593 SENP2 Non-synonymous 4 6q22.32 rs4549631 C/T 0.481 0.28 (0.04) 5.56 × 10 −11 85,593 PRELID1P1 Intergenic 5 8q24.13 rs16898691 G/C 0.040 −0.92 (0.11) 5.71 × 10 −16 79,976 KLHL38 Non-synonymous 6 12q13.3 rs2926743 A/G 0.257 −0.32 (0.05) 9.40 × 10 −11 85,593 NACA Non-synonymous 7 15q26.3 rs4966020 G/A 0.387 −0.27 (0.04) 2.99 × 10 −9 85,593 IGF1R Intronic 8 20p12.3 rs961253 A/C 0.357 0.30 (0.04) 1.20 × 10 −11 85,593 CASC20 Intergenic Previously identified loci 9 1p32.3 rs11588271 A/G 0.333 −0.34 (0.05) 7.59 × 10 −14 85,593 CDKN2C Intergenic 10 1p13.1 rs4074536 C/T 0.305 −0.29 (0.05) 8.27 × 10 −10 85,593 CASQ2 Non-synonymous 11 2p22.2 rs7562790 G/T 0.424 0.37 (0.04) 4.34 × 10 −17 85,593 CRIM1 Intronic 12 2p22.2 rs17020136 C/T 0.185 0.38 (0.07) 1.02 × 10 −8 59,876 HEATR5B Intronic 13 3p22.2 rs6795970 A/G 0.371 0.80 (0.05) 9.19 × 10 −70 85,593 SCN10A Non-synonymous 14 3p21.1 rs4687718 A/G 0.164 −0.36 (0.06) 1.19 × 10 −8 83,134 TKT Intronic 15 5q33.2 rs13165478 A/G 0.377 −0.68 (0.04) 6.74 × 10 −52 85,593 HAND1 Intergenic 16 6p21.2 rs9470361 A/G 0.249 0.84 (0.05) 1.21 × 10 −63 85,593 CDKN1A Intergenic 17 6q22.31 rs11153730 C/T 0.475 0.56 (0.04) 1.99 × 10 −38 85,593 SLC35F1 Intergenic 18 7p14.2 rs1362212 A/G 0.144 0.55 (0.06) 1.22 × 10 −18 85,593 TBX20 Intergenic 19 7p12.3 rs7784776 G/A 0.397 0.27 (0.04) 1.18 × 10 −9 85,593 IGFBP3 Intergenic 20 7q31.2 rs3807989 A/G 0.427 0.40 (0.04) 2.14 × 10 −19 85,593 CAV1 Intronic 21 12q24.21 rs3825214 G/A 0.200 0.46 (0.05) 1.10 × 10 −17 85,593 TBX5 Intronic 22 12q24.21 rs7966651 T/C 0.270 −0.38 (0.05) 6.74 × 10 −15 85,593 TBX3 Intergenic 23 13q22.1 rs1886512 A/T 0.380 −0.36 (0.05) 3.17 × 10 −13 70,887 KLF12 Intronic 24 14q24.2 rs11848785 G/A 0.237 −0.44 (0.05) 5.59 × 10 −18 85,593 SIPA1L1 Intronic 25 17q21.32 rs17608766 C/T 0.127 0.70 (0.07) 9.81 × 10 −27 85,593 GOSR2 UTR3 26 17q24.2 rs9912468 G/C 0.416 0.43 (0.05) 2.34 × 10 −21 79,976 PRKCA Intronic 27 18q12.3 rs663651 G/A 0.446 −0.44 (0.05) 6.59 × 10 −18 61,604 SETBP1 Non-synonymous 28 20q11.22 rs3746435 C/G 0.190 −0.36 (0.06) 2.67 × 10 −10 79,976 MYH7B Non-synonymous Top panel: novel loci; bottom panel: previously identified loci Locus index number for each independent locus, Band cytogenetic band in which the lead SNP for the locus resides, dbSNPID dbSNP rs-number of the lead SNP of the locus, A1/A2 coded/non-coded alleles, cMAF cumulative minor allele frequency, beta(se) effect size (standard error) in ms, PPvalue, ntotal number of individuals analyzed for this variant, Nearest gene (nearest) gene, Annotation variant function (protein coding) Prins et al. Genome Biology (2018) 19:87 Page 3 of 17 significant difference consistent with an interaction with sex (P=2.3×10 −5 ). Additionally, no further evidence was observed in the replication analyses alone (P male =7.95× 10 −4 ,β(SE) = −0.30(0.09), N males = 50,457), (P female = 3.55 × 10 −2 ,β(SE) = −0.17(0.08), N females = 61,417). Association of coding and non-coding variants with QRS duration Among the eight newly identified loci in the sex-combined analysis, five had lead variants that were non-synonymous: CCDC141 (Coiled-Coil Domain Containing 141); KLHL38 (Kelch Like Family Member 38); DLEC1 (Deleted in Lung and Esophageal Cancer 1); NACA (Nascent Polypeptide-Associated Complex Alpha subunit); and SENP2 (SUMO1/Sentrin/SMT3 Specific Protease 2). Suggestive evidence for association of the same non-synonymous variant in CCDC141 (rs17362588; P=4.75×10 −7 ) and an intronic variant in KLHL38 (rs11991744; P=1.25×10 −7 ) with QRS duration was shown in two earlier GWAS [24,25]. DLEC1 has recently been suggested to have a possible role as a tumor suppressor [26], and while specific roles for KLHL38 and CCDC141 (a centrosome associated protein) have not yet been elucidated, they show the highest expression in skeletal and/or cardiac tissue, respectively, among the tissues examined in the Genotype-Tissue Expression (GTEx) Portal database (http://www.gtexportal.org)[27]. Two of the novel loci, NACA and SENP2, have established roles in cardiac development and dysfunction. NACA produces the isoform skNAC (skeletal NACA) and acts as a skeletal muscleand heart-specific transcription factor and is critical for ventricular cardiomyocyte expansion [28]. Cardiac-specific knockdown of skNAC in a Drosophila Hand4.2-Gal4 driver cell-line results in severe cardiac defects [19]. Cardiac-specific overexpression of SENP2,a SUMO-specific protease, leads to congenital heart defects and cardiac dysfunction [29]. In the sex-stratified analysis, the association with ZHX3 (Zinc Fingers and Homeoboxes 3) was also driven by an amino acid changing variant. ZHX3 encodes a transcriptional repressor whose functions are largely unknown. However, the sex-specific association might be explained by hormonal changes that have previously been hypothesized to explain a variety of sex-specific differences observed in ECG measures and conduction disorders [30, 31]. A sex-specific association of ZHX3 has also been previously shown for total cholesterol levels (the effect is only significant in men) [32]. We further identified an intronic variant in the IGF1R (Insulin Like Growth Factor 1 Receptor) locus and two intergenic variants: rs4549631 at locus 6q22.32 and rs961253 at locus 20p12.3. Interestingly, when queried against results from the GTEx project portal [27] for blood and eight tissues (including adipose [subcutaneous], artery [aorta, coronary, tibial], heart [atrium, appendage, left ventricle], lung, muscle [skeletal], nerve [tibial], skin [sun exposed], and thyroid), the lead intronic variant in IGF1R (rs4966020; MAF EA/AA 0.36/0.63) is a left ventricle tissue-specific cis-eQTL (P=2.4×10 −7 ). The variant is also in strong linkage disequilibrium with the strongest cis-eQTL for this tissue (rs4966021, P=5×10 −8 ). IGF1R promotes physiological hypertrophy but protects against cardiac fibrosis [33]; the signaling pathways induced by its binding partner, IGF1, regulate contractility, metabolism, hypertrophy, autophagy, senescence, and apoptosis in the heart [34]. The nearest genes for the two intergenic variants are PRELID1P1 (PRELI Domain Containing 1 Pseudogene 1 [locus 6q22.32]) and CASC20 (Cancer Susceptibility Candidate 20 [non-protein-coding]; locus 20p12.3)—the former a pseudogene and the latter a non-protein-coding gene, both with currently uncharacterized function. Rare ADAMTS6 variants are associated with QRS duration By collapsing rare variants in genes as functional units and jointly testing these for association, substantial statistical power-gains can be achieved [35]. We, therefore, performed gene-based analyses using both the Sequence Kernel Association Test (SKAT) (Additional file 1:Table S4) and burden test (T1) (Additional file 1:TableS5), because these tests have optimal power under different scenarios. Analyses were restricted to variants with MAF < 1% in a total of 16,085 genes. One gene-based significant association (P<5.18×10 −7 ) was identified in ADAMTS6 (A Disintegrin-Like And Metalloproteinase with Thrombospondin Type 1 Motif 6; P SKAT =8.18×10 −8 ,Table2), when including only variants classified as damaging (see “Methods”). Four additional genes showed suggestive evidence of association (P<1×10 −4 )(Table2). The ADAMTS6 gene-based signal is driven by two rare non-synonymous variants: rs61736454 (p.Ser90Leu) and rs114007286 (p.Arg603Trp), which have allele frequencies of 0.0018 and 0.0021, respectively (Additional file 1:Table S6). Notably, a look-up in the independent deCODE QRS duration analysis showed that rs61736454 was highly significant, however not exome-wide ([P=2.65×10 −7 , β(SE) = 3.01(0.58)], MAF = 0.002, N= 59,903), and was extremely well imputed (info score = 0.995). Importantly, after meta-analysis with discovery exome summary statistics, the signal reached exome-wide significance ([P=8.96×10 −13 ,β(SE) = 2.75(0.38)], N= 145,496), underscoring the robustness of our initial discovery signal driver. Data for rs114007286 were not available. ADAMTS6 is a highly constrained gene, with a probability of loss of function intolerance score of 1.0 (pLI = 1.0) (Exome Aggregation Consortium [ExAC], Cambridge, MA, USA; http://exac.broadinstitute.org/). The p.Ser90Leu variant lies within the ADAMTS6 propeptide, which is predicted to be important for Prins et al. Genome Biology (2018) 19:87 Page 4 of 17 initiation of folding, because the homologous ADAMTS9 propeptide is an intramolecular chaperone essential for its secretion [36]. The second variant, p.Arg603Trp, is located in the N-terminal-most TSR domain (TSR1) of ADAMTS6. This domain is the target of protein-O-fucosylation, which is a QC signal that prevents secretion of ADAMTS proteins that are improperly folded [37]. ADAMTS6 is necessary for cardiac development and expression of gap junction protein Cx43 ADAMTS6 belongs to a family of metalloproteases that mediates extracellular proteolytic processing of extracellular matrix (ECM) components and other secreted molecules. ADAMTS6 is closely related to ADAMTS10, which interacts with and accelerates assembly of fibrillin-1, mutations in which cause Marfan syndrome [38]. This suggests that ADAMTS6 could regulate cardiac ECM. While no specific ADAMTS6 substrates have been unequivocally identified, it was reported to regulate focal adhesions, epithelial cell–cell interactions, and microfibril assembly in cultured cells [39]. We show by RNA in situ hybridization that Adamts6 is expressed in the atrioventricular and septal cushions and myocardium of the embryonic heart, with expression persisting into adult ventricular, trabecular, and septal myocardium (Fig. 1a–d). Mice with recessive Adamts6 mutations were recovered in a forward genetic screen [40] (Fig. 1e and f). One mutation (p.Met1Ile) affects the start codon and is predicted null. The second mutation (p.Ser149Arg) lies in the propeptide. Both mutations cause prenatal/neonatal lethality with identical congenital heart defect phenotypes (Additional file 1: Table S7), comprising double outlet right ventricle (Fig. 1j, Additional file 3: Video S1), atrioventricular septal defect (Fig. 1k), and ventricular hypertrophy (Fig. 1j and l). Ventricular conduction relies on cardiomyocyte coupling through gap junctions, with connexin 43 (Cx43) being the predominant myocardial gap junction protein in the human and mouse myocardium. Gja1 (encoding Cx43) knockout mice exhibit slow conduction, QRS prolongation, and increased susceptibility to ventricular arrhythmias [41–43], consistent with its role in mediating electrical coupling required for efficient propagation of ventricular depolarization. While Adamts6 heterozygous (Adamts6 m/+ ) adult mice are viable and without structural heart defects (Additional file 2: Figure S9), their ventricular myocardium shows reduced Cx43 staining (Fig. 2a and b). Western blot shows reduction of Cx43 protein in the adult Adamts6 m/+ myocardium (Fig. 2c and d). Interestingly, parallel quantitative real-time polymerase chain reaction (qRT-PCR) shows unchanged Gja1 messenger RNA (mRNA) expression (Fig. 2e), suggesting post-transcriptional regulation. Analysis of embryonic day 14.5 homozygote Adamts6 m/m mutants shows Cx43 is completely absent in the ventricular myocardium (Fig. 2a and b). Thus, whereas Adamts6 m/m mice have severe structural heart defects and Cx43 deficiency, Adamts6 m/+ hemizygosity leads to reductioninCx43expressionintheventricleswithout defects in cardiac morphogenesis. Together these findings suggest the QRS prolongation in individuals with rare pathogenic ADAMTS6 variants could arise from impaired myocardial connectivity due to Cx43 reduction. Rare ADAMTS6 coding variants lead to impaired ADAMTS6 secretion To determine the functional consequences of the two predicted pathogenic human ADAMTS6 coding variants from the exome-chip analysis (p.Ser90Leu and p.Arg603Trp), myc-tagged ADAMTS6 constructs with the variants introduced by site-directed mutagenesis were expressed in HEK293F cells. Western blotting was used to compare the levels of mutant and wild type (WT) myc-tagged ADAMTS6 in the transfected cell lysates and medium. As positive and negative controls, respectively, we transfected the known pathogenic murine variant (p.Ser149Arg) and two rare non-synonymous human ADAMTS6 variants predicted to be benign (p.Ser210Leu and p.Met752Val). Western blotting confirmed that the Adamts6 p.Ser149Arg variant was not secreted (Fig. 3a). The predicted human pathogenic variants show much reduced secretion comparedtotheWTandbenignvariants(Fig.3b–d). Table 2 Gene-based test association results (for genes with variants classified as damaging) Gene N SNPs cMAF beta(se) T1-Burden P T1-Burden P SKAT Protein function Cardiac-specific involvement ADAMTS6 12 0.0097 −0.72 (0.23) 1.48 × 10 −3 8.18 × 10 −8 Zinc-dependent protease – CSRP3 3 0.0048 1.38 (0.31) 9.65 × 10 −6 9.10 × 10 −6 Regulator of myogenesis Myocyte cytoarchitecture maintenance FHOD3 17 0.0171 0.00 (0.17) 9.86 × 10 −1 1.82 × 10 −5 Actin filament assembly Myofibril development and repair ISM1 5 0.0037 1.47 (0.36) 5.05 × 10 −5 5.88 × 10 −5 Angiogenesis inhibitor – TBX5 8 0.0171 −0.32 (0.17) 5.21 × 10 −2 7.80 × 10 −5 T-box transcription factor Cardiac development and cell cycle control Displayed are the top five genes that have the lowest Pvalues in the SKAT test (for genes with damaging variants) Gene gene in which variants were collapsed, N SNPs number of variants used in the collapsed variant test, cMAF cumulative minor allele frequency of variants in the test, beta(se) T1burden effect size (standard error) in ms, P T1-Burden Pvalue of T1-burden test, P SKAT Pvalue of SKAT test, Protein function function of the protein encoded by respective gene, Cardiac-specific involvement, literature support for physiological involvement of the protein in the heart Prins et al. Genome Biology (2018) 19:87 Page 5 of 17 Significantly, the molecular masses of the secreted p.Ser90Leu and p.Arg603Trp variants observed in cell lysate are comparable to that of the WT protein, indicating normal glycosylation and propeptide excision, which are essential for ADAMTS zymogen conversion to their mature forms [44]. These results suggest that heterozygous individuals have a reduction of secreted ADAMTS6 to 50% of normal, implying reduced proteolytic activity. The resulting disruption of proteolytic remodeling could potentially affect cell– cell and cell–matrix interactions essential for efficient Cx43 gap junction assembly. However, the rs61736454 (p.Ser90Leu) and rs114007286 (p.Arg603Trp) variants were associated with longer and shorter QRS duration, respectively. The reduced secretion observed was more profound for the rs61736454 variant compared to rs114007286, and the assay does not predict what impact a small amount of secreted protein may have, nor how it interacts in the presence of other modifier genes/variants carried by the same individual. Additionally, the two variants might affect overall protein function and interaction with binding partners in different ways. Conclusions In a meta-analysis of data from 77,898 participants of European ancestry and 7695 of African descent in our discovery cohort participating in the Cohorts for Heart and Aging Research in Genomic Epidemiology (CHARGE) Exome-Chip ECG consortium, we identified 28 loci associated with QRS duration. With the addition of 111,874 individuals of European ancestry from the UK Biobank and deCODE cohorts, all 34 variants across the 28 loci passed the exome-chip-wide significance threshold, indicating our results are robust. Furthermore, effect size directions between discovery and replication remained consistent and Pvalues of non-replicating variants in the replication analysis alone were all below nominal significance (P< 0.05). Novel loci include genes involved in cardiac development and dysfunction, some of which are highly expressed in Fig. 1 Adamts6 cardiac expression, sequence conservation, and cardiac anomalies in Adamts6-deficient mice. a–dAdamts6 (red punctate signal)is expressed in the outflow tract (a,blue arrowhead), heart valves (a,yellow arrowhead), atria (a,green arrowhead), and ventricular myocardium (a,orange arrowhead,b-d). e,fDiagram of the two Adamts6 mutant alleles recovered: Met1Ile and Ser149Arg. The sequence alignment shows conservation of the Ser149 residue in ADAMTS6 across species. g–lCongenital heart defects observed in Adamts6 Ser149Arg (Adamts6 m/m ) mutant embryos. A WT mouse heart with normal atrial, ventricular, and outflow tract anatomy (g), an intact atrioventricular septum (d), and normal ventricular myocardium (i). Homozygous Adamts6 Ser149Arg mutants (Adamts6 m/m ) exhibit a spectrum of congenital heart defects, such as a double outlet right ventricle (j, in which the aorta and pulmonary artery both arise from the right ventricle; see Additional file 3: Video S1) or an atrioventricular septal defect (AVSD) (k, in which the atrial and ventricular septa fail to form). Thickening of the ventricular wall is commonly observed, indicating ventricular hypertrophy (l). These mutant hearts (j–l) are shown at embryonic day (E)16.5 but their development is delayed, giving an appearance similar to WT hearts at E14.5 (as shown in (g–i)). Ao aorta, AVSD atrioventricular septal defect, LA left atrium, LV left ventricle, Pa pulmonary artery, RA right atrium, RV right ventricle. Scale bar: (a) 500 μm; (b–d)50μm; (g–l)1mm Prins et al. Genome Biology (2018) 19:87 Page 6 of 17 skeletal and/or cardiac tissue. To establish further evidence for these novel loci and mechanisms underlying each association, future functional experiments are essential. The present study also highlights the efficacy of large-scale population-based exome-chip analysis for discovery of non-synonymous coding variants with significant functional effects. In gene-based tests, we identified an association between ventricular depolarization and rare non-synonymous variants in ADAMTS6, a gene not previously implicated in cardiac conduction. We chose to focus on this novel locus and seek functional validation as the association was driven by multiple rare coding variants that were predicted to be damaging by in silico tools. The coding variants driving the association in the population study and the mutations identified in the mouse forward genetic screen all impair ADAMTS6 secretion, indicating reduction/ loss of function. Significantly, although heterozygosity of the variants in mice is not associated with structural heart defects, we detected reduction of Cx43 gap junctions in the ventricular myocardium. Homozygous Adamts6 mutantsshowcompletelossofCx43 gap junctions as well as structural heart defects, implying a dosage effect. Together, these findings indicate that ADAMTS6 has a novel role in regulating gap junction-mediated ventricular depolarization, with quantitative reduction in ADAMTS6 causing cardiac conduction perturbation. While our study focuses on cardiac conduction, the findings support the potential broad utility of large-scale exome-chip analysis for interrogating coding variants associated with other physiological or clinical parameters. Methods Discovery association analyses Study cohorts All participating studies formed the CHARGE EKG exome-chip consortium, including those belonging to the CHARGE consortium and external studies to investigate the role of functional variation in electrocardiographic traits. Twenty-two cohorts participated in the QRS duration analysis effort representing a maximum total sample size of 85,593 samples, consisting of 77,898 participants of European ancestry (91%) and 7695 of African descent. Individual study details and characteristics are summarized in Additional file 1:TableS1. Fig. 2 Reduction of Cx43 intercalated disk gap junction staining in Adamts6-deficient mice. a,bCx43 staining (green)(a) is reduced throughout ventricular myocardium in embryonic day (E) 14.5 Adamts6 m/m embryos and 6-week and 12-month Adamts6 m/+ mice and quantified in (b). DAPI (blue) was used to visualize cell nuclei. c,dRepresentative western blot (c) and quantification (d) shows reduced Cx43 in three pairs of 6-week Adamts6 m/+ and WT myocardium controls. Gapdh was used as a loading control. eNo change in Gja1 RNA level in 6-week and 12-month Adamts6 m/+ myocardium as compared to control. Scale bar: 50 μm. *P≤0.01. E embryonic, W weeks, M months Prins et al. Genome Biology (2018) 19:87 Page 7 of 17 Phenotype measurements We analyzed QRS duration measured in milliseconds. In each study, individuals were excluded from the analyses if these had a QRS duration of > 120 ms, atrial fibrillation (AF) on baseline electrocardiogram, a history of myocardial infarction or heart failure, had Wolff–Parkinson– White syndrome (WPW), a pacemaker, or used Class I and class III blocking medications (those medications with prefix C01B* according to the Anatomical Therapeutic Chemical (ATC) Classification System, http:// www.whocc.no/atcddd/)[45]. For cohorts that were disease case-control studies, we included only the control subjects in our analyses irrespective of the nature of the case disease. Genotyping and quality control Each participating study performed genotyping using the Illumina HumanExome BeadChip / HumanCoreExome platforms. Owing to the difficulty of accurately detecting and assign genotype calls for rare variants (MAF < 1%), an initial core set of CHARGE cohorts, comprising approximately 62,000 samples, assembled intensity data into a single project for a joint improved calling. The quality of the joint calling was assessed through investigating the concordance of genotypes in samples having both exome-chip and exome-sequence data, described extensively elsewhere [46,47]. Using the curated clustering files from the CHARGE central calling effort, several cohorts within our study re-called their genotypes. The remainder of participating studies used either Gencall [48] or zCall [49], or a combination of both. Full details concerning the genotyping and quality control for each cohort are summarized in Additional file 1: Table S1. Individual studies performed sample-level genotype QC filtering for call rate, removing autosomal heterozygosity outliers, gender mismatches, duplicates as established by identity by descent (IBD) analysis, and removed ethnic outliers as determined by multidimensional scaling. Poorly called variants were typically removed by filtering for Hardy-Weinberg equilibrium test Pvalue (pHWE), call rate, and filtering removing poorly clustering variants. Each study aligned their data reference strand to the Illumina forward strand using a central SNP allele reference and annotation file (SNP info file) [46] for the Illumina Exome Chip. Variants were all mapped to GRCh37/hg19. Only variants present within the SNP info file were initially considered for analyses, 247,871 in total. Next, we filtered out 9252 variants that Fig. 3 A mouse Adamts6 ENU mutant and predicted damaging ADAMTS6 variants have impaired secretion. a,bRepresentative western blots using anti-Myc antibody show a major molecular species of 150 kDa in HEK293F cell lysates, corresponding to the ADAMTS6 zymogen (Z). In contrast, the culture medium of cells transfected with WT ADAMTS6 shows a 130 kDa species, corresponding to mature (M, i.e. furin-processed) ADAMTS6. aThe p.Ser149Arg murine variant is not secreted into the culture medium. bThe predicted damaging human variants, p.Ser90Leu and p.Arg603Trp, have reduced secretion, whereas the predicted benign variants, p.Ser210Leu and p.Met752Val, are secreted normally. Lysate and medium of HEK293F cells transfected with an empty vector (EV) lack immunoreactivity. The membrane was subsequently re-blotted using an anti-GAPDH monoclonal antibody to demonstrate comparable sample loading. c,dDensitometry of ADAMTS6 signal in lysates (c) and medium (d) shows reduced secretion of p.Ser90Leu and p.Arg603Trp variants and normal secretion of p.Ser210Leu and p.Met752Val into the medium, relative to the WT control (*P≤0.01 for n = 3 transfections of each vector) Prins et al. Genome Biology (2018) 19:87 Page 8 of 17 failed QC in the joint calling effort, as well as 6591 variants with inconsistent reference alleles across studies (a total of 11,392 unique SNPs), and considered furthermore only autosomal and chromosome X variants, and only those that were polymorphic in our study, leaving an initial set of 228,164 variants for analysis. For our single variant analyses, we only included variants with MAF > 0.012% (equal to a minor allele count [MAC] of 10), 162,199 in total. Statistical methods All association analyses were carried out using the R-package seqMeta [50]. Each study ran the “prepScores”function and adjusted their analyses for age, gender, body mass index (BMI), height, principal components, and study-specific covariates when appropriate (details in Additional file 1: Table S1). The output of this function is an R “list”object (“a prepScores object”), stored in an .RData file, where each element corresponds to a gene, and contains the scores and MAFs for variants, as well as a matrix of the covariance between the scores at all pairs of SNPs within a gene. All studies performed both gender combined and separated analyses, in addition to separation by ancestry. Using the prepScores objects from each study, we performed meta-analyses using the “singlesnpMeta()”for single variant meta-analyses, and the “burdenMeta”and “skatMeta()”functions of SeqMeta. Coefficients and standard errors from seqMeta can be interpreted as a “one-step” approximation to the maximum likelihood estimates. Ancestry groups were analyzed both separate and combined at the meta-analysis level. For single variant meta-analyses, we included all variants with a MAC ≥10 in order to have well-calibrated type I error rates [51]. Statistical significance was defined using Bonferroni corrections. For single variants, maximally 162,199 variants were included in five separate analyses after filtering for MAC: European and African ancestry separated and combined (n = 3); and sex-stratified analyses (n = 2), resulting in a Bonferroni corrected Pvalue of α=0.05 / 162,199 variants / 5 analyses = 6.17 × 10 −8 . Suggestive sexually dimorphic associations were identified by performing sex-stratified meta-analyses, totaling 39,907 women and 31,702 men, including only from cohorts that had both male and female samples. Variants were deemed to be suggestive sex-specific when reaching below a Pvalue threshold of exome-wide significance (P<6.17×10 −8 ) in one sex and above nominal significance in the other (P> 0.05). For gene-based tests, also performed using seqMeta using the “prepScores”objects from individual cohorts, we assigned variants to genes by annotating all variants on the Exome Chip using ANNOVAR [52] following RefSeq [53] gene definitions mapped to human genome build 37 (hg19). In the collapsed variant tests, we included only variants with MAF < 1% and included only genes for which two or more variants were present (n = 16,085). We performed both SKAT [54] and T1 burden [55] tests, for three different functional sets of variants limited to the following: (I) all variants; (II) missense, nonsense, splice, and indel variants; (III) “damaging”: the same variants as in group II, except for missense only including those that are predicted to be damaging by at least two out of four functional prediction algorithms (Polyphen2 [56], SIFT [57], Mutation Taster [58], and LRT [59]). For the gene-based tests, we used a Bonferroni corrected Pvalue significance threshold of α=0.05 / 16,085 genes / 2 different tests / 3 functional variant classes = 5.18 × 10 −7 . Wedefineaphysicallyindependentlocusasthegenomic region that contains variants within 250 kb on either side of LD-independent lead SNPs (exome-wide significant variants with r2 < 0.1), where LD calculations were based on European ancestry. Following this definition, in certain cases LD-independent lead variants are present in overlapping regions, complicating the definition and reporting of associated genetic loci and harbored genes. Therefore, we annealed loci if LD-independent exome-wide significant variants were < 250 kb from each other. Where lead SNPs from previous analyses were not contained in these regions, we considered these as novel. LD calculations were performedontheIlluminaExomeChipgenotypedatafrom the TwinsUK cohort [60] (n = 1194), using PLINK 1.9 [61]. Replication association analyses Study cohort: UK biobank (UKB) UK Biobank (www.ukbiobank.ac.uk) is a prospective study of 500,000 volunteers, comprising relatively even numbers of men and women aged 40–69 years old at recruitment, with extensive baseline, and follow-up clinical, biochemical, genetic, and outcome measures. Approximately 95,000 individuals were recruited for a Cardio test using a stationary bicycle in conjunction with a four-lead electrocardiograph device at the initial assessment (2006–2008) and ~ 20,000 individuals performed the test again (the first repeat assessment: 2011–2013). The Cardio test, thereafter known as the exercise test, started with 15 s of rest (pre-test), followed by 6 min of exercise (cycling) with an increasing workload, and a 1-min recovery period without exercise. To improve accuracy, we calculated an average QRS waveform by aligning all QRS complexes present in a window of 15 s from the resting stage. Ectopic beats and artifacts were removed. Then, we calculated the correlation between each individual QRS complex and the average QRS waveform and removed those with a correlation coefficient < 0.8. Finally, we repeated the calculation of the average QRS waveform by only considering those highly correlated individual QRS complexes. The QRS width was measured from the average QRS waveform as the interval Prins et al. Genome Biology (2018) 19:87 Page 9 of 17 13. Nolte IM, Jansweijer JA, Riese H, Asselbergs FW, van der Harst P, Spector TD, et al. A comparison of heritability estimates by classical twin modeling and based on genome-wide genetic relatedness for cardiac conduction traits. Twin Res Hum Genet. 2017;20:489–98. 14. Holm H, Gudbjartsson DF, Arnar DO, Thorleifsson G, Thorgeirsson G, Stefansdottir H, et al. Several common variants modulate heart rate, PR interval and QRS duration. Nat Genet. 2010;42:117–22. 15. Li J, Huo Y, Zhang Y, Fang Z, Yang J, Zang T, et al. Familial aggregation and heritability of electrocardiographic intervals and heart rate in a rural Chinese population. Ann Noninvasive Electrocardiol. 2009;14:147–52. 16. Mutikainen S, Ortega-Alonso A, Alén M, Kaprio J, Karjalainen J, Rantanen T, et al. Genetic influences on resting electrocardiographic variables in older women: a twin study. Ann Noninvasive Electrocardiol. 2009;14:57–64. 17. Sotoodehnia N, Isaacs A, de Bakker PIW, Dörr M, Newton-Cheh C, Nolte IM, et al. Common variants in 22 loci are associated with QRS duration and cardiac ventricular conduction. Nat Genet. 2010;42:1068–76. 18. Ritchie MD, Denny JC, Zuvich RL, Crawford DC, Schildcrout JS, Bastarache L, et al. Genomeand phenome-wide analyses of cardiac conduction identifies markers of arrhythmia risk. Circulation. 2013;127:1377–85. 19. van der Harst P, van Setten J, Verweij N, Vogler G, Franke L, Maurano MT, et al. 52 genetic loci influencing myocardial mass. J Am Coll Cardiol. 2016;68:1435–48. 20. Sudlow C, Gallacher J, Allen N, Beral V, Burton P, Danesh J, et al. UK biobank: an open access resource for identifying the causes of a wide range of complex diseases of middle and old age. PLoS Med. 2015;12:e1001779. 21. Gudbjartsson DF, Helgason H, Gudjonsson SA, Zink F, Oddson A, Gylfason A, et al. Large-scale whole-genome sequencing of the Icelandic population. Nat Genet. 2015;47:435–44. 22. Macfarlane PW, McLaughlin SC, Devine B, Yang TF. Effects of age, sex, and race on ECG interval measurements. J Electrocardiol. 1994;27(Suppl):14–9. 23. Okin PM, Roman MJ, Devereux RB, Kligfield P. Gender differences and the electrocardiogram in left ventricular hypertrophy. Hypertension. 1995;25:242–9. 24. den Hoed M, Eijgelsheim M, Esko T, Brundel BJJM, Peal DS, Evans DM, et al. Identification of heart rate-associated loci and their effects on cardiac conduction and rhythm disorders. Nat Genet. 2013;45:621–31. 25. Sano M, Kamitsuji S, Kamatani N, Hong K-W, Han B-G, Kim Y, et al. Genomewide association study of electrocardiographic parameters identifies a new association for PR interval and confirms previously reported associations. Hum Mol Genet. 2014;23:6668–76. 26. Wang Z, Li L, Su X, Gao Z, Srivastava G, Murray PG, et al. Epigenetic silencing of the 3p22 tumor suppressor DLEC1 by promoter CpG methylation in nonHodgkin and Hodgkin lymphomas. J Transl Med. 2012;10:209. 27. Consortium GTE. Human genomics. The genotype-tissue expression (GTEx) pilot analysis: multitissue gene regulation in humans. Science. 2015;348:648–60. 28. Park CY, Pierce SA, von Drehle M, Ivey KN, Morgan JA, Blau HM, et al. skNAC, a Smyd1-interacting transcription factor, is involved in cardiac development and skeletal muscle growth and regeneration. Proc Natl Acad Sci U S A. 2010;107:20750–5. 29. Kim EY, Chen L, Ma Y, Yu W, Chang J, Moskowitz IP, et al. Enhanced desumoylation in murine hearts by overexpressed SENP2 leads to congenital heart defects and cardiac dysfunction. J Mol Cell Cardiol. 2012;52:638–49. 30. James AF, Choisy SCM, Hancox JC. Recent advances in understanding sex differences in cardiac repolarization. Prog Biophys Mol Biol. 2007;94:265–319. 31. Yang P-C, Clancy CE. Gender-based differences in cardiac diseases. J Biomed Res. 2011;25:81–9. 32. Teslovich TM, Musunuru K, Smith AV, Edmondson AC, Stylianou IM, Koseki M, et al. Biological, clinical and population relevance of 95 loci for blood lipids. Nature. 2010;466:707–13. 33. Huynh K, McMullen JR, Julius TL, Tan JW, Love JE, Cemerlang N, et al. Cardiac-specific IGF-1 receptor transgenic expression protects against cardiac fibrosis and diastolic dysfunction in a mouse model of diabetic cardiomyopathy. Diabetes. 2010;59:1512–20. 34. Troncoso R, Ibarra C, Vicencio JM, Jaimovich E, Lavandero S. New insights into IGF-1 signaling in the heart. Trends Endocrinol Metab. 2014;25:128–37. 35. Lee S, Abecasis GR, Boehnke M, Lin X. Rare-variant association analysis: study designs and statistical tests. Am J Hum Genet. 2014;95:5–23. 36. Koo B-H, Longpré J-M, Somerville RPT, Alexander JP, Leduc R, Apte SS. Regulation of ADAMTS9 secretion and enzymatic activity by its propeptide. J Biol Chem. 2007;282:16146–54. 37. Wang LW, Dlugosz M, Somerville RPT, Raed M, Haltiwanger RS, Apte SS. Ofucosylation of thrombospondin type 1 repeats in ADAMTS-like-1/punctin-1 regulates secretion: implications for the ADAMTS superfamily. J Biol Chem. 2007;282:17024–31. 38. Kutz WE, Wang LW, Bader HL, Majors AK, Iwata K, Traboulsi EI, et al. ADAMTS10 protein interacts with fibrillin-1 and promotes its deposition in extracellular matrix of cultured fibroblasts. J Biol Chem. 2011;286:17156–67. 39. Cain SA, Mularczyk EJ, Singh M, Massam-Wu T, Kielty CM. ADAMTS-10 and -6 differentially regulate cell-cell junctions and focal adhesions. Sci Rep. 2016;6:35956. 40. Li Y, Klena NT, Gabriel GC, Liu X, Kim AJ, Lemke K, et al. Global genetic analysis in mice unveils central role for cilia in congenital heart disease. Nature. 2015;521:520–4. 41. Thomas SA, Schuessler RB, Berul CI, Beardslee MA, Beyer EC, Mendelsohn ME, et al. Disparate effects of deficient expression of connexin43 on atrial and ventricular conduction: evidence for chamber-specific molecular determinants of conduction. Circulation. 1998;97:686–91. 42. Gutstein DE, Morley GE, Tamaddon H, Vaidya D, Schneider MD, Chen J, et al. Conduction slowing and sudden arrhythmic death in mice with cardiacrestricted inactivation of connexin43. Circ Res. 2001;88:333–9. 43. Danik SB, Liu F, Zhang J, Suk HJ, Morley GE, Fishman GI, et al. Modulation of cardiac gap junction expression and arrhythmic susceptibility. Circ Res. 2004;95:1035–41. 44. Longpré J-M, McCulloch DR, Koo B-H, Alexander JP, Apte SS, Leduc R. Characterization of proADAMTS5 processing by proprotein convertases. Int J Biochem Cell Biol. 2009;41:1116–26. 45. World Health Organization. WHO | The Anatomical Therapeutic Chemical Classification System with Defined Daily Doses (ATC/DDD). http://www.who. int/classifications/atcddd/en/. Accessed 12 Dec 2017. 46. Grove ML, Yu B, Cochran BJ, Haritunians T, Bis JC, Taylor KD, et al. Best practices and joint calling of the HumanExome BeadChip: the CHARGE consortium. PLoS One. 2013;8:e68095. 47. Wessel J, Chu AY, Willems SM, Wang S, Yaghootkar H, Brody JA, et al. Lowfrequency and rare exome chip variants associate with fasting glucose and type 2 diabetes susceptibility. Nat Commun. 2015;6:5897. 48. Illumina Inc. Illumina GenCall Data Analysis Software. GenCall software algorithms for clustering, calling, and scoring genotypes. San Diego: Technology Spotlight. 2005. http://www.illumina.com/Documents/products/ technotes/technote_gencall_data_analysis_software.pdf. 49. Goldstein JI, Crenshaw A, Carey J, Grant GB, Maguire J, Fromer M, et al. zCall: a rare variant caller for array-based genotyping: genetics and population analysis. Bioinformatics. 2012;28:2543–5. 50. Voorman A, Brody J, Chen H, Lumley T, Davis B. seqMeta: Meta-Analysis of Region-Based Tests of Rare DNA Variants. 2017. https://cran.r-project.org/ web/packages/seqMeta/index.html 51. Ma C. Statistical Methods for Low-frequency and Rare Genetic Variants. 2014. https://deepblue.lib.umich.edu/handle/2027.42/110435. Accessed 12 Dec 2017. 52. Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res. 2010;38:e164. 53. Pruitt KD, Tatusova T, Maglott DR. NCBI Reference sequence (RefSeq): a curated non-redundant sequence database of genomes, transcripts and proteins. Nucleic Acids Res. 2005;33:D501–4. 54. Wu MC, Lee S, Cai T, Li Y, Boehnke M, Lin X. Rare-variant association testing for sequencing data with the sequence kernel association test. Am J Hum Genet. 2011;89:82–93. 55. Li B, Leal SM. Methods for detecting associations with rare variants for common diseases: application to analysis of sequence data. Am J Hum Genet. 2008;83:311–21. 56. Adzhubei IA, Schmidt S, Peshkin L, Ramensky VE, Gerasimova A, Bork P, et al. A method and server for predicting damaging missense mutations. Nat Methods. 2010;7:248–9. 57. Ng PC, Henikoff S. SIFT: predicting amino acid changes that affect protein function. Nucleic Acids Res. 2003;31:3812–4. 58. Schwarz JM, Cooper DN, Schuelke M, Seelow D. MutationTaster2: mutation prediction for the deep-sequencing age. Nat Methods. 2014;11:361–2. 59. Chun S, Fay JC. Identification of deleterious mutations within three human genomes. Genome Res. 2009;19:1553–61. 60. Spector TD, Williams FMK. The UK adult twin registry (TwinsUK). Twin Res Hum Genet. 2006;9:899–906. 61. Chang CC, Chow CC, Tellier LC, Vattikuti S, Purcell SM, Lee JJ. Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience. 2015;4:7. 62. Loh P-R, Tucker G, Bulik-Sullivan BK, Vilhjálmsson BJ, Finucane HK, Salem RM, et al. Efficient Bayesian mixed-model analysis increases association power in large cohorts. Nat Genet. 2015;47:284–90. Prins et al. Genome Biology (2018) 19:87 Page 16 of 17 63. Bulik-Sullivan BK, Loh P-R, Finucane HK, Ripke S, Yang J, Schizophrenia Working Group of the Psychiatric Genomics Consortium, et al. LD score regression distinguishes confounding from polygenicity in genome-wide association studies. Nat Genet. 2015;47:291–5. 64. Mägi R, Morris AP. GWAMA: software for genome-wide association metaanalysis. BMC Bioinformatics. 2010;11:288. 65. Staley JR, Blackshaw J, Kamat MA, Ellis S, Surendran P, Sun BB, et al. PhenoScanner: a database of human genotype-phenotype associations. Bioinformatics. 2016;32:3207–9. 66. Prins BP, Mead TJ, Brody JA, Sveinbjornsson G, Ntalla I, Bihlmeyer NA, et al. Exome-chip meta-analysis identifies novel loci associated with cardiac conduction, including ADAMTS6, Data sets. dbGAP. https://www.ncbi.nlm. nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000287.v6.p1 Prins et al. Genome Biology (2018) 19:87 Page 17 of 17