The electrocardiogram (ECG) is an important tool for diagnosing, monitoring and evaluating risk in patients with cardiovascular disease (CVD; Lin et al., 2013; Pelto et al., 2013). ECG measurements, such as PR interval, QRS complex duration, and QT interval, are used for the diagnosis and prediction of cardiac arrhythmias and sudden cardiac death (SCD; Kolder et al., 2012). Myocardial depolarization and repolarization time are measured by the QT interval: the time between the onset of the QRS complex and the end of the T wave. QT shortening or prolongation has been associated with an increased risk for arrhythmias and SCD (Newton-Cheh and Shah, 2007). PR interval and QRS duration are measures of cardiac conduction time; QRS duration reflects conduction through the ventricular myocardium, while PR interval measures atrial and atrioventricular conduction from the sinoatrial node to the ventricular myocardium, primarily through the atrioventricular node (Cheng et al., 2015; Mozos and Caraba, 2015).

There are significant genetic contributions to ECG measurements; genome-wide association studies (GWAS) identified at least 71 common variants associated with their variability (Arking et al., 2006, 2014; Newton-Cheh et al., 2007, 2009; Pfeufer et al., 2009, 2010; Holm et al., 2010; Sotoodehnia et al., 2010). A number of these associations were established in loci containing genes that encode proteins with previously known roles in heart development and function, such as cardiac transcription factors; sodium, calcium, and potassium ion channels; genes with a role in myocardial electrophysiology; and others involved in the conduction of electrical impulses (Kolder et al., 2012). These include ARHGAP24, SETBP1, LRIG1, CREBBP, MEIS1. TBX20, and TBX5. Some ion channel encoding genes, such as SCN5A, HERG, KCNE1, and KCNE2, have been associated with long QT syndrome (LQTS; Tristani-Firouzi et al., 2001), atrial fibrillation (AF) and Brugada Syndrome (Hedley et al., 2011). Collectively, however, these loci explain only modest proportions of phenotypic variability; GWAS SNPs specific for each trait account for limited trait heritability (17% for QRS, 4% for QT, and 2% for PR) (Silva et al., 2015).

Genome-wide association studies generally interrogate only common variants, typically of small effect. Families, in addition to being robust against population stratification, may be enriched for less frequent variants, which can potentially be identified by linkage and fine mapping. The aim of this study, therefore, was to discover less frequent variants using linkage analysis in a large family-based study, the Erasmus Rucphen Family Study (ERF).

Table 1 shows the characteristics of the participants included in the discovery linkage analyses and exome sequencing, as well as the exome chip replication sample. There were no significant differences between the largely overlapping linkage and exome sequence groups. The replication sample was considerably older, and was characterized by increased frequency of hypertension (and BP differences), increased PR interval and decreased QT interval compared to the discovery samples. The three ECG traits studied (the QT, QRS, and PR intervals) demonstrated only modest pair-wise correlations in the discovery dataset (Supplementary Table 1).

Supplementary Table 2 shows the linkage results for the ECG traits, which yielded a total of five regions with suggestive LOD scores (LOD > 1.9). QT was suggestively linked to two regions, on chromosome 1 (LOD = 2.63) and on chromosome 2 (LOD = 2.05). A suggestive LOD score for QRS was observed on chromosome 1 (LOD = 2.52) and, for PR, two suggestive regions were located on chromosomes 9 and 14 with LOD scores of 2.20 and 2.29, respectively (Supplementary Table 2). Plots of the linked regions are shown in Figure 1.

Our analysis of coding variants in these linked regions revealed 55,050 variants in coding regions of genes under the peaks, as described in Table 2. Of these mutations, 1334 had a frequency less than or equal to 5%, 437 were predicted to be damaging by at least two of the prediction software packages used, and six were nonsense variants. By linkage peak, there were 207 missense damaging mutations and two nonsense mutations on 14 and 113 missense damaging mutations and two nonsense mutations on 2q32 for QT; 51 missense mutations and one nonsense mutation on 1p36 for QRS; and 29 missense mutations on 9q21 and 31 missense mutations and one nonsense mutation on 14q12 for PR. In total, 21 variants had nominal regression P-values less than 0.05 (the smallest P-values under each linkage peak were P = 2.8 × 10^-4 for QT on chromosome 1, P = 2.3 × 10^-2 for QT on chromosome 2, P = 2.6 × 10^-2 for QRS on chromosome 1, P = 1.9 × 10^-2 for PR on chromosome 9, and P = 1.9 × 10^-2 for PR on chromosome 14) without reaching the significance levels needed to account for multiple comparisons (Supplementary Table 3). Looking for known genes under the linkage peaks (Supplementary Table 4), we found two variants previously related to heart failure, TTN (rs72648923; P = 5.5 × 10^-2, MAF = 1.4 × 10^-2) and HSD3B1 (P = 3.9 × 10^-2 MAF = 1.1 × 10^-2). Neither achieved statistical significance after Bonferroni correction, although both genes were marginally associated with QT. Only a single variant, a C > G (Ser > Cys) variant in FCRL2 (rs74608430; P = 2.8 × 10^-4, MAF = 1.9 × 10^-2), approached the Bonferroni threshold for multiple-testing (P = 1.5 × 10^-4). This variant, under the linkage peak on chromosome 1q23.1 for QT, is highly conserved (scorePhastCons = 0.998) and also predicted by PolyPhen-2 to be damaging (0.999). In the whole ERF population, rs74608430 explained 2.42% of the heritability of QT (reducing the LOD to 1.1; h² = 0.87%; P = 0.02). This finding was not replicated in the RS (P = 0.12, β = 0.14). A sequence kernel association test analysis of the gene also failed to achieve significance in the replication sample (P = 0.44).

Not much is known about the function of FCRL2. Among the functions predicted by Gene Network are the regulation of cytosolic Ca²⁺ levels (P = 3.3 × 10^-3) and AMPK stimulated fatty acid oxidation in muscle (P = 4.1 × 10^-3). In the GEO database, FCRL2 expression was higher in AF (Barth et al., 2005; Yeh et al., 2013). Supplementary Figure 1A shows the genes co-expressed with FCRL2, according to Gene Network. Two genes that have been associated with ECG outcomes by GWAS emerge: ARHGAP24, associated with PR, and SETBP1, associated with QRS (Holm et al., 2010; Pfeufer et al., 2010; Sotoodehnia et al., 2010). In the chromosome 1 region linked to QT, looking for co-expression, we found correlations between DMRTA2. CEP350, and MPL with genes previously associated with ECG traits: DMRTA2 is co-expressed with LRIG1, a QRS associated gene (Supplementary Figure 1B); MPL is in a module with MEIS1, associated with PR (Supplementary Figure 1C); and CEP350 interacts with CREBBP, associated with QT (Supplementary Figure 1D). These three genes are not in linkage disequilibrium with each other. At the chromosome 2q34 locus linked with QT, a heart failure gene, TTN, was under the linkage peak. According to Gene Network analysis, expression of TTN is related to expression of three previously known QT genes (ATP1B, TCEA3, and PLN) and two QRS and PR associated genes (TBX20 and TBX5) (Supplementary Figure 1E) (Holm et al., 2010; Pfeufer et al., 2010; Sotoodehnia et al., 2010; Arking et al., 2014). Additionally, SPHKAP, on chromosome 2 under the QT linkage peak, is co-expressed with TBX5 (Supplementary Figure 1F).

Linkage analysis is an important tool for the identification of genomic regions influencing trait variability. The role of TPM1 mutations with sudden death is a clear example of a locus discovered by linkage analysis (Ott et al., 2011; Mango et al., 2016). The advantages of family studies include control of heterogeneity and population stratification (Reich et al., 2005; Ott et al., 2011). We performed a linkage study on ECG measurements and identified five suggestive regions (1p35.1, 1q24.2, 2q34, 9p22.2, 14q11.2). Rare variant analysis in these regions uncovered two genes related to heart failure, TTN (P = 5.5 × 10^-2) and HSD3B1 (P = 3.9 × 10^-2) and one gene with unknown cardiac function FCRL2 (P = 2.8 × 10^-4). None of them reaches statistical significance level after correction for multiple comparisons.

This study was conducted in a large, well-characterized family-based cohort, ascertained on the basis of genealogy and not phenotype. Multiple levels of genetic data, including a linkage panel and exome sequence data, provided a powerful dataset for identifying variants that may not be easily discovered with GWAS. Unfortunately, exome data was not available in the whole cohort, which could limit our ability to identify causal variants. Additionally, the sequence data did not include extra-genic or intronic variants that may be responsible for the observed linkage peaks.

Our analysis of rare coding variants in these linkage regions revealed 55,050 variants in coding regions. One thousand three hundred and thirty-four of these mutations had a frequency less than or equal to 5% and 437 were predicted to be damaging; none reached the significance threshold accounting for multiple comparisons. These variants spanned genes, including TTN and HSD3B1, which have been previously related to CVDs. HSD3B, a gene on chromosome 1 (1p13.1), has two isoforms (HSD3B1 and HSD3B2) that were found to be associated with an increase in plasma aldosterone (Shimodaira et al., 2010). Changes in circulating aldosterone levels can modulate BP and hypertrophy (HT). A genome wide linkage analysis revealed that HSD3B1 is a locus for BP variation (Shimodaira et al., 2010).

Another interesting gene covered by these variants was TTN; this gene encodes a sarcomeric protein named Titin, with a crucial role in sarcomeric structural integrity and muscle elasticity. Mutations in TTN have been shown to cause heart failure in humans. Additionally, mouse models with TTN mutations exhibit weak heart contractility and heart failure (Gerull et al., 2002; Xu et al., 2002; May et al., 2004) and hearts of mutant embryos displayed weak spontaneous contraction (May et al., 2004). Additionally, the TTN network includes three QT associated genes, ATP1B, TCEA3, and PLN. TBX320, a QRS associated gene; and TBX5 (a QRS and QT associated gene).

We also identified a less frequent C > G missense variant (rs74608430) in the FCRL2 gene under the linkage peak on chromosome 1p23.1. This variant explains 2.42% (h² = 0.87%, P = 0.02) of the total genetic variance of QT (h² = 36%) in the ERF population. FCRL2 has not been previously described with respect to cardiac function. Bioinformatics resources, however, showed that FCRL2 expression is associated with ARHGAP24 and SETBP1 expression, two genes implicated in ECG variability by GWAS. This suggests that FCRL2 may be relevant for heart function. FCRL2 is expressed mostly in liver, heart, testis and kidney⁸. Gene Network predicts that it may be relevant for cytosolic Ca²⁺ levels and AMPK stimulated fatty acid oxidation in muscle. These are plausible pathways for QT function. This finding for rs74608430, however, was not replicated in the RS, in which the MAF was 2.9 × 10^-2. The absence of replication could be related to environmental differences influencing complex gene-environment interactions between these two study groups (Teo et al., 2009). Another plausible explanation is that, due to longer stretches of linkage disequilibrium in the family-based ERF sample, rs74608430 is tagging another variant in ERF and this is not the case in the general population.

Further, Ingenuity analysis revealed that FCRL2 is correlated with some microRNAs (such as miR-1263, miR337-5p, miR-4699-3p, miR518e-3p, miR-507, miR3689a-5p, miR-507, miR-3622a-5p, miR-450b-5p, miR-4720-3p, and miR-1253). Among these, miR-337-5p is known to be differentially expressed in patients with valvular heart disease and patients with chronic AF (Cooley et al., 2012). This is consistent with the GEO database at NCBI⁹, which suggests that FCRL2 is upregulated in patients with AF and dilated cardiomyopathy. In summary, the bioinformatics data available for this gene supports the hypothesis that FCRL2 may be involved in heart function, and, specifically, related to ECG variability.

Additional interesting genes have been uncovered under the linkage peaks. First, the PR linkage peak on chromosome 14 contains damaging variants in the alpha and beta subunits of cardiac myosin MYH6 and MYH7. Previous studies showed that genetic variants in these two genes have been found in hypertrophic cardiomyopathy (Anan et al., 1994; Lankford et al., 1995; Rayment et al., 1995; Niimura et al., 2002; Carniel et al., 2005), dilated cardiomyopathy (Kamisago et al., 2000; Carniel et al., 2005) and atrial septal defect (Ching et al., 2005). Second, we found TNNT2 under the linkage peak on chromosome 1 for QT, which harbors known mutations underlying hypertrophic cardiomyopathy (Thierfelder et al., 1994) and familial dilated cardiomyopathy (Kamisago et al., 2000).

No explanatory variants were found for the other loci, for which there are a number of potential explanations. Linkage peaks are not precise in highlighting the location of the causal variant; even the region of interest cannot be easily pinpointed. Additionally, we did not take into account alternative forms of genetic variation, such as structural and copy number variations (CNVs) or repeats in the linkage regions. Lastly, causal rare variants may be located outside the coding sequence, which we did not include in our sequencing analyses.

Although the combination of linkage and exome sequencing did not lead to the identification of a causal variant, suggestive linkage regions contain a number of plausible candidate genes, including FCRL2. TTN, MYH6, MYH7, TNNT2, and HSD321. Further analysis will need to be performed to demonstrate the involvement of these proteins in ECG measurements. We could not explain these with exonic sequence variants, so they will require more extensive follow-up, but provide potentially important indicators of the location of variation influencing ECG.

CS: Formal analysis, writing – original draft preparation; IZ: Formal analysis; NA: Formal analysis; AD: Formal analysis; EvL: Formal analysis; JK: Formal analysis, investigation, software; MvB: Formal analysis; BS: Investigation, resources; AU: Investigation, resources; AK: Formal analysis, software; JW: Investigation, resources; RW: Writing – original draft preparation, supervision; BO: Investigation, resources; TA: Formal analysis, supervision; CvD: Conceptualization, formal analysis, investigation, resources, writing – original draft preparation, supervision; AI: Conceptualization, formal analysis, writing – original draft preparation, supervision.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.