Translating GWAS-identified loci for cardiac rhythm and rate using an in vivo image- and CRISPR/Cas9-based approach

A meta-analysis of genome-wide association studies (GWAS) identified eight loci that are associated with heart rate variability (HRV), but candidate genes in these loci remain uncharacterized. We developed an image- and CRISPR/Cas9-based pipeline to systematically characterize candidate genes for HRV in live zebrafish embryos. Nine zebrafish orthologues of six human candidate genes were targeted simultaneously in eggs from fish that transgenically express GFP on smooth muscle cells (Tg[acta2:GFP]), to visualize the beating heart. An automated analysis of repeated 30 s recordings of beating atria in 381 live, intact zebrafish embryos at 2 and 5 days post-fertilization highlighted genes that influence HRV (hcn4 and si:dkey-65j6.2 [KIAA1755]); heart rate (rgs6 and hcn4); and the risk of sinoatrial pauses and arrests (hcn4). Exposure to 10 or 25 µM ivabradine—an open channel blocker of HCNs—for 24 h resulted in a dose-dependent higher HRV and lower heart rate at 5 days post-fertilization. Hence, our screen confirmed the role of established genes for heart rate and rhythm (RGS6 and HCN4); showed that ivabradine reduces heart rate and increases HRV in zebrafish embryos, as it does in humans; and highlighted a novel gene that plays a role in HRV (KIAA1755).

). Two embryos with missed calls in more than two targeted sites were excluded from the analysis. In the remaining 381 embryos, sequencing call rate was > 98% for all but one targeted site, and missed calls were imputed to the mean. Seventy-eight embryos had a missed call for neo1b (Supplementary Table S4). Since these 78 embryos showed a similar distribution for all outcomes compared with the remaining embryos (not shown), and since the mutant allele frequency of neo1b was 93.4% for the successfully sequenced embryos, missing neo1b calls were also imputed to the mean. We observed a normal distribution of the number of mutated alleles across the nine targeted orthologues (range 2-14 mutated alleles, Supplementary Fig. S6). Mutant allele frequencies-based on previously unreported variants located within ± 30 bp of the CRISPR/Cas9 cut sites-ranged from 4.1% for hcn4l to 93.4% for neo1b (Supplementary Table S4). We did not identify any mutations at three predicted off-target sites in the 381 successfully sequenced embryos (Supplementary Table S2).

Effects of mutations in candidate genes on sinoatrial pauses and arrests.
At 2dpf, each additional CRISPR/Cas9 mutated allele in the first exon of hcn4 resulted in a > 2.5-fold higher odds of sinoatrial pauses or arrests ( Table 2, Supplementary Tables S5, S6). Embryos showing pauses or arrests during positioning and orienting in the microscope's field of view, but not during image acquisition were not included in the statistical analysis. Still, it is worth noting that pauses were observed during positioning in four of the nine embryos that later turned out to be compound heterozygous for CRISPR/Cas9-induced nonsense mutations in hcn4; and in four of the eight embryos carrying a mutated allele in hcn4 as well as in hcn4l. Importantly, zebrafish embryos can survive without a functionally beating heart at this stage of development, thanks to adequate tissue oxygenation by diffusion 23 . This allowed us to observe genetically driven sinoatrial arrests that would have been lethal in embryos of most other species.
Of the 39 embryos with a sinoatrial pause during the acquisition at 2dpf; two embryos died before imaging at 5dpf, and only one also showed a pause or arrest at 5dpf. Since only six new pauses were observed at 5dpf, the statistical power to detect genetic effects on sinoatrial pauses was low at 5dpf. Earlier we performed a pilot experiment in 406 2dpf offspring of an in-cross of heterozygous carriers of the sa11188 mutation 24 , a nonsense mutation in exon 2 of hcn4. Of these 406 embryos, 95, 206 and 105 carried nonsense mutations in 0, 1 and 2 hcn4 alleles at 2dpf, respectively. Of the 406 embryos, 349 passed quality control at 5dpf, or which 81, 177 and 91 embryos carried 0, 1 and 2 mutated alleles, respectively. The statistical power to detect effects of nonsense mutations in hcn4 on sinoatrial pauses was thus substantially higher in the pilot study. Combining data from the CRISPR/Cas9 and sa11188 experiments provided us with 51 and 43 sinoatrial pauses at 2 and 5dpf, and confirmed an > twofold higher odds of sinoatrial pauses for each additional mutated hcn4 allele at 2 and at 5dpf ( Table 2). Embryos carrying nonsense mutations in both hcn4 alleles even had ninefold higher odds of sinoatrial pauses at 5dpf than embryos free from CRISPR/Cas9-induced and sa11188 mutations in hcn4 (P = 1.4 × 10 -5 , Table 2). Genotype distributions did not deviate from Hardy Weinberg equilibrium in either the CRISPR/Cas9 or sa11188 experiment, confirming that nonsense mutations in hcn4 were not lethal between 2 and 5dpf.
Of the remaining candidate genes examined in the multiplexed CRISPR/Cas9 experiment, only mutations in syt10 showed a trend for higher odds of sinoatrial pauses at 2dpf (Supplementary Table S6). This is likely a transient effect, since the six embryos with nonsense mutations in both syt10 alleles that showed a sinoatrial pause at 2dpf did not die between 2 and 5dpf and did not show a sinoatrial pause at 5dpf. In fact, each additional mutated allele in syt10 tended to protect from sinoatrial pauses at 5dpf; none of the 41 embryos with nonsense mutations in both syt10 alleles showed a sinoatrial pause at 5dpf (Supplementary Table S5). Effects of mutations in candidate genes on heart rate variability and heart rate. We next examined the effect of CRISPR/Cas9-induced mutations on HRV and heart rate in embryos free from sinoatrial pauses and arrests. Embryos with CRISPR/Cas9-induced mutations in hcn4 tended to have a higher HRV and a www.nature.com/scientificreports/ higher heart rate at 2dpf (Fig. 3, Supplementary Fig. S3, Supplementary Tables S7, S8). From 2 to 5dpf, embryos with CRISPR/Cas9-induced mutations in hcn4 on average had a larger decrease in HRV and a larger increase in heart rate, resulting in a lower HRV and a higher heart rate at 5dpf (Figs. 3, 4, 5, Supplementary Fig. S3, Supplementary Tables S7, S8). At 5dpf, standardized effect sizes of mutations in hcn4 on HRV and heart rate were of similar magnitude but opposite direction (Fig. 5, Supplementary Tables S7, S8). For heart rate, effect estimates were directionally consistent across the CRISPR/Cas9 and sa11188 experiments, but were more conservative in the latter, both at 2 and at 5dpf (Fig. 5). A lower frame rate in the sa11188 experiment (i.e. 20 frames/s) meant we could not examine the effect of sa11188 mutations in hcn4 on HRV in the pilot study. Across the remaining candidate genes, the 17 embryos with one CRISPR/Cas9 mutated rgs6 allele on average had a nearly 0.5 SD units lower heart rate. The effect size for HRV was ~ twofold lower and did not reach significance (Fig. 3, Supplementary Fig. S3, Supplementary Tables S4, S7). Furthermore, embryos with mutations in si:dkey-65j6.2 (i.e. KIAA1755) tended to have: (1) a higher HRV at 2dpf; (2) a larger decrease in HRV from 2 to 5dpf; and (3) a higher HRV at 5dpf (Fig. 3, Supplementary Fig. S3, Supplementary Table S7). While we did not observe an effect of mutations in si:dkey-65j6.2 or quo on heart rate at 2 or 5dpf, embryos with mutations in si:dkey-65j6.2 tended to have a smaller increase in heart rate from 2 to 5dpf. An opposite trend was observed in  Table S7). Similarly, embryos with nonsense mutations in both neo1a alleles showed a larger increase in heart rate from 2 to 5dpf than embryos free from CRISPR/ Cas9-induced mutations, without showing a difference in heart rate at 2 or 5dpf (Fig. 3, Supplementary Table S8). This suggests that 5dpf may have been too early to detect true genetic effects of mutations in KIAA1755 and/or NEO1 orthologues on heart rate.

Effect of mutations in candidate genes on body size.
Besides having a lower heart rate at 2dpf, embryos with CRISPR/Cas9-induced mutations in rgs6 were shorter at 2 and 5 dpf, and had a smaller dorsal body surface area normalized for body length at 5dpf (Supplementary Figs. S4, S7, Supplementary Table S9). Embryos with mutations in si:dkey-65j6.2 tended to be leaner at both 2 and 5dpf, while embryos with mutations in quo were larger at both time points (Supplementary Figs. S4, S7, Supplementary Tables S9, S10). Embryos with mutations in neo1a tended to be larger at 2dpf but not at 5dpf (Supplementary Figs. S4, S7, Supplementary  Tables S9, S10).

Transcriptomic analyses.
Mutagenesis by targeting a proximal site in the protein-coding sequence using CRISPR/Cas9 has been shown to trigger a compensatory upregulation of the expression of transcripts with sequence similarity 25 . To explore if such compensation may have influenced our results, we next targeted hcn4 and/or hcn4l using CRISPR/Cas9 at the single cell stage, pooled five embryos per sample at 5dpf, and performed a qRT-PCR analysis to examine compensatory effects on transcripts with at least 75% sequence similarity to the main zebrafish hcn4 transcript (Supplementary Tables S11, S12). Compared with control-injected embryos, targeting hcn4 or hcn4 and hcn4l simultaneously-using the same gRNAs and protocols used to generate CRISPR/ Cas9 founders-resulted in a lower expression of hcn4 at the CRISPR/Cas9-targeted site at 5dpf, confirming ontarget activity for the hcn4 gRNA (Fig. 6). We observed no effect of targeting hcn4 on expression of the last exon of hcn4, which may reflect the low proportion of embryos with a nonsense mutation in both hcn4 alleles in the phenotypic screen (i.e. 2.4%, Supplementary Table S4). However, we did observe what may be a compensatory increase in the expression of zebrafish orthologues of HCN1 and HCN2 in samples of hcn4-targeted embryos (Fig. 6). While the effect of targeting hcn4 on the expression of CABZ01086574.1 (likely an orthologue of HCN2) and the trend for an effect on hcn1 expression were driven by one sample of hcn4-targeted embryos, the compensatory effect on hcn2b expression was more robust (Supplementary Table S13). Possible compensatory effects were observed when targeting hcn4, but not when targeting hcn4 and hcn4l simultaneously.
Nanopore off-target sequencing. We next used in vitro Nanopore off-target sequencing (Nano-OTS) 26 in DNA from an adult Tg(acta2:GFP) positive fish used to generate CRISPR/Cas9 founders to explore if offtarget mutagenic activity may have influenced the results for hcn4; and if it may have played a role in the relatively large number of missed sequencing calls for neo1b (Supplementary Table S4). No off-target mutations were Table 2. Effect of mutations in hcn4 on sinoatrial pauses in data from the CRISPR/Cas9 F 1 and sa11188 experiments combined. Associations between the presence of sinoatrial pauses during the 30 s recording and the number of mutated alleles in hcn4 weighted by the predicted effect of those mutations on protein function (additive); or in embryos with nonsense mutations in both hcn4 alleles vs. embryos free from CRISPR/Cas9induced or sa11188 mutations in hcn4 (2 vs. 0) at 2 days post-fertilization (dpf), 5dpf, or at either time point (any). In the CRISPR/Cas9 experiment, associations were examined using logistic regression, adjusting for time of day, batch and for the weighted number of mutated alleles in the other genes as fixed factors; in the sa11188 experiment and combined analysis, associations were examined using mixed models (xtmelogit in Stata), with embryos nested in batches and experiments and adjusting for time of day as a fixed factor. At 2dpf, 26 and 94 unaffected larvae were dropped from the additive and '2 vs. 0' analyses in the CRISPR/Cas9 F1 larvae due to multicollinearity. These observations were included in the combined analysis. www.nature.com/scientificreports/ identified for the hcn4 and neo1b gRNAs. In spite of five mismatches, we did observe one off-target mutation for the hcn4l gRNA, inside the coding region of the glutamine transporter slc38a3a ( Supplementary Fig. S8). This off-target is unlikely to have influenced our results, since we did not observe effects of CRISPR/Cas9-induced mutations in hcn4l, and since the off-target mutagenic activity was even lower than the on-target activity for this gRNA ( Supplementary Fig. S8). Furthermore, in vitro off-target activity does not necessarily imply in vivo mutagenic activity. One off-target-in spite of four mismatches-was observed for the neo1a gRNA. Since this off-target was also characterized by low mutagenic activity and was located > 30 kb away from the nearest gene ( Supplementary  Fig. S8), large indel mutations at the neo1b target site due to neo1a or neo1b gRNA off-target activity are highly unlikely. Hence, it seems safe to conclude that either the 78 missing sequencing calls for neo1b were missing at random-e.g. due to inherently present variants interfering with primer binding-or that they were caused by large indel mutations in neo1b resulting from on-target neo1b gRNA activity that did not subsequently influence outcomes of interest.
Potential druggability. Three of the six human candidate genes taken forward for experimental follow-up showed evidence for a role in cardiac rhythm and/or rate in zebrafish embryos (i.e. HCN4, RGS6, KIAA1755). Only a trend for an effect on sinoatrial pauses was observed for the SYT10 orthologue, albeit with an opposite direction of effect at 2 and 5dpf. We next examined whether these four genes are already targeted by existing, FDA-approved medication using the drug-gene interaction (DGI) database (www.dgidb .org). Only HCN4 is currently targeted 27 , i.e. using ivabradine, an open channel blocker of I f channels 28 (see below). We next identified predicted interaction partners of the proteins encoded by the four genes, using GeneMania 29 and STRING 30 . Some partners are targeted by FDA-approved medication (Supplementary Table S14), amongst which are several Heart rate variability Heart rate Ivabradine. Given the role of mutations in hcn4 on HRV and heart rate, we next examined the effect of 24 h of exposure to 0, 10 or 25 µM ivabradine in DMSO on these outcomes at 5dpf, in embryos free from CRISPR/Cas9-induced mutations. One embryo with a sinoatrial pause was excluded from the analysis, as were   . Association of heart rate variability and heart rate at 2-and 5-days post-fertilization (dpf) and effect of nonsense mutations in both hcn4 alleles. Top: the association of heart rate variability and heart rate at 2 and 5dpf, with embryos carrying nonsense mutations in both hcn4 alleles shown as pink diamonds. Bottom: the mean ± standard deviation of heart rate variability and heart rate at 2dpf and 5dpf for embryos with nonsense mutations in both hcn4 alleles (n = 5) vs. embryos without CRISPR/Cas9-induced mutations in hcn4 (n = 147).
Scientific RepoRtS | (2020) 10:11831 | https://doi.org/10.1038/s41598-020-68567-1 www.nature.com/scientificreports/ 13 embryos with suboptimal image or image quantification quality. In the remaining 118 embryos, ivabradine treatment resulted in a dose dependent higher HRV and lower heart rate, i.e. directionally opposite to the effects of mutations in hcn4 (Fig. 5, Supplementary Fig. S9, Supplementary Table S15). In fact, the effect on heart rate was of similar size-but opposite direction-for 10 µM ivabradine when compared with nonsense mutations in both hcn4 alleles, in data from the CRISPR/Cas9 and sa11188 experiments combined (Fig. 5). Interestingly, the effect of ivabradine on heart rate was 1.5 to 1.6-fold higher than its effect on HRV (Fig. 5, Supplementary  Table S15).

Discussion
Large-scale, in vivo follow-up studies of candidate genes in GWAS-identified loci remain sparse. Here we present an objective, image-based pipeline to systematically characterize candidate genes for cardiac rhythm, rate, and conduction-related disorders. Using a zebrafish model system, we confirmed a role for genes previously implicated in heart rate and rhythm (rgs6 and hcn4); show that effects of ivabradine on heart rate and rhythm are opposite when compared with the early-stage effects of mutations in one of the drug's molecular targets (hcn4); identified a previously unanticipated gene influencing heart rate variability (si:dkey-65j6.2, i.e. KIAA1755); and observed effects of previously unanticipated genes on early growth and development (rgs6, quo, si:dkey-65j6. 2, neo1a). In addition, we confirmed that mutations in hcn4 increase the odds of sinoatrial pauses and arrests 15 . The latter adds weight to the notion that GWAS-identified common variants for complex traits can flag genes for which rare, detrimental mutations cause severe, early-onset disorders 7,31,32 . We show here that an image-and CRISPR/Cas9-based zebrafish model system can be used for systematic characterization of candidate genes in GWAS-identified loci for complex traits. In the near future, integration of evidence across a range of species and approaches should close the existing gap between association and functional understanding. This will no doubt yield new targets that can be translated into efficient new medication for prevention and treatment of complex diseases. The locus harboring HCN4 has been identified in GWAS for HRV 3 , heart rate 4 and atrial fibrillation 33 . HCN4 belongs to the family of I f or "funny" channels, aptly named for being activated upon hyperpolarization and a non-selective permeability for Na + and K + . It is expressed in the sinoatrial node 34,35 and plays an important role in cardiac pace making 36 . The heart rate lowering agent ivabradine 37 is an open channel blocker of I f channels, and has been shown to reduce heart rate and increase HRV in humans 38,39 . It also reduces cardiovascular and all-cause mortality in heart failure patients 40 . In line with findings in humans, one day of treatment with ivabradine resulted in a dose dependent lower heart rate and higher HRV in 5dpf zebrafish embryos. Since the drug had a 1.5 to 1.6-fold larger effect on heart rate than on HRV, its protective effect on cardiovascular mortality in heart failure patients may indeed be driven by its heart rate lowering effect 28 . The larger effect size for heart rate than for HRV also suggests that ivabradine may influence HRV at least in part through its non-vagal effects on heart rate, supporting an intrinsic effect of heart rate on HRV, as has been postulated earlier 41 . Exercise traininginduced bradycardia has previously been shown to result from HCN4 downregulation 42 . Our findings suggest that downregulation of HCN4 may be at least partly responsible for the higher HRV in exercisers, in contrast or addition to a higher vagal tone in exercisers 43 . With only one affected embryo in the ivabradine experiment, we could not examine the role of the drug in prevention of sinoatrial pauses. www.nature.com/scientificreports/ Across two separate experiments, we observed that embryos with nonsense mutations in both hcn4 alleles have ninefold higher odds of sinoatrial pauses at 5dpf, as well as a higher heart rate and lower HRV in embryos free from sinoatrial pauses. Hence, nonselective open I f channel blocking using ivabradine and nonsense mutations in hcn4 have directionally opposite effects on heart rate and HRV in 5dpf zebrafish embryos. Others previously showed that HCN4 +/− humans are typically characterized by bradycardia 44 ; Hcn4 −/− mice die prenatally, without arrhythmias and with lower heart rates when compared with Hcn4 +/− and Hcn4 +/+ mice 36 ; and 2dpf zebrafish embryos with experimentally downregulated hcn4 expression have a higher odds of sinoatrial arrests and lower heart rate compared with un-injected controls 15 . In this light, a higher heart rate in 5dpf zebrafish embryos with nonsense mutations in both hcn4 alleles could be considered unexpected. However, results from our qRT-PCR experiment suggest that this higher heart rate is likely driven by a compensatory upregulation of the expression of HCN1 and HCN2 orthologues. The potential compensatory increase in hcn1 and CABZ01086574.1 expression is driven by only one sample, but since technical triplicates show coherent results, this may reflect the low mutagenic efficiency of the hcn4 gRNA. The likely compensatory increase in hcn2b expression upon CRISPR/cas9 targeting of hcn4 was more robust. Results from a Nanopore-OTS screen showed that off-target activity of the hcn4 gRNA is highly unlikely to explain our results. Thus, it seems likely that a higher heart rate in embryos with nonsense mutations in hcn4-likely facilitated by a higher expression of HCN1 and HCN2 orthologues-serves as an attempt to prevent sinoatrial arrests and sudden cardiac death. This physiological attempt to compensate a genetic defect will no doubt be followed by the bradycardia observed in other species once the heart can no longer cope with the persistently elevated workload. Besides from an absence of the compensatory response when downregulating gene expression using morpholino oligonucleotides 45 , the lower heart rate previously reported in 2dpf zebrafish embryos with downregulated hcn4 15 could also have resulted from RNA toxicity, off-target effects, or a developmental delay caused by the microinjection itself 15,46 .
Regulator of G protein signaling 6 (RGS6) plays a role in the parasympathetic regulation of heart rate 47 and is a negative regulator of muscarinic signaling, thus decreasing HRV to prevent bradycardia. Common variants www.nature.com/scientificreports/ near RGS6 have been identified in GWAS for resting heart rate 6,48,49 , heart rate recovery after exercise 49,50 and HRV 3 . An eQTL analysis using GTEx data showed that the minor T-allele in rs4899412-associated with lower HRV-is also associated with a higher expression of RGS6 in whole blood and tibial nerve. This is directionally consistent with humans with loss-of-function variants in RGS6 showing higher HRV 51 . Also in line with this, Rgs6 −/− mice were previously characterized by lower heart rate, higher HRV and higher susceptibility to bradycardia and atrioventricular block 52 . In our study, zebrafish embryos that were heterozygous for CRISPR/ Cas9-induced mutations in rgs6 on average had a nearly 0.5 SD units lower heart rate at 2dpf. Hence, our results for heart rate at 2dpf are in line with studies in mice. The effect on heart rate was no longer apparent at 5dpf, in spite of a slightly better statistical power. We did not detect an effect of mutations in rgs6 on HRV or risk of sinoatrial pauses or arrests at 2 or 5dpf. Non-synonymous SNPs in KIAA1755 have been identified in GWAS for HRV 3 and heart rate 5 . In our analysis, mutations in si:dkey-65j6.2 tend to result in a higher HRV at 2 and 5dpf, but do not affect heart rate. In humans, the minor C-allele of the HRV-associated rs6123471 in the non-coding 3′ UTR of KIAA1755 tends to be associated with a lower expression of KIAA1755 in the left ventricle and atrial appendage, amongst other tissues (GTEx 53 ); as well as with a higher HRV 3 and a lower heart rate 4 . Hence, a higher HRV at 2 and 5dpf for each additional mutated allele in si:dkey-65j6.2 is directionally consistent with results in humans. Our results suggest that the GWAS-identified association in the KIAA1755 locus-and possibly other loci-with heart rate may have been driven by HRV. This would explain why the eleven HRV-associated loci that showed evidence of an association with heart rate all did so in the expected (i.e. opposite) direction from a phenotypic point of view 3 . KIAA1755 is a previously uncharacterized gene that shows a broad expression pattern, including different brain regions and the left atrial appendage (GTEx 53 ). Future mechanistic studies are required to distill how KIAA1755 influences heart rhythm.
In addition to a comparison of specimens with nonsense mutations in both alleles vs. zero mutated alleles, we examined the role of mutations in each gene using an additive model. The much larger sample size of this alternative approach helped put results from the more traditional approach into perspective. In our pipeline, we decided to raise founders to adulthood, and to phenotypically characterize and sequence F 1 embryos, with stable genotypes. Screening the F 0 or F 2 generation would have each had their own advantages and disadvantages. Wu and colleagues described a powerful approach to efficiently disrupt and analyze F 0 embryos 54 . However, this approach requires drawing conclusions from results in injected, mosaic founders, and is not compatible with multiplexing of multiple gRNAs and target sites. Screening the F 2 generation on the other hand requires hand picking of F 1 fish with suitable mutations, which limits the throughput.
The large differences in mutant allele frequency across the targeted genes can be attributed to several factors. First, mosaic founders were in-crossed six times, using random mating. Hence, different founder fish may have produced the screened F 1 offspring in each crossing. Secondly, while all gRNAs were pre-tested for mutagenic activity, mutated alleles detected in the test injections that didn't affect the germline were not included in our screen. Thirdly, mutations that are embryonic lethal prior to 2dpf did not make it into the screen.
Of the three genes for which we show effects of mutations on cardiac rhythm and/or rate, only HCN4 is currently targeted by FDA-approved drugs 55 . Although the other two genes are not highlighted as being druggable by small molecules, targeting via antisense oligonucleotides, antibodies, or other approaches can be explored. Furthermore, expanding our search revealed several druggable interaction partners of putative causal genes that are already targeted by FDA-approved anti-hypertensive and neuromuscular blocking agents, amongst others. SNAP25, a proposed interaction partner of SYT10, is the target of botulinum toxin A. Injection of botulinum toxin A into the neural ganglia of CAD/atrial fibrillation patients was recently shown to reduce the occurrence of post-operational atrial fibrillation 56 . For existing drugs that target interaction partners of putative causal genes, it is worthwhile examining the effect on cardiac rhythm, rate and development, since repurposing FDA-approved drugs would imply the quickest and safest route to the clinic. Quantifying possibly unknown beneficial or adverse side effects of these drugs related to cardiac rhythm would also be informative.
Four potential limitations of our approach should be discussed. Firstly, acquiring 30 s recordings implies that false negatives for sinoatrial pauses or arrests are inevitable. We decided to exclude the embryos with a sinoatrial pause or arrest during positioning under the microscope from the analysis, because the case status of these embryos cannot be confirmed objectively, but they are not appropriate controls either. This limitation will at most have resulted in conservative effect estimates. Secondly, CRISPR/Cas9 gRNAs with predicted off-target effects free from mismatches were avoided. However, two of the selected targets-i.e. for hcn4 and si:dkey-65j6.2-had predicted off-target activity with three mismatches in galnt10 and dclk1, respectively, at the time we designed the gRNAs (Supplementary Table S2). Human orthologues of these two genes have previously been associated with heart rate variability-related traits 57 (dclk1), as well as with carotid intima-media thickness 58 and body mass index 59-61 (galnt10). However, none of the 381 embryos carried CRISPR/Cas9-induced mutations at the three predicted off-target sites (Supplementary Table S2). Furthermore, in vitro genome-wide Nano-OTS only revealed two off-target sites across the four gRNAs of HCN4 and NEO1 orthologues, both with low mutagenic activity. Taken together, this implies that gRNA off-target mutagenic activity is extremely unlikely to have influenced our results. Thirdly, we in-crossed mosaic founders (F 0 ) and phenotypically screened and sequenced the F 1 generation. For some genes, this yielded a very small number of embryos with 0 or 2 mutated alleles (i.e. for rgs6, hcn4, and hcn4l), resulting in a low statistical power to detect a true role for mutations in these genes. In spite of this limitation, we still detected significant effects of mutations in rgs6 and hcn4, and confirmed our results for mutations in hcn4 on heart rate and risk of sinoatrial pauses in an independent experiment with more statistical power. Unfortunately, the frame rate of image acquisition in the pilot experiment was too low to also replicate effects of mutations in hcn4 on HRV. Finally, we recorded the atrium only, to enable a higher frame rate, a higher resolution in time for HRV quantification, and a higher statistical power to detect small genetic effects on HRV. As a result, any ventricular abnormalities that may have occurred were not registered, Scientific RepoRtS | (2020) 10:11831 | https://doi.org/10.1038/s41598-020-68567-1 www.nature.com/scientificreports/ and uncontrolled atrial contractions may thus reflect atrial fibrillation, premature atrial contractions, high atrial rate, or atrial tachycardia. Strengths of our study include its repeated measures design, which enabled us to capture genetic effects at different stages of early development in zebrafish, as well as genetic effects on phenotypic changes over time. Furthermore, the throughput of the setup allowed us to examine the effect of mutations in multiple genes simultaneously, in a larger than usual sample for in vivo genetic screens. Our results demonstrate that a large sample size is paramount to robustly detect genetic effects on complex traits in zebrafish embryos when screening the F 1 generation, even for nonsense mutations. Identifying CRISPR/Cas9-induced mutations allele-specifically using a custom-written algorithm helped us distinguish between heterozygous and compound heterozygous embryos, which in turn helped pinpoint the effect of mutations in these genes. Also, our study is based on a high-throughput imaging approach with objective, automated quantification, as compared with manual counting and annotation of heart rate in most 4 but not all 12,62 earlier studies. Finally, for all genes that showed an effect on HRV, observed effects were directionally consistent with eQTL associations in humans. This further emphasizes the strength of our model system and the robustness of our findings.
In conclusion, our large-scale imaging approach shows that zebrafish embryos can be used for rapid and comprehensive follow-up of GWAS-prioritized candidate genes for HRV and heart rate. This will likely increase our understanding of the underlying biology of cardiac rhythm and rate, and may yield novel drug targets to prevent cardiac death.

Methods
Candidate gene selection. Candidate genes in GWAS-identified loci for HRV were identified as described in detail in Nolte et al. 3 . Of the 18 identified candidate genes, six were selected for experimental follow-up. This selection was based on overlap with findings from GWAS for heart rate (KIAA1755, SYT10, HCN4, GNG11) 4 , as well as with results from eQTL analyses in sinoatrial node and brain (RGS6). Additional candidate genes from the same or nearby loci were also selected for experimental follow-up, i.e. NEO1, which resides next to HCN4. Zebrafish orthologues of the human genes were identified using Ensembl, as well as using a comprehensive synteny search using Genomicus 63 (Supplementary Table S1). Of the selected genes, GNG11, SYT10 and RGS6 have one orthologue in zebrafish, and HCN4, NEO1 and KIAA1755 each have two orthologues, resulting in a total of nine zebrafish orthologues for six human candidate genes (Table 1). CRISPR/Cas9-based mutagenesis. All nine zebrafish genes were targeted together using a multiplexed CRISPR/Cas9 approach 20 . Briefly, guide-RNAs (gRNAs) were selected using ChopChop 64 and CRISPRscan 65 (Supplementary Table S2), based on their predicted efficiency, a moderate to high GC-content, proximal location in the protein-coding sequence, and absence of predicted off-target effects without mismatches. Oligonucleotides were designed as described 21 , consisting of a T7 or SP6 promoter sequence (for gRNAs starting with 'GG' or 'GA' , respectively), a gene-specific gRNA-target sequence, and an overlap sequence to a generic gRNA. The gene-specific oligonucleotides were annealed to a generic 80 bp long oligonucleotide at 98 °C for 2 min, 50 °C for 10 min, and 72 °C for 10 min. The products were checked for correct length on a 2% agarose gel. The oligonucleotides were subsequently transcribed in vitro using the manufacturer's instructions (TranscriptAid T7 high yield transcription kit/MEGAscript SP6 transcription kit, both ThermoFisher Scientific, Waltham, USA). The gRNAs were purified, after which the integrity of the purified gRNAs was examined on a 2% agarose gel. The zebrafish codon-optimized plasmid pT3TS-nls-zCas9-nls was used as a template to produce Cas9 mRNA 66 . The plasmid was linearized with Xba1, and then purified using the Qiaprep Spin Miniprep kit (Qiagen, Hilden, Germany). The DNA was transcribed using the mMESSAGE mMACHINE T3 Transcription Kit (ThermoFisher Scientific, Waltham, USA), followed by LiCl precipitation. The quality of the RNA was confirmed on a 1% agarose gel.

Husbandry and microinjections. A zebrafish line with GFP-labelled α-smooth muscle cells
Tg(acta2:GFP) 21 was used to visualize the beating heart. To this end, eggs from an in-cross of Tg(acta2:GFP) fish were co-injected with a mix of Cas9 mRNA (final concentration 150 ng/μL) and all nine gRNAs (final concentration 25 ng/μL each) in a total volume of 2nL, at the single-cell stage. CRISPR/Cas9 injected embryos were optically screened for the presence of Tg(acta2:GFP) at 2 days post fertilization (dpf), using an automated fluorescence microscope (EVOS FL Cell imaging system, ThermoFisher Scientific, Waltham, USA). Tg(acta2:GFP) carriers were retained and raised to adulthood in systems with circulating, filtered and temperature controlled water (Aquaneering, Inc, San Diego, CA). All procedures and husbandry were conducted in accordance with Swedish and European regulations, and have been approved by the Uppsala University Ethical Committee for Animal Research (C142/13 and C14/16).
Experimental procedure imaging. The mosaic founders (F 0 generation) were only used to yield offspring (F 1 embryos) for experiments. To reach the experimental sample size, founders were in-crossed six times using random mating. Eggs were collected after founders were allowed to reproduce for 45 min, to minimize variation in developmental stage. Fertilized eggs were placed in an incubator at 28.5 °C. At 1dpf, embryos were dechorionated using pronase (Roche Diagnostics, Mannheim, Germany).
At 2dpf, embryos were removed from the incubator and allowed to adapt to controlled room temperature (21.5 °C) for 20 min. Individual embryos were exposed to 100 μg/mL Tricaine (MS-222, Sigma-Aldrich, Darmstadt, Germany) for 1 min before being aspirated, positioned in the field of view of a fluorescence microscope, and oriented dorsally using a Vertebrate Automated Screening Technology (VAST) BioImager (Union Biometrica Inc., Geel, Belgium). We subsequently acquired twelve whole-body images, one image every 30 degrees of rotation, using the camera of the VAST BioImager, to quantify body length, dorsal and lateral surface area and Scientific RepoRtS | (2020) 10:11831 | https://doi.org/10.1038/s41598-020-68567-1 www.nature.com/scientificreports/ volume, as well as the presence or absence of cardiac edema. The VAST BioImager then positioned and oriented the embryo to visualize the beating atrium and triggered the upright Leica DM6000B fluorescence microscope to start imaging using an HCX APO L 40X/0.80W objective and L5 ET, k filter system (Micromedic AB, Stockholm, Sweden). Images of the beating atrium were acquired for 30 s at a frame rate of 152 frames/s using a DFC365 FX high-speed CCD camera (Micromedic AB, Stockholm, Sweden). After acquisition, the embryos were dispensed into a 96-well plate, rinsed from tricaine, and placed back into the incubator. The procedure was repeated at 5dpf, to allow capturing of genetic effects that influence HRV and heart rate differently at different stages of development 22 . After imaging at 5dpf, the embryos were once again dispensed into 96-well plates, sacrificed, and stored at − 80 °C for further processing.

Quantification of cardiac traits and body size.
A custom-written MATLAB script was used to convert the images acquired by the CCD camera into quantitative traits. To acquire the heart rate, each frame of the sequence was correlated with a template frame. The repeating pattern yields a periodic graph from the correlation values and by detecting the peaks in the graph we can assess the heart rate. The template frame should represent one of the extreme states in the cardiac cycle, i.e. end-systole or end-diastole. To detect these frames, we examined the correlation between the first 100 frames. The combination of frames that showed the lowest correlation corresponded to the heart being in opposite states. One of these frames was chosen as the template 67 . This numeric information was subsequently used to quantify: (1) heart rate as the inverse of RR-interval; (2) the standard deviation of the normal-to-normal RR interval (SDNN); and (3) the root mean square of successive heart beat interval differences (RMSSD). Finally, a graph of pixel changes over time was generated across the 30 s recording to help annotate the script's performance. Inter-beat-intervals were used to objectively quantify sinoatrial pauses (i.e. the atrium stops contracting for longer than 3 × the median inter-beat-interval of the embryo, Fig. 2, Supplementary Recording S1) and sinoatrial arrests (i.e. the atrium stops contracting for longer than 2 s, Fig. 2, Supplementary Recording S1). The graphs of pixel changes over time were also used to identify embryos with other abnormalities in cardiac rhythm. Such abnormalities were annotated as: uncontrolled atrial contractions (Fig. 2, Supplementary Recording S2); abnormal cardiac morphology (i.e. a tube-like atrium, Fig. 2, Supplementary Recording S3); or impaired cardiac contractility (i.e. a vibrating rather than a contracting atrium, Supplementary Recording S3). These phenotypes were annotated independently by two investigators (BvdH and MdH), resulting in an initial concordance rate > 90%. Discrepancies in annotation were discussed and reevaluated to reach consensus. Bright-field images of the embryos were used to assess body length, dorsal and lateral surface area, and body volume. Images were automatically segmented and quantified using a custom-written CellProfiler 68 pipeline, followed by manual annotation for segmentation quality. Embryos with suboptimal segmentation quality due to the presence of an air-bubble on the capillary, a bent body, an incomplete rotation during imaging, partial capturing of the embryo, or an over-estimation of size were replaced by images with a 180° difference in rotation, or excluded from the analysis for that outcome if the second image was also sub-optimally segmented. The embryo was excluded from the analysis for body volume if more than four of the 12 images had a bad segmentation. Imaging, image quantification and image quality control were all performed blinded to the sequencing results.
Quality control of phenotype data. A series of quality control steps was performed to ensure only highquality data was included in the genetic association analysis (Supplementary Fig. S2). First, graphs indicating that one or more true beats were missed by the script were removed from the analysis a priori ( Supplementary  Fig. S2). Second, embryos showing uncontrolled atrial contractions, sinoatrial pauses or arrests, edema, abnormal morphology and/or reduced contractility, and embryos showing a sinoatrial pause or arrest during positioning under the microscope were excluded from the analyses for HRV and heart rate ( Supplementary Fig. S2). The latter were also excluded from the analysis for sinoatrial pauses and arrests, since we cannot ascertain case status for sinoatrial pauses and arrests in the same rigorous manner for such embryos, but they are not appropriate controls either. Third, genetic effects were only examined for cardiac outcomes with at least ten cases, i.e. for sinoatrial pauses and arrests only.
Sample preparation for sequencing. After imaging at 5dpf, embryos were sacrificed and DNA was Scientific RepoRtS | (2020) 10:11831 | https://doi.org/10.1038/s41598-020-68567-1 www.nature.com/scientificreports/ Processing of sequencing data. A custom-written bioinformatics pipeline was developed in collaboration with the National Bioinformatics Infrastructure Sweden, to prepare .fastq files for analysis. First, a customwritten script was used to de-multiplex the .fastq files by gene and well. PEAR 70 was then used to merge pairedend reads, followed by removal of low-quality reads using FastX 71 . The reads were then mapped to the wildtype zebrafish genome (Zv11) using STAR 72 . Next, we converted files from .sam to .bam format using samtools 73 , after which variants-mostly indels and SNVs-were called allele specifically using a custom-written variant calling algorithm in R (Danio rerio Identification of Variants by Haplotype-DIVaH). A summary of all unique sequences identified at each targeted site is shown in Supplementary Fig. S5. All unique variants (Supplementary  Table S3) located within ± 30bps of the CRISPR/Cas9-targeted sites that were identified across the two alleles were subsequently pooled, and used for functional annotation using Ensembl's VEP 74 . Naturally occurring variants based on sequencing of un-injected Tg(acta2:GFP) embryos or Ensembl were excluded. In absence of a continuous score provided by variant effect prediction algorithms, we attributed weights of 0.33, 0.66 and 1 for variants with a predicted low, moderate and high impact on protein function. Pilot experiments suggested that this was a more powerful approach than not weighing. Transcript-specific dosage scores were then calculated by retaining the variant with the highest predicted impact on protein function for each allele, target site, and embryo, followed by summing the scores across the two alleles at each target site and embryo. Since all transcripts within a target site were affected virtually identically, we only used the main transcript of each target site for the genetic association analysis. Most embryos had successfully called sequences in all nine CRISPR/Cas9-targeted sites. Embryos with missing calls at more than two of the nine targeted sites were excluded from the genetic association analysis. For the remaining 381 embryos, missing calls were imputed to the mean dosage of the transcript. For neo1b, calling failed in 78 embryos (Supplementary Table S4). The imputed mean dosage for the main transcript of neo1b was still included in the genetic association analysis, since: (1) the mutant allele frequency in successfully called embryos was very high (i.e. 0.934) and an imputed call thus likely closely resembles the truth; and (2) the distribution of embryos with a missed call was similar for all outcomes when compared with the remaining embryos. Hence, using an imputed dosage was concluded to influence the results less than to either exclude the gene from the analysis while it had been targeted; or to discard the 78 embryos that had been successfully phenotyped and sequenced for the other targeted genes. However, the results for neo1b should still be interpreted in light of its call rate. The mutant allele frequency was low for hcn4, hcn4l and rgs6. No embryos carried CRISPR/Cas9-induced mutations within ± 30 bp of any of the three putative off-targets regions (i.e. dclk1b and both galnt10 orthologues, Supplementary Table S2). qPCR to evaluate targeting hcn4 and hcn4l. We performed a qRT-PCR experiment to examine if targeting hcn4 and/or hcn4l resulted in a compensatory upregulation of the expression of transcripts with > 75% sequence similarity to the main transcript of hcn4. To this end, fertilized eggs of Tg(acta2:GFP) positive fish were either: (1) left un-injected; or injected at the single cell stage with: (2) Cas9 mRNA only; (3) hcn4 and hcn4l gRNA only; (4) hcn4 gRNA and Cas9 mRNA; (5) hcn4l gRNA and Cas9 mRNA; or (6) hcn4 and hcn4l gRNAs and Cas9 mRNA. The same protocols, quantities and gRNAs were used as described above for generating founders (see "Husbandry and microinjections" section). Embryos were raised to 5dpf in an incubator at 28.5 °C. On the morning of day 5, within each of the six conditions, multiple pools of five embryos per sample were flash frozen in liquid nitrogen. The experiment was performed twice, with all conditions generated on both occasions, to generate a total of 69 samples for conditions 1 to 6 described above. Samples were homogenized in Trizol using a 20-gauge needle and a 1 mL syringe. RNA was extracted using the TRIzol™ Plus RNA Purification Kit and Phase-maker™ Tubes Complete System (Cat.No: A33254, Invitrogen, Waltham, MA). The quality and concentration of the extracted RNA was measured with the Agilent RNA 6000 nano kit (Cat-No: 5067-1511, Agilent, Santa Clara, CA) on an Agilent 2100 Bioanalyzer. In vitro transcription of 100 ng of RNA per sample was performed using the SuperScript IN VILO Master Mix (Cat.No: 11755500, Invitrogen), and quantitative PCR was performed in triplicate using the PowerUp SYBR Green Master Mix on either an AriaMx (Agilent) or a StepOnePlus (Applied Biosystems, Waltham, MA) Real-Time PCR System. Pooled cDNA samples were used to optimize the primer concentrations and generate standard curves for the calculation of primer efficiencies. Primers with an efficiency close to 100% and a single peak in melt curve analysis were selected (Supplementary Table S11). We used 1 ng of cDNA as input for the rest of the reactions with the following cycle conditions: 50 °C for 2 min (1 cycle), 95 °C for 2 min (1 cycle), 95 °C for 3 s and 60 °C for 30 s (40 cycles). Differences between the three experimental conditions (conditions 4 to 6) and the two injected controls combined (conditions 2 and 3) on gene expression at the hcn4 and hcn4l CRISPR/Cas9-targeted site and at the last exon, as well as for CABZ01086574.1 (HCN2), hcn2b, hcn3 and hcn1 (Supplementary Table S11) were examined using Pfaffl's method 75 , taking into account primer efficiencies (Supplementary Table S12); expression levels in non-injected controls as calibrator; and expression of mob4 as the reference gene 76 .
Nanopore off-target sequencing (OTS). We next used Nanopore off-target sequencing (Nano-OTS) to explore whether CRISPR/Cas9 off-target activity may have influenced the results for hcn4, and whether off-targets may have caused the relatively large number of missed calls for neo1b. To this end, DNA was extracted using the MagAttract HMW kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions for extraction from tissue. Genomic DNA was fragmented to 20 kb using Megaruptor 2 (Diagenode, Liege, Belgium) and size selected with a 10 kb cut-off using the Blue Pippin system (Sage Science, Beverly, MA). Libraries for Nano-OTS were prepared as described recently by Höijer et al. 26 . In short, ribonucleoproteins (RNPs) were prepared as a pool using the hcn4, hcn4l, neo1a and neo1b gRNAs. Next, 3 µg of fragmented and size-selected DNA was dephosphorylated and digested by Cas9 using the RNPs, and DNA ends were dA-tailed. Finally, sequencing Scientific RepoRtS | (2020) 10:11831 | https://doi.org/10.1038/s41598-020-68567-1 www.nature.com/scientificreports/ adapters from the SQK-LSK109 kit (Oxford Nanopore Technologies, Oxford, UK) were ligated to the Cas9cleaved ends. Sequencing was performed on one R9.4.1 flow cell using the MinION system (Oxford Nanopore Technologies). Guppy v3.4 was used for base calling. Analysis and calling of OTS sites were performed as described by Höijer et al. 26 , using the GRCz11 reference genome.
Examining druggability of putative causal genes. All human orthologues of the zebrafish genes for which we observed (a trend for) an effect were explored in the drug-gene interaction database (DGIdb) v3.0.2 27 . Possible interaction partners of the human proteins and genes (Supplementary Table S14) were distilled from STRING v10.5 30 and GeneMania 29 , respectively. The GeneMania search was limited to physical interactions and pathway data.
Characterizing the effect of treatment with ivabradine. To examine the effect of exposure to ivabradine, the procedures described in "Experimental procedure imaging" section were repeated in embryos free from CRISPR/Cas9-induced mutations. On the morning of day 4, embryos were placed in 0, 10 or 25 µM ivabradine in DMSO, for 24 h (Sigma-Aldrich, St Louis, MO). On the morning of day 5, conditions were blinded and embryos were imaged as described in the "Experimental procedure imaging" section. Embryos from the three blinded conditions were imaged in an alternating manner, to prevent differences in developmental stage from influencing the results. The experiment was performed four times to generate a total of 44 embryos per condition. One embryo with a sinoatrial pause and 13 embryos with suboptimal image or image quantification quality were excluded from the analysis, leaving 118 embryos for the analysis of effects on HRV and heart rate.
Statistical analysis. The standard deviation of NN-intervals (SDNN) and the root mean square of successive differences (RMSSD) were strongly correlated at 2dpf (r 2 = 0.78) and at 5dpf (r 2 = 0.83), so a composite endpoint 'HRV' was calculated as the average of SDNN and RMSSD. In embryos free from sinoatrial pauses and arrests; abnormal cardiac morphology; impaired cardiac contractility; and edema, we inverse-normally transformed HRV and heart rate at 2dpf (n = 279) and 5dpf (n = 293) to a mean of 0 and a standard deviation of 1, so effect sizes and their 95% confidence intervals can be interpreted as z-scores. This approach also ensured a normal distribution of all quantitative outcomes, and allowed a comparison of effect sizes across traits. For body size analyses, we first normalized dorsal surface area, lateral surface area and volume for body length, and subsequently inverse normalized these outcomes as well as body length for the statistical analysis. Effects of CRISPR/Cas9-induced mutations in candidate genes on sinoatrial pauses and arrests at 2dpf, and on sinoatrial pauses at 5dpf were examined using a logistic regression analysis. In embryos free from sinoatrial pauses or arrests during the image acquisition, effects of mutations in candidate genes and of ivabradine on HRV and heart rate were examined using hierarchical linear models at 2 and 5dpf separately (Stata's xtmixed), adjusted for the time of imaging (fixed factor), and with embryos nested in batches (random factor with fixed slope). If a sinoatrial pause or arrest was observed before but not during imaging, embryos were not included in the statistical analysis. In a sensitivity analysis, we mutually adjusted associations for HRV and heart rate for the other outcome (i.e. for heart rate and HRV) by adding the outcome as an independent variable to the model. Effects of CRISPR/Cas9-induced mutations in candidate genes on body size were also examined using hierarchical linear models, as described above.
For dichotomous outcomes and continuous outcomes alike, genetic effects were examined using an additive model, with dosage scores for all nine CRISPR/Cas9 targeted sites as independent exposures in the same model, i.e. mutually adjusted for effects of mutations in the other targeted sites. Additionally, for the six genes where at least five embryos carried CRISPR/Cas9-induced nonsense mutations in both alleles, these embryos were compared with embryos free from CRISPR/Cas9-induced mutations at that site, adjusting for dosage scores in the eight remaining genes. For each outcome, associations were examined for the main transcript of the orthologue.
In the qRT-PCR experiment, the effect of targeting hcn4, hcn4l, or hcn4 & hcn4l on gene expression was examined using multiple linear regression analysis, adjusting for batch.
P-values < 0.05 were considered to reflect statistically significant effects. All statistical analyses were performed using Stata MP version 14.2 (StataCorp, College Station, TX).