Associations of carotid intima media thickness with gene expression in whole blood and genetically predicted gene expression across 48 tissues

Abstract

Carotid intima media thickness (cIMT) is a biomarker of subclinical atherosclerosis and a predictor of future cardiovascular events. Identifying associations between gene expression levels and cIMT may provide insight to atherosclerosis etiology. Here, we use two approaches to identify associations between mRNA levels and cIMT: differential gene expression analysis in whole blood and S-PrediXcan. We used microarrays to measure genome-wide whole blood mRNA levels of 5647 European individuals from four studies. We examined the association of mRNA levels with cIMT adjusted for various potential confounders. Significant associations were tested for replication in three studies totaling 3943 participants. Next, we applied S-PrediXcan to summary statistics from a cIMT genome-wide association study (GWAS) of 71 128 individuals to estimate the association between genetically determined mRNA levels and cIMT and replicated these analyses using S-PrediXcan on an independent GWAS on cIMT that included 22 179 individuals from the UK Biobank. mRNA levels of TNFAIP3, CEBPD and METRNL were inversely associated with cIMT, but these associations were not significant in the replication analysis. S-PrediXcan identified associations between cIMT and genetically determined mRNA levels for 36 genes, of which six were significant in the replication analysis, including TLN2, which had not been previously reported for cIMT. There was weak correlation between our results using differential gene expression analysis and S-PrediXcan. Differential expression analysis and S-PrediXcan represent complementary approaches for the discovery of associations between phenotypes and gene expression. Using these approaches, we prioritize TNFAIP3, CEBPD, METRNL and TLN2 as new candidate genes whose differential expression might modulate cIMT.

Introduction

As a marker of atherosclerosis, carotid intima media thickness (cIMT) is a predictor of coronary heart disease (CHD) and stroke (1,2). cIMT can be used to evaluate the progression of subclinical atherosclerosis before symptoms become clinically evident. Like CHD and stroke, cIMT has a moderate heritability (3–7), and numerous genetic loci for cIMT have been identified through genome-wide association studies [GWAS (8–17)]. However, the genetic variants at these loci collectively explain only a small fraction of the heritability of cIMT (8–17) and the ability of these variants to predict incident cardiovascular disease remains limited (18–21). Furthermore, most variance in cIMT is not explained by traditional cardiovascular risk factors, signaling that cIMT may represent an alternative modulator of cardiovascular disease that should be investigated (22–25). Besides genetic association studies, alternative approaches harnessing genomic data may yield new loci associated with atherosclerosis.

One such approach is transcriptome-wide differential gene expression analysis, in which phenotypes are tested for association with gene expression levels in terms of mRNA levels, as measured by microarrays or RNA sequencing, in a hypothesis-free manner. No large-scale transcriptome-wide differential expression study has been performed on cIMT, although several differential expression studies in whole blood have already identified genes whose expression is associated with cardiovascular disease (26–32). Emerging techniques for imputing the association of genetically determined gene expression values from GWAS summary statistics, such as S-PrediXcan (33), hold promise to improve discovery of gene expression-phenotype associations by studying regulation of gene expression in a broader set of tissues and detecting the sharing of regulatory mechanisms across multiple tissues (34–39). To date, neither of these two approaches has been used to study genetic associations with cIMT.

In this study, we sought to compare the findings of each approach and their implications for the etiology of carotid intimal thickening and atherosclerosis. To this end, we comprehensively examined the transcriptome for associations with cIMT, by performing an ex vivo study of directly measured gene expression in whole blood of 5647 individuals and an in silico study of predicted gene expression in 48 tissues using S-PrediXcan.

Results

Baseline characteristics

Our discovery analysis included 5647 individuals of European descent from four different cohorts. The replication analysis included a total of 3953 individuals of European descent across three cohorts. All studies measured cIMT using similar methods; detailed descriptions are given in the Supplemental Methods. Baseline characteristics of the studies included in the discovery and replication differential expression analyses are shown in Table 1. The mean age of participants across the studies ranged from 45.3 to 70.2 years. The percentage of women ranged from 35.1 to 65.0%, and the mean body-mass index (BMI) ranged from 25.2 to 29.9 kg/m². The percentage of participants with cardiovascular disease ranged from 0.6 to 39.2%.

Table 1

Baseline characteristics of the participants of the seven studies included in the discovery, and replication analyses

	Discovery				Replication
	KORA F4	LIFE-Adult	Rotterdam study	SHIP-TREND	FHS	LIFE-Heart	NESDA
Sample size	836	2973	856	982	1958	1752	243
Age	70.2 (5.3)	57.6 (12.5)	59.7 (8.0)	50.1 (13.7)	66.6 (8.9)	62.7 (10.9)	45.3 (12.6)
Sex (% women)	50.5	48.1	53.4	56.0	54.2	35.1	65.0
BMI (kg/m2)	29.0 (4.5)	27.4 (4.6)	27.7 (4.6)	27.3 (4.5)	28.4 (5.3)	29.9 (5.0)	25.2 (3.3)
HDL cholesterol (mmol/l)	1.4 (0.4)	1.6 (0.5)	1.4 (0.4)	1.5 (0.4)	1.5 (0.5)	1.3 (0.4)	1.7 (0.4)
Total cholesterol (mmol/l)	5.7 (1.0)	5.6 (1.1)	5.5 (1.1)	5.5 (1.1)	4.8 (1.0)	5.4 (1.2)	5.1 (1.0)
Lipid-lowering medication use (% yes)	24.3	15.2	27.0	7.3	49.0	35.6	NA
Systolic blood pressure (mmHg)	128.5 (19.1)	129.0(16.8)	134.5 (20.1)	124.3 (16.9)	128.8 (17.5)	138.3 (19.4)	134.0 (19.8)
Diastolic blood pressure mmHg)	73.9 (9.8)	75.5 (9.9)	82.9 (11.6)	76.5 (9.7)	73.4 (10.2)	83.8 (11.5)	81.5 (11.4)
Antihypertensive medication use (% yes)	56.8	44.5	27.3	29.3	53.6	81.3	14.0
Type 2 diabetes (% yes)	13.9	14.5	9.23	0.2	13.9	31.6	3.0
Current smoking (% yes)	6.2	20.9	27.1	18.4	7.7	16.8	30.0
Prevalent cardiovascular disease (% yes)	5.4	4.8	6.0	0.6	15.7	39.2	NA
cIMT	0.97 (0.13)	0.75 (0.15)	0.96 (0.19)	0.73 (0.17)	0.71 (0.18)	0.79 (0.15)	0.61 (0.12)

Values are mean (SD) or percentages.

Abbreviations: BMI refers to body-mass index. HDL refers to high-density lipoprotein. cIMT refers to carotid intima media thickness.

Discovery analysis

All four studies using Illumina arrays used a linear mixed model where cIMT was transformed and adjusted for various confounders. We ran a separate model for each gene expression probe with cIMT and gene mRNA levels as the independent and dependent variables, respectively. Results were meta-analyzed using an inverse-variance fixed effects model implemented in METAL and adjusted for multiple testing (40). A quantile-quantile (QQ) plot of the discovery analysis is provided in Supplementary Material, Fig. S1. As shown in Table 2, after Bonferroni correction (P-value < 9.2 × 10⁻⁷) four probes were significantly associated with cIMT: ILMN_1780861 and ILMN_1688775 mapping to METRNL (P-value = 4.8 × 10⁻¹⁰ and P-value = 2.8 × 10⁻⁸), ILMN_1702691 mapping to TNFAIP3 (P-value = 1.2 × 10⁻⁷) and ILMN_1782050 mapping to CEBPD (P-value = 2.8 × 10⁻⁷). Forest plots illustrating the association of the four probes in each of the discovery cohorts are provided in Supplementary Material, Fig. S2. The association of ILMN_1782050 (CEBPD) with cIMT showed strong heterogeneity across the included studies (I² = 0.81, P-value = 0.001), whereas the two probes mapping to METRNL, ILMN_1780861 (I² = 0.27, P-value = 0.2) and ILMN_1688775 (I² = 0.24, P-value = 0.3), showed moderate heterogeneity, and ILMN_1702691 (TNFAIP3) showed no heterogeneity (I² = 0, P-value = 0.7) Intensities of all four probes were inversely associated with cIMT (Supplementary Material, Fig. S3). There was widespread correlation between the probes, with pairwise correlations ranging from low to high (Pearson’s r: 0.05–0.89), indicating partially independent effects (Fig. 1).

Table 2

Association of cIMT with significant probes in the differential gene expression discovery (N = 5647) and replication (N = 3953) analyses

			Discovery		Replication
Probe ID	Locus	Gene	Beta	P-value	Z score	P-value
ILMN_1702691	6q23.3	TNFAIP3	−0.46	1.2 × 10−7	−1.70	0.090
ILMN_1782050	8q11.21	CEBPD	−0.39	2.8 × 10−7	−2.16	0.031
ILMN_1780861	17q25.3	METRNL	−0.57	4.8 × 10−10	−1.63	0.10
ILMN_1688775	17q25.3	METRNL	−0.49	2.8 × 10−8	−1.70	0.089

Figure 1

Correlation (Pearson’s r) between the four probes that were significantly associated with cIMT in (A) the Rotterdam Study and (B) LIFE-Adult. ^*Indicates that the correlation is statistically significant.

Replication

For the replication phase, probe-level results from FHS were combined with gene-level results from LIFE-Heart and NESDA, and meta-analyzed in METAL. None of the four probes indicating significantly associated mRNA levels in the discovery analysis replicated after adjusting for multiple testing (P-value < 0.0125). However, the direction of effects were consistent with the discovery analysis.

Additional adjustments

For studies using Illumina arrays, we performed additional analyzes adjusting for other atherosclerosis-relevant covariates. As shown in Supplementary Material, Fig. S4A, most effect sizes did not change substantially when we additionally adjusted for total/high-density lipoprotein (HDL) cholesterol ratio, systolic blood pressure, BMI, prevalent type 2 diabetes, lipid-lowering medication and antihypertensive medication (correlation r² = 0.93). As shown in Supplementary Material, Fig. S4B, effect sizes also remained stable when we repeated the analysis excluding participants with prevalent CHD (correlation r² = 0.96). For the four significant probes in particular, effect sizes decreased when adjusted for additional covariates, though all probes remained nominally significant (Supplementary Material, Table S1A). When participants with prevalent CHD were excluded, effect sizes remained stable or slightly increased. Of the four probes, the probe mapping to TNFAIP3 was the most stable, with effect estimates changing by <10% after adjustment (Supplementary Material, Table S1A).

Pathway analysis

We performed an analysis using ConsensusPathDB to identify pathways with an overrepresentation of genes mapping to the 73 probes (see Supplementary Material, Table S1B) that were suggestively associated [false discovery rate (FDR) < 0.25] with cIMT across the five studies that used Illumina arrays (41). The associated genes were overrepresented in 43 pathways, many of which were related to inflammation and immunity (Supplementary Material, Table S2). The five most significant pathways were ‘Natural killer cell mediated cytotoxicity’, ‘Graft-versus-host disease’, ‘Allograft rejection’, ‘Type I diabetes mellitus’ and ‘Exercise-induced Circadian Regulation’.

Associated genetic variants

To identify genetic variants associated with the expression levels of genes found in whole blood, we used an existing database of cis- and trans-expression quantitative trait locus (eQTL) associated with complex diseases (42). The association of genetic variants with mRNA levels of METRNL, TNFAIP3 and CEBPD is shown in Supplementary Material, Table S4A. We identified variants associated with the expression levels of METRNL and TNFAIP3. For METRNL only nearby genetic variants (cis) were associated, whereas for TNFAIP3 both nearby variants as well as a distant (trans) variant were associated. None of the three variants was associated with cIMT in a large published GWAS, which was based on data from 71 128 individuals [Supplementary Material, Table S4B; (13)].

S-PrediXcan analysis

S-PrediXcan uses summary statistics from GWAS studies and models developed using gene expression data from the genotype tissue expression (GTEx) project to identify associations between predicted gene expression and a phenotype of interest. We used S-PrediXcan on summary statistics from a CHARGE GWAS of cIMT with 71 128 European individuals to compute the association between cIMT and genetically predicted gene expression. The number of genes that passed quality control and were included in the S-PrediXcan analysis ranged from 1633 in kidney cortex tissue to 9980 in tibial nerve tissue (Supplementary Material, Table S5), and the median predicted performance r² ranged from 0.064 (interquartile range (IQR): 0.03–0.142) for skeletal muscle tissue to 0.135 (IQR: 0.094–0.213) in kidney cortex tissue (Supplementary Material, Table S5). QQ-plots for each of the tested tissues are shown in Supplementary Material, Fig. S5A. After Benjamini–Hochberg correction, 86 tissue-gene pairs were significantly (FDR < 0.05) associated with cIMT, involving 36 different genes and 33 different tissues (Supplementary Material, Table S6). Limiting these results to a single lead tissue-gene pair at each locus, corresponded to 22 loci that were significantly associated with cIMT. Descriptions of the associated loci are provided in Supplementary Material, Table S6. Of the 86 tissue-gene pairs found to be significantly associated with cIMT, 21% (18/86) involved artery-related tissues, more than what would be expected by chance (P-value: 5.24 × 10⁻⁶). In addition, 12 out of the 22 genes have not, to our knowledge, been previously associated with cIMT (Supplementary Material, Table S6).

We replicated these results by applying S-PrediXcan to an independent GWAS of cIMT in the UK Biobank. Of the 86 tissue-gene pairs, 34 replicated after Bonferroni correction (P-value < 5.8E−4; Supplementary Material, Table S7). Comparing z-scores for the 86 significant gene pairs in the discovery and replication S-PrediXcan results, we observed high concordance, with a pairwise correlation of 0.78. After limiting the results to a single lead tissue-gene pair, six loci successfully replicated, including the MCPH1-AS1, CBFA2T3, BCAR1, KIAA1462, CDH13 and TLN2 loci (Table 3). Forest plots showing S-PrediXcan effect estimates for the six replicating genes across all tissues assessed are shown in Supplementary Material, Fig. S6. All six replicating genes showed a consistent direction of effect. Of the replicating loci, one gene, TLN2, has not previously been associated with cIMT. Regional plots for the associations of genetic variants at these six loci in the CHARGE and UK Biobank GWAS are shown in Supplementary Material, Fig. S7. The single variant association results for variants in the S-PrediXcan models for the six replicating loci are shown in Supplementary Material, Table S8A and B, for the CHARGE and UKBB GWAS studies, respectively.

Table 3

Replication analysis of CHARGE significant S-PrediXcan genes in UK Biobank

		CHARGE			UK Biobank
Gene	Tissue	Effect size	P-value	FDR	Effect size	P-value	FDR	Novel
MCPH-AS1	Nerve Tibial	0.008	6.73E−06	0.0278	0.023	2.62E−11	7.26E−06	No
CBFA2T3	Testis	−0.155	1.41E−08	0.0003	−0.241	3.40E−07	0.006	No
BCAR1	Esophagus Mucosa	−0.013	4.35E−07	0.0039	−0.027	3.44E−07	0.006	No
KIAA1462	Artery Tibial	0.039	2.85E−09	0.0002	0.061	1.83E−06	0.022	No
CDH13	Artery Aorta	0.012	1.53E−05	0.0494	0.020	1.32E−04	0.249	No
TLN2	Adipose Subcutaneous	−0.084	7.16E−06	0.0291	−0.138	1.97E−04	0.280	Yes

Comparison of differential expression and S-PrediXcan analyses

We examined the correlation of the regression coefficients in the CHARGE whole blood differential gene expression analysis with those from the whole blood component of our S-PrediXcan study. The overall pairwise correlation between beta coefficients was low (Pearson’s r: 0.024). When stratified by quartiles of model fit r², the beta coefficients had a mixed upward trend (Pearson’s r: Q1: 0.024; Q2: 0.034; Q3:0.011 and Q4:0.059) with a dip in correlation observed in the third quantile. There was modest enrichment of concordance between effect size for genes that were more significantly associated with trait in one or both analyses. For genes where the association P-value from both analyses was ≤ 0.05, the correlation of the beta coefficients was 0.167. When including all genes that were nominally significant (P-value ≤ 0.05) in each method regardless of significance in the other analysis, the correlation was 0.087 and 0.073 for S-PrediXcan and differential gene expression analysis, respectively (Supplementary Material, Fig. S8A and B).

We then explored the six replicating genes found in the S-PrediXcan analyses to determine how their association with cIMT compared with our findings in the CHARGE differential gene expression analysis study. None of these genes identified in S-PrediXcan were significant (P < 0.0071; 0.05/7) in the differential gene expression results (Supplementary Material, Table S9A). When exploring the three significant genes found in the CHARGE differential gene expression analysis in our S-PrediXcan analysis, none of the genes were significant at P < 0.0042; (0.05/12; Supplementary Material, Table S9B). We correct for multiple testing using a Bonferroni correction in the DE analysis and FDR in the S-PrediXcan analysis. To test the sensitivity of these comparisons with the multiple testing correction approach used, we performed an additional Benjamini–Hochberg FDR-based multiple testing correction for the differential gene expression analysis. Using this approach, we find two additional probes that were significantly associated with cIMT: ILMN_1661599 (P-value = 1.6 × 10⁻⁶) and ILMN_1876266 (P-value = 4.4 × 10⁻⁶), mapping to DDIT4 and GJA3, respectively (Supplementary Material Table S2). These genes were not significantly associated using S-PrediXcan and neither gene has previously been associated with cIMT in GWAS.

We also compared the proportion of genes previously identified in GWAS in each approach, with the hypothesis that our S-PrediXcan results are enriched for genes previously implicated in GWAS, and we find that for nominally significant genes (i.e. P < 0.05 in our results), there is no significant enrichment of GWAS-identified loci, P = 0.22. When we restrict to more significant genes (i.e. P < 0.005 in our results), we see a statistically significant enrichment, P = 0.042, of GWAS loci for S-PrediXcan, relative to DE (Supplementary Material, Table S10).

Discussion

We performed the first large-scale differential gene expression study meta-analysis of cIMT based on whole-blood mRNA levels of over 5600 participants. We complemented this ex vivo study of differential gene expression in whole blood with an in silico study of predicted gene expression using S-PrediXcan across 48 tissues. In our ex vivo study of differential gene expression, we identified four oligonucleotide probes indicating differential mRNA levels associated with cIMT that mapped to three genes, namely: TNFAIP3, CEBPD and METRNL, although these associations did not convincingly replicate in an independent sample. In our in silico study of predicted gene expression, we were able to identify five known loci and one novel locus that were significantly associated with cIMT and that replicated in an independent dataset. These results lend further support to previous findings that predicted gene expression association methods like S-PrediXcan can both identify novel loci through increased power to aggregate smaller effects and also replicate previous GWAS results (43).

Results from our differential expression analysis have compelling functional evidence of involvement in inflammation. TNFAIP3 is involved in the negative feedback regulation of NF-kappaB (44), but it may also inhibit IFNγ/STAT1 signaling (45). Thus, it represents an anti-inflammatory protein, and low expression of TNFAIP3 has been associated with inflammatory disorders such as rheumatoid arthritis (46,47). Furthermore, a study in mice found that Tnfaip3 slowed the progression of atherosclerosis by reducing NF-kappaB activity (48). The proposed protective effects of TNFAIP3 are in line with our study, in which expression was inversely associated with cIMT. CEBPD encodes CCAAT/Enhancer Binding Protein Delta (C/EBP-Delta), a transcription factor regulating several inflammatory genes (49). Previous studies suggest that C/EBP-Delta can be both pro-inflammatory and anti-inflammatory: C/EBP-Delta may amplify the NF-kappaB response (50,51), however, C/EBP-Delta has been shown to have an anti-inflammatory role in pancreatic β-cells and brain pericytes (52,53), inhibiting the accumulation of amyloid plaques in Alzheimer’s disease (54). In our study, increased expression of CEBPD in blood is associated with less atherosclerosis as measured by cIMT. METRNL encodes meteorin-like protein, which increases thermogenesis in brown and beige adipocytes, and increases the expression of genes encoding anti-inflammatory proteins (55). Recent studies have also shown that serum METRNL concentrations are negatively correlated with TNF-a and IL-6 in patients with type 2 diabetes mellitus and coronary artery disease (56,57), as well as negatively associated with cIMT (57), which is consistent with our observations. The potential effects on both adiposity and inflammation could explain the inverse association of METRNL expression with cIMT in our study (52,53).

All three genes identified in the differential gene expression analysis thus appear to be related to inflammation. This is not surprising, given the importance of inflammation in atherosclerosis, and the fact that expression levels were measured in the blood, where we expect most mRNA to originate from white blood cells. None of the three genes was reported to be significantly associated with interleukin-6 levels in a whole blood differential gene expression study (58). There has been no previous large-scale differential gene expression study of cIMT, but several studies of CHD have been conducted. None of the three genes we report were significant in these previous studies (26–32).

The genes identified in the ex vivo differential expression analysis also did not include genes that had been previously identified by GWAS of cIMT, CHD or stroke (13,59,60). This lack of overlap may have several possible explanations. First, the genetic background of atherosclerosis and CHD may be differentially reflected through genetic variation and gene expression levels. In a large-scale transcriptome-wide association study of blood pressure, only 2 out of 34 genes were previously reported to be linked to hypertension, and none were identified through GWAS (61). This pattern may hold for other traits, including cIMT. Second, although blood is a relevant tissue for atherosclerosis, it may not be the tissue in which the genes identified by GWAS are primarily expressed.

In contrast, five of the six associated loci that were significantly associated with cIMT using S-PrediXcan had been previously identified by GWAS on cIMT (13,14,17), supporting the validity of our in silico approach and revealing that variation in gene expression because of genetic variants is likely involved in disease pathogenesis (62). The one gene that had not been previously reported by GWAS was TLN2, encoding Talin-2. Although the function of Talin-2 is less clearly understood than that of Talin-1, Talin-1 and Talin-2 are closely related isoforms that are part of a family of cytoskeletal proteins that are important for organ and vascular tissue development (63–65). Talins help integrin proteins to connect the cytoskeleton with the extracellular matrix and are thought to provide tension transmission between the contractile apparatus of the actin cytoskeleton and the extracellular matrix (66–68). Talins are major constituents of the focal adhesion molecules present in the endothelial cells that line blood vessels (68,69). Altered expression of Talins, including TLN2, may impair the ability of adhesion molecules in endothelial cells to withstand variations in blood pressure, speed and force (68). A number of genes encoding proteins involved in cell mechanostability, including adhesion molecules and TLN2, have been linked with the progression of atherosclerosis (70–73). The reduction of TLN2 in endothelial cells may impair the intercellular gap junctions and allow macrophage accumulation and consequent arterial wall thickening and plaque formation (73–75), consistent with our observed negative association between cIMT and genetically determined expression of TLN2. We examined results across all assessed tissues for TLN2 and each of the other five replicating genes from the S-PrediXcan analysis (Supplementary Material, Fig. S5A–F). We found that results for available artery tissue models for all replicating genes showed evidence of effect, as the 95% confidence interval for the effect size did not include a null effect, even if they did not survive multiple testing corrections. The effect direction was also the same across arterial tissues for each gene.

To our knowledge, this is the first report of a direct comparison between the results of S-PrediXcan and transcriptome-wide differential gene expression analysis in the same tissue. We note that comparison of the two approaches is complicated by factors outside of differences in the methods themselves, including differences in the number of tissues assessed and sample sizes available for each, however, we did not observe a strong correlation between results for these different analysis methods. One possible explanation is that effect estimates from genetically predicted expression models capture the intrinsic effects of genetic variation on gene expression with cIMT, whereas the effect estimates from differential gene expression analysis of measured gene expression levels also represent changes that may be the result of the tissue-specific environment, in this case, atherosclerotic vascular disease. As a result, associations with gene expression levels using the ex vivo approach are likely to be affected by confounding and reverse causation. Additionally, the genetic component of gene expression as computed by S-PrediXcan is often quite minor (median r² = 0.088), meaning that differences in expression of some genes may not be captured using this approach. Indeed, total cis- and trans-heritability of gene expression is quite variable (76), depending on the properties of regulation of expression, each approach may be better suited to examining the impact of changes in expression for different genes.

Strengths of our differential expression analysis include the large sample sizes for both discovery and replication, the hypothesis-free approach and the strict correction for multiple testing. Although we consider blood to be a relevant tissue for the expression of genes associated with atherosclerosis, not having access to microarray or RNA sequencing data for other, less easily accessible tissues is a limitation of this study. Microarray transcriptomic data also have known limitations for detection of low abundance transcripts and a narrower dynamic range compared with RNA sequencing (77), which may have reduced power in our differential expression analysis. Use of RNA sequencing in future studies would reduce these limitations. Furthermore, the interpretation of the results is challenging because it is difficult to distinguish the direction of causality: genes whose expression influences atherosclerosis and genes whose expression is influenced by atherosclerosis may both be differentially expressed in this type of study. Likewise, probes that are associated with cIMT may not be captured because of the stochastic nature of gene expression. A host of intrinsic and extrinsic factors, including environmental stimuli, host factors and batch effects, could lead to underestimation of true effects at some probes (76–78). Although this documented low signal-to-noise ratio of gene expression data are not expected to lead to spurious associations, it is likely to reduce statistical power and prevent some probes from reaching significance thresholds. Comparatively, genetically determined gene expression is definitionally unaffected by questions of causality and other similar confounders; this may explain the low level of concordance between S-PrediXcan and differential gene expression analysis. Finally, although we adjusted for batch effects, cell types, and, in an additional analysis, traditional cardiovascular risk factors, the associations described in our differential expression study may be affected by residual confounding, environmental variation and experimental variation.

We attempted to mitigate some of the limitations of our differential expression analysis by using S-PrediXcan, which examines multiple tissues and is not affected by reverse causation or confounding. A limitation of using S-PrediXcan is that the statistical power to detect associations between cIMT and gene expression is constrained by the degree to which gene expression of each gene is genetically determined. This approach will miss associations with genes for which gene expression is determined mostly by environmental factors, even those where gene expression changes precede development of disease and may in fact be causal. Considering blood tissue procurement in GTEx was pre-mortem for some donors and post-mortem for others, the donor’s cause of death and timing of tissue procurement can influence the quality of collected tissues, like whole blood (78). This could help explain the low correlation observed between the differential gene expression analysis and S-PrediXcan results.

We performed a transcriptome-wide differential gene expression analysis and for the first time detected three genes whose mRNA levels in whole blood were associated with atherosclerosis as measured by cIMT. The proteins encoded by all three genes are reported to be involved in inflammation, with TNFAIP3 and METRNL having well described anti-inflammatory properties. We also applied S-PrediXcan to existing GWAS studies of cIMT in CHARGE and the UK Biobank, implicating TLN2 as a novel gene involved in intimal thickening. Notably, we identified a lack of concordance between results from our ex vivo differential gene expression analysis and our in-silico S-PrediXcan analysis. There was only a very modest enrichment of concordance of effect for genes that were more significantly associated with cIMT in one or both approaches, suggesting that these approaches capture different components of gene expression and its association with cIMT.

Materials and methods

Whole blood differential expression study population

This study was conducted within the framework of the Cohorts for Heart and Aging Research in Genomic Epidemiology (CHARGE) consortium (79). Individuals from four population-based cohort studies were included in the discovery analysis: 836 from the Cooperative Heath Research in the Region of Augsburg (KORA F4) study (79), 2973 from LIFE-Adult (80,81), 856 from the Rotterdam Study (82) and 982 from the Study of Health in Pomerania [SHIP-TREND (83)]. The total discovery sample size was 5647. The replication analysis was done in 1958 individuals from the Framingham Heart Study [FHS (84)], 1752 individuals from the LIFE-Heart study (85) and 243 individuals from the Netherlands Study of Depression and Anxiety (NESDA) study (86). Further details are provided in the Supplemental Methods. All studies were approved by appropriate research ethics committees and all participants signed informed consent prior to participation.

Measurement of cIMT

cIMT of the common carotid artery was measured with high-resolution B-mode ultrasonography. cIMT was calculated as the mean of the maximum cIMT of the near and far walls of the right and left common carotid arteries. When the intima media thickness of the near walls was unavailable, only the far walls were used. Further details are provided in the Supplemental Methods.

Measurement of gene expression levels

Genome-wide gene expression levels in whole blood, detected by up to 54 124 oligonucleotide probes that detected gene-specific mRNA levels of more than 25 000 genes were determined using the Illumina HumanHT-12 Gene Expression BeadChip v3.0 or v4.0, except in FHS and NESDA, where the Affymetrix Exon Array ST 1.0 and Affymetrix Human Genome U219 Array were used. In LIFE-Heart peripheral blood mononuclear cells were used. Gene expression levels were measured in blood that was drawn at the same time as cIMT was measured. See Supplemental Methods for additional details.

Differential expression analysis

Analyses were conducted separately in each study using a standardized analysis plan. cIMT was natural log transformed, and each study used a linear mixed model adjusting for batch effects (examples: array ID and position on array) as random effects and adjusting for further technical covariates (examples: RNA quality and storage time between sampling and RNA isolation), cell types (examples: granulocytes, lymphocytes and monocytes), age, sex, fasting state and smoking status as fixed effects (Model 1). In addition, FHS included familial relatedness as a random effect. We ran a separate model in R for each gene expression probe, using cIMT and the covariates as independent variables, and gene mRNA levels as the dependent variable. The study-specific results were meta-analyzed using inverse-variance fixed effects meta-analysis implemented in METAL (40). We used a Bonferroni correction to adjust for multiple testing.

In the replication, probe-level results in LIFE-Heart were combined with gene-level results from FHS and NESDA using meta-analysis of P-values (taking sample size and direction of effect into account), as implemented in METAL (40). We used a Bonferroni correction to correct the significance threshold for the number of probes included in the replication.

Among the five studies that used Illumina arrays, we performed additional analyses including further covariates relevant to atherosclerosis: total/HDL cholesterol ratio, systolic blood pressure, BMI, prevalent type 2 diabetes, lipid-lowering medication and antihypertensive medication (Model 2). We also repeated the original model in only those individuals with data available on all the additional covariates (Model 1). Finally, we reran the full model excluding individuals with prevalent CHD (Model 3). All meta-analyses were performed using inverse-variance fixed effects meta-analysis implemented in METAL (40).

We computed pairwise correlations (Pearson’s r) between any significantly associated probes in the participants from the Rotterdam Study and LIFE-Adult. We used ConsensusPathDB to identify pathways that showed an overrepresentation of genes mapping to probes that were suggestively associated (FDR < 0.25) with cIMT across the five cohorts that used Illumina arrays (41,87). Pathways with an FDR < 0.05 were considered significant.

Finally, we used an existing database (http://genenetwork.nl/bloodeqtlbrowser/) to identify genetic variants associated with expression levels of significant genes in whole blood of over 5000 individuals (42). Both genetic variants closer than (cis) and further than (trans) 250 kb were included, but only suggestive associations (FDR < 0.05) are included in the database. We then queried variants identified in a published GWAS of cIMT (8).

S-PrediXcan

We applied S-PrediXcan to summary statistics from a CHARGE GWAS of cIMT that included 71 128 individuals of European ancestry originating from 31 different studies (13,79,88). In the CHARGE GWAS, cIMT was defined as the mean of the maximal measurements from the near and far walls of the internal carotid arteries for both the left and right sides (8). For the S-PrediXcan replication analysis we used summary statistics from an independent GWAS of cIMT in the UK Biobank, including 22 179 individuals of British ancestry (14). In the UK Biobank GWAS, cIMT was defined as the largest of the four maximum measurements of the far wall from both the left and right of the distal common carotid artery (14). Summary statistics for the CHARGE GWAS study have been made publicly available and can be accessed here https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000930.v9.p1.

We used S-PrediXcan to estimate the association of genetically determined mRNA levels with cIMT using S-PrediXcan models built in tissue-specific gene expression data from the GTEx project (33). The GTEx v8 models include gene expression samples of 48 different tissues with at least 65 samples (89). S-PrediXcan model-building methods have been previously described (90). These models assign weights to each cis-eQTL variant, and include measures of variances and covariances of genetic markers in the GTEx data to account for linkage disequilibrium between SNPs. All model files are stored in a publicly available resource (http://predictdb.org/). We applied S-PrediXcan models for 48 different tissues to results from GWAS on cIMT in the CHARGE (14,79). The tissue-specific results were filtered to include genes with a cross-validation correlation >0.1. Genetically determined mRNA levels in a given tissue were considered to be significantly associated with cIMT and were taken forward for replication if the FDR was < 0.05. A Bonferroni correction was used to adjust the significance threshold for the number of statistical tests performed in the replication analysis.

We tested for enrichment of GWAS-identified genes in the S-PrediXcan results relative to the DE results by first creating a list of genome-wide significant single variants associated with IMT in Strawbridge et al. (14) and Franceschini et al. (13) We then created regions ±500 kb surrounding the reported single variants and queried our artery tissue (aortic, coronary and tibial) S-PrediXcan results and DE results for genes or probes that overlapped those regions with P < 0.05 and again with P < 0.005. We performed a one-sided Fisher’s exact test to test that the proportion of genes that overlapped with GWAS loci was higher in S-PrediXcan than in DE.

Acknowledgements

We thank all study participants, volunteers and study personnel that made this work possible. We thank Lesca Holdt very much for helping to create the LIFE-Heart gene-expression data, and thank Kay Olischer and Annegret Unger very much for technical assistance. We thank Kerstin Wirkner very much for running the LIFE study center, and thank Sylvia Henger very much for LIFE-Adult data quality control. We thank Ms Mila Jhamai, Ms Sarah Higgins, Mr Marijn Verkerk and Jeroen van Rooij, MSc for their help in creating the RNA array expression database for the Rotterdam Study.

Conflict of interest statement. The authors declare no competing interests.

Funding

This work was supported by American Heart Association Grant [grant numbers 17POST33350042, 18CDA34110116 to P.V.]; National Heart, Lung, and Blood Institute [grant numbers R01HL146860 to P.V, R01HL142302 to J.B. and L.P., R01HL105756]; National Institute of Aging [grant number R01AG061351 to J.B. and L.P]. The LIFE-Adult study, LIFE-Heart study, and the Leipzig Research Centre for Civilization Diseases was supported by the European Union; the European Regional Development Fund; and the Free State of Saxony [project numbers 713-241202, 713-241202, 14505/2470, 14575/2470]. The Framingham Heart Study of the National Heart, Lung, and Blood Institute of the US National Institutes of Health is supported by the Boston University School of Medicine; the Systems Approach to Biomarker Research in Cardiovascular Disease (SABRe) project; the National Heart, Lung and Blood Institute Division of Intramural Research; and the United States National Institute of Health [contract N01-HC-25195]. The Cooperative Health Research in the Region of Augsburg (KORA F4) study and the Helmholtz Zentrum München – German Research Center for Environmental Health were supported by the German Federal Ministry of Education and Research; the State of Bavaria; the Munich Center of Health Sciences, Ludwig-Maximilians-Universität LMUinnovativ grant; the German Federal Ministry of Health (BMG); the Ministry of Culture and Science of the State North Rhine-Westphalia; and the German Center for Diabetes Research which is funded by the and the European Union Seventh Framework Programme [grant numbers n°603288, n°313010, n°602736]. The Netherlands Study of Depression and Anxiety is supported by VU University Medical Center; GGZ ingest; Arkin; Leiden University Medical Center; GGZ Rivierduinen; University Medical Center Groningen; Lentis; GGZ Friesland; GGZ Drenthe; Scientific Institute for Quality of Healthcare Netherlands Institute for Health Services Research; Netherlands Institute of Mental Health and Addiction (Trimbos Institute); the Netherlands Heart Foundation [grant number 2006B258]; the US National Institute of Mental Health [grant number RC2 MH089951]; and the Geestkracht program of the Netherlands Organization for Health Research and Development [grant number 10-000-1002]. The Rotterdam Study is supported by the Human Genotyping Facility of the Genetic Laboratory of the Department of Internal Medicine, Erasmus Medical Center; Erasmus Medical Center and Erasmus University Rotterdam; the Netherlands Organization for the Health Research and Development (ZonMw); the Research Institute for Diseases in the Elderly (RIDE); the Ministry of Education, Culture, and Science; the Ministry for Health, Welfare and Sports; the European Commission (DG XII); and the Municipality of Rotterdam. The Study of Health in Pomerania – TREND is funded by Siemens Healthcare; the Federal State of Mecklenburg, West Pomerania; the `Center of Knowledge Interchange' program of the Siemens AG; and the Federal Ministry of Education and Research [grant numbers 01ZZ9603, 01ZZ0103, 01ZZ0403, 03IS2061A, 03ZIK012].

Author notes

Andy B. Castaneda, Lauren E. Petty, Jennifer E. Below, Abbas Dehghan and Paul S. de Vries contributed equally as junior or senior author.