A Meta-analysis of Gene Expression Signatures of Blood Pressure and Hypertension

Genome-wide association studies (GWAS) have uncovered numerous genetic variants (SNPs) that are associated with blood pressure (BP). Genetic variants may lead to BP changes by acting on intermediate molecular phenotypes such as coded protein sequence or gene expression, which in turn affect BP variability. Therefore, characterizing genes whose expression is associated with BP may reveal cellular processes involved in BP regulation and uncover how transcripts mediate genetic and environmental effects on BP variability. A meta-analysis of results from six studies of global gene expression profiles of BP and hypertension in whole blood was performed in 7017 individuals who were not receiving antihypertensive drug treatment. We identified 34 genes that were differentially expressed in relation to BP (Bonferroni-corrected p<0.05). Among these genes, FOS and PTGS2 have been previously reported to be involved in BP-related processes; the others are novel. The top BP signature genes in aggregate explain 5%–9% of inter-individual variance in BP. Of note, rs3184504 in SH2B3, which was also reported in GWAS to be associated with BP, was found to be a trans regulator of the expression of 6 of the transcripts we found to be associated with BP (FOS, MYADM, PP1R15A, TAGAP, S100A10, and FGBP2). Gene set enrichment analysis suggested that the BP-related global gene expression changes include genes involved in inflammatory response and apoptosis pathways. Our study provides new insights into molecular mechanisms underlying BP regulation, and suggests novel transcriptomic markers for the treatment and prevention of hypertension.


Author Summary
The focus of blood pressure (BP) GWAS has been the identification of common DNA sequence variants associated with the phenotype; this approach provides only one dimension of molecular information about BP. While it is a critical dimension, analyzing DNA variation alone is not sufficient for achieving an understanding of the multidimensional complexity of BP physiology. The top loci identified by GWAS explain only about 1 percent of inter-individual BP variability. In this study, we performed a meta-analysis of gene expression profiles in relation to BP and hypertension in 7017 individuals from six studies. We identified 34 differentially expressed genes for BP, and discovered that the top BP signature genes explain 5%-9% of BP variability. We further linked BP gene expression signature genes with BP GWAS results by integrating expression associated SNPs (eSNPs) and discovered that one of the top BP loci from GWAS, rs3184504 in SH2B3, is a trans regulator of expression of 6 of the top 34 BP signature genes. Our study, in conjunction Introduction Systolic and diastolic blood pressure (SBP and DBP) are complex physiological traits that are affected by the interplay of multiple genetic and environmental factors. Hypertension (HTN) is a critical risk factor for stroke, renal failure, heart failure, and coronary heart disease [1]. Genome-wide association studies (GWAS) have identified numerous loci associated with BP traits [2,3]. These loci, however, only explain a small proportion of inter-individual BP variability. In aggregate the 29 loci reported by the International Consortium of Blood Pressure (ICBP) consortium GWAS account for about one percent of BP variation in the general population [3]. Most genes near BP GWAS loci are not known to be mechanistically associated with BP regulation [3]. Therefore, further studies are needed to determine whether the genes implicated in GWAS demonstrate functional relations to BP physiology and to uncover the molecular actions and interactions of genetic and environmental factors involved in BP regulation.
Alterations in gene expression may mediate the effects of genetic variants on phenotype variability. We hypothesized that characterizing gene expression signatures of BP would reveal cellular processes involved in BP regulation and uncover how transcripts mediate genetic and environmental effects on BP variability. We additionally hypothesized that by integrating gene expression profiling with genetic variants associated with altered gene expression (eSNPs or eQTLs) and with BP GWAS results, we would be able to characterize the genetic architecture of gene expression effects on BP regulation.
Several previous studies have examined the association of global gene expression with BP [4,5] or HTN [6,7]. Most of these studies, however, were based on small sample sizes and lacked replication [4,5,6,7]. To address this challenge, we conducted an association study of global gene expression levels in whole blood with BP traits (SBP, DBP, and HTN) in six independent studies. In order to avoid the possibility that the differentially expressed genes we identified reflect drug treatment effects, we excluded individuals receiving anti-hypertensive treatment. The eligible study sample included 7017 individuals: 3679 from the Framingham Heart Study (FHS), 972 from the Estonian Biobank (EGCUT), 604 from the Rotterdam Study (RS) [8], 597 from the InCHIANTI Study, 565 from the Cooperative Health Research in the Region of Augsburg [KORA F4] Study [9], and 600 from the Study of Health in Pomerania [SHIP-TREND] [10]. We first identified differentially expressed BP genes in the FHS (n = 3679) followed by external replication in the other five studies (n = 3338). Subsequently, we performed a meta-analysis of all 7017 individuals from the six studies, and identified 34 differentially expressed genes associated with BP traits using a stringent statistical threshold based on Bonferroni correction for multiple testing of 7717 unique genes. The differentially expressed genes for BP (BP signature genes) were further integrated with eQTLs and with BP GWAS results in an effort to differentiate downstream transcriptomic changes due to BP from putatively causal pathways involved in BP regulation.

Clinical characteristics
After excluding individuals receiving anti-hypertensive treatment, the eligible sample size was 7017 (FHS, n = 3679; EGCUT, n = 972; RS, n = 604; InCHIANTI, n = 597; KORA F4, n = 565 SHIP is part of the Community Medicine Research net of the University of Greifswald, Germany, which is funded by BMBF (grants no. 01ZZ9603, 01ZZ0103, and 01ZZ0403), the Ministry of Cultural Affairs as well as the Social Ministry of the Federal State of Mecklenburg-West Pomerania, and the network 'Greifswald Approach to Individualized Medicine (GANI_MED)' funded by the BMBF (grant 03IS2061A). Generation of genome-wide data has been supported by the BMBF (grant no. 03ZIK012). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. and SHIP-TREND, n = 600). Clinical characteristics of participants from the four studies are presented in Table 1. The mean age varied across the cohorts (FHS = 51, EGCUT = 36, RS = 58, InCHIANTI = 71, KORA F4 = 72 and SHIP-TREND = 46 years) as did the proportion of individuals with hypertension (11% in FHS, 19% in EGCUT, 35% in RS, 45% in InCHIANTI, 26% in KORA, and 12% in SHIP).

Identification and replication of differentially expressed BP signature genes
At a Bonferroni corrected p<0.05, we identified 73, 31, and 8 genes that were differentially expressed in relation to SBP, DBP, and HTN, respectively in the FHS, which used an Affymetrix array for expression profiling, and 6, 1, and 1 genes in the meta-analysis of the 5 cohorts that used an Illumina array (Illumina cohorts): EGCUT, RS, InCHIANTI, KORA F4 and SHIP-TREND (S1 Table). For each differentially expressed BP gene in the FHS or in the Illumina cohorts, we attempted replication in the other group. At a replication p<0.05 (Bonferroni corrected), 13 unique genes that were identified in the FHS were replicated in the Illumina cohorts, including 10 for SBP (CD97, TAGAP, DUSP1, FOS, MCL1, MYADM, PPP1R15A, SLC31A2, TAGLN2, and TIPARP), 5 for DBP (CD97, BHLHE40, PRF1, CLC, and MYADM), and 2 for HTN (GZMB and MYADM) ( Table 2). Each of the unique BP signature genes in the Illumina cohorts, 6 for SBP (TAGLN2, BHLHE40, MYADM, SLC31A2, DUSP1, and MCL1), 1 for DBP (BHLHE40) and 1 for HTN (SLC31A2), replicated in the FHS. All 6 Illumina cohorts BP signature genes that replicated in the FHS were among the 13 FHS BP signature genes that replicated in the Illumina cohorts. The BP signature genes identified in the FHS showed enrichment in the Illumina cohorts at pi1 = 0.88, 0.75, and 0.99 for SBP, DBP, and HTN respectively (pi1 value indicates the proportion of significant signals among the tested associations [11]; see details in the Methods section). Fig. 1 shows that the mean gene expression levels of the top BP signature genes were consistent with the BP phenotypic changes observed in the FHS and the Illumina cohorts.
The 73 SBP signature genes in the FHS (55 of these 73 genes were measured in the Illumina cohorts) at a Bonferroni corrected p<0.05 in aggregate explained 9.4% of SBP phenotypic variance in the Illumina cohorts, and the 31 DBP signature genes from the FHS (22 of these 31 genes were measured in the Illumina cohorts) in aggregate explained 5.3% of DBP phenotypic variance in the Illumina cohorts. These results suggest that in contrast to common genetic variants identified by BP GWAS, which explain in aggregate only about 1% of inter-individual BP variation [3], changes in gene expression levels explains a considerably larger proportion of phenotypic variance in BP.  Meta-analysis of the six cohorts identifies differentially expressed BP signature genes A meta-analysis of differential expression across all six cohorts revealed 34 differentially expressed BP genes at p<0.05 (Bonferroni corrected for 7717 genes that were measured and passed quality control in the FHS and Illumina cohorts), including 21 for SBP, 20 for DBP, and 5 for HTN ( Table 2 and S2 Fig.). All of the 34 differentially expressed BP signature genes showed directional consistency in the FHS and the Illumina cohorts ( Table 2). The 34 BP signature genes included all 13 genes that were cross-validated between the FHS and the Illumina cohorts. Of the 34 BP signature genes, 27 were positively correlated with BP and only 7 genes were negatively correlated. MYADM and SLC31A2 were top signature genes for SBP, DBP, and HTN. At FDR<0.2, 224 unique genes were differentially expressed in relation BP phenotypes including 142 genes for SBP, 137 for DBP, and 45 for HTN (details are reported in the S1-S2 Text, and S3-S5 Table).

Functional analysis of differentially expressed BP signature genes
We used gene set enrichment analysis (GSEA) to identify the biological process and pathways associated with gene expression changes in relation to SBP, DBP, and HTN in order to better understand the biological themes within the data. As shown in Table 3, the GSEA of genes whose expression was positively associated with BP showed enrichment for antigen processing and presentation (p<0.0001), apoptotic program (p<0.0001), inflammatory response (p<0.0001), and oxidative phosphorylation (p = 0.0018). The negatively associated genes showed enrichment for nucleotide metabolic process (p<0.0001), positive regulation of cellular metabolic process (p<0.0001), and positive regulation of DNA dependent transcription (p = 0.0021).

Genetic effects on expression of BP signature genes
Among the 34 BP signatures genes from the meta-analysis of all 6 studies, 33 were found to have cis-eQTLs and 26 had trans-eQTLs ( Fig. 2A and S2 Table) based on whole blood profiling [12,13]. Of these, six master trans-eQTLs mapped to either five or six BP signature genes (no master cis-eQTL was identified). Five master trans-eQTLs (rs653178, rs3184504, rs10774625, rs11065987, and rs17696736) were located on chromosome 12q24 within the same linkage disequilibrium (LD) block (r 2 >0.8, Fig. 2B). We retrieved a peak cis-and trans-eQTL for each BP signature gene. The peak cis-eQTL explained 0.2-20% of the variance in the corresponding transcript levels, in contrast, the peak trans-eQTL accounted for very little (0.02-2%) of the corresponding transcript variance. Westra et al. also reported a similar small proportion of variance in transcript levels explained by trans-eQTLs [12]. We then linked the cis-and trans-eQTLs of the 34 BP signature genes with BP GWAS results from the ICBP Consortium [3] and the NHGRI GWAS Catalog [14] (Fig. 2 and S2 Table). We did not find any cis-eQTLs for the top BP signature genes that also were associated with BP in the ICBP GWAS [3]. However, the 6 master trans-eQTLs were all associated with BP at p<5e-8 in the ICBP GWAS [3] and were associated with multiple complex diseases or traits ( Table 4). For example, rs3184504, a nonsynonymous SNP in SH2B3 that was associated in GWAS with BP, coronary heart disease, hypothyroidism, rheumatoid arthritis, and type 1 diabetes [12], is a trans-eQTL for 6 of our 34 BP signature genes from the meta-analysis (FOS, MYADM, PP1R15A, TAGAP, S100A10, and FGBP2; Fig. 2A-B and Table 4). These 6 genes are all highly expressed in neutrophils, and their expression levels are correlated significantly (average r 2 = 0.04, p<1e-16). rs653178, intronic to ATXN2 and in perfect LD with rs3184504 (r 2 = 1), also is associated with BP and multiple other diseases in the NHGRI GWAS Catalog [14]. It also is a trans-eQTL for the same 6 BP signature genes ( Table 4). These two SNPs are cis-eQTLs for expression SH2B3 in whole blood (FDR<0.05), but not for ATXN2 (FDR = 0.4). We found that the expression of SH2B3 is associated with expression of MYADM, PP1R15A, and TAGAP (at Bonferroni corrected p<0.05), but not with FOS, S100A10, or FGBP2. The expression of ATXN2 was associated with expression of 5 of the 6 genes (PP1R15A was not associated). S3 Fig. shows the coexpression levels of the eight genes that were cis-or transassociated with rs3184504 and rs653178 genotypes. These results suggest that there may be a pathway or gene co-regulatory mechanism underling BP regulation involving these genes that is driven by this common genetic variant (rs3184504; minor allele frequency 0.47) or its proxy SNPs. We further checked whether the cis-or trans-eQTLs for the top 34 BP signature genes are associated with other diseases or traits in the NHGRI GWAS catalog [14]. We identified 12 cis-eQTLs (for 8 genes) and 6 trans-eQTLs (for 6 genes) that are associated with other diseases or traits in the NHGRI GWAS catalog [14] ( Table 4).

Discussion
Our meta-analysis of gene expression data from 7017 individuals from six studies identified and characterized whole blood gene expression signatures associated with BP traits. Thirtyfour BP signature genes were identified at Bonferroni corrected p<0.05 (224 genes were identified at FDR<0.2, reported in the S1 Text). Thirteen BP signature genes replicated between the FHS and Illumina cohorts. The top BP signature genes identified in the FHS (55 genes for SBP  [3] are highlighted with blue triangles. eQTL-transcript pair genes that are BP signature genes from analysis of differential gene expression in relation to BP are depicted by red circles. B) Regional association plots for rs3184504 proxy QTLs that showing association with BP in ICBP GWAS [3]. −log10(p) indicated the −log10 transformed DBP association p values in ICBP GWAS [3]. Color coding indicates the strength (measured by r 2 ) of LD of each SNP with the top SNP (rs3184504). Five master trans-eQTLs (also BP GWAS SNPs) for BP signature genes are labeled in the figure. This figure was drawn by LocusZoom [32]. Among the 34 BP signature genes (at Bonferroni corrected p<0.05), only FOS [15] and PTGS2 [16] have been previously implicated in hypertension. We did not find literature * rs653178, intronic to ATXN2 and in tight linkage disequilibrium with rs3184504 (r 2 = 1), was also associated with BP in ICBP GWAS and all the 6 genes; + A proxy SNP rs4698412 at LD r 2 = 1 associated with the same trait; $ A proxy SNP rs4389526 at LD r 2 = 1 associated with the same trait; § indicated eQTL were identified from [12]. support for a direct role of the remaining signature genes in BP regulation. However, we found several genes involved in biological functions or processes that are highly related to BP, such as cardiovascular disease (GZMB, ANXA1, TMEM43, FOS, KCNJ2, PTGS2, and MCL1), angiogenesis (VIM and TIPARP), and ion channels (CD97, ANXA1, S100A10, PRF1, ANTXR2, SLC31A2, TIPARP, and KCNJ2). We speculate that these genes may be important for BP regulation, but further experimental validation is needed. Seven of the 34 signature genes, including KCNJ2, showed negative correlation of expression with BP. KCNJ2 is a member of the potassium inwardly-rectifying channel subfamily; it encodes the inward rectifier K+ channel Kir2.1, and is found in cardiac, skeletal muscle, and nervous tissue [17]. Most outward potassium channels are positively correlated with BP. Loss-offunction mutations in ROMK (KCNJ1, the outward potassium channel) are associated with Bartter's syndrome, and ROMK inhibitors are used in the treatment of hypertension [18,19]. Previous studies reported that greater potassium intake is associated with lower blood pressure [20,21,22,23]. These data suggest that KCNJ2 up-regulation may be a means of lowering BP.
By linking the BP signature genes with eQTLs and with BP GWAS results, we found several SNPs that are associated with BP in GWAS and that also are trans associated with several of our top BP signature genes. For example, rs3184504, a non-synonymous SNP located in exon 3 of SH2B3, is associated in GWAS with BP, coronary heart disease, hypothyroidism, rheumatoid arthritis, and type I diabetes [12]. rs3184504 is a common genetic variant with a minor allele frequency of approximately 0.47; the rs3184504-T allele is associated with an increment of 0.58 mm Hg in SBP and of 0.48 mm Hg in DBP [2]. rs3184504 is a cis-eQTL for SH2B3, expression of this gene was not associated with BP or hypertension in our data. However, rs3184504 also is a trans-eQTL for 6 of our 34 BP signature genes: FOS, MYADM, PP1R15A, TAGAP, S100A10, and FGBP2. These 6 genes are highly expressed in neutrophils [12], and are coexpressed. Prior studies have suggested an important role of neutrophils in BP regulation [24]. We speculate that these 6 BP signature genes, all driven by the same BP-associated eQTL, point to a critical and previously unrecognized mechanism involved in BP regulation. Further experimental validation is needed.
One limitation of our study is the use of whole blood derived RNA for transcriptomic profiling. GSEA showed that the top enriched biological processes for the differentially expressed BP genes include inflammatory response. Numerous studies have shown links between inflammation and hypertension [25,26,27]. The top ranked genes in inflammatory response categories provide a guide for further experimental work to recognize the contributions of inflammation to alterations in BP regulation. We speculate that using similar approaches in other tissues might identify additional differentially expressed BP signature genes.
In conclusion, we conducted a meta-analysis of global gene expression profiles in relation to BP and identified a number of credible gene signatures of BP and hypertension. Our integrative analysis of GWAS and gene expression in relation to BP can help to uncover the genetic and genomic architecture of BP regulation; the BP signature genes we identified may represent an early step toward improvements in the detection of susceptibility, and in the prevention and treatment of hypertension.

Study population and ethics statement
This investigation included six studies (the Framingham Heart Study (FHS), the Estonian Biobank (EGCUT), the Rotterdam Study (RS) [8], the InCHIANTI Study, the Cooperative Health Research in the Region of Augsburg (KORA F4) Study [9], and the Study of Health in Pomerania (SHIP-TREND) [10], each of which conducted genome-wide genotyping, mRNA expression profiling, and had extensive BP phenotype data. Each of the six studies followed the recommendations of the Declaration of Helsinki. The FHS: Systems Approach to Biomarker Research (SABRe) in cardiovascular disease is approved under the Boston University Medical Center's protocol H-27984. Ethical approval of EGCUT was granted by the Research Ethics Committee of the University of Tartu (UT REC). Ethical approval of the InCHIANTI study was granted by the Instituto Nazionale Riposo e Cura Anziani institutional review board in Italy. Ethical approval of RS was granted by the medical ethics committee of the Erasmus Medical Center. The study protocol of SHIP-TREND was approved by the medical ethics committee of the University of Greifswald. KORA F4 is a population-based survey in the region of Augsburg in Southern Germany which was performed between 2006 and 2008. KORA F4 was approved by the local ethical committees. Informed consent was obtained from each study participant.
Hypertension (HTN) was defined as SBP !140 mm Hg or DBP !90 mm Hg. We excluded individuals receiving anti-hypertensive treatment because of the possibility that some of the differentially expressed genes we identified would reflect treatment effects. The eligible study sample included 7017 individuals: 3679 from FHS, 972 from EGCUT, 604 from RS, 597 from InCHIANTI, 565 from KORA F4, and 600 from SHIP-TREND.

Identification and replication of differentially expressed genes associated with BP
The association of gene expression with BP was analyzed separately in each of the six studies (Equation 1). A linear mixed model was used in the FHS in order to account for family structure. Linear regression models were used in the other five studies. In each study, gene expression level, denoted by geneExp, was included as the dependent variable, and explanatory variables included blood pressure phenotypes (SBP, DBP, and HTN), and covariates included age, sex, body mass index (BMI), cell counts, and technical covariates. A separate regression model was fitted for each gene. The general formula is shown below, and the details of analyses for each study are provided in the S2 Text and S6 Table. geneExp The overall analysis framework is provided in S1 Fig. We first identified differentially expressed genes associated with BP (BP signature genes) in the FHS samples (Set 1) and attempted replication in the meta-analysis results from the Illumina cohorts (Set 2, see Methods, Meta-analysis). We next identified BP signature genes in the Illumina cohorts (Set 2), and then attempted replication in the FHS samples (Set 1). The significance threshold for pre-selecting BP signature genes in discovery was at Bonferroni corrected p = 0.05 (in FHS, corrected for 17,318 measured genes [17,873 transcripts], and in illumina cohorts, corrected for 12,010 measured genes [14,222 transcripts] that passed quality control). Replication was established at Bonferroni corrected p = 0.05, correcting for the number of pre-selected BP signatures genes in the discovery set. We computed the pi1 value to estimate the enrichment of significant p values in the replication set (the Illumina cohorts) for BP signatures identified in the discovery set (the FHS) by utilizing the R package Qvalue [11]. Pi1 is defined as 1-pi0. Pi0 value provided by the Qvalue package, represents overall probability that the null hypothesis is true. Therefore, pi1 value represents the proportion of significant results. For genes passing Bonferroni corrected p<0.05 in the discovery set for SBP, DBP and HTN, we calculated pi1 values for each gene set in the replication set.

Meta-analysis
We performed meta-analysis of the five Illumina cohorts (for discovery and replication purposes), and then performed meta-analysis of all six cohorts. An inverse variance weighted meta-analysis was conducted using fixed-effects or random-effects models by the metagen() function in the R package Meta (http://cran.r-project.org/web/packages/meta/index.html). At first, we tested heterogeneity for each gene using Cochran's Q statistic. If the heterogeneity p value is significant (p<0.05), we will use a random-effects model for the meta-analysis, otherwise use a fixed-effects model. The Benjamini-Hochberg (BH) method [28] was used to calculate FDR for differentially expressed genes in relation to BP following the meta-analysis of all six cohorts. We also used a more stringent threshold to define BP signature genes by utilizing p<6.5e-6 (Bonferroni correction for 7717 unique genes [7810 transcript] based on the overlap of FHS and illumina cohort interrogated gene sets).

Estimating the proportion of variance in BP attributable to BP signature genes
To estimate the proportion of variances in SBP or DBP explained by a group of differentially expressed BP signature genes (gene 1, gene 2, . . ., gene n), we used the following two models: Full model: covariates Null model: The proportion of variance in BP attributable to the group of differentially expressed BP signature genes (h 2 BP sig ) was calculated as: The proportion of the variance in BP phenotypes attributable to the FHS BP signature genes was estimated in the five Illumina cohorts, respectively, and then the average proportion values were reported. In turn, the proportion of the variance in BP phenotypes attributable to the Illumina BP signature genes was estimated in the FHS.
Identifying eQTLs and estimating the proportion of variance in gene expression attributable to single cis-or trans-eQTLs SNPs associated with altered gene expression (i.e. eQTLs) were identified using genome-wide genotype and gene expression data in all available FHS samples (n = 5257) at FDR<0.1 (Joehanes R, submitted, 2014, and a brief summary of methods and results are provided in the S2 Text). A cis-eQTL was defined as an eQTL within 1 megabase (MB) flanking the gene. Other eQTLs were defined as trans-eQTLs. We combined the eQTL list generated in the FHS with the eQTLs generated by meta-analysis of seven other studies (n = 5300) that were also based on whole blood expression [12].
For every BP signature gene, we estimated the proportion of variance in the transcript attributable to the corresponding cis-or trans-eQTLs (h 2 eQTL ) using the formula:

Functional category enrichment analysis
In order to understand the biological themes within the global gene expression changes in relation to BP, we performed gene set enrichment analysis [29] to test for enrichment of any gene ontology (GO) biology process [30] or KEGG pathways [31]. "Metric for ranking gene" parameters were configured to the beta value of the meta-analysis, in order to look at the top enriched functions for BP associated up-regulated and down-regulated gene expression changes respectively. One thousand random permutations were conducted and the significance level was set at FDR 0.25 to allow for exploratory discovery [29].

Members of International Consortium for Blood Pressure GWAS (ICBP)
Steering And in turn, we replicated the BP signature genes identified in Illumina cohorts in FHS cohort. Fourth, we conducted a meta-analysis in the six cohorts and reported the BP signature genes passing Bonferroni corrected p<0.05 (corrected for 7717 genes). And finally, we cross-analyzed the BP signature genes with blood eQTLs as well as with BP GWAS results to identify the BP signature genes having BP GWAS eQTLs.