Novel Genetic Locus Influencing Retinal Venular Tortuosity Is Also Associated With Risk of Coronary Artery Disease

Supplemental Digital Content is available in the text.

(ODradius), and establishment of the standard retina coordinates and zones (Figure I), the 6 thickest arterioles and 6 thickest venules appearing in a zone extending from the optic disc boundary to 2 optic disc diameters out was sampled. This was used to calculate the median (TortA) and maximum (TortAmax) arteriolar tortuosity and the median (TortV) and maximum (TortVmax) venular tortuosity. Central retinal arteriolar equivalent (CRAE), central retinal vein equivalent (CRVE) and Arteriole-to-Venule ratio (AVR) qualify vessel calibers and were measured in a zone 2 to 3 optic disc radii from the center of the optic disc. Tortuosity is a quantitative measure of how much blood vessels twist and turn. Smaller values indicate straighter vessels. It has been implicated with various diseases, e.g. retinopathy of prematurity and other systemic diseases 6,7 . Several algorithms have been proposed to quantify tortuosity consistently with clinical grading; the algorithm used in VAMPIRE 3.1 was chosen based on the study by Lisowska et al. 8 , which compared five classes of algorithms and provides an introduction to the topic. Results indicated that curvature-based algorithms performed best over a range of image characteristics, assuming high sampling rates of the vessel centre lines are available, which is the case in our work. Briefly, given the centre line of a vessel (defined as a segment of vein or artery between two consecutive junctions or the end of a retinal zone) the algorithm computes the total squared curvature divided by the segment's arc-length in pixels. Robust curvature estimates are provided by Annunziata's algorithm 9 . Tortuosity is dimensionless and therefore does not have units; the range of its values and scale depends on the algorithm adopted for its calculation. For comparisons between studies where different algorithms may have been used relative values (e.g. rank, as opposed to absolute values) or suitable normalizations are performed. In our case, venular tortuosity ranges from 2.51×10 -6 to 1.27×10 -3 and arterial tortuosity ranges from 1.45×10 -6 to 1.55×10 -3 . Figure Ib. shows an example of high venular and low arteriolar tortuosity vessels, with corresponding values computed by the VAMPIRE algorithm. Replication Cohorts. Retinal fundus images using digital fundus camera from 1091 individuals from LBC1936 were collected at the recruitment stage and three years later, retinal traits were measured at a subsequent wave of testing using SIVA v3.1 10 , at a mean age of 72.5 years (SD 0.7). Similarly, 897 and 976 individual's retinal fundus images centred between the macula and optic disc from Croatia-Korčula and Croatia-Split cohorts were collected using digital fundus camera and retinal traits were quantified using SIVA v3.1 10 . SIVA is a semi-automated software which can be used to measure the retinal vascular parameters including retinal vascular tortuosity and vascular caliber from retinal fundus images. After automatic detection of the optic disc, it placed a grid with reference to the center of the optic disc. Then the tortuous vessels were identified and tortuosity traits including TortA, and TortV were measured using the standard grading protocol by the software; this process was monitored by trained graders and adjusted manually if necessary.
Genotyping, quality control and imputation Discovery cohorts. Quality control (QC) assessment were performed using PLINK1.9 11 . The poor quality variants, samples were excluded based on the quality control (QC) criteria included the following: SNPs call rate < 95%, Hardy-Weinberg equilibrium (HWE) P value < 10 -6 , sample call rate < 95%, sample relatedness (IBD >0. 8), and mismatch between reported and genotypic gender information. QC'd genotype data were imputed using IMPUTE2 12 on the basis of 1000 Genome Projects reference panel for all population. Finally, ancestry information of the individuals was derived using EIGENSTRAT 13 and first three principal components (PCs) were used for the association analyses to adjust the population stratification. ORCADES samples were genotyped with either the Illumina HumanHap300 bead chip (n=890) or the Illumina Omni1M (n=304) or Illumina Omni Express bead chips (n=1073). Alleles were called in Bead Studio/Genome Studio (Hap300/Omni) using Illumina cluster files. Subjects were excluded if they fulfilled any of the following criteria: genotypic call rate <98%, mismatch between reported and genotypic sex, unexpectedly low genomic sharing with first or second degree relatives, excess autosomal heterozygosity, and outliers identified by IBS clustering analysis. We excluded SNPs on the basis of minor allele frequency (<0.01/monomorphism), HWE (P<10 -6 ) and call rate (<97%). Given the very high overlap in SNPs between the two Omni chips, the intersection of QC'd SNPs was used to impute and phase individuals' genotyped on the Omni arrays together, whilst the Hap300 individuals were phased using SHAPEIT version 2 14 and imputed, separately. Imputation was carried out using IMPUTE2 12 and the 1,000 genomes reference panel. All ancestries phase1 integrated v3 reference panel, with a secondary reference panel of local exome sequences, sequenced using the Agilent Sure Select All Exon Kit v2.0. Illumina 100 bp paired end reads (average 30x depth), derived from 90 ORCADES subjects were chosen to optimally represent the haplotypes present. Imputations for the Hap300 and Omni subjects were then combined to form a combined panel of 37.5m SNPs for 2222 subjects 15 . Replication Cohorts. For LBC1936 individuals were excluded based on unresolved gender discrepancy, relatedness, call rate (≤0.95), and evidence of non-Caucasian descent. SNPs were included if they met the following conditions: call rate ≥ 0.98, minor allele frequency ≥ 0.01, and HWE test with P ≥ 0.001. Imputation to the 1000 Genomes Phase I v3 (March 2012 release) reference set was performed using minimac software 16 . For Croatia-Korčula and Croatia-Split samples and markers were excluded based on the following QC metrics; SNPs call rate < 98%, HWE with P value < 10 -6 , sample call rate < 97%, MAF < 1%, outliers identified by IBS clustering analysis and unresolved gender discrepancy. Pre-phasing was performed using SHAPEIT version 2 14 . Imputation was carried out using IMPUTE2 12 software and 1000G Phase I v3 (March 14, 2012) reference panel.

Power calculation
The statistical power of detecting SNPs associations with the quantitative traits in two stage GWAS was calculated using the GWASPower/QT 17 .

In-silico functional annotation
Top SNPs were queried in the HaploReg v4.1 18 to catalogue the all SNPs near noncoding variants with r 2 > 0.8, and RegulomeDB 19 and GWAS catalog 20 databases used to explore the known and predicted regulatory elements and relevant genetic association studies. Functional effects of the top genes were predicted using the Encyclopedia of DNA Elements (ENCODE) project 21 and Roadmap Epigenomics projects, HaploReg, UCSC Genome Browser 22 , and RegulomeDB. We used the expression Quantitative Trait Loci (eQTL) browser database in Genotype-Tissue Expression (GTEx) 23 to examine the cis-eQTLs for the top retinal traits associated SNPs mapped to the gene within the genomic region. Co-localization analysis was performed using eCAVIAR to investigate whether the same significant variant is causal in both GWAS of retinal venular tortuosity and gene expression in the relevant tissues and to identify causal gene 24 . We tested genes with significant cis-eQTL association (ACTN4 and CAPN12) for the GWAS locus associated with TortV by analysing the lead variants using summary data from the present study and eQTL from GTEx v6. The co-localization posterior probability (CLPP) score was estimated for each variant in a GWAS locus for a given gene. Gene Visible web database from Genevestigator, open-access version which integrates manually curated gene expression data from microarray and RNAseq experiments. It was used to find the expression level of the genes in human tissues, cell lines, cancers or perturbations, associated with tortuosity traits. It is a large and growing reference database which integrates systematically annotated and manually curated gene expression data from high-quality microarray data (Affymetrix array) from several organisms 25 .

In-silico look-ups of the novel variants for clinical outcomes
We performed in-silico look-ups of variants of interest for cardiovascular related outcomes including coronary artery disease, myocardial infarction, hypertension, heart rate, HDLC, LDLC, Atrial fibrillation (AF) and triglycerides using summary results from different largescale genome-wide studies on cardiovascular risk factors. The Coronary Artery Disease (CARDIoGRAMplus C4D) consortium 26 (CARDIoGRAMplusC4D) 1000 Genomes-based meta-analysis data comprised of 60,801 CAD cases and 123,504 controls from European, South Asian, and East Asian descent. In the Global Lipid Genetics Consortium, genetic data from 188,577 individuals of European, East Asian, South Asian, and African ancestry were used to examine the genetic loci associated with blood lipids levels. The International consortium for blood pressure 27 (ICBP) GWAS investigated the genetic loci associated with systolic and diastolic blood pressure traits in 200,000 individuals of European descent. A recent large-scale meta-analysis of GWAS includes 65,446 participants of multi-ethnic cohorts with AF identified several variants associated with AF.
UK Biobank data comprised of 112,008 participants who had a measure of pulse rate at the main interview and had genotype data 28 . We extracted the imputed genotypes for these SNPs from the interim release data set of the UK Biobank and performed multiple linear regressions including covariates of age, gender, and the first ten principal components obtained using EIGENSTRAT. , and non-tortuous arteries (2.41e -5 ), with values computed by the algorithm; Zone B represents the ring 0.5-1 optic disc parameters away from the centre; Zone C represents the ring extending from optic disc boundary to 2 optic disc diameters away.

Figure VII. eQTL annotation of top GWAS variants for quantitative retinal vessel traits. Y-axis represents tissue-specific gene expression data which is normalized by rank normalization method while x-axis shows the GWAS lead variant genotypes.
A) TortA associated SNP, rs7991229 is correlated with COL4A2 expression in heart left ventricle; B) TortA associated SNP, rs9515212 is correlated with COL4A2 expression in heart left ventricle tissue; C) TortV associated SNP, rs1808382 is correlated with CAPN12 expression in artery aorta; D) TortV associated SNP, rs1808382 is correlated with ACTN4 expression in artery aorta.                       Table XIII. Co-localization of eQTLs with TortV-associated GWAS SNPs in heart and blood vessel tissues using eCAVIAR.  TortV associated SNPs are strong LD (r 2 >=0.8) with rs11083475 (Bold text), previously reported SNP for elevated resting heart rate.