An epigenome-wide association analysis of cardiac autonomic responses among a population of welders.

DNA methylation is one of the potential epigenetic mechanisms associated with various adverse cardiovascular effects; however, its association with cardiac autonomic dysfunction, in particular, is unknown. In the current study, we aimed to identify epigenetic variants associated with alterations in cardiac autonomic responses. Cardiac autonomic responses were measured with two novel markers: acceleration capacity (AC) and deceleration capacity (DC). We examined DNA methylation levels at more than 472,506 CpG probes through the Illumina Infinium HumanMethylation450 BeadChip assay. We conducted separate linear mixed models to examine associations of DNA methylation levels at each CpG with AC and DC. One CpG (cg26829071) located in the GPR133 gene was negatively associated with DC values after multiple testing corrections through false discovery rate. Our study suggests the potential functional importance of methylation in cardiac autonomic responses. Findings from the current study need to be replicated in future studies in a larger population.

Introduction

Despite the fact that numerous therapies have been developed in cardiovascular medicine, cardiovascular disease (CVD) remains one of the leading causes of morbidity and mortality worldwide. Understanding the biological mechanisms and predictors of CVD may have important implications in prevention and treatment of this public health concern. The cardiac autonomic nervous system consists of two major branches—the sympathetic branch and the parasympathetic branch—and it is essential in modulation of cardiac electrophysiology. Failure of parasympathetic control, which is often measured by conventional heart rate variability (HRV) parameters, has been a strong predictor of cardiovascular mortality in both high-risk and low-risk populations. Meanwhile, the associations of reduced HRV with various environmental and occupational pollutants, such as particulate matter air pollution with aerodynamic diameter <2.5 µm (PM2.5), ozone, and heavy metals, have been consistently reported. However, the biological mechanisms by which these pollutants trigger alterations in cardiac autonomic responses remain unclear.

Previous genome-wide association studies (GWAS) have identified several genetic variants associated with cardiac autonomic responses. More recent studies suggested that epigenetic regulation, such as DNA methylation, may also play an important role in alterations in cardiac autonomic responses. PM2.5 exposure has been associated with both gene-specific methylation (e.g., iNOS gene), and HRV changes in previous air pollution studies. To date, no epigenome-wide association study (EWAS) has been conducted to systematically examine DNA methylation and cardiac autonomic responses. Understanding this relationship may help elucidate the biological pathway of cardiac autonomic responses, which may further have implications for the management and development of therapies for cardiovascular disease.

Welders have higher exposure to welding fumes than other trades involved with welding. Welding fumes often contain a variety of hazards, such as metals, gases, and chemicals, which can cause various cardiovascular disorders. We previously observed in several cross-sectional studies that welders had decreased HRV parameters following a 6-hour work shift. Further, exposure to metal-rich welding fumes across the work shift was strongly associated with decreased HRV parameters. In this study, we conducted an EWAS to identify epigenetic variants that were associated with exposure-induced cardiac autonomic responses. Two novel markers—acceleration capacity (AC) and deceleration capacity (DC)—were used as indices of cardiac autonomic responses.

Results

Characteristics of the study population are shown in Table 1. The study population was comprised of 75 male welders. The majority of them were Caucasian (93.3%) with an age average of 41.6 years (range: 21.6–71.2). Close to 40% were current smokers, and average BMI was 28.7 kg/m2. The baseline DC and AC values were 7.64 [Standard deviation (SD) = 3.64] ms and −6.42 (SD = 4.05) ms, respectively. Compared with baseline, the estimated post-work natural killer (NK) cell, monocyte, and granulocyte proportions were significantly higher, whereas post-work CD4+ T cell, CD4+ T cell, and B cell proportions were not significantly different (Table 2).

Table 1.

Demographics of study population (n = 75).

Characteristic
Male	75 (100)a
Current smoker	28 (37.3)a
Average age (years)	41.6 ± 12.8b
Age range (years)	21.6–71.2
Average BMI (kg/m2)	28.7 ± 5.3b
Acceleration capacity (ms)	−6.42 ± 4.05b
Deceleration capacity (ms)	7.64 ± 3.64b

n (%)

mean ± SD

Table 2.

Estimated cell type composition (%) by time.

Cell type	Post-worka	Baselinea	P valueb
CD8+T cell	4.27 (4.66)	4.08 (4.59)	0.39
CD4+T cell	13.89 (4.67)	13.73 (4.21)	0.81
NK cell	4.47 (4.61)	2.93 (3.22)	<0.01
B cell	3.68 (2.5)	3.92 (2.35)	0.08
Monocyte	7.64 (2.45)	7.08 (2.5)	0.04
Granulocyte	64.15 (9.23)	66.27 (8.13)	<0.01

Unadjusted mean percentage (SD); n = 208

Comparison between prior and post exposure using linear mixed regression models.

We examined the associations of methylation levels at 472,506 CpG probes in the whole-genome with AC and DC, while adjusting for age, BMI, smoking status, work-shift, time of sample collection, and cell type composition for whole blood. For DC analyses, only one CpG (cg26829071), located in the GPR133 gene on chromosome 12 (Fig. 2), reached genome-wide significance using a threshold of false discovery rate (FDR) q-value ≤ 0.1(Table 3). There was a significant negative association between methylation levels at this CpG and DC (β = −0.32, SE = 0.06; P = 7.75E-08). In addition, we observed suggestive evidence of associations at cg12991522 (β = 0.47, SE = 0.09; P = 1.82E-07), located in the PPL gene on chromosome 16, with DC. For AC analyses, one CpG (cg15273468) was genome-wide significant with FDR q values ≤ 0.1 (Table 4). However, this CpG was not linked with any known functional genes.

Figure 2.

Manhattan Plots for AC and DC analysis. The Manhattan plot denoting the P-values for the association of DNA methylation levels at autosomal sites (X-axis).

Table 3.

List of top-ranking CpGs associated with DC.

Name	Chromosome	Strand	Gene name	Location	β	SE	P value	FDR q value
cg26829071	12	F	GPR133	Body	−0.32	0.06	7.75E-08	0.04
cg12991522	16	F	PPL	1stExon	0.47	0.09	1.82E-07	0.04
cg12071328	11	F	NELL1	Body	0.39	0.08	2.81E-06	0.44
cg20928238	14	R			−0.24	0.05	6.11E-06	0.64
cg01520463	2	F	SDC1	Body	0.18	0.04	8.03E-06	0.64
cg25048985	12	F			−0.29	0.06	1.05E-05	0.64
cg12454595	22	R	FAM19A5	Body	0.68	0.15	1.14E-05	0.64
cg12479444	7	F			0.62	0.13	1.21E-05	0.64
cg12044338	1	F	AQP10	TSS1500	0.41	0.09	1.27E-05	0.64
cg06716283	19	F	RELB	TSS1500	0.39	0.09	1.42E-05	0.64

Table 4.

List of top-ranking CpGs associated with AC.

Name	Chromosome	Strand	Gene name	Location	β	SE	P value	FDR q value
cg15273468	10	R			−0.31	0.05	6.60E-08	0.03
cg19458608	3	R	LOC729375	TSS1500	−0.92	0.19	5.67E-06	1.00
cg14004557	11	F			−0.28	0.06	1.81E-05	1.00
cg26905268	19	R	KANK2	Body	−0.54	0.12	2.00E-05	1.00
cg05434952	17	R	CCDC137	1stExon	−0.95	0.22	2.45E-05	1.00
cg10414946	11	F	MS4A2	Body	−0.42	0.10	2.64E-05	1.00
cg02069944	17	R	AZI1	Body	1.28	0.29	2.83E-05	1.00
cg27445005	18	R			0.49	0.11	2.98E-05	1.00
cg11945929	7	F			−0.15	0.03	3.05E-05	1.00
cg24078451	6	F	PSMB8	1stExon	−0.36	0.08	3.83E-05	1.00

The Quantile-Quantile (Q-Q) plots suggested sizable inflation in AC analyses and better calibrated test statistic in DC analyses. To quantify this observation, we calculated the genomic inflation factor (λ) as the median of the observed distribution of the test statistic divided by the expected median. We observed λ = 0.81 for AC analyses and λ = 1.02 for DC analyses (Fig. 1). Taken together, the observed associations for DC analyses were unlikely due to population stratification.

Figure 1.

Quantile-Quantile (Q-Q) plot. The observed P-values (Y-axis) were plotted against the expected P-values under the null hypothesis (X-axis). The red diagonal line denotes the pattern under null hypothesis.

Discussion

In the current study, we aimed to identify epigenetic locations at which the methylation levels were associated with cardiac autonomic changes among a population of welders. There was a significantly negative association between methylation at cg26829071 and DC values. Annotation analysis indicated that this significant CpG was located in the gene body of the GPR133 gene.

GPR133 encodes an adhesion G-protein-coupled receptor (aGPCR), which is commonly characterized by long extracellular N termini that are composed of a seven transmembrane spanning domain. It is primarily expressed in the central nervous system as well as heart (ventricles, atria, and septal tissues). Bohnekamp et al. have reported a concentration-dependent relationship between GPR133 and intracellular cyclic adenosine 3′,5′-cyclic adenosine monophosphate (cAMP) levels, suggesting that GPR133 may be coupled to the Gs protein, activate adenylate cyclase, and stimulate G protein cascades. In the heart, the activation of adenylate cyclase is associated with various cardiovascular effects including modulation of heart rate. Adenylyl cyclase may activate the production of intracellular cAMP, which serves as the second messenger that further binds to protein kinase A and modulates cardiac contractility. Meanwhile, intracellular cAMP is also essential for the generation of action potential in the sinoatrial node. However, further research is needed to replicate previous findings and to explain the role of GRP133 in cardiac autonomic dysfunction.

Both toxicological studies and GWAS have documented associations of several subtypes of GPCRs with adverse cardiovascular effects, such as hyperproliferative vascular malformations (GPR124), myocardial wall thinning (GPR126), stroke (CELSR1), and myocardial infarction (CELSR2). To date, several GWAS studies have been conducted to explore genetic contributions to cardiac autonomic responses. Arking et al. identified the NOS1AP (CAPON) gene associated with the electrocardiographic (ECG) QT interval variations among participants from the KORA cohort in Germany. Marroni et al. observed similar findings and they additionally identified the associations of variants in GPR133 gene with the ECG RR interval alterations. Newton-Cheh et al. examined 70,987 common genetic variants and six conventional heart rate variability (HRV) parameters among 1345 participants from the Framingham Heart Study Original and Offspring cohort and found there was no genomic hit that yielded a genome-wide significance. Our study, however, expands the literature with EWAS and suggests that epigenetic regulation of the GPR133 gene may also play a role in heart rate modulation.

The role of DNA methylation in regulation of gene-expression may vary depending upon different genomic contexts. DNA methylation in the promoter sequences is known to downregulate gene expression, whereas methylation in the gene body is often positively associated with gene expression. Decreased DC, which reflects impaired cardiac autonomic function, has been identified as a strong predictor of cardiovascular mortality among patients after myocardial infarction. Our study suggests the important role of DNA methylation in cardiac autonomic dysfunction. However, the underlying biological mechanisms of whether and how methylation at cg26829071 may affect gene expression and induce decreased DC remain to be elucidated.

The second-ranked CpG associated with DC is located in the PPL gene (encoding periplakin). The PPL gene is a protein-coding gene responsible for keratinocyte intercellular adhesion as well as tissues integrity. Decreased periplakin expression has been associated with urothelial carcinoma of the urinary bladder (UCB), and serum periplakin has been suggested as a biomarker for UCB. However, there was no evidence suggesting a role for periplakin in cardiovascular disease. Its association with cardiac autonomic responses may warrant further research.

To the best of our knowledge, this is the first EWAS investigating epigenetic variants associated with cardiac autonomic responses. GPCRs are one of most studied receptor families that are used as pharmacological targets. In the heart, the adrenergic GPCR signaling pathways are within the major targets of pharmaceuticals for the treatment of cardiovascular disease. For example, beta-blockers have been extensively used for the management of arrhythmia, hypertension, and chronic heart failure. Our study indicated a potential role of GPR133 in cardiac autonomic dysfunction, primarily through affecting the deceleration capacity of heart rate. Future studies may investigate the association of the GPR133 gene with various cardiovascular diseases, such as abnormal heart rhythms and heart failure.

Our study is limited by the relatively small sample size and lack of independent replication. Therefore, our findings need to be confirmed in future studies in a larger population. In addition, previous studies have reported strong associations of acute and chronic PM2.5 exposure with DNA methylation in several genes as well as cardiac autonomic responses. Both short-term and long-term particulate matter exposure from welding fumes were significantly associated with AC and DC as well as conventional HRV parameters changes among this welder cohort. PM2.5 exposure has also been negatively associated with gene-specific (iNOS) methylation in a population of boilermaker construction workers as well as elderly men from the VA Normative Aging Study. Hence, we cannot preclude the potential confounding effects of PM2.5 due to lack of data in this study. Meanwhile, there might be other potential confounders in the occupational settings, such as heavy metals and psychological factors, in addition to PM2.5.

In conclusion, results of our study indicate that methylation at the GPR133 gene may play a role in cardiac autonomic dysfunction. However, we cannot preclude chance or bias due to lack of replication, small sample size, and potential confounding by pollutants exposure. Hence, our preliminary findings need to be confirmed in future research.

Materials and methods

Study population

The Harvard Boilermakers Study is a prospective cohort study conducted among a population of welders from a local boilermaker union in Quincy, MA. These welders primarily assemble or weld boilers in large power plants. The study population included 75 welders recruited from 6 sampling occasions in January 2003, January 2004, June 2010, January 2011, June 2011, and June 2012. Workers were allowed to participate multiple times within the 6 sampling occasions. The inclusion criteria were: a) male welders ≥ 18 years of age; b) unionized welders including both apprentice and journeyman; c) contributed blood and ECG samples at least once. The exclusion criteria were self-reported physician diagnosed cardiovascular disease prior participation.

During each sampling occasion, participants were monitored at a union welding school on an approximate 6-hour workday. Their major occupational activities included welding, grinding, and cutting tasks. We collected urine, blood, and resting ECG recordings at both prior- (baseline) and post-work. We also collected self-administrated questionnaires including age, height, weight, current smoking status and/or smoking history, medical history, and medication use during the previous six months from each participant at baseline, with this information being updated at each sampling occasion. The Harvard T. H. Chan School of Public Health Institutional Review Board reviewed and approved the study protocol and we obtained a written informed consent from each participant at each sampling occasion.

Electrocardiographic recording and sample analysis

During each sampling occasion, participants were fitted with a 7-lead ambulatory ECG Holter monitor. We collected 12-minute resting ECG recordings from each participant at both prior- and post-work. During the 12-minute resting period, participants were asked to remain seated and quiet; walking, talking, or eating was not allowed. The ECGs samples were then sent to the Cardiovascular Epidemiology Research Unit of Beth Israel Deaconess Medical Center (Boston, Massachusetts, USA), where trained technicians blinded to exposure status processed and analyzed these samples. The methods of processing ECGs samples for AC and DC analysis have been discussed in our previous study. Briefly, to compute DC values, RR intervals longer than the immediate preceding interval were defined as anchors; to compute AC values, RR intervals shorter than the immediate preceding interval were defined as anchors. Segments of interval data that had the same size around these anchors were identified and aligned at those anchors. Upon alignment of all segments, RR intervals at all defined anchors (X0), immediate preceding (X1) and following the anchors (X-1) were averaged separately. The quantities of AC or DC were obtained by computing the difference between the sum of X0 and X1 and the sum of (X-1) and (X-2). For ECG sample analysis, the first two minutes of the ECG recordings were discarded to allow for acclimation. AC and DC quantities were summarized every 5 minutes through the PRSA method, as described by Bauer et al. A total of 208 5-minute AC and DC values were obtained and repeated measurements of AC and DC values were modeled separately for each subject.

DNA methylation profiling and data quality control

We collected whole blood samples (n = 208) from each participant at both prior- and post-work through venous phlebotomy in EDTA tubes. Plasma was extracted and blood pellets were bisulfite-converted. DNA methylation levels of the entire genome that covers more than 480,000 CG dinucleotide (CpG) probes from both pre- and post-shift blood samples were determined through Infinium HumanMethylation450 BeadChip assay (Illumina, Inc.) following the Infinium HD Methylation Assay protocol guide. The BeadChips were scanned using the Illumina iScan and raw data was imported into GenomeStudio where image intensities were extracted. A methylation score (β value) was quantified on a scale from 0 to 1 to represent the percentage of methylated signal at each CpG, such that a value of 0 represents fully unmethylated signal and a value of 1 represents fully methylated signal.

Following background subtraction and dye bias adjustment, both sample-level and probe-level quality control procedures were performed. Samples with detection P-value > 0.05 in more than 5% probes were excluded; probes with detection P-value > 0.05 in more than 5% samples or probes with very low variation (CV < 5%) were omitted. Sex chromosomes and single nucleotide polymorphism-associated probes were also excluded. In addition, Houseman's algorithm in minfi was applied to estimate the cell type composition in blood samples. Beta-mixture quantile normalization was conducted for probe design bias correction. The β values were normalized with functional normalization and batch effect adjustment was performed with ComBat. In total, there were 472,506 CpG probes included in the final association analyses.

Statistical analysis

In the current study, we compared the post-work blood cell type composition with baseline composition through linear mixed effects models. We fitted separate linear mixed effects models with random intercepts to examine whether DNA methylation at each CpG was associated with AC and DC. A number of covariates including age (continuous in years), body mass index (BMI; continuous weight in kilograms divided by height squared in meters), current smoking status (non-smoker, former smoker, or current smoker), time (baseline or post-work) and cell type composition were adjusted for in all models as they have been suggested as potential confounders. In addition, all models were adjusted for the variable “time in day,” which reflects the time when blood and ECG samples were collected, to account for the potential circadian variations of AC and DC.

To correct for multiple comparisons, we computed the false discovery rate (FDR) adjusted P- value (q value). We considered an FDR q value ≤ 0.1 as statistically significant.

We fitted the QQ plot to visualize the expected distribution of test statistics of the association analyses across the CpG probes (X-axis) vs. the observed values (Y-axis). We also calculated the genomic inflation factor (lambda) to assess the bulk inflation and false positive rate. All analyses were performed with R 3.2.2 (R Core Team 2015).