In vivo characterization of candidate genes for heart rate variability identifies culprits for sinoatrial pauses and arrests

A meta-analysis of genome-wide association studies (GWAS) recently identified eight loci that are associated with heart rate variability (HRV) in data from 53,174 individuals. However, functional follow-up experiments - aiming to identify and characterize causal genes in these loci - have not yet been performed. We developed an image‐ and CRISPR-Cas9-based pipeline to systematically characterize candidate genes for HRV in live zebrafish embryos and larvae. Nine zebrafish orthologues of six human candidate genes were targeted simultaneously in fertilized eggs from fish that transgenically express GFP on smooth muscle cells (Tg(acta2:GFP)), to visualize the beating heart. An automated analysis of 30s recordings of 384 live zebrafish atria at 2 and 5 days post-fertilization helped identify genes that influence HRV (kiaa1755 and gngt1); heart rate (kiaa1755); sinoatrial pauses and arrests (syt10, hcn4 and kiaa1755); and early cardiac development (gngt1, neo1a). Hence, comprehensively characterizing candidate genes in GWAS-identified loci for HRV in vivo helped us identify previously unanticipated culprits for life-threatening cardiac arrhythmias.


INTRODUCTION
Heart rate variability (HRV) reflects the inter-beat variation of heart rate. HRV is 52 controlled by the sinoatrial node, which receives input from the autonomic nervous 53 system. Autonomic imbalance, which has been associated with work stress and other 54 modifiable or non-modifiable risk factors (Thayer et al. 2010), is reflected in lower 55 HRV. Lower HRV has been associated with higher cardiac morbidity and mortality 56 (de Bruyne et al. 1999), and with a higher risk of all-cause mortality (Dekker et al. 57 1997). HRV can be quantified non-invasively using the RR interval of a 12-lead ECG, 58 making HRV a useful clinical marker for perturbations of the autonomic nervous 59 system. However, the genetic basis of HRV remains largely elusive. 60 Recently, we and others identified the first loci that are robustly associated with 61 HRV, using a meta-analysis of genome-wide association studies (GWAS) with data 62 from 53,174 participants (Nolte et al. 2017). Eleven of the 17 associated single 63 nucleotide polymorphisms (SNPs) overlapped with loci that we previously identified 64 as being associated with resting heart rate (den Hoed et al. 2013). The heart rate 65 associated loci in turn had on aggregate been associated with altered cardiac 66 conduction and risk of sick sinus syndrome (den Hoed et al. 2013). In silico 67 functional annotation of the five loci that are associated with both HRV (Nolte et al. 68 2017) and heart rate (den Hoed et al. 2013) previously resulted in the prioritization of 69 six candidate genes that are anticipated to be causal (Nolte et al. 2017). Functional 70 follow-up experiments -ideally in vivo -are required to conclude if these genes are 71 indeed causal, and examine if they influence HRV, heart rate, and/or cardiac 72 conduction. 73 Mouse models are commonly used in cardiac research, but mice show substantial 74 differences in cardiac rate and electrophysiology compared with humans (Poon and 75 Brand 2013). Inherently, these differences complicate extrapolation of results from 76 mouse models to humans (Poon and Brand 2013). Additionally, rodents are not 77 suitable for high-throughput, in vivo screens of cardiac rhythm and rate. Such screens 78 are essential to systematically characterize positional candidate genes in the large 79 number of loci that have now been identified by GWAS for cardiac rhythm, rate 80 (Eppinga et al. 2016) and conduction (Arking et al. 2014). Hence, novel model 81 systems that facilitate systematic, in vivo characterization of a large number of 82 candidate genes are desirable. 83 In recent years, the zebrafish has become an important model system for genetic 84 and drug screens for human disease (Lieschke and Currie 2007; MacRae and Peterson 85 2015). Amongst the many advantages of zebrafish larvae as a model system are their 86 short generation time, large number of offspring produced regularly, relatively low 87 maintenance costs, rapid early development, and optical transparency. Zebrafish have 88 a fully functional, beating heart at ~24 hours post-fertilization, and the ECG of the 89 two-chambered zebrafish heart is comparable to that of humans (Liu et al. 2016). 90 Fluorescently labeled transgenes facilitate visualization of cell types and tissues of 91 interest, which can now be accomplished in high-throughput thanks to advances in 92 automated positioning of non-embedded, live zebrafish embryos  2013). In addition, the zebrafish has a well-annotated genome, with orthologues of at 94 least 71.4% of human genes (Howe et al. 2013). These genes can be targeted 95 efficiently and in a multiplexed manner (i.e. in multiple genes simultaneously) thanks 96 to recent advances in Clustered, Regulatory Interspaced, Short Palindromic Repeats 97 (CRISPR) and CRISPR-associated systems (Cas) (Varshney et al. 2016). All 98 characteristics combined make zebrafish embryos an attractive model system to 99 systematically characterize candidate genes for cardiac rhythm, rate and conduction. 100 The aim of this study was to objectively characterize the most promising candidate 101 genes in GWAS-identified loci for HRV and heart rate for a role in cardiac rhythm, 102 rate and function using a large-scale, image-based screen in zebrafish embryos and 103 larvae.

106
Experimental pipeline 107 Our experimental pipeline is outlined in Figure 1. Founder mutants (F 0 generation) 108 were generated in which all nine zebrafish orthologues of six human candidate genes 109 were targeted using a multiplexed CRISPR-Cas9 approach (Table 1, Supplemental   110   Table S1) (Varshney et al. 2016). Generating mutants in a background with a 111 transgenically expressed fluorescent label on smooth muscle cells (Tg(acta2:gfp)) 112 allowed us to visualize the beating atrium with minimal background noise. F 0 mutants 113 were in-crossed, and the atria of 384 F 1 embryos were recorded for 30s at 2 days post 114 fertilization (dpf). In addition, we captured bright field images of the embryos to 115 quantify body length, surface area and volume. Of the 384 imaged embryos, 326 were 116 imaged again at 5dpf (Supplemental Fig. S1), to capture genetic effects on cardiac 117 outcomes and body size at two key stages of early cardiac development (Hou et al.  rsID refers to the lead SNP(s) of the GWAS-identified loci (Nolte et. al, 2017). Multiple independent SNPs per locus can be present. For simplicity, quo is referred to as kiaa1755_1, and si:dkey-65j6.2 as kiaa1755_2. 0.9%); cardiac edema (n=14, 4.2%); uncontrolled atrial contractions (n=9, 2.7%); or 128 an abnormal cardiac morphology and impaired cardiac contractility (n=6, 1.8%, 129 Supplemental Movie S4). In the 279 embryos without cardiac abnormalities that 130 survived quality control at 2dpf (Supplemental Fig. S1), HRV was 15.5±10 ms 131 (mean±SD) and heart rate 208±24 bpm. In the 291 larvae included in the analysis at 132 5dpf, HRV was 13.18±17.06 ms and heart rate 265±35 bpm (Table 2).  Table S2). A total of 164 unique mutations were 140 identified across the nine CRISPR-targeted sites, ranging from three unique mutations 141 in hcn4l to 37 in kiaa1755_2 (Supplemental Table S3). Frameshift mutations were 142 most common (48.5%), followed by missense variants (24.5%), in-frame deletions 143 (14.7%), and synonymous variants (5.5%). Eighty-five, sixty-nine and nine variants 144 were predicted to have a high, moderate, and low impact on protein function,  Information based on data from 279 embryos at 2 days post-fertilization (dpf) and 291 larvae at 5 dpf. HRV was calculated by averaging SDNN and RMSSD. HRV, SDNN and RMSSD are shown in ms, whereas HR is shown in beats per min. 1 st and 3 rd refers to the respective quartile. HRV: heart rate variability; SDNN: standard deviation of the normal-to-normal interval; RMSSD: root mean square of successive heart beat interval differences; HR: heart rate.
respectively (Supplemental Table S3). Mutant allele frequencies ranged from 4.1% 146 for hcn4l to 92.0% for neo1b (Supplemental Table S4). Each additional mutated allele in the main transcript of kiaa1755_1 was associated 162 with a higher HRV and a trend for a higher heart rate at 2dpf (Figure 2, Supplemental   163   Tables S6-S7). Mutations in kiaa1755_1 were also associated with a larger body 164 volume at 5dpf (Supplemental Table S8). Embryos with mutations in kiaa1755_2 165 tended to be shorter and smaller at 2dpf (Supplemental Table S8), and had higher 166 odds of sinoatrial pauses and a higher heart rate at 5dpf (Figure 2, Supplemental   167   Tables S5-S7). 168 Each additional mutated allele in hcn4 was associated with higher odds of 169 sinoatrial pauses and arrests at 2dpf (Supplemental Table S5). In larvae free from 170 severe cardiac abnormalities, each additional mutated allele in hcn4 was associated 171 with a trend for higher heart rate at 5dpf (Figure 2, Supplemental Tables S6-S7),   172 without affecting body size (Supplemental Table S8). Mutations in hcn4l were 173 positively associated with lateral body surface area at 5dpf (Supplemental Table S8 Table   181 S5), but not with body size.

182
Each additional mutated allele in syt10 was associated with higher odds of a 183 sinoatrial pause at 2dpf (Supplemental Table S5), but did not affect HRV, heart rate 184 or body size (Supplemental Tables S6-S8) Fig. 2: Effect of mutations in candidate genes on heart rhythm and rate at 2 (n=279) and 5dpf (n=291). Dots and whiskers show the effect size and its 95% confidence interval for each additional mutated allele, weighted by its predicted effect on protein function. Effects were adjusted for the number of mutated alleles in the other genes, as well as for heart rate or heart rate variability and time of day (fixed factor), with larvae nested in batches (random factor). Genes are ordered by chromosomal position in humans.      2.43E-01 Associations of cardiac outcomes with the number of mutated alleles across the nine orthologues, weighted by their predicted effect on protein function. At 2 and 5 days post fertilization (dpf), associations were analyzed using logistic regression for outcomes that were observed in at least five larvae. Associations were adjusted for time of day and for the weighted number of mutated alleles in the other genes as fixed factors. '-' indicates that parameters could not be estimated. For simplicity, kiaa1755_1 refers to quo, and kiaa1755_2 to si:dkey-65j6. 740 -Associations of heart rate variability (HRV) and heart rate (mutually adjusted) with the number of mutated alleles across the nine orthologues, weighted by their predicted effect on protein function. Associations were examined for the main transcript of each targeted zebrafish orthologue using hierarchical linear models (xtmixed in Stata), using inverse-normally transformed outcomes, so effect sizes and SEs can be interpreted as z-score units. Associations were adjusted for time of day and the weighted number of mutated alleles in the other genes as fixed factors, with larvae nested in batches (random factor). For simplicity, kiaa1755_1 refers to quo, and kiaa1755_2 to si:dkey-65j6.2.

5 291
Heart rate Supplemental  Associations of heart rate variability (HRV) and heart rate (mutually adjusted) with the number of mutated alleles across the nine orthologues, weighted by their predicted effect on protein function. Associations were examined for the main transcript of each targeted zebrafish orthologue using hierarchical linear models (xtmixed in Stata), using inverse-normally transformed outcomes, so effect sizes and SEs can be interpreted as z-score units. Larvae with a sinoatrial pause before or after the recording were excluded from the analysis in this sensitivity analysis. Analyses were adjusted for time of day and the weighted number of mutated alleles in the other genes as fixed factors, with larvae nested in batches (random factor). For simplicity, kiaa1755_1 refers to quo, and kiaa1755_2 to si:dkey-65j6.2.

Body volume
Associations were analyzed for the main transcript of each targeted zebrafish orthologue using hierarchical linear models (xtmixed in Stata) using inverse-normally transformed outcomes, so effect sizes and SEs can be interpreted as z-score units. Associations were adjusted for time of day and for the weighted number of mutated alleles in the other genes as fixed factors, with larvae nested in batches (random factor). For simplicity, kiaa1755_1 refers to quo, and kiaa1755_2 to si:dkey-65j6.2.