Common and rare genetic markers of lipid variation in subjects with type 2 diabetes from the ACCORD clinical trial

Study population

The ACCORD trial (NCT00000620; https://clinicaltrials.gov/ct2/show/NCT00000620) had a double 2 × 2 factorial design, consisting of 10,251 recruited subjects with type 2 diabetes mellitus and either a history of cardiovascular disease (CVD) or at least two known risk factors for cardiovascular disease, such as documented atherosclerosis, albuminuria, dyslipidemia, hypertension, smoking, or obesity (Buse, 2007). Subjects were randomized to either intensive or standard glycemia treatment strategies (targeting HbA1c < 6.0 vs. HbA1c between 7.0 and 7.9). A subset of 4,733 subjects were further randomized to intensive versus standard blood pressure management (targeting systolic blood pressure of <120 mm Hg versus <140), and the remaining 5,518 subjects were randomized to intensive versus standard lipid management (fenofibrate versus placebo, with all subjects on simvastatin) (ACCORD Study Group, 2010b). The age range for subjects with a history of CVD was 40–79, and for those with no prior CVD history, 55–79. Body mass index (BMI) was limited to a maximum of 45, and serum creatinine to 1.5 mg per deciliter. Median length of follow-up was 4.7 years, and the primary outcome was the first occurrence of nonfatal myocardial infarction or nonfatal stroke or death from cardiovascular causes.

Participants in ACCORD were given an option to provide a blood sample for future genetic studies, and over 80% of participants agreed to do so. After genotyping and quality control (see Supplemental Information 1) of DNA extracted from these samples, the population for the current study included 7,844 subjects. A comparison of baseline characteristics of this subset of ACCORD to the entire ACCORD patient population is given in Table 1. The subjects included in the current study were similar to the overall ACCORD participants regarding all baseline demographic and clinical factors, and in the use of lipid lowering drugs prior to the trial. Phenotypes consisted of baseline measurements for TC (mg/dL), LDL (mg/dL), HDL (mg/dL) and TG (mg/dL). Data for TC, LDL, HDL, and TG were extracted from the publicly released data files and then log transformed as shown in Fig. 1 to meet parametric assumptions of the statistical models described below. The Protocol Review Committee, appointed by the National Heart, Lung, and Blood Institute (NHLBI), approved the study protocol. Each ACCORD participant provided written informed consent using procedures reviewed and approved by each clinical site’s local institutional review board and based on a template provided by the study group that was approved and subsequently centrally monitored by the Coordinating Center and the NHLBI (IRB: FWA00003429). The portion of the informed consent document describing the genetics component of ACCORD uses the multilevel approach recommended by the NHLBI (National Heart Institute et al., 1997).

Genotyping

DNA extracted from stored white blood cells from ACCORD subjects who consented to genetic research were obtained from the ACCORD Central Laboratory at the University of Washington. DNA were quantitated using PicoGreen, and 200 ng of each DNA was aliquoted into 96-well plates. Three duplicate ACCORD samples and two DNAs from HapMap trios (obtained from Coriell Cell Repositories, Camden, NJ) were included as controls on each plate. Amplification, fragmentation, and hybridization of DNA to Axiom Biobank 1 genotyping arrays (Affymetrix, Santa Clara, CA) were carried out using Axiom Reagent Kits, and hybridization, ligation, washing, staining and scanning of the arrays were carried out on the GeneTitan MC instrument (Affymetrix). Initial plate QC was performed using Affymetrix Genotyping Console Software and genotype calling was done using Affymetrix Power Tools (v1.15.0) with the Axiom GT1 algorithm, which is a modified version of the BRLMM-P algorithm that adapts generic prior cluster positions to the data using an EM algorithm.

The Axiom Biobank Genotyping Array from Affymetrix includes ∼246,000 common genetic variants selected to optimize genome-wide coverage for association in addition to ∼335,000 variants focused on the coding portion of the genome, and including non-synonymous, frameshift, indel and splice variants as well as other known loss-of-function and disease causing mutations. Most of this second group of variants was identified by exome sequencing in up to 26,000 individuals, and the majority of them are uncommon or rare. Additional variants to cover known eQTLs and pharmacogenomics markers are also included on the array. A total of 628,679 probes were genotyped. After genotype quality control, genotypes were available for 583,613 variants, of which 89,212 were monomorphic. After imputation based on 1,000 genomes haplotype data, an additional 26,862,499 imputed variants with an “info” metric > 0.5 were retained for association testing. The genotype information is available in dbGaP (http://www.ncbi.nlm.nih.gov/dbgap).

Data processing

Quality control

Genotypes were subjected to rigorous quality control based on genotyping quality metrics, duplicate concordance, Mendelian segregation (in HapMap trios included on the genotyping plates), Hardy–Weinberg Equilibrium, and predicted gender (see Supplemental Material). Cryptic relatedness was identified using KING (v1.4), and one member of each pair with a kinship coefficient >  ${\left(\frac{1}{2}\right)}^{\frac{5}{2}}=0.1768$ was removed from the analysis data set (random seed: 1,485) (Manichaikul et al., 2010). Principal components (PCs) based on the genotype data were computed using EIGENSTRAT (v4.2) and were used to control population stratification (Price et al., 2006). The first two PCs are shown in Fig. 2, where the marker colors and labels represent the self-reported ethnic background for each sample. Genotype imputation was accomplished using a two-step approach where the genotype calls were first pre-phased using SHAPEIT2 (v2.r778) and then imputation was conducted using IMPUTE2 (v2.3.0) (Howie, Donnelly & Marchini, 2009; Howie, Marchini & Stephens, 2011; Delaneau, Marchini & Zagury, 2012; Delaneau, Zagury & Marchini, 2013). Both steps used the 1,000 Genomes Phase1 integrated haplotypes reference panel (release date Dec 2013) from the IMPUTE2 website. Probes significantly deviating from HWE (χ² > 19.51, p-value < 10⁻⁵) in at least two of the four main ethnic subgroups were excluded from the imputation process. Values of the cleaned probe set for white samples are shown in Fig. S14.

All samples were pre-phased together since this approach was previously found to improve phasing in a cohort with diverse ancestries compared to phasing groups separately (Delaneau, Zagury & Marchini, 2013). The imputation process began by splitting the pre-phased haplotype data into 5 Mbp non-over-lapping segments for each chromosome, using the same reference panel from the pre-phasing step. Only those variants with an “info” metric > 0.5 were retained for association testing, resulting in a total of 26,862,499 imputed variants (∼71.7% of total imputed variants). Additional information regarding the pre-phasing and imputation process can be found in the File S1.

Covariate selection

Here, we take a combined approach to variable selection to address potential confounding variables. Some covariates are forced into the model based on the results of previous studies or on expert knowledge related to the phenotype, while variable selection, via backwards selection, is subsequently performed on candidate variables to identify covariates specific to the ACCORD dataset. Covariate names and descriptions can be found in Tables S3 and S4. A substantial proportion of the cohort was taking lipid lowering medications at the time baseline lipid measurments were taken (e.g., 63% were on a statin). Lipid lowering concomitant medications were forced into the model to prevent confounding of the baseline lipid measurements. The percentage of the cohort on each of these lipid lowering medications is in Table 1 and additional information is in Table S3. Years since diabetes diagnosis and years of hyperlipidemia diagnosis were mean imputed due to large numbers of patients missing those records. Only samples with complete phenotype and covariate data were retained. A correlation matrix for all covariates is shown in Fig. S18. Three pairs of covariates were flagged for high correlation, and one of each pair was removed prior to analysis: glomerular filtration rate and serum creatinine were negatively correlated (r =  − 0.77), systolic and diastolic blood pressure were positively correlated (r = 0.53) and BMI and waist circumference were positively correlated (r = 0.63). See Supplemental Material for additional information regarding covariate selection.

Heritability approximation

The phenotypic variation explained by genome-wide variants is estimated using the software tool genome-wide complex trait analysis (GCTA; v1.22) (Yang et al., 2011). An estimate of a phenotype’s heritability can help to determine if the results from association tests seem plausible. Additional information regarding the heritability approximation is in the Supplemental Material.

Common variant analysis

Association between a phenotype and single common variant is tested using the linear regression model ${\displaystyle y={\alpha }_0+X\alpha +{\beta }_{g}g+ϵ,}$ where y is the phenotype, α₀ is the intercept, X are the covariates, α are the covariate regression parameters, β_g is the regression parameter for the variant, g is the additively coded genotype and ϵ is the error term. Genotyped variants are tested using PLINK, where $g_i\in \left\{0,1,2\right\}$ is the number of minor alleles for the ith individual. Imputed variants are tested using a linear regression model in the statistical programming language, R, where $g_i=p_i\left(Aa\right)+2p_i\left(aa\right)$ is the dosage score computed from the posterior probabilities for genotypes Aa and aa (R Development Core Team, 2014).

The resulting test statistics from the common variant analysis were adjusted for genomic inflation (Devlin & Roeder, 1999). Genomic inflation values were used as a guide for selecting the MAF threshold of 3% to separate common and rare variant analyses. Additional information about the MAF selection criteria can be found in the Supplemental Material. The results from the common variant tests were considered statistically significant based on a p < 1 × 10⁻⁸.

Rare variant analysis

Current approaches can be commonly divided into burden and non-burden approaches. Burden tests collapse a set of rare variants from a region of interest (e.g., a gene) into a single variable, which is then tested for association with a phenotype. However, a major limitation of simple burden tests is that they cannot account for the possible direction (positive or negative association) of a rare variant effect (Wu et al., 2011). One non-burden rare variant test that allows for different directions and magnitudes of effects for each variant is the sequence kernel association test (SKAT) (Wu et al., 2011). The approach used here is to apply a suite of tests that are representative of commonly used methods in the literature that cover both burden and non-burden approaches.

Simple collapsing methods use indicator, proportion and weighted approaches chosen a priori that can easily be computed from genotype data. More sophisticated approaches can incorporate variant-specific information and variable thresholds for collapsing. Three burden tests are used here and are based on simple collapsing approaches. These methods first create a collapse score, c, and then test for association using the linear regression model ${\displaystyle y={\alpha }_0+X\alpha +{\beta }_{c}c+ϵ,}$ where y is the phenotype, α₀ is the intercept, X are covariates, α are the covariate regression parameters, β_c is the regression parameter for the collapse score and ϵ is the error term. Association p-values are computed under the null hypothesis H₀:β_c = 0. The selected burden tests differ in how c is computed.

The first two tests are based on RVT1 and RVT2 originally proposed by Morris & Zeggini (2010). We adapted a slight change to RVT1 where $c_i=\left({\Sigma }_j{g}_{ij}\right)∕2n_i$, which counts the total number of rare alleles rather than the number of variants with rare alleles. In RVT2, c_i is simply an indicator function $c_i=I\left(r_i>0\right)$. The third test uses a weighting scheme similar to the one proposed in the weighted sum statistic (WSS) case-control method (Madsen & Browning, 2009). While burden tests use simple linear regression, SKAT uses a linear mixed model and tests for association with a variance-component score test (Wu et al., 2011). While SKAT has been found to be more powerful than burden tests when the variants have different directions of effect, it is less powerful when all variant effects are in the same direction. The balance between SKAT and burden tests was addressed by the optimal test, SKAT-O, where a combination of the two approaches is optimized (Lee, Wu & Lin, 2012). Additional information regarding the rare variant testing implemented here is in the Supplemental Material.

Here we use both SKAT and SKAT-O in addition to the three burden tests mentioned above. Only those variants <3% MAF were included in the rare variant approaches. The non-burden approaches were implemented using the R package “SKAT” (v0.95) (Lee, Miropolsky & Wu, 2013).

The use of multiple rare variant tests for each gene compounds the problem of multiple-hypothesis testing. Application of a correction for multiple testing, such as a false discovery rate (FDR) (Storey, 2002) approach is not straightforward due to each gene having a set of p values, one for each rare variant test. This complication can be resolved by combining the set of p-values into a single p-value for each gene, while accounting for the dependence of p-values using the “correlated Lancaster procedure” described by Dai et al. to combine the p-values from each test into a single p-value for each gene (Dai, Leeder & Cui, 2014). Combining the set of p-values for each rare variant test into a single p-value for each gene allows straightforward application of FDR controlling procedures. Computing q-values was done using the R package, q value (v1.36.0) (Dabney, Storey & Warnes, 2004). Additional information about combining the rare-variant p-values is in the Supplemental Material.

Common variant analysis

All common genotyped and imputed variants were tested individually for association using PLINK or a linear model in the statistical programing language, R, respectively. A set of LD pruned (VIF < 1.5), genotyped probes was used to compute genomic inflation values for different MAF thresholds, as shown in Fig. S20, and a MAF threshold of 3% was subsequently selected. The inflation values at MAF = 3% are 1.00, 1.00, 1.02 and 1.00 for TC, LDL, HDL and TG, respectively. A total of 292,816 genotyped and 7,812,348 imputed variants had MAF > 3% and were included in the common variant analysis. Quantile–quantile (QQ) plots for the common variants are shown in Fig. S21, along with lambda values. The location of significantly associated loci can be seen on a Manhattan plot, where the negative log transformed, genomic control adjusted p-values versus genomic location for each of the phenotypes are shown in Figs. 3A–3D for TC, LDL, HDL, and TG, respectively.

Three peaks reaching genome-wide levels of statistical significance (p < 1 × 10⁻⁸) were observed in the TC association. The lead markers in these regions are located 109.82 Mb on Chr1, 126.5 Mb on Chr8, and 45.41 Mb on Chr 19 (Fig. 3A). These SNPs are located in genes CELSR2, TRIB1, and APOE, respectively (Table 2). In the LDL common variant analysis, two significant peaks were observed at 109.82 Mb on Chr1 and 45.41 Mb on Chr 19, located in CELSR2 and APOE genes, respectively (Fig. 3B and Table 2). Five peaks represent statistically significant associations with HDL (Fig. 3C). The lead SNPs in these peaks were located 19.82 Mb on Chr 8, 58.68 Mb on Chr 15, 56.99 Mb on Chr 16, 47.17 Mb on Chr 18, and 45.41 Mb on Chr 19, although the lead SNP in the peak on Chr 19 (rs429358) is not in HWE (p < 1 × 10⁻⁵) in all populations combined. These lead SNPs are in the LPL, LIPC, CETP, LIPG, and APOE genes, respectively (Table 2). Lastly, seven peaks were observed with SNPs significantly associated with TG levels. The lead SNPs in these peaks are located 63.07 Mb on Chr1, 27.73 Mb on Chr 2, 73.02 Mb on Chr 7, 19.82 Mb and 125.58 Mb on Chr 8, 116.65 Mb on Chr 11, and 45.41Mb on Chr 19 (Fig. 3D). The SNP (rs6982502) on Chr 8 was significantly out of HWE (p < 1 × 10⁻⁵) in all populations combined. These SNPs are located in the ANGPTL3, GCKR, MLXIPL, LPL, TRIB1, ZNF259, and APOE genes, respectively (Table 2).

Rare variant analysis

All functional variants (annotated based on GRCh37p13 build and retained those with accession numbers beginning with “NM” and having codes for stop-gain, missense, stop-loss, frameshift, cds-indel, splice-3 and splice-5) and a MAF < 3% were considered for incorporation into the gene-based, rare variant analysis. After QC, a total of 16,480 genes contained at least two variants below the MAF threshold. In total, 146,689 genotyped and 73,295 imputed variants went into the analysis, with median of nine variants per gene. Mean minor allele frequencies for all the rare variants in an individual gene ranged from 6.37e–05 to 2.53e–02. The suite of three burden and two nonburden tests were applied and the resulting p-values were combined using the correlated Lancaster procedure, as described above. QQ plots for each of the phenotypes are shown in Fig. 4, where the marker colors indicate q-value significance thresholds.

A total of 11 genes were significantly associated with TC, LDL, HDL, or TG in the combined rare variant analysis (q < 0.1), with the number of rare variants in these genes ranging from two to 22 (Table 3). Variants in PCSK9 were associated with TC and LDL. LPL, CNOT2, CETP, HNF1B, ANGPTL4, and HPN-AS1 were associated with HDL, and APOC3, PAFAH1B2, ANGPTL4, SIRPD, and UBE2L3 were significantly associated with TG (Table 3).