Figures
Abstract
Heart rate variability is an important risk factor for cardiovascular disease and all-cause mortality. The acetylcholine pathway plays a key role in explaining heart rate variability in humans. We assessed whether 443 genotyped and imputed common genetic variants in eight key genes (CHAT, SLC18A3, SLC5A7, CHRNB4, CHRNA3, CHRNA, CHRM2 and ACHE) of the acetylcholine pathway were associated with variation in an established measure of heart rate variability reflecting parasympathetic control of the heart rhythm, the root mean square of successive differences (RMSSD) of normal RR intervals. The association was studied in a two stage design in individuals of European descent. First, analyses were performed in a discovery sample of four cohorts (n = 3429, discovery stage). Second, findings were replicated in three independent cohorts (n = 3311, replication stage), and finally the two stages were combined in a meta-analysis (n = 6740). RMSSD data were obtained under resting conditions. After correction for multiple testing, none of the SNPs showed an association with RMSSD. In conclusion, no common genetic variants for heart rate variability were identified in the largest and most comprehensive candidate gene study on the acetylcholine pathway to date. Future gene finding efforts for RMSSD may want to focus on hypothesis free approaches such as the genome-wide association study.
Citation: Riese H, Muñoz LM, Hartman CA, Ding X, Su S, Oldehinkel AJ, et al. (2014) Identifying Genetic Variants for Heart Rate Variability in the Acetylcholine Pathway. PLoS ONE 9(11): e112476. https://doi.org/10.1371/journal.pone.0112476
Editor: Marie-Pierre Dubé, Universite de Montreal, Canada
Received: February 13, 2014; Accepted: October 15, 2014; Published: November 10, 2014
Copyright: © 2014 Riese et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: This work was supported by a Marie Curie International Re-Integration Grant from the EU to HS (IRG grant 44949) and a grant from the National Heart Lung and Blood Institute to XW (HL0836530). Participating centers of the TRacking Adolescents’ Individual Lives Survey (TRAILS) include various departments of the University Medical Center and University of Groningen, the Erasmus University Medical Center Rotterdam, the University of Utrecht, the Radboud Medical Center Nijmegen, and the Parnassia Bavo group, all in the Netherlands. TRAILS has been financially supported by various grants from the Netherlands Organization for Scientific Research NWO (Medical Research Council program grant GB-MW 940-38-011; ZonMW Brainpower grant 100-001-004; ZonMw Risk Behavior and Dependence grants 60-60600-97-118; ZonMw Culture and Health grant 261-98-710; Social Sciences Council medium-sized investment grants GB-MaGW 480-01-006 and GB-MaGW 480-07-001; Social Sciences Council project grants GB-MaGW 452-04-314 and GB-MaGW 452-06-004; NWO large-sized investment grant 175.010.2003.005; NWO Longitudinal Survey and Panel Funding 481-08-013), the Dutch Ministry of Justice (WODC), the European Science Foundation (EuroSTRESS project FP-006), Biobanking and Biomolecular Resources Research Infrastructure BBMRI-NL (CP 32), the participating universities, and Accare Center for Child and Adolescent Psychiatry. The authors are grateful to all adolescents, their parents and teachers who participated in this research and to everyone who worked on this project and made it possible. TWINS was supported by the Netherlands Organisation for Health Research and Development (ZonMw 904-57-130) and the UK–Netherlands Partnership Program in Science (BR 56–481 and BR 96–229), which is jointly run and financed by the British Council and the Netherlands Organization for Scientific Research (NWO). The BP Stress Study and Georgia Cardiovascular Twin Study were supported by the National Institute of Health (HL086530 and HL69999, respectively). The genotyping conducted in this study was supported by the National Institute of Health (R21 HL086530) and American Heart Association (SDG0730156N). NESDA was supported by the Geestkracht program of ZonMW [grant 10-000-1002]; matching funds from universities and mental health care institutes involved in NESDA (GGZ Buitenamstel-Geestgronden, Rivierduinen, University Medical Center Groningen, GGZ Lentis, GGZ Friesland, GGZ Drenthe). NTR was supported by the Genetic-epidemiology of anxious depression: analyses of twin-family data in the framework of the NEtherlands Twin-family Study of Anxious Depression (NETSAD; ZonMW 940-37-024) and the Centre for Medical Systems Biology (CMSB) 1, WP Genetic epidemiology & gene discovery for cardiovascular risk factors (Netherlands Genomics Initiative/Netherlands Organisation for Scientific Research NWO). Genotyping for NESDA and NTR was funded by the Genetic Association Information Network (GAIN) of the Foundation for the US National Institutes of Health, and analysis was supported by grants from GAIN and the NIMH (MH081802). Genotype data were obtained from dbGaP (www.ncbi.nlm.nih.gov/dbgap, accession number phs000020.v1.p1). PREVEND genetics is supported by the Dutch Kidney Foundation (Grant E033), the National Institutes of Health (grant LM010098), The Netherlands organisation for health research and development (NWO VENI grant 916.761.70), and the Dutch Inter University Cardiology Institute Netherlands. Part of the statistical analyses were carried out on the Genetic Cluster Computer (http://www.geneticcluster.org) hosted by SURFsara and financially supported by the Netherlands Scientific Organization (NWO 480-05-003 PI: Posthuma) along with a supplement from the Dutch Brain Foundation and the VU University Amsterdam. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Heart rate variability (HRV) is the variation over time of the interval between consecutive heart beats. HRV in the respiratory frequency range is a specific marker of parasympathetic control of the heart rhythm [1], and can be reliably assessed by the root mean square of successive differences (RMSSD) of normal RR intervals [2]. Prior research has produced an extensive list of HRV correlates including a broad range of somatic and mental health problems [3]. High HRV is a sign of good adaptability, often used as a marker of well-functioning cardiac autonomic control mechanisms. Reduced HRV is a predictor of hypertension [4], all-cause mortality, arrhythmic events, and sudden death after acute myocardial infarction as well as in the general population [4]–[7].
The parasympathetic (vagal) branch of the autonomic nervous system has an inhibitory influence on the pace-making activity of the sinoatrial node of the heart. Parasympathetic regulation of the heart is mediated by acetylcholine neurotransmission. Acetylcholine activates mainly two types of receptors, the muscarinic and nicotinic receptors. Muscarinic receptors are found on all effector cells that are stimulated by the postganglionic cholinergic neurons of the parasympathetic nervous system. Nicotinic receptors are found on the postganglionic neurons of the autonomic ganglia. The two types of receptors have different functions, and specific drugs can be used to stimulate or block one or the other type [8].
Although there is consistent evidence for the influence of genetic factors on HRV from twin studies showing heritabilities up to 51% [9]–[11], very few studies have tried to identify the genetic polymorphisms responsible for this heritability [12]. Genes involved in the regulatory pathways of acetylcholine, the neurotransmitter of the parasympathetic nervous system, seem plausible candidates to harbor such polymorphisms. In the present study we examined eight key genes (listed in Table 1) involved in biosynthesis, transport, breakdown, and receptor binding of acetylcholine. With the exception of the choline transporter [13], these genes have not been investigated in genetic studies of HRV before.
Neumann et al. [13] genotyped one single nucleotide polymorphism (SNP, rs3731683) in SLC5A7, and observed that this polymorphism was significantly associated with HRV in 413 white adults. The ACHE, CHAT, and SLC18A13 genes have been investigated as functional candidates for Alzheimer’s disease [14], [15]. The fact that patients with Alzheimer’s disease have reduced HRV [16], [17] and this autonomic dysfunction can be normalized by acetylcholinesterase inhibitor treatment [16] suggests that these three genes are also candidate genes for HRV. SLC18A3 lies within the first intron of CHAT [18].
The genes for the β4, α3, and α5 subunits of nicotinic receptors (CHRNB4, CHRNA3, and CHRNA5) are clustered on chromosome 15q24 and physically linked [19]. Interestingly, the only genome-wide linkage scan conducted so far for HRV [20] yielded a multipoint LOD score of 1.84 on chromosome 15 close to this gene cluster. Among the five subtypes of muscarinic receptors (M1-M5), the mammalian heart predominantly expresses M2 receptors [21]. In CHRM2 knockout mice, bradycardia caused by vagal stimulation was completely abolished [22].
In the current study, we assessed whether common genetic variants in these eight key genes of the acetylcholine pathway were associated with RMSSD in a total of 3429 individuals of European descent from four cohorts. Findings were replicated in 3155 subjects from three further cohorts of European descent. The key objective of this study was to identify genetic variants in the acetylcholine pathway contributing to variation in RMSSD.
Methods
Subjects
For the discovery stage, we included subjects from four cohorts: the Georgia cohort (GA) [23], the TRacking Adolescents’ Individual Lives Survey population cohort (TRAILS-P) [24], the Twin Interdisciplinary Neuroticism Study (TWINS) [25], and the Netherlands Study of Depression and Anxiety (NESDA) cohort [26]. As replication samples The Prevention of Renal and Vascular End-stage Disease study (PREVEND) [27], the TRacking Adolescents’ Individual Lives Survey clinical cohort (TRAILS-CC) [24], and the Netherlands Twin Registry (NTR) [28] were used. A short description of the different cohorts is given below; the general characteristics are listed in Table 2. The study protocol was approved by the Research Ethics Board at The Institutional Review Board at the Medical College of Georgia, USA (GA), Central Committee on Research Involving Human Subjects (CCMO; TRAILS-P, TRAILS-CC), Medical Ethics Committee of the University Medical Center Groningen, Groningen (TWINS, PREVEND), and the Ethical Review Board of the VU University Medical Center (NTR, NESDA). Participants were treated in accordance with the guidelines of the declaration of Helsinki, and all measurements were carried out with their adequate understanding and written consent provided by all subjects and by parents if subjects were <18 years.
The GA cohort included subjects from two studies: the BP Stress Study and the Georgia Cardiovascular Twin Study. Both are ongoing longitudinal studies, which have followed the subjects more than 10 years including roughly equal numbers of Caucasians and African Americans and men and women. The BP Stress Study was established in 1989 with subjects aged 7–16 years at baseline [29], while the Georgia Cardiovascular Twin Study was established in 1996 with subjects aged 7–25 years at baseline [30]. For the twin study, the twin pairs were reared together and zygosity was determined using five standard microsatellite markers. Subjects in both studies were recruited from the southeastern United States and were overtly healthy and free of any acute or chronic illness based on parental report. Study design, selection criteria and the criteria to classify subjects as Caucasians or African Americans for these two studies have been described before [31], [32]. A total of 682 Caucasians were available for this study of which i) 261 subjects (mean±SD age, 23.1±2.88 years; range, 15.6 to 31.8 years) from the BP stress cohort who had RMSSD measured on visit 13 and, ii) 421 twins (mean±SD age, 17.9±3.95 years; range, 11.9 to 32.9 years; 45 singletons and 188 twin pairs) from the Georgia Cardiovascular Twin Study who had RMSSD data available from the third visit.
The TRAILS-P cohort included subjects from the ongoing longitudinal population cohort, which started in 2001. Sample selection procedures and methods of study have been described before [24], [33], [34]. For the current study, data of 1029 10-to-13-year-old Dutch preadolescents (47.6% boys, mean±SD age, 11.5±0.53 years) were used for whom RMSSD data were available from the first measurement wave. A more detailed description of the group that participated in the cardiovascular measurements has been given elsewhere [35].
The TWINS cohort was drawn from the Groningen Twin Register, a 3-wave study including approximately 800 twin pairs from the Northern part of the Netherlands [25]. At the second wave a selection of 125 female twin pairs between 18 and 30 years old were invited to visit a psychophysiological laboratory to perform various rest and mental stress tasks. For the current study, valid RMSSD data of 216 subjects were available [11]. Zygosity of this target group was determined using 10 microsatellite markers.
The NESDA cohort was drawn from a larger 8-year longitudinal study conducted among 652 persons without depression or anxiety disorder and 2329 with a remitted or current diagnosis of depression or an anxiety disorder [26]. For the current study, 1502 NESDA subjects were included from which RMSSD data were collected on the baseline measurement session of the study and who had been genotyped on a genome-wide platform.
The PREVEND cohort was drawn from a larger community-based prospective cohort which was originally initiated to investigate the natural course of increased albuminuria levels and its association to renal and cardiovascular disease [27], [36]. PREVEND subjects completed the baseline survey from 1997 to 1998. At the second screening (2001 to 2003) RMSSD data were obtained from 2793 subjects.
The TRAILS-CC cohort included subjects from the ongoing longitudinal cohort of clinically referred subjects following the same protocol as the TRAILS-P cohort described above. Sample selection procedures and methods of study have been described before [24], [33], [37]. For the current study data of 307 TRAILS-CC 10-to-12-year-old Dutch preadolescents (69.3% boys, mean±SD age, 11.11±0.48 years) were used for whom RMSSD data were obtained at the first measurement wave.
The NTR cohort that was established in the late 1980s, consisted of Dutch adult twins (mean age ±SD, 30.83±14.42 years) that had participated in one of various studies of cardiovascular regulation [38], [39] and had been part of the GAIN genotyping effort [40]. RMSSD data obtained during a resting condition were used.
In all cohorts subjects with heart disease and those taking antidepressant or anticholinergic medications were excluded from the analyses.
Heart rate variability measurement and processing
For the GA cohorts, the continuous RR intervals were recorded with BioZ (CardioDynamics, San Diego, CA) impedance monitors. Four dual sensors connected to the BioZ were placed on the subject’s neck and thorax, which formed four electrocardiogram (ECG) vectors. These ECG vectors can be detected by the BioZ. The BioZ then converted the RR interval into beat-to-beat heart rate with a precision of 2 decimal places, including a record of the real time of each beat. At the 10th, 12th, and 14th minute of the 15-minute resting period, a 30-second continuous period of RR intervals was recorded. Prior to calculation of HRV parameters, the raw RR interval data were firstly processed for artificial readings using the following two criteria: (1) RR intervals were between 300 and 2000 ms; (2) the successive RR interval ratios were between 0.8 and 1.2. Dedicated HRV software (Kubios HRV analysis) [41] was used to compute the RMSSD from the RR intervals. RMSSD was calculated as the average of the three 30-second segments of continuous RR intervals.
In the TRAILS-P and TRAILS-CC cohorts, RMSSD was computed on the ECG-derived RR interval time series. ECG measurements took place in a quiet room at school (TRAILS-P) or at the hospital research center (TRAILS-CC). First, subjects lied down and were encouraged to relax and not to move or speak. Next, the procedure was explained to them and data acquisition did not commence until the signal was considerably stable (i.e., as visually observed through the levelling out of the signals after initial fluctuations), usually within five minutes. RR intervals were registered for four minutes in the supine position. To obtain RR intervals with sufficient resolution for HRV determination, a special interpolation algorithm was used, increasing the time resolution for R-peak detection by a factor of 2.5. After initial visual inspection and exclusion of unusable signals, stationarity of the time series was checked and artifacts were corrected using the CARSPAN program (IECProgamma, Rijks Universiteit Groningen, Groningen, The Netherlands) [42]. A more detailed description of cardiovascular data assessment, analysis, attrition, and internal reliability has been given before [35].
In the TWINS cohort RMSSD was computed on the ECG-derived RR interval time series, which were assessed in a quiet laboratory with a three-lead ECG during four standardized conditions: a rest condition, two stress conditions, and a second rest condition. Cardiovascular measurements started after the subjects had relaxed in sitting position for at least 10 minutes. Each condition lasted for about five minutes. The average of the two rest conditions was used in the current study. The same data cleaning procedure as described above for the TRAILS-P and TRAILS-CC studies was used. Further details have been described before [43], [44].
In the NESDA cohort RMSSD was computed on the ECG-derived RR interval time series measured with the Vrije Universiteit Ambulatory Monitoring System (VU-AMS) [45], [46]. Subjects wore the VU-AMS device during a large part of the baseline assessment while participating in other assessments (medical examination, an interview and a computer task). Further details on the assessment have been described before [47].
In the PREVEND cohort RMSSD was computed on the RR interval time series based on the beat-to-beat blood pressure signals measured with the Portapres device (FMS Finapres Medical systems BV, Amsterdam, The Netherlands; [48]). Subjects were measured in the supine resting condition, after resting for five minutes, during 15 minutes, while holding their hand with the Portapres cuff at heart level. They were not allowed to talk or to move during the measurement. In short, from the 15 minute recording an artifact free period of at least 100 seconds but no longer than five minutes was selected for RMSSD calculations. The same data cleaning procedure as described above for the TRAILS-P, TRAILS-CC, and TWINS studies was used. Further details have been described before [49].
In the NTR cohort RMSSD was computed on the ECG-derived RR interval time series measured with the VU-AMS device during a 24-h ambulatory monitoring procedure which has been described before [38], [39]. Prompted by the VU-AMS device, the subjects filled out a diary every 30 (±10) minutes in which they logged their physical activities and bodily postures since the previous prompt. This diary information was combined with simultaneously recorded vertical accelerometer data to divide the entire 24-h recording into periods of homogeneous activities/postures (median duration 26.2 minutes). From the total set of ambulatory RMSSD measurements only the evening measurements were selected where subjects were seated calmly for at least five minutes but no longer than 60 minutes. The average of these measurements was taken as their RMSSD resting value.
SNP selection
For the GA cohorts, tagging SNPs were selected from HapMap phase II data (July, 2007) for seven genes including CHAT, SLC18A3, SLC5A7, CHRNB4, CHRNA3, CHRNA5 and CHRM2 using Haploview pairwise tagging with a r2≥0.8. All the common SNPs with a minor allele frequency (MAF) ≥5% within the gene as well as −10 kb upstream and +10 kb downstream of the gene were covered. Since the HapMap had a limited number of SNPs available in the ACHE gene at the time when tagging SNPs were selected (July, 2007), the SNP data from the SeattleSNPs Variation Discovery Resource (www.pga.gs.washington.edu), which had completely re-sequenced ACHE in 23 white Americans, was used for tagging SNPs selection. Table 1 lists the physical position, number of the tagging SNPs selected, the genotyped and imputed SNPs for the eight key genes. Since SLC18A3 is within the first intron of the CHAT gene and CHRNB4, CHRNA3, and CHRNA5 are within one gene cluster, these genes were listed together as two loci in Table 1.
For the TRAILS-P, TRAILS-CC and TWINS cohorts, the same tagging SNPs as selected for the Caucasian GA cohorts described above were selected. For the NESDA cohort, which had genome wide genotyping data, SNPs within the same physical regions were selected for this study.
DNA extraction and genotyping
In the GA cohorts, genotyping of blood derived DNA was performed by BioServe Biotechnologies (Laurel, MD) using Sequenom’s high-throughput matrix-assisted laser desorption/ionization time-of-flight (MALDI-TOF) mass spectrometry. All genotyping results were reviewed manually for quality control. Two ‘no template’ controls and two positive controls were included on each plate. In addition, 107 sets of blinded controls were distributed throughout the plates for quality control purposes. There was excellent observer agreement in the randomly selected duplicates that were included for quality control purposes. The kappa statistic ranged from 0.89 to 1, with an average of 0.974.
In a subsample of the TRAILS-P, TRAILS-CC and TWINS cohorts, DNA was extracted from blood samples or (in a few cases) buccal swabs (Cytobrush) using a manual salting out procedure [50]. Genotyping was performed on the Golden Gate Illumina BeadStation 500 platform (Illumina Inc., San Diego, CA, USA), following the manufacturers protocol. Genotyping data and clustering was performed in BeadStudio 3.0 (Illumina Inc., San Diego, CA, USA). Concordance between DNA replicates showed a genotyping accuracy of 100%. Data cleaning was in line with standard procedures [51]. The call rate in the combined set of TRAILS-P, TRAILS-CC, and TWINS after data cleaning was 99.8%. In addition to this candidate gene genotyping, genome-wide genotyping of TRAILS-CC was performed using the Illumina HumanCytoSNP12 v2 beadchip assay (Illumina, Inc; San Diego, CA, USA). Genotypes were called with the Illumina GenomeStudio software package (Illumina, Inc). Details on genotyping and quality control have been given elsewhere [52].
In the NESDA and NTR cohorts, genome-wide genotyping of a subsample of the subjects was performed as part of the GAIN initiative to detect genetic variation associated with major depression [40]. Genotyping was performed on Affymetrix 500 K GeneChip arrays. A detailed description of genotyping and quality control criteria has been given previously [53].
For PREVEND genome-wide genotyping of a subsample of the cohort was performed using the Illumina HumanCytoSNP12 v2 beadchip assay (Illumina, Inc; San Diego, CA, USA). Genotypes were called with the Illumina GenomeStudio software package (Illumina, Inc). Details on genotyping and quality control have been described before [54].
Imputation
Non-genotyped SNPs within the eight genes, as well as −10 kb upstream and +10 kb downstream of the genes, were imputed. For the GA, TRAILS-P and TWINS cohorts, for seven genes (CHAT, SLC18A3, SLC5A7, CHRNB4, CHRNA3, CHRNA5 and CHRM2), imputation was performed using Phase II CEU HapMap data (release 24, build 36) as the reference database using IMPUTE version 0.3.2[55], while for the ACHE gene, imputation was performed using Phase III CEU HapMap data (release 2, build 36) because coverage of this gene in HapMap Phase II (release 24) was minimal. For the NESDA cohort imputation was done with IMPUTE version 2.1.2 [56] using the ALL panel of the 1000 Genomes phase I integrated variant set of March 2012. For the NTR and TRAILS-CC cohorts imputation was done with IMPUTE version 0.3.2 [55] using Phase II CEU HapMap data (release 22, build 36) reference population. For the PREVEND cohort imputation was done with Beagle version 3.1.0 [57] using the Phase II CEU HapMap (release 23a, build 36) reference population.
Statistical analysis
Hardy-Weinberg equilibrium (HWE) was tested in each cohort by a χ2 test with 1 df. To prevent inflated significance, these tests were performed in data including only one of the twins, chosen at random, in the GA and TWINS cohorts. For the GA cohort, SNPs were excluded if they violated HWE (p<0.01) and at the same time the kappa statistic for the randomly selected duplicates on the genotype of these SNPs was less than 0.90. For the TRAILS-P and TWINS cohorts, SNPs were excluded if they violated HWE (p<0.01). For the NESDA, NTR, PREVEND, and TRAILS-CC cohorts SNPs with a HWE p-value below10−5 were excluded.
The genotyped and imputed SNPs were examined for their associations with RMSSD under an additive model of inheritance with adjustment for age and sex. For imputed SNPs dosages were calculated and used in the analyses to properly account for uncertainty in the imputed genotypes. RMSSD was log-transformed in the analyses. For the GA and TWINS cohorts, generalized estimating equations were used with exchangeable correlation structure to account for familial correlations. A fixed effects meta-analysis on these four cohorts was conducted using the inverse variance weighted method. The Q and I2 statistics were used to estimate between study heterogeneity. To help interpret I2’s magnitude we used guidelines from Higgins and colleagues [58] in which I2≤25% marks low, I2 = 50% medium, and I2≥75% high heterogeneity.
To optimize the chances of discovery we used a relatively lenient p<0.05 cut-off in the discovery stage. That is, all SNPs with p<0.05 in the discovery stage were taken forward for replication in three independent samples: PREVEND, TRAILS-CC and NTR. For all three replication cohorts results for these SNPs were extracted from genome-wide analyses of log-transformed RMSSD, who analyzed it by means of SNPTEST [55] using an additive model with age and sex as covariates and taking the uncertainty in the imputed genotypes into account. Results from discovery and replication phase were subsequently combined in a meta-analysis using METAL [59]. When not mentioned analyses were performed in StataSE software (version 12, StataCorp, College Station, Texas, USA).
Results
The characteristics of the subjects included in the analysis in each of the cohorts are given in Table 2. Of the 95 SNPs selected in the design stage, 87 were successfully genotyped in the GA cohort and 90 in the TRAILS-P and TWINS cohorts (Table 1). No further SNPs were excluded as all met our criteria for inclusion (see methods). All imputed SNPs in the discovery cohorts had at least an imputation quality (info score) of 0.5. In total, 443 genotyped and imputed SNPs were included in the association analysis with RMSSD.
Discovery stage
In the discovery phase the association between the candidate SNPs and RMSSD were tested in four cohorts (n = 3429; GA, TRAILS-P, TWINS and NESDA). For 25 SNPs an association with RMSSD (p<0.05) was found (see Table 3).
When testing the genetic heterogeneity of the SNPs in this phase, the point estimate for I2 was significant for only one SNP (rs17269293; I2 = 75.8%). Heterogeneity for all other SNPs was either completely absent or low and non-significant, confirming that performing a fixed effects meta-analysis on these data is justified (Table 3).
Replication stage and the overall meta-analysis
To validate the association found in the discovery phase, the 25 lead SNPs (Table 3), were taken forward to the replication phase in three independent cohorts (n = 3311 PREVEND, TRAILS-CC and NTR). Imputation quality was very high for NTR (all 25 SNPs: info score >0.72; 21 SNPs: info score >0.90) and TRAILS-CC (all 25 SNPs: info score >0.78; 24 SNPs: info score >0.94) and somewhat lower for PREVEND, which had 18 SNPs of high quality imputation (r2>0.73), but 7 SNPs with lower imputation quality (r2<0.35). For one of these SNPs (rs333214 in the SLC5A7 gene) data had to be excluded from the PREVEND cohort as the imputation quality was too low (r2<0.10) resulting in lower replication sample size (n = 516). None of the SNPs were replicated with a nominal p<0.05, except for one SNP. However, this SNP (rs3731683) showed an effect in the opposite direction as observed in the discovery stage (see the opposite beta coefficients for this SNP in Table 3). Thus, none of the SNPs were replicated.
The overall meta-analysis combining discovery and replication stage (n = 6740) gave no significant Bonferroni corrected (p<0.0005, conservatively assuming 100 tests for the almost 100 independent tagging SNPs from Table 1) findings for any of the 25 SNPs (Table 3).
Discussion
The current study comprehensively tested the association between all common (MAF≥5%) variants in eight key genes of the acetylcholine pathway and an established measure of HRV, the RMSSD, in subjects of European descent. The eight key genes were carefully selected in this two-stage candidate gene study on the basis of their involvement in biosynthesis, transport, breakdown, and receptor binding of acetylcholine. However, none of the common genetic variants in these eight genes of the acetylcholine pathway were significantly associated with RMSSD.
Neumann and colleagues [13] genotyped a single SNP (rs333229) in the choline transporter gene SLC5A7 in their study of 413 individuals of European ancestry (aged 30 to 54 years). They found an association with the power in the high frequency (HF) band (HF power of Fourier transformed inter-beat-interval time series of RR intervals of the ECG, is another measure of HRV comparable to RMSSD) obtained during five minutes rest: compared with GG homozygotes, T allele carriers had higher HF power. In a more recent study on air pollution in 61 glucose intolerant subjects, T allele carriers of this SNP (rs333229, which was one of the 139 tested SNPs) had lower SDNN in response to air pollution [60]. Interestingly, the exact same SNP in SLC5A7 (rs333229) was genotyped in the current study but did not survive our discovery stage (i.e., was not significantly associated with RMSSD). As such, the significant results in the above two underpowered studies should probably be interpreted as chance findings, because our discovery stage alone (n = 3429) already had more than eight times the sample size of the largest of those two previous studies (n = 413).
Our study has a number of strengths compared to previous candidate gene studies that investigated only a single SNP in a single gene and/or had (very) small sample sizes. We comprehensively tested 443 genotyped and imputed common SNPs in eight key genes of the acetylcholine pathway in which our selected tagging SNPs covered all SNPs with MAF≥5% and r2≥0.8 ranging from −10 kb upstream to +10 kb downstream of those genes. Another strong point of our study is that we aimed to replicate the initially discovered associations in independent samples. As replication failed, our initial findings should be labeled as false-positive findings, reinforcing the idea that independent replication is essential in candidate gene studies. Our two stage association study with large combined sample size of n = 6740 had more than 80% power with a Bonferroni corrected alpha of 0.0005 for 100 tests to detect SNPs explaining as little as 0.32% of the variance in RMSSD. As such it is highly unlikely that common SNPs in the investigated genes explain a considerable part of individual differences in RMSSD in rest.
A potential limitation of our study is that the RMSSD measures in the different cohorts were assessed in both sitting and supine positions, which might have introduced heterogeneity in our RMSSD data. Using these heterogeneous phenotypes might have decreased our power to detect associations somewhat. Not all SNPs were genotyped in our study, which may have limited our power to find or confirm significant associations. However, common SNPs in the discovery cohorts were well imputed (info score >0.5), because we selected tagging SNPs with good coverage of the genes and used 1000 Genomes imputation in NESDA. For the replication cohorts the imputation quality was dependent on the coverage of the GWAS chip and the imputation software used. Imputation quality was very high for NTR and TRAILS-CC, but somewhat lower for 7 SNPs in PREVEND reducing the effective sample size and power for these SNPs. We guarded against potential biases caused by the imputation process such as inflated significance, through properly accounting for uncertainty in these imputed genotypes in our analyses. Furthermore, the findings and conclusions of the current study are not generalizable to individuals of non-European descent.
Another, perhaps more compelling, reason for our null-finding is that our approach was based on currently existing knowledge about the acetylcholine pathway. Using a hypothesis-driven design is a strength, but has also clear downsides. First of all, not all known elements of the acetylcholinergic signaling cascade were considered. Second, current knowledge about the acetylcholine pathways may be vastly incomplete. Third, RMSSD is known to be a complex phenotype influenced by multiple genes [61], physiological and psychological systems, and their interactions, many of which may be outside of the acetylcholine pathways.
In conclusion, no common genetic variants for HRV were identified in the largest and most comprehensive candidate gene study on the acetylcholine pathway to date. Future gene finding efforts for RMSSD may want to focus on hypothesis free approaches such as the genome-wide association study.
Author Contributions
Conceived and designed the experiments: HR CAH AJO RTG PVDH DIB GW BWJHP EJCDG XW HS. Performed the experiments: HR XD SS AMVR JL NV CMML JJH GW EJCDG. Analyzed the data: LMM XD SS IMN PJVDM XW. Wrote the paper: HR LMM CAH XD SS AJO AMVR IMN PJVDM JL RTG PVDH NV CMML DIB JJH GW BWJHP EJCDG XW HS.
References
- 1. Task Force of the European Society of Cardiology and the North American Society of pacing and Electrophysiology (1996) Heart rate variability: standards of measurement, physiological interpretation and clinical use. Circulation 93: 1043–1065.
- 2. Goedhart AD, Van der Sluis S, Houtveen JH, Willemsen G, De Geus EJC (2007) Comparison of time and frequency domain measures of RSA in ambulatory recordings. Psychophysiology 44: 203–215.
- 3. Thayer JF, Yamamoto SS, Brosschot JF (2010) The relationship of autonomic imbalance, heart rate variability and cardiovascular disease risk factors. Int J Cardiol 141: 122–131.
- 4. Singh JP, Larson MG, Tsuij H, Evans JC, O’Donnell CJ, et al. (1998) Reduced heart rate variability and new-onset hypertension: insights into pathogenesis of hypertension: the Framingham Heart Study. Hypertension 32: 293–297.
- 5. Thayer JF, Lane RD (2007) The role of vagal function in the risk for cardiovascular disease and mortality. Biol Psychol 74: 224–242.
- 6. Buccelletti F, Gilardi E, Scaini E, Galiuto L, Persiani R, et al. (2009) Heart rate variability and myocardial infarction: systematic literature review and metanalysis. Eur Rev Med Pharmacol Sci 13: 299–307.
- 7. Dekker JM, Crow RS, Folsom AR, Hannan PJ, Liao D, et al. (2000) Low heart rate variability in a 2-minute rhythm strip predicts risk of coronary heart disease and mortality from several causes - The ARIC study. Circulation 102: 1239–1244.
- 8.
Guyton AC, Hall JE (2006) Textbook of medical physiology. Philadelphia: Elsevier Saunders.
- 9. Snieder H, Boomsma DI, Van Doornen LJ, de Geus EJ (1997) Heritability of respiratory sinus arrhythmia: dependency on task and respiration rate. Psychophysiology 34: 317–328.
- 10. Wang X, Thayer JF, Treiber F, Snieder H (2005) Ethnic differences and heritability of heart rate variability in African- and European American youth. Am J Cardiol 96: 1166–1172.
- 11. Riese H, Rosmalen JGM, Ormel J, van Roon AM, Oldehinkel AJ, et al. (2007) The genetic relationship between neuroticism and autonomic function in female twins. Psychol Med 37: 257–267.
- 12. Imumorin IG, Dong Y, Zhu H, Poole JC, Harshfield GA, et al. (2005) A gene-environment interaction model of stress-induced hypertension. Cardiovasc Toxicol 5: 109–132.
- 13. Neumann SA, Lawrence EC, Jennings JR, Ferrell RE, Manuck SB (2005) Heart rate variability is associated with polymorphic variation in the choline transporter gene. Psychosom Med 67: 168–171.
- 14. Cook LJ, Ho LW, Wang L, Terrenoire E, Brayne C, et al. (2005) Candidate gene association studies of genes involved in neuronal cholinergic transmission in Alzheimer’s disease suggests choline acetyltransferase as a candidate deserving further study. Am J Med Gen Part B-Neuropsychiatr Gen 132B: 5–8.
- 15. Mubumbila V, Sutter A, Ptok U, Heun R, Quirin-Stricker C (2002) Identification of a single nucleotide polymorphism in the choline acetyltransferase gene associated with Alzheimer’s disease. Neurosci Lett 333: 9–12.
- 16. Giubilei F, Strano S, Imbimbo BP, Tisei P, Calcagnini G, et al. (1998) Cardiac autonomic dysfunction in patients with Alzheimer disease: Possible pathogenetic mechanisms. Alzheimer Dis Assoc Disord 12: 356–361.
- 17. Aharonperetz J, Harel T, Revach M, Benhaim SA (1992) Increased sympathetic and decreased parasympathetic cardiac innervation in patients with Alzheimers disease. Arch Neurol 49: 919–922.
- 18. Erickson JD, Varoqui H, Schafer MKH, Modi W, Diebler MF, et al. (1994) Functional identification of a vesicular acetylcholine transporter and its expression from a cholinergic gene locus. J Biol Chem 269: 21929–21932.
- 19. Duga S, Solda G, Asselta R, Bonati MT, Dalpra L, et al. (2001) Characterization of the genomic structure of the human neuronal nicotinic acetylcholine receptor CHRNA5/A3/B4 gene cluster and identification of novel intragenic polymorphisms. J Hum Genet 46: 640–648.
- 20. Singh JP, Larson MG, O’Donnell CJ, Tsuji H, Corey D, et al. (2002) Genome scan linkage results for heart rate variability (the Framingham Heart Study). Am J Cardiol 90: 1290–1293.
- 21. Brodde OE, Bruck H, Leineweber K, Seyfarth T (2001) Presence, distribution and physiological function of adrenergic and muscarinic receptor subtypes in the human heart. Basic Res Cardiol 96: 528–538.
- 22. Fisher JT, Vincent SG, Gomeza J, Yamada M, Wess J (2004) Loss of vagally mediated bradycardia and bronchoconstriction in mice lacking M-2 or M-3 muscarinic acetylcholine receptors. Faseb Journal 18: 711–739.
- 23. Ge DL, Dong YB, Wang XL, Treiber FA, Snieder H (2006) The Georgia Cardiovascular Twin Study: Influence of genetic predisposition and chronic stress on risk for cardiovascular disease and type 2 diabetes. Twin Res Hum Genet 9: 965–970.
- 24. Ormel J, Oldehinkel AJ, Sijtsema J, van Oort F, Raven D, et al. (2012) The TRacking Adolescents’ Individual Lives Survey (TRAILS): Design, current status, and selected findings. J Am Acad Child Adolesc Psychiatry 51: 1020–1036.
- 25. Riese H, Rijsdijk FV, Snieder H, Ormel J (2013) Twin Interdisciplinary Neuroticism Study (TWINS). Twin Res Hum Genet 16: 268–270.
- 26. Penninx BWJH, Beekman ATF, Smit JH, Zitman FG, Nolen WA, et al. (2008) The Netherlands Study of Depression and Anxiety (NESDA): rationale, objectives and methods. Int J Methods Psychiatr Res 17: 121–140.
- 27. Pinto-Sietsma SJ, Janssen WM, Hillege HL, Navis G, de Zeeuw D, et al. (2000) Urinary albumin excretion is associated with renal functional abnormalities in a nondiabetic population. J Am Soc Nephrol 11: 1882–1888.
- 28. Willemsen G, Vink JM, Abdellaoui A, den Braber A, van Beek JHDA, et al. (2013) The adult Netherlands Twin Register: Twenty-five years of survey and biological data collection. Twin Res Hum Genet 16: 271–281.
- 29. Li ZB, Snieder H, Su SY, Ding XH, Thayer JF, et al. (2009) A longitudinal study in youth of heart rate variability at rest and in response to stress. Int J Psychophysiol 73: 212–217.
- 30. Wang XL, Ding XH, Su SY, Li ZB, Riese H, et al. (2009) Genetic influences on heart rate variability at rest and during stress. Psychophysiology 46: 458–465.
- 31. Dekkers JC, Snieder H, van den Oord EJ, Treiber FA (2002) Moderators of blood pressure development from childhood to adulthood: a 10-year longitudinal study. J Pediatr 141: 770–779.
- 32. Snieder H, Dong Y, Barbeau P, Harshfield GA, Dalageogou C, et al. (2002) Beta2-adrenergic receptor gene and resting hemodynamics in European and African American youth. Am J Hypertens 15: 973–979.
- 33. Huisman M, Oldehinkel AJ, de WA, Minderaa RB, de Bildt A, et al. (2008) Cohort profile: the Dutch ‘TRacking Adolescents’ Individual Lives ‘Survey’; TRAILS. Int J Epidemiol 37: 1227–1235.
- 34. de Winter A, Oldehinkel AJ, Veenstra R, Brunnekreef JA, Verhulst FC, et al. (2005) Evaluation of non-response bias in mental health determinants and outcomes in a large sample of pre-adolescents. Eur J Epidemiol 20: 173–181.
- 35. Dietrich A, Riese H, van Roon AM, Van Engelen K, Ormel J, et al. (2006) Spontaneous baroreflex sensitivity in (pre)adolescents. J Hypertens 24: 345–352.
- 36. Joosten H, Izaks GJ, Slaets JPJ, de Jong PE, Visser ST, et al. (2011) Association of cognitive function with albuminuria and eGFR in the general population. Clin J Am Soc Nephrol 6: 1400–1409.
- 37. de Winter A, Oldehinkel AJ, Veenstra R, Brunnekreef JA, Verhulst FC, et al. (2005) Evaluation of non-response bias in mental health determinants and outcomes in a large sample of pre-adolescents. Eur J Epidemiol 20: 173–181.
- 38. Kupper NH, Willemsen G, van den Berg M, de Boer D, Posthuma D, et al. (2004) Heritability of ambulatory heart rate variability. Circulation 110: 2792–2796.
- 39. Kupper N, Willemsen G, Posthuma D, de Boer D, Boomsma DI, et al. (2005) A genetic analysis of ambulatory cardiorespiratory coupling. Psychophysiology 42: 202–212.
- 40. Boomsma DI, Willemsen G, Sullivan PF, Heutink P, Meijer P, et al. (2008) Genome-wide association of major depression: description of samples for the GAIN Major Depressive Disorder Study: NTR and NESDA biobank projects. Eur J Hum Genet 16: 335–342.
- 41. Niskanen JP, Tarvainen MP, Ranta-Aho PO, Karjalainen PA (2004) Software for advanced HRV analysis. Comput Methods Programs Biomed 76: 73–81.
- 42.
Mulder LJM, van Dellen HJ, van der Muelen P, Opheikens B (1988) CARSPAN: a spectral analysis program for cardiovascular time series. In: Maarse FJ, Mulder L, Akkerman A, editors. Computers in psychology: methods, instrumentation and psychodiagnostics. Amsterdam: Swets and Zeitlinger. pp. 39–47.
- 43. Riese H, Rijsdijk FV, Ormel J, van Roon AM, Neeleman J, et al. (2006) Genetic influences on baroreflex sensitivity at rest and during mental stress. J Hypertens 24: 1779–1786.
- 44. Rijsdijk FV, Riese H, Tops M, Snieder H, Brouwer WH, et al. (2009) Neuroticism, recall bias and attention bias for valenced probes: a twin study. Psychol Med 39: 45–54.
- 45. De Geus EJC, Willemsen GHM, Klaver CHAM, Van Doornen LJP (1995) Ambulatory measurement of respiratory sinus arrhytmia and respiration rate. Biol Psychol 41: 205–227.
- 46. Willemsen GHM, De Geus EJC, Klaver CHAM, Van Doornen LJP, Carroll D (1996) Ambulatory monitoring of the impedance cardiogram. Psychophysiology 33: 184–193.
- 47. Licht CM, de Geus EJ, Zitman FG, Hoogendijk WJ, van DR, et al. (2008) Association between major depressive disorder and heart rate variability in the Netherlands Study of Depression and Anxiety (NESDA). Arch Gen Psychiatry 65: 1358–1367.
- 48.
Wesseling KH (1990) Finapres, continuous noninvasive finger arterial pressure based on the method of Peñáz. In: Meyer-Sabellek W, Gotzen R, Anlauf M, Steingard RJ, editors. Blood Pressure Measurements. Frankfurt am Main: Steunkopff. pp. 161–172.
- 49.
Muñoz ML, van Roon AM, Riese H, Oostenbroek E, Westrik I, et al. (2014) Validity of (ultra-)short heart rate variability recordings. unpublished data.
- 50. Miller SA, Dykes DD, Polesky HF (1988) A simple salting out procedure for extracting DNA from human nucleated cells. Nucleic Acids Research 16: 1215.
- 51.
Nolte IM, McCaffery JM, Snieder H (2010) Candidate gene and genome-wide association studies in behavioral medicine. In: Steptoe A, editors. Handbook of behavioral medicine: methods and applications. New York: Springer.
- 52. Ganesh SK, Tragante V, Guo W, Guo YR, Lanktree MB, et al. (2013) Loci influencing blood pressure identified using a cardiovascular gene-centric array (vol 22, pg 1663, 2013). Hum Molec Genet 22: 3394–3395.
- 53. Sullivan PF, De Geus EJC, Willemsen G, James MR, Smit JH, et al. (2009) Genome-wide association for major depressive disorder: a possible role for the presynaptic protein piccolo. Mol Psychiatry 14: 359–375.
- 54. Verweij N, Mahmud H, Leach IM, de Boer RA, Brouwers FP, et al. (2013) Genome-wide association study on plasma levels of midregional-proadrenomedullin and C-terminal-pro-endothelin-1. Hypertension 61: 602–608.
- 55. Marchini J, Howie B, Myers S, McVean G, Donnelly P (2007) A new multipoint method for genome-wide association studies by imputation of genotypes. Nat Genet 39: 906–913.
- 56.
Howie BN, Donnelly P, Marchini J (2009) A flexible and accurate genotype imputation method for the next generation of genome-wide association studies. Plos Genetics 5.
- 57. Browning BL, Browning SR (2009) A unified approach to genotype imputation and haplotype-phase inference for large data sets of trios and unrelated individuals. Am J Hum Genet 84: 210–223.
- 58. Higgins JPT, Thompson SG, Deeks JJ, Altman DG (2003) Measuring inconsistency in meta-analyses. BMJ 327: 557–560.
- 59. Willer CJ, Li Y, Abecasis GR (2010) METAL: fast and efficient meta-analysis of genomewide association scans. Bioinformatics 26: 2190–2191.
- 60. Hampel R, Breitner S, Schneider A, Zareba W, Kraus U, et al. (2012) Acute air pollution effects on heart rate variability are modified by SNPs involved in cardiac rhythm in individuals with diabetes or impaired glucose tolerance. Environ Res 112: 177–185.
- 61. Neumann SA, Tingley WG, Conklin BR, Shrader CJ, Peet E, et al. (2009) AKAP10 (I646 V) functional polymorphism predicts heart rate and heart rate variability in apparently healthy, middle-aged European-Americans. Psychophysiology 46: 466–472.