Dissecting the association of autophagy-related genes with cardiovascular diseases and intermediate vascular traits: A population-based approach

Autophagy is involved in cellular homeostasis and maintenance and may play a role in cardiometabolic health. We aimed to elucidate the role of autophagy in cardiometabolic traits by investigating genetic variants and DNA methylation in autophagy-related genes in relation to cardiovascular diseases and related traits. To address this research question, we implemented a multidirectional approach using several molecular epidemiology tools, including genetic association analysis with genome wide association studies data and exome sequencing data and differential DNA methylation analysis. We investigated the 21 autophagy-related genes in relation to coronary artery disease and a number of cardiometabolic traits (blood lipids, blood pressure, glycemic traits, type 2 diabetes). We used data from the largest genome wide association studies as well as DNA methylation and exome sequencing data from the Rotterdam Study. Single-nucleotide polymorphism rs110389913 in AMBRA1 (p-value = 4.9×10−18) was associated with blood proinsulin levels, whereas rs6587988 in ATG4C and rs10439163 in ATG4D with lipid traits (ATG4C: p-value = 2.5×10−15 for total cholesterol and p-value = 3.1×10−18 for triglycerides, ATG4D: p-value = 9.9×10−12 for LDL and p-value = 1.3×10−10 for total cholesterol). Moreover, rs7635838 in ATG7 was associated with HDL (p-value = 1.9×10−9). Rs2447607 located in ATG7 showed association with systolic blood pressure and pulse pressure. Rs2424994 in MAP1LC3A was associated with coronary artery disease (p-value = 5.8×10−6). Furthermore, we identified association of an exonic variant located in ATG3 with diastolic blood pressure (p-value = 6.75×10−6). Using DNA methylation data, two CpGs located in ULK1 (p-values = 4.5×10−7 and 1×10−6) and two located in ATG4B (2×10−13 and 1.48×10−7) were significantly associated with both systolic and diastolic blood pressure. In addition one CpG in ATG4D was associated with HDL (p-value = 3.21×10−5). Our findings provide support for the role of autophagy in glucose and lipid metabolism, as well as blood pressure regulation.

Introduction Autophagy is a highly conserved catabolic process involved in the degradation of long-lived proteins and dysfunctional organelles [1]. Under conditions of oxidative stress, hypoxia and nutrient deprivation, autophagy is activated as a key mechanism of cell survival by degradation and recycling of damaged organelles and protein aggregates [2]. A disturbance in removal of these non-functional cellular components is a general impairment in housekeeping of the cell, and could lead to important phenotypic changes at cell and tissue levels [3]. Several studies, mainly based on animal models, have provided evidence about the role of autophagy in the progression of ageing and in particular in atherosclerosis, inflammation, and cardiovascular disease [2,[4][5][6][7][8]. A number of autophagy-related mechanisms have been studied for cell survival and in some instances the role of genetic variants involved in molecular mechanisms of autophagy were investigated in relation to cardiometabolic health in animal models [9]. Systemic knockout of autophagy-related genes (ATG) in mice has shown the role of dysfunctional autophagy in hyperglycemia [10], hypoinsulinemia and increased basal Ca2+ concentrations in vascular smooth muscle cells [11,12]. However, this association has not been studied at population level. Genetic and epigenetic variations in ATG could affect the autophagy process in human cells, modify certain metabolic traits and eventually cause susceptibility to cardiometabolic disorders [13,14]. In this study, by using Genome-Wide Association Studies (GWAS) data and data from the Rotterdam Study, a population-based prospective cohort study, we aimed to examine the association of genetic and epigenetic variation in ATG with intermediate vascular traits and cardiovascular outcomes.

Materials and methods
The Rotterdam Study has been approved by the Medical Ethics Committee of the Erasmus MC (registration number MEC 02.1015) and by the Dutch Ministry of Health, Welfare and Sport (Population Screening Act WBO, license number 1071272-159521-PG). The Rotterdam Study has been entered into the Netherlands National Trial Register (NTR; www.trialregister. nl) and into the WHO International Clinical Trials Registry Platform (ICTRP; www.who.int/ ictrp/network/primary/en/) under shared catalogue number NTR6831. All participants provided written informed consent to participate in the study and to have their information obtained from treating physicians.

Identification of autophagy-related genes
To select the set of autophagy-related genes to be studied, a literature search was performed and was updated by August 2018 in PubMed using the MeSH terms: Autophagy, genomewide association studies, candidate gene studies, and knockout experimental models. We identified a total of 30 genes related to autophagy and associated pathways. We further included the genes in the pathway analysis done by QIAGEN's Ingenuity Pathway Analysis Software (IPA, http://www.ingenuity.com, Fig 1) to determine genes enriched for autophagy. Core Analysis implemented by IPA was used to interpret the role of these genes in the context of biological processes, pathways and molecular networks. IPA uses right-tailed Fisher exact test to identify enriched canonical pathways and diseases associated to these genes. A pvalue < 0.05 was considered significant. A set of 21 genes showed significant enrichment for autophagy (p-value = 1.33×10 −40 ) (Table 1). We studied genetic variants within the genes and those located within 10 MB upstream and downstream of the gene. We excluded SNPs with a reported minor allele frequency (MAF) less than 1% in NCBI.

Genetic association
For each gene, we studied the association of the genetic variants intermediate vascular traits and cardiovascular outcomes. We extracted the summary statistics for the association of the selected SNPs with traits and diseases from the most recent GWAS meta-analysis of 12 traits: fasting glucose, fasting insulin, fasting proinsulin [15][16][17][18], type 2 diabetes [19], systolic blood pressure (SBP), diastolic blood pressure (DBP) [20], total cholesterol (TC), triglycerides (TG), HDL cholesterol (high density lipoprotein), LDL cholesterol (low density lipoprotein) [21] and coronary artery disease (CAD) [22]. Details of the consortia included are provided in S1 Note and number of subjects included in each study is shown in Table 2.

Exonic variants
We studied the association of exonic variants and DNA methylation status in ATG with lipid, glycemic and intermediate vascular traits using data from the Rotterdam Study. The Rotterdam Study is an ongoing prospective, population-based cohort study started in 1990 including subjects from Ommoord district in the city of Rotterdam (The Netherlands). The objectives and details of the study are described elsewhere [24]. The Rotterdam Study comprises a total of 14,926 subjects aged 45 years and over who were recruited in three sub cohorts in 1989-1993, 2000-2001, and 2006 [24]. Genomic DNA was extracted from peripheral blood mononuclear cells. Whole-Exome sequencing (WES) was performed at the Rotterdam Genomics Autophagy activation is controlled by several factors, including amino acid and nutrient deprivation, which stimulates mTORrelated pathways. Upon mTOR inhibition by starvation or rapamycin, ATG13/ULK1/FIP200 form a stable complex, essential for autophagy activity. Highly coordinated autophagy-related proteins are recruited to the phagosphore to participate in the autophagosome formation. ATG7 is one of the crucial regulators of this process. Lipidation of LC3-I known as LC3-II is attached to the autophagosome membrane and is widely used to monitor autophagy induction. When autophagosome formation is completed, the fusion between autophagosome and lysosome is mediated by lysosomal membrane protein LAMP-2. ATG12 conjugation to ATG5 is a crucial step in autophagosome formation. ATG7 regulates ubiquitin-like reaction that mediates this mechanism and subsequent processes in the autophagy machinery. Core in Erasmus Medical Centre. WES was conducted on 2,998 paired-end sequenced samples using the Illumina hiSeq2000 (2×100bp reads). Indels and single nucleotide variants were filtered out and evaluated using GAKs Variant Evaluation; variants with a call rate <0.97, duplicate samples, duplicate variants and >5% of missing genotypes were considered as exclusion criteria. A total of 2628 samples passed through all technical quality control and GATKs Haplotype Caller was used to call SNPs and indels simultaneously. Annovar software tool was used to functionally annotate each genetic variant. Each variant was coded as 0, 1 and 2 representing two reference alleles, one reference allele and one mutated allele and two mutated alleles, respectively. Based on the functional annotation, we identified a total of 13 loss of function (LoF) variants defined as changes in DNA sequence predicted to completely disrupt the formation and/or function of a protein.

DNA methylation data collection and normalization
DNA was extracted from whole peripheral blood (stored in EDTA tubes) by standardized salting out methods. The Infinium HumanMethylation450 BeadChip array was employed to determine the gene methylation status using DNA from whole blood samples of 468 individuals from third visit of the second sub-cohort, 731 individuals from first visit of the third subcohort, and 251 individuals from second visit of the third sub-cohort (no-overlap

Clinical traits and outcomes
Lipid measurements were carried out using venous blood samples obtained from all participants of the Rotterdam Study. Total cholesterol, high-density lipoprotein cholesterol, and triglycerides were measured on the COBAS 8000 Modular Analyzer (Roche Diagnostics GmbH). Low density lipoprotein cholesterol levels were estimated indirectly from measurements of TC, HDL, and TG by means of the Friedewald equation [25]. The corresponding interassay coefficients of variation was <2.1%. During each visit, blood pressure was measured twice in the right arm using a random-zero sphygmomanometer, in sitting position, after a resting period of 5 minutes. After 2006 Omron M6 Comfort and Omron M7 devices were used. The measurements were taken in duplicate. The average of the 2 measurements was used in the analyses. A qualified physician at the research center collected data on indication for use of BP-lowering medication during the interview.
Glucose and insulin were measured from venous blood samples. Fasting insulin and glucose were measured on the COBAS 8000 Modular Analyzer (Roche Diagnostics GmbH). The interassay coefficients of variations are <8% and <1.4% for insulin and glucose respectively. Type 2 diabetes was defined according to recent WHO guidelines, as a fasting blood glucose � 7.0 mmol/L, a non-fasting blood glucose � 11.1 mmol/L (when fasting samples were absent), or the use of blood glucose lowering medication [26]. Information regarding the use of diabetes medication was collected from both structured home interviews and linkage to pharmacy records [27].

Statistical analysis
SNP pruning and Genome-wide association studies look up. We found 5398 SNPs in the 21 identified autophagy genes. Linkage disequilibrium (LD)-based SNP pruning implemented in PLINK software was applied using a genetic correlation threshold of 0.5 to calculate the number of independent SNP per gene [28]. Using LD pruning, we identified 604 independent SNPs in 21 genes directly involved in autophagy pathway. We applied Bonferroni correction for 604 SNPs and 4 traits. The study-wide significance threshold was set at 2.07×10 −5 (0.05/604 SNPs×4 traits).
Exome sequencing. We used Seqmeta, an R package to perform region-based tests of rare DNA variants, to conduct multiple regression models in order to determine association between LoF variants and intermediate vascular traits including age and sex as a covariates. The package allows to obtain the scores and MAF for each gene with the function prepScores and to calculate both gene-based and single variant-based associations based on the prepScores output [29,30]. The significance level was based on Bonferroni correction to adjust for multiple comparisons according to the number of variants and traits being tested. Here, we set a pvalue threshold as 9.62×10 −4 (0.05/(4 traits×13 LoF)) Epigenome-wide association studies. We conducted epigenome wide association studies (EWAs) on blood pressure, glucose, insulin, HDL, LDL, triglycerides and total cholesterol. Beta-values to quantify the DNA methylation levels were used in a linear mixed-effect regression model with the outcomes of interest as a dependent variable. We fit the primary regression model adjusting by age and sex, cell counts, technical covariates and smoking. Further adjustments including, hypertensive or lipid lowering treatment, prevalent diabetes mellitus were employed for blood pressure and lipids levels respectively. We combined the results of the three sub-cohorts by conducting a fixed effect meta-analysis using the inverse variance method implemented in METAL [31]. We established the target set of ATG-CpGs considering the CpG islands found in upstream, downstream and promoter regions of each gene. CpG sites per gene were obtained using the UCSC genome browser (http://genome.ucsc.edu/). Statistically significant association was set at 5.63×10 −5 based on Bonferroni correction for 296 CpG sites and three trait groups (glycemic traits, lipids and blood pressure). Next to the identification of epigenetic marks in ATG genes, we extended our search of methylation quantitative trait loci (meQTLs) by the identification of both cis and trans associated SNPs using a large meQTLs dataset from mQTLdb resource (http://www.mqtldb.org/) [32].

Genetic association
The results of genetic association of the SNPs and cardiometabolic traits are presented in S1 Table. The lead significant variants are presented at Table 3. Six SNPs at three genes, AMBRA1 (Lead SNP: rs11038913), ATG13 (rs8914), and ATG16L1 (rs4944804) passed the significant threshold for association with pro-insulin serum levels. Total cholesterol level was significantly associated with 34 SNPs at ATG4D (rs10439163) and one SNP at ATG4C (rs6587988). HDLcholesterol was associated with 41 variants at ATG7 (rs7635838). One SNP at ATG4D (rs10439163) was associated with LDL-cholesterol and the one at ATG4C (rs6587988) was further associated with triglyceride levels. Moreover, 56 SNPs and 38 SNPs located in ATG7 were associated with SBP and pulse pressure, respectively. Two SNPs (rs2424994 and rs6088521) located at MAP1LC3A showed an association with CAD.
In addition to this significant association, we also observed borderline significant association between SNPs at ATG7 with SBP, ATG4B with fasting proinsulin, and MAP1LC3A with HDL-cholesterol. No associations were found between ATG and type 2 diabetes.

Loss of function variants
We used data from 2,628 participants of the first sub-cohort of the Rotterdam Study. Baseline characteristics of the individuals are presented in S2 Table. We identified 13 LoF stop-gain or splicing mutations in ATG4C, ULK1, GABARAPL2, GABARAP, ATG4D, ATG3, ATG10 and ATG5 (Table 4). We used single variant-based analysis implemented in SeqMeta to assess the association of LoF mutations with intermediate vascular traits. We identified one exonic variant in ATG3 associated with DBP and fasting glucose (p-value for SBP: 0.054). The C/A variant found in this gene is considered a nonsense mutation and is related with stop-gain function of the protein. We also found a nominally significant association between an exonic variant in ULK1 and triglycerides (p-value = 0.0012). The results of association between these variants and studied traits and diseases are mentioned in S3 Table. DNA methylation at autophagy genes In total, 1450 individuals with a mean±SD age of 60.6 ± 5.3 years were included in the analysis. The association of the studied CpGs and all traits and diseases are presented in S4 Table. In a fully adjusted model, we found differential methylation at ATG4B and ULK1 were associated with blood pressure (Table 5). Hypermethylation in these CpGs were associated with lower SBP/DBP. On the other hand, no significant associations between DNA methylation and glycemic and lipid traits were found in the study population. Using mQTLdb database, we identified five trans-meQTLs associated with cg02710553. The variants were located in an intergenic region close to MAP3K7. Furthermore, 23 trans-meQTLs associated cg06006530 were identified close to CDH18 gene. meQTLs of cg02710553 were found to be associated with CAD at nominal significance level (p-value = 0.02). Association between autophagy and cardiometabolic disorders: The Rotterdam Study

Discussion
In this study we used genetic and epi-genetic factors to assess the role of autophagy genes in cardiometabolic traits and disorders at population level. To this end, we utilized a number of genetic variants, LoF variants and DNA methylation in autophagy related genes. All approaches indicated associations between autophagy genes and certain cardiometabolic traits, mainly blood pressure, lipid levels and proinsulin levels but not coronary artery disease. Our findings on the association between genetic variants in AMBRA1, ATG13 and ATG16L1 and proinsulin levels are in agreement with previous evidence reporting autophagy deficiency as an important determinant in the pathogenesis of insulin resistance and diabetes [33]. AMBRA1 participates in the activation of beclin-1-regulated autophagy and favors the autophagosome core complex with the participation of other ATG proteins such as ATG13 and ATG16L1 [34]. The upstream autophagy-signaling network controlled by AMBRA1 is crucial for the metabolic response to a vast number of stress stimuli, ranging from starvation to hypoxia or DNA damage. Dysfunctional autophagy has been suggested as a key process related to an impaired proinsulin/ insulin homeostasis observed in pancreatic beta cells from experimental models [35]. Conditional knockout of ATG7 in high-fat-fed C57BL/6 mice has resulted in declined insulin secretion, impaired glucose tolerance and degeneration of pancreatic islets [36]. In autophagydeficient cells, the insulin secretion is restrained enabling the proinsulin accumulation in secretory granules and increased secretion in response to stimuli [37]. At population level, in line with our findings, an exome array analysis conducted in 8229 individuals showed association between genetic variants in AMBRA1 and ATG13 with fasting proinsulin concentrations [38]. The role of impaired autophagy in proinsulin degradation has highlighted the importance of the autophagy pathway on the design of novel therapeutic strategies, aimed to manipulate proinsulin clearance as means to increase the insulin secretion in diabetic population.
Our findings on common genetic variants in ATG and blood lipids support the hypothesis that autophagy may play a role in dyslipidemia. The role of impaired autophagy and lipid metabolism has been established from the identification of significant increased levels of hepatic triglycerides and cholesterol content in hepatocyte-specific knockout mice of ATG5 [39]. Moreover ATG7 knockout mice models have displayed severe morphological abnormalities in the structure of white adipocytes, as well as an aggravated insulin resistance with increased lipid content and inflammatory changes [40,41]. In summary, upregulation of autophagy leads to a decrease of triglycerides and cholesterol in plasma, reduced lipid store as well as LDL oxidation and free fatty acid B-oxidation and an increase of folding and traffic proteins [42]. Further investigations are needed to clarify this association and potential pathways linking autophagy with blood lipids levels. Recent findings of a critical role for macroautophagy in the metabolism and storage of cellular lipids have now suggested that alterations in autophagy may mediate human disorders marked by excessive cellular lipid stores. In contrast, mice with endothelial specific deletion of ATG7 have shown normal blood pressure and normal vessel architecture compared to wild types [43]. We further identified a suggestive association between LoF variants annotated in ATG3 and glucose levels. It has been shown that glycogen autophagy in newborns serves as a mechanism of glucose homeostasis [44]. In adult animals, the administration of glycogen autophagy-inhibiting insulin triggers a reduced rate of breakdown of liver glycogen by autophagy [45]. Impaired mechanisms related to autophagosomal glucose production and the influence of gluconeogenesis may lead to a dissociation of gluconeogenic glucose production from blood glucose levels [46]. On the other hand, we found no association between changes in DNA methylation patterns and glycemic traits. This might be explained by the fact that methylation is a tissue-specific process. Therefore, SNPs could be operating independently from methylation patterns.
We studied DNA methylation at autophagy-related genes and found differential patterns at ULK1 and ATG4B associated with blood pressure levels. Increased blood pressure is determined by a complex machinery regulated, among others, by the renin-angiotensin system (RAS) [47]. The interaction between autophagy and RAS has been previously examined. Porrello et al, provided the first evidence for an interplay between these two mechanisms in cardiomyocytes. This study reported that rat cardiomyocytes developed and augmented (via the angiotensin II receptor type 1) or inhibited (via the angiotensin II receptor type 2) autophagic response on stimulation by angiotensin II [48]. The association between autophagy and blood pressure found in our study could be explained by its interaction with RAS and its role in hypertension. Given the fact that autophagy and RAS are both involved in many pathophysiological processes, further investigation is warranted to better understand the molecular mechanism behind this interaction and its role in pathological conditions. In addition, we found an association between DNA methylation in ATG4D and HDL levels. ATG4D is known to participate in the delipidation of GABARAP-L1, whereas the silencing of ATG4D abrogates GABAR-AP-L1 autophagosome formation [49]. Experimental evidence of the role of this gene on lipid transport/metabolism is currently lacking.
Epi-genetic associations are subject to confounding and reverse causation. One approach to overcome these biases in molecular epidemiology is to use genetic instruments. In this study, we found that meQTL at cg02710553 was associated with CAD. This might indicate that the association of cg2710553 and CAD is causal.
We used both genomics and epi-genomics data to examine the role of autophagy related genes in cardiometabolic traits and diseases. Although we have several significant association, we did not find a consistent association with a certain gene or a certain trait across all approaches. It should be noted that genomic and epi-genomic approaches are not always pointing to the same gene. Moreover, it is not a surprise if a different gene in a certain pathway has more or less importance for a certain trait. The high complexity of the autophagy mechanism may be a major contributor of the heterogeneity of the findings in our results [50].
Our study has several strengths. First, we used data from largest GWAS, which has provided the best statistical power that could be achieved. Second, we used both genetic and epi-genetic approaches towards the research question. Third, we used exome sequence data, which in contrast to GWAS, is not based on a studying a proxy variant. Our study also had several limitations. First, the sample size for exome sequence analysis was much smaller than the GWAS, thus the statistical power was significantly lower. Loss of function variants are mostly found in low frequency, indicating that they are enriched for mildly deleterious polymorphisms suppressed by negative natural selection [51]. Therefore, smaller populations might also provide enough statistical power to detect their effect, however, this could not be ruled out that lack of association in our analysis might be a results of small sample size.
Second, DNA methylation pattern is tissue/cell line specific. As it is common in epidemiologic studies we used DNA methylation in whole blood, a cell type mixture, which might not be relevant for some of the traits. In our study DNA methylation was captured from the assessment of leukocytes. Although leukocytes are not relevant tissue for some of the traits that we have assessed in our analysis, they are the main tissue available in large scale in epidemiologic studies and evidence has demonstrated that methylation patterns might correlate between blood and the relevant tissues, suggesting that the use of blood tissue could yet be informative [52][53][54][55]. However, such findings should be validated in subsequent studies.
In conclusion, this study is the first to examine the role of autophagy-related genes in intermediate vascular traits using a population-based approach. We have characterized the role of common and rare genetic variants as well as epigenetic variations of autophagy-related genes on several traits. Despite the heterogeneity of our findings across approaches employed and traits evaluated, we found many associations between autophagy genes and cardiometabolic traits and diseases. The integral approach covered by this study could contribute to further analysis evaluating the role of autophagy in other human diseases and traits, as well as the design of experimental studies targeting other autophagy-related genes and/or associated pathways.
Supporting information S1 Note. Brief description of the consortia included in this study. (DOCX) S1