Original Research ARTICLE
Integrative Genomics Analysis Unravels Tissue-Specific Pathways, Networks, and Key Regulators of Blood Pressure Regulation
- 1Department of Integrative Biology and Physiology, University of California, Los Angeles, Los Angeles, CA, United States
- 2Department of Genetics, Department of Developmental Biology, Stanford University Medical Center, Stanford, CA, United States
- 3The National Heart Lung and Blood Institute's Framingham Heart Study, Framingham, MA, United States
- 4The Population Sciences Branch and the Division of Intramural Research, National Heart, Lung and Blood Institute, Bethesda, MD, United States
Blood pressure (BP) is a highly heritable trait and a major cardiovascular disease risk factor. Genome wide association studies (GWAS) have implicated a number of susceptibility loci for systolic (SBP) and diastolic (DBP) blood pressure. However, a large portion of the heritability cannot be explained by the top GWAS loci and a comprehensive understanding of the underlying molecular mechanisms is still lacking. Here, we utilized an integrative genomics approach that leveraged multiple genetic and genomic datasets including (a) GWAS for SBP and DBP from the International Consortium for Blood Pressure (ICBP), (b) expression quantitative trait loci (eQTLs) from genetics of gene expression studies of human tissues related to BP, (c) knowledge-driven biological pathways, and (d) data-driven tissue-specific regulatory gene networks. Integration of these multidimensional datasets revealed tens of pathways and gene subnetworks in vascular tissues, liver, adipose, blood, and brain functionally associated with DBP and SBP. Diverse processes such as platelet production, insulin secretion/signaling, protein catabolism, cell adhesion and junction, immune and inflammation, and cardiac/smooth muscle contraction, were shared between DBP and SBP. Furthermore, “Wnt signaling” and “mammalian target of rapamycin (mTOR) signaling” pathways were found to be unique to SBP, while “cytokine network”, and “tryptophan catabolism” to DBP. Incorporation of gene regulatory networks in our analysis informed on key regulator genes that orchestrate tissue-specific subnetworks of genes whose variants together explain ~20% of BP heritability. Our results shed light on the complex mechanisms underlying BP regulation and highlight potential novel targets and pathways for hypertension and cardiovascular diseases.
Hypertension or elevated blood pressure (BP) is among the most prevalent, treatable risk factors for coronary artery disease (CAD), heart failure, and stroke. Systolic (SBP) and diastolic blood pressure (DBP) traits are highly heritable, with the total heritability for European/African ancestry individuals estimated to be 20/27% and 39/50% for SBP and DBP, respectively (1). Large-scale GWAS have successfully established ~800 genetic loci for SBP, DBP, and hypertension in multiple ethnic groups (2) and the number will continue to rise as the sample size increases. Some of the loci contain genes previously known or suspected to regulate BP (such as ADM and NPPA) and the remaining loci are novel.
Despite the successful identification of novel genetic loci associated with BP regulation through GWAS, pinpointing the causal genes and underlying mechanisms mediating the effects of these loci is not straightforward (3). Integration of genetic information with functional data, such as genetic variants associated with altered gene expression, or expression quantitative trait loci (eQTLs), is of critical importance to pinpoint the causal genes and their associated pathogenic mechanisms (4). Another critical gap is that the known ~800 GWAS BP loci with significance together only explains a small fraction of total BP heritability (2), a phenomenon called the missing heritability or the dark matter. Recent evidence suggests that the missing heritability can be explained by numerous additional genetic loci with moderate or subtle effects that are well under the genome-wide significance level in GWAS as well as the interactions between multiple genetic loci (5). Such multigenic interactions can be captured using gene regulatory networks, and our recent modeling of blood pressure networks using whole blood transcriptome data has highlighted the importance of inflammatory pathways in blood pressure regulation (6). However, BP regulation most likely involves many more genes functioning in numerous biological processes in diverse tissues such as kidney, heart, liver, and the vasculature (2, 7). Systematic modeling of multidimensional omics data that capture tissue-specific blood pressure networks informed by genetic loci of strong, moderate, to subtle effects will likely provide a more comprehensive understanding of BP mechanisms.
Here we employed an integrative genomics strategy leveraging multiple genetic and genomic datasets. Previous applications of this strategy have successfully identified novel mechanisms of CAD and other complex diseases (4, 8, 9). Once integrated, these multidimensional datasets have the potential to delineate the genes, pathways, and epistatic interaction subnetworks associated with BP that are informed by genetic signals with a wide range of effect sizes.
Overall Analysis Design
Figure 1 shows the general flowchart of the study. First, we utilized the human genetic association data (i.e., GWAS) from the ICBP, which provides the full spectrum of statistical associations between SNPs and clinically measured SBP and DBP (not limiting to the top significant loci). Second, we curated eQTLs (SNPs under eQTLs are defined as eSNPs) from diverse tissues, which have been confirmed to be enriched for complex disease loci (10) and provide functional support for tissue-specific connections between SNPs and genes in a data-driven manner. To further enrich functional annotation, we incorporated information from the Encyclopedia of DNA Elements (ENCODE) studies (11). Third, to provide a holistic view of the organization of genes and reveal the most important regulatory hubs in a given tissue, we included knowledge-based metabolic and signaling pathways and data-driven gene networks from various tissues to improve the detection of multigenic disease processes. The general analytical pipeline, Mergeomics, has been implemented as an open-access web-server (http://mergeomics.research.idre.ucla.edu/) (12) as well as an open-access R Bioconductor package (https://bioconductor.org/packages/release/bioc/html/Mergeomics.html) (13). Ethical standards and procedures were followed throughout the study.
Figure 1. The integrative genomics framework for identifying genetically informed biological processes and networks and for prioritizing key drivers of BP regulation.
GWAS of SBP, DBP, Hypertension, and CAD
The summary statistics of GWAS for SBP, DBP, and hypertension was obtained from the ICBP, which was formed by two consortia, the CHARGE-BP consortium (Cohorts for Heart and Aging Research in Genomic Epidemiology—blood pressure) and the GBPGEN consortium (Global Blood Pressure Genetics Consortium) (14). (dbGaP accession: phs000585.v2.p1). The study is comprised of 200,000 individuals of European descent in a multi-stage design from 29 studies. More than 906,600 SNPs were genotyped using Affymetrix Genome-Wide Human SNP Array 6.0. Imputation was further carried out to obtain information for up to 2.6 million SNPs using the HapMap CEU (Utah residents with ancestry from northern and western Europe) panel. SNPs with minor allele frequency (MAF) <1% were removed. Finally, a total of ~2.5 million SNPs tested for association with systolic and diastolic blood pressure were used in our study. The 1,000 Genomes-based CAD GWAS dataset was retrieved from CARDIoGRAMplusC4D Consortium (http://www.cardiogramplusc4d.org/media/cardiogramplusc4d-consortium/data-downloads/cad.additive.Oct2015.pub.zip) (15). The hypertension and CAD GWAS were used to connect the SBP/DBP related findings to disease conditions. All statistical association p-values for all SNPs, regardless of significance level, were used in our downstream analysis.
Mapping SNPs to Genes and Removal of SNPs in Linkage Disequilibrium (LD)
We used different mapping methods that are based on (a) chromosomal distance, (b) eQTLs, or (c) ENCODE to link GWAS SNPs to their potential target genes.
(a) We used a standard distance-based approach where a SNP was mapped to a gene if within 50 kb of the respective gene region.
(b) The expression levels of genes can be seen also as quantitative traits in GWAS. Hence, it is possible to determine eQTLs and the expression SNPs (eSNPs) within the eQTLs that provide a functionally motivated mapping from SNPs to genes. Moreover, the eSNPs within the eQTL are specific to the tissue where the gene expression was measured and can therefore provide mechanistic clues regarding the tissue of action when intersected with BP-associated SNPs. Results from eQTL studies in human adipose tissue, artery, liver, brain, and blood (16– 28) were combined with the eQTLs from the same tissues in the GTEx database (29). In addition, we included eQTLs from 44 additional tissues from studies with smaller sample sizes including the GTEx (30) (https://storage.googleapis.com/gtex_analysis_v7/single_tissue_eqtl_data/GTEx_Analysis_v7_eQTL_all_associations.tar.gz) and eQTLs in kidney at FDR <5% (Kim et al., unpublished data) but the statistical power of these eQTL sets was limited. We included both cis-eSNPs (within 1 Mb distance from gene region) and trans-eSNPs (beyond 1 Mb from gene region), at a false discovery rate (FDR) <5%. However, for these tissues, the number of eSNPs after LD trimming was small and lacked power to detect BP-associated signals. To improve statistical power, we included eSNPs at P < 1.0E-5 from these 44 tissues as “suggestive” eQTL sets.
(c) In addition to eQTLs and distance-based SNP-gene mapping approaches, we integrated functional data sets from the Regulome database (11) which annotates SNPs in regulatory elements in the Homo sapiens genome based on the results from the ENCODE studies (31).
Using the above mapping approaches, the following sets of SNP-gene mappings: eSNP adipose, eSNP artery, eSNP liver, eSNP blood, eSNP brain, eSNP all (i.e., combing all the tissue-specific eSNPs above), Distance (chromosomal distance-based mapping), Regulome (ENCODE-based mapping), Combined (combing all the above methods), and 44 suggestive eQTL sets.
We observed a high degree of LD in the eQTL, Regulome, and distance-based SNPs, and this LD structure may cause artifacts and biases in the downstream analysis. For this reason, we devised an algorithm to remove SNPs in LD while preferentially keeping those with a strong statistical association with SBP/DBP. We chose a LD cutoff (r2 < 0.7) to remove redundant SNPs in high LD.
Knowledge-Based Biological Pathways
We included 1,827 canonical pathways from the Reactome (Version 45), Biocarta, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) databases (32, 33). In addition to the curated pathways, we constructed two positive control BP pathways based on known BP loci (p < 1.0E-5) and candidate genes from the GWAS Catalog (GWAS p < 5.0E-8) (34) for SBP and DBP separately. We also curated hypertension/CAD positive control gene sets based on GWAS Catalog (p < 1.0E-5). In addition, the CAD positive control genes were complemented with the CADgene V2.0 database, which contains 583 CAD related genes and detailed CAD association information from about 5,000 publications. These gene sets serve as positive controls to validate our computational method.
Data-Driven Modules of Co-expressed Genes
Beside the canonical pathways, we used co-expression modules that were derived from a collection of genomics studies of liver, adipose tissue, aortic endothelial cells, brain, blood, kidney, and muscle (GEO accession numbers: GSE7965, GSE25506, GSE9588, GSE24335, GSE20142, GSE20332, GSE22070, GSE2814, GSE3086, GSE2814, GSE3086, GSE3087, and GSE3088, and GSE30169) (16–19, 21, 22, 35–38). For each dataset, we extracted the normalized gene expression profile and reconstructed co-expression networks using the established WGCNA R package (39). Modules with size smaller than 10 genes were excluded to avoid statistical artifacts, yielding a total of 2,705 co-expression modules in this study. We included these tissue-specific co-expression networks to confirm whether known tissue types for BP could be objectively detected and whether any additional tissue types are also important for BP regulation.
These data-driven modules along with the knowledge-driven pathways in the previous section were used together to capture gene sets containing functionally related genes in a wide variety of tissue and functional settings.
Marker Set Enrichment Analysis (MSEA)
We applied MSEA (13) to identify pathways/co-expression modules that demonstrate enrichment for genetic association with SBP, DBP, hypertension, or CAD using the same parameters. MSEA employs a chi-square like statistic with multiple quantile thresholds to assess whether a pathway or co-expression module shows enrichment of disease SNPs compared to random chance based on the full spectrum of association statistics for each GWAS dataset. For each pathway or co-expression module, 10,000 permuted gene sets were generated, and enrichment P-values were determined from a Gaussian distribution approximated using the enrichment statistics from the permutations as detailed in Shu et al. (13). Benjamini-Hochberg FDR was estimated across all pathways and modules tested for each GWAS. Pathways or co-expression modules with FDR < 5% in at least one SNP-gene mapping were considered statistically significant.
Construction of Independent Supersets and Confirmation of BP Association
Because the significant pathways or co-expression modules were from multiple sources, there were overlapping or nested structures among the gene sets. To make the results more meaningful, we constructed independent supersets that captured the core genes from groups of redundant pathways and co-expression modules (Figure 1). We merged the 42 common pathways associated with SBP/DBP using a merging algorithm in Mergeomics (13). After merging, we annotated each superset based on function enrichment analysis of the known pathways from the Gene Ontology and KEGG databases (Bonferroni-corrected P < 0.05 in Fisher's exact test). The supersets were given a second round of MSEA to confirm their significant association with BP using Bonferroni corrected P < 0.05 as the cutoff.
Key Driver Analysis (KDA)
We used a KDA algorithm (40) to identify potential key driver (KD) genes of the BP-associated supersets. KDA overlays BP-associated gene sets that were discovered by MSEA onto graphical network models detailing molecular interactions among genes to see if a particular subnetwork was significantly enriched for disease genes, using a chi-square like statistic analogous to the one used for MSEA. Statistical significance of KDs was estimated by permuting the gene labels in the network for 10,000 times and estimating the P-value based on the null distribution. To control for multiple testing, stringent adjustment (FDR < 1%) was used to focus on the top robust KDs. Graphical networks used sources including tissue-specific Bayesian regulatory networks available for seven tissues (including cardiac muscle, artery, adipose, blood, liver, brain, and kidney, described previously (7) and a protein-protein interaction network from Human Protein Reference Database (HPRD) database (41).
To cross-validate the top ranked KDs in silico, we used literature-mining methods (PolySearch, COREMINE, and T-HOD), searched mouse phenome database (http://www.informatics.jax.org/), and examined their association with BP in the latest GWAS to assess supporting evidence for their role in BP regulation. We also retrieved gene essentiality information (42) from an exome sequencing study of 60,706 humans (Supplementary Table 5) (42) to evaluate the biological importance of the predicted KDs. This study considered genes with a probability of being loss-of-function (LoF) intolerant >= 0.9 to be essential genes, where counts of the observed and expected variants were utilized to predict whether a given gene is significantly intolerant to LoF. Genes were classified into three groups within the context of tolerance to LoF, null/complete toleration, recessive/heterozygous toleration, and haploinsufficient/intolerant.
Heritability Estimates of GWAS Hits and KD Subnetworks
We used the Heritability Estimator from Summary Statistics (HESS) (43) to estimate the total genome-wide SNP heritability (h2) as well as the portion explained by the KD subnetworks of SBP/DBP. We define GWAS hits as SNPs with p < 5.0E-8. To correct for any potential bias in the large numbers of KD subnetwork genes for SBP/DBP, we implemented a permutation strategy by generating 10,000 random gene sets of matching sizes for the top KD subnetworks of SBP/DBP to derive permutation-based significant P-values for the heritability estimates, where PKD = (number of > )/10000. The reference heritability of SBP/DBP was reported by a recent study (44), which studied 2,889 twin pairs not on any BP-lowering therapy from the Twins UK cohort, and estimated the additive genetic variance for baseline BP, long-term average BP, BP trajectory (rate of change of BP in mmHg/year) and BP variability (coefficient of variation and average real variability over an average of 3.2 visits).
Identification of BP-Associated Pathways and Co-expression Network Modules
Out of the 4,532 gene sets (1,827 canonical pathways and 2,705 data-driven co-expression modules), we identified 96 and 78 that were significantly associated with SBP and DBP, respectively, in at least one SNP-to-gene mapping approach (FDR < 5%; Supplementary Table 1). As expected, the predefined positive controls based on GWAS catalog for SBP and DBP were among the top signals for their corresponding traits.
Among the significant signals, 42 gene sets were shared between SBP and DBP (Table 1). These included the previously documented hypertension processes, such as angiotensin II (45), cell-cell junction organization (46), and muscle contraction (47). Other plausible pathways included NOTCH signaling, platelet production, insulin secretion/signaling, immune and inflammation, corticosteroids, and cell cycle pathways (Table 1). Among the gene sets unique to either SBP or DBP were “Wnt signaling” and “mammalian target of rapamycin (mTOR) signaling” for SBP, and “cytokine network” and “tryptophan catabolism” for DBP (Supplementary Table 1).
The use of tissue-specific eQTLs allowed us to implicate the potential tissues where the BP-associated processes may function (Supplementary Table 1). Out of the 44 tissue-specific eQTL sets, those from the adipose, liver, blood, and brain tissues appeared to be more informative, although this could be due to the higher statistical power resulting from the abundance of eQTLs from these well-studied tissues. Adipose eQTLs were shown to be informative for the majority of the pathways or co-expression modules identified, and insulin secretion/signaling, autophagy, cell cycle and protein modification processes were mainly identified when adipose eQTLs were used. Use of smaller eQTLs datasets such as those from heart and kidney tissues, which were previously implicated in BP (7), did not yield significant BP-associated pathways. This is likely due to the limited statistical power instead of lacking biological relevance. Indeed, by lowering the stringency of eQTLs from FDR < 5% to P < 1.0E-5 to include more eQTLs from these tissues, we found suggestive pathways (e.g., response to wounding, cell-cell adhesion, and G-protein coupled receptor signaling) informed by eQTLs from artery, heart, pancreas, and testis tissues (FDR < 5%; Supplementary Table 2).
Construction of Supersets for SBP/DBP and Functional Annotation
We next focused on the 42 significant gene sets (pathways/co-expression modules) shared between SBP and DBP as they reflect reproducible signals for BP regulation. To minimize redundancy in the biological processes captured in these gene sets, we merged 35 overlapping gene sets into 12 independent supersets and kept the other 7 non-overlapping gene sets intact (Table 1). The resulting 19 non-overlapping gene sets represent a diverse range of molecular pathways, including NOTCH signaling, protein catabolism, cell adhesion and junction, cardiac and smooth muscle contraction, phosphatidylinositol signaling, insulin secretion/signaling, immune/inflammation, integrin signaling, gene expression regulation, cell cycle, angiotensin II-induced JNK activation, platelet, autophagy, and corticosteroids and cardioprotection. We confirmed that these merged supersets still captured the BP-relevance demonstrated by their subcomponents by performing a second round of MSEA (Supplementary Table 3).
Identification and Prioritization of key Regulators in the BP-Associated Supersets
To identify and prioritize the central regulatory components (termed as key drivers, KD) among the large number of genes in the BP supersets, we performed a key driver analysis (KDA; details in Methods) (13) on the 19 supersets shared between SBP and DBP using 8 graphical networks including 7 tissue-specific Bayesian networks [cardiac muscle, artery, adipose, blood, liver, brain, and kidney, as these tissues have been implicated in BP by previous (7) and our studies] and a PPI network. For each superset in each network, we focused on the top five ranked KDs satisfying FDR < 1%, yielding 295 unique KD genes, among which 29 were shared by >=2 supersets and 6 (COX5B, FN1, COL4A2, COX4I1, NDUFS3, and HLA-C) were shared in >=3 networks (Supplementary Table 4). The low number of shared KDs for the BP-associated gene sets between tissue-specific networks suggests tissue-specific regulation of BP genes and pathways.
We cross-checked these top KD genes for their role in BP and functional significance through a comprehensive in silico analysis using multiple literature mining tools, databases of gene knockout and mutation mouse models, and recent GWAS findings (Supplementary Table 5). Our search revealed both known KDs showing connections to hypertension or relevant conditions in one or more types of studies (99 KDs; such as GNAS, SLC2A3, IRS1, ADM, and SERPINE1) and relatively novel KDs in BP studies (such as SPTBN1 and GNB1). Moreover, the top KDs were significantly enriched with essential human genes (41) (Supplementary Table 5; Fisher Exact test P = 5.12e-7).
Additionally, a number of suggestive genes arose from the vascular tissue data, most notably those implicated in extracellular matrix (ECM) homeostasis and smooth muscle contraction with emphasis on angiotensin signaling. The genes highlighted as having specific effects on the ECM are TIMP2, PSAP, FN1, VCAN, and LOX. Those implicated in smooth muscle contraction include CALD1, DUSP5, and TAGLN.
Top KD Subnetworks of BP Regulation
We retrieved the top KD subnetworks in BP-relevant tissues including cardiac muscle, artery, adipose, blood, liver, brain, and kidney, as well as in the PPI network. As shown in Supplementary Figure 1, BP-associated processes are closely orchestrated by KDs in the adipose, artery, blood, and liver subnetworks. The top KD subnetworks demonstrated high tissue specificity, for example, subnetworks involved in insulin secretion/signaling were only overrepresented in adipose and liver tissues (Supplementary Figure 1). Similar results were found for the BP subnetworks from the other tissues (the detailed interactions in the eight KD networks are listed in Supplementary Table 6).
BP Heritability Explained by the KD Subnetworks
We estimated the total genome-wide SNP heritability (h2) and the fractions of heritability explained by the top KD subnetworks using HESS (43) based on the ICBP GWAS summary statistics. We found that while the significant GWAS SNPs (14) can only explain 4.82% and 4.46% of trait variance in SBP and DBP, respectively, which agreed with the estimates from a recent BP GWAS study (2), KDs and genes in the top KD subnetworks in the various tissues explained much higher proportions of heritability (Table 2). For example, the top 54 KDs and their subnetwork genes in the adipose tissue (Supplementary Figure 1) can explain 19.8% and 21.8% of the heritability of SBP and DBP, respectively. This increase in heritability estimates for the KD subnetworks over GWAS top hits is less likely to be driven by the number of network genes, as random subnetworks containing matching numbers of genes explained much lower proportions of heritability (P < 0.0001).
Overlap of BP Pathways/Networks With Those of Hypertension and CAD
We explored the relationship between BP-associated pathways/networks with those for hypertension and CAD. Out of the 42 SBP/DBP common pathways (Table 1), five showed significant enrichment for hypertension GWAS (14) signals (FDR < 0.05), including “Antigen processing and presentation,” “Insulin signaling pathway,” and “Integration of energy metabolism,” and the two positive control sets for SBP and DBP. Six of the SBP/DBP common pathways also showed significant enrichment for CAD GWAS (15) signals, including “ALK in cardiac myocytes,” “Factors involved in megakaryocyte development and platelet production,” “Integrin cell surface interactions,” “Meiotic recombination,” and “Antigen processing and presentation”. In addition, 19 out of the top BP KDs (Supplementary Table 5; e.g., SHC1, FN1, APOB, COL4A1, RELA, and ADM) were CAD susceptibility genes based on GWAS studies (8, 48). Notably, 10 KD genes highlighted by our analysis ( e.g., TGFBR2, LAMC1, KIF15, and RBPMS) have only just been implicated as causal genes by the most recent blood pressure GWAS study (2). Moreover, many hypertension/CAD genes (Supplementary Table 7) are in the subnetworks of the top KDs (Supplementary Table 6). These molecular level overlaps support the strong mechanistic connections between BP and incident diseases.
While elevated BP has a significant genomic contribution, this phenotype has been historically multidimensional in its physiological induction. Previous high-throughput GWAS and transcriptome studies have revealed a plethora of genomic changes contributing to BP in different populations. However, an integrative systems analysis fully utilizing the complementarity of diverse omics data has not been conducted to capture a comprehensive, tissue-specific view of BP regulation. Due to the highly interconnected nature of the genome and molecular regulatory systems, finding the gene sets most relevant in pathogenesis through the inherent noise of each expressed gene is challenging. Based on the predictions made via the “omnigenic” model, the vast majority of heritability arises from the multitude of genes which have indirect effects on disease through their association with central regulators (49). Our previous studies on other complex diseases indeed support that GWAS genes are mostly peripheral genes in gene networks that are regulated by key driver genes (8, 9).
Here, to investigate the gene regulatory networks and pathways governing BP regulation, we integrated the full association spectrum of BP GWAS (not limited to top GWAS hits, but included moderate and subtle signals as well), functional genomics information (eQTLs and ENCODE), knowledge-driven pathways, and data-driven networks to uncover biological processes and tissue-specific networks mediating the actions of BP GWAS signals. Our study revealed a diverse set of pathways significantly associated with SBP/DBP based on genetic evidence. These pathways are connected in tissue-specific networks via central regulators, both known (e.g., GNAS, SLC2A3, and IRS1) and novel (e.g., SPTBN1 and GNB1), revealing critical interactions and regulatory cascades underlying BP maintenance. The KDs and subnetworks collectively explain ~20% heritability for SBP/DBP, substantiating the much improved power of our network approach in capturing key BP genes and processes compared to conventional approaches. Furthermore, many of the BP pathways and networks were shared with CAD, thus shedding light on the mechanistic connections between CAD and its clinical risk factor hypertension.
Compared to several previous integrative network studies of BP involving whole blood transcriptome and BP GWAS, (20, 33) our study incorporated eQTL and network information from a comprehensive list of tissue types that capture pathways or networks that are either systemic or tissue-specific, revealing complex, systems-level regulation of BP (summarized in Figure 2). The detected processes can be categorized into a number of general classes, such as muscle contraction, immune and inflammation, cell adhesion and junction, protein catabolism, and angiotensin pathway. Therefore, our omics-driven systems biology study depicts a much more comprehensive landscape for BP regulation that includes both known processes targeted by antihypertensive drugs (50) (angiotensin II receptor blockers targeting the angiotensin pathway) and additional processes that may lead to novel therapeutic strategies for hypertension and CAD.
Figure 2. Pathways associated with systolic and/or diastolic blood pressure with FDR determined by MSEA analysis from various tissues: (A) Adipose. (B) Blood. (C) Brain. (D) Liver (E) Artery.
The use of tissue-specific functional genomics information allowed us to detect potential tissue specificity in the BP processes identified in our study. For example, autophagy regulation and insulin signaling showed genetic association with BP traits when adipose tissue eQTLs were used (Table 1); our network models also suggest that the insulin signaling pathway in adipose tissue closely connects with other critical BP processes, such as cell-cell adhesion and cell junction organization (Figure 2, Supplementary Figure 1). The link of adipose tissue biology and insulin resistance to blood pressure homeostasis has been noted before (51, 52). Our results suggest that autophagy and insulin signaling are adipose pathways contributing to BP regulation. The connection between autophagy and insulin signaling is supported by previous evidence that adipose-specific deletion of autophagy modulators (e.g., ATG7) reduces white adipose mass, increases the proportion of brown adipocytes and subsequently promotes the oxidation of free fatty acids, leading to enhanced insulin sensitivity (53). A recent study unraveled that spermidine, an anti-aging molecule, reduced SBP and prevented cardiac hypertrophy and a decline in diastolic function through autophagy-related protein ATG5 (54). Another study identified autophagy modulators as potential therapeutic targets of FDA-approved antihypertensive drugs, such as the imidazoline receptors and l-type Ca2+ channels (55). These lines of evidence substantiate that autophagy and insulin resistance processes in adipose may contribute to BP regulation.
Beside the common pathways between SBP and DBP discussed above, the trait-specific processes are also of value to help understand the differential regulatory mechanisms between the two traits. For instance, we found that Wnt signaling is more associated with SBP while tryptophan catabolism tends to be more involved in DBP regulation (Supplementary Table 1). In the central nervous system of spontaneously hypertensive rats, Wnt signaling regulates SBP by downregulating a glycogen synthase kinase 3β-mediated pathway to enhance insulin signaling (56). As for DBP-specific tryptophan catabolism, the central enzyme indoleamine 2,3-dioxygenase (IDO) in the pathway metabolizes the essential amino acid tryptophan to kynurenine. There was an inverse association between IDO activity and DBP but no association with SBP (57). Further investigations of these trait-specific pathways are warranted.
In addition to retrieving BP-relevant processes and pathways in a tissue-specific fashion, our network modeling detected potential KDs of these processes. These KDs include both known BP genes (e.g., GNAS, SLC2A3, and IRS1) and novel genes (e.g., SPTBN1 and GNB1). The KD SPTBN1 in “Positive Control,” is associated with multiple BP GWAS susceptibility genes (including CPEB4, FBXL19, NPR3, and HIPK2) in the adipose network (Supplementary Figure 1). It encodes βII-spectrin, which has been proven to play an essential role in the regulation of bone mineral density (58). On the other hand, bone mineral density was significantly lower in hypertensive subjects compared with normotensive subjects and was inversely correlated with SBP (59). Another novel KD GNB1 in “Insulin regulation” encodes a G protein β subunit, multiple mutations of which were found to affect the protein interface that binds Gα subunits (GNAS), which has been genetically and clinically associated with hypertension (60). These lines of evidence support the potential importance of the novel KDs in BP regulation.
Several interesting observations emerged from relationship between the top KDs we found and BP GWAS signals. First, we observed a lack of GWAS signals in the predicted KDs for direct genetic association with BP, with only 22 (e.g., PTPN11, INSR, SLC2A4, and GNAS) out of the 295 predicted KDs to be candidate BP GWAS genes (2, 61). This may be a result of negative selection pressure because of the critical roles of the KDs in tissue-specific networks. In support of this hypothesis, the top KDs we found were significantly overrepresented within biologically essential genes which demonstrates few protein-truncating variants due to the critical nature of their functions (42). This helps explain why KDs are commonly missed in GWAS and the power of our network approach in uncovering potential missing disease genes and processes. The lack of KD signals in GWAS does not undermine their critical regulatory role in disease processes, as the KD subnetworks collectively explain significantly more heritability than the subnetworks of random genes as well the top BP GWAS loci for both SBP and DBP (Table 2). These results substantiate the improved ability of network approach in capturing key BP genes and processes compared to conventional approaches (Table 2). Additionally, these data indicate that the KDs regulate both the significant BP GWAS genes which likely have stronger influence on pathophysiological outcomes and the genes with more subtle effects which account for the missing heritability. Furthermore, our KD subnetworks retrieved based on a previous BP GWAS meta-analysis with a smaller sample size, was able to predict causal BP genes such as KIF15, LAMC1, TGFBR2, and RBPMS which hadn't been highlighted until the latest BP GWAS (2) with a larger sample size and power. This strongly supports the accuracy and predictive power of our computational pipeline. These genes are linked to pathways potentially associated with mechanisms modulating blood pressure. Particularly, LAMC1 coding for Laminin Subunit Gamma 1, has been implicated in the regulation of vascular lumen development (62). Also, proteomic analysis suggests this subunit is involved in ECM remodeling during venous hypertension in varicose veins (63). Furthermore, TGFBR2 and RBPMS show a high degree of expression in smooth muscle cells, the former heavily implicated with predisposition for arterial aneurysms which suggests dysfunctions of the arterial wall, while the latter is involved in regulation of smooth muscle plasticity during development (64, 65).
Moreover, we had a number of suggestive genes arise from our vascular tissue data, highlighted as having specific effects on the ECM including TIMP2, PSAP, FN1, VCAN, and LOX, each of which have been previously implicated in hypertension (63, 66–68). With the ECM being crucial for maintaining the structural integrity of the vessel wall, it seems intuitive that genes affecting ECM homeostasis may contribute to hypertension. One such gene is TIMP2, coding for tissue inhibitor of metalloproteinases 2, which is part of a peptidase family involved in the degradation of the ECM, linked with resistant hypertension and arterial stiffness (66). Additionally, our vascular tissue data suggested genes implicated in smooth muscle contraction such as CALD1, DUSP5, and TAGLN (69–72).
Despite the comprehensive data integration and discoveries discussed above, there are limitations in our study. First, mapping GWAS SNPs to candidate genes is not straightforward. Chromosome location-based mapping lacks direct evidence for the functions of the reported genes, whereas functional data-supported mapping suffers from incomplete coverage of tissue and lack of power in identifying weak cis-association and trans-regulation. To address these issues, we incorporated different SNP-gene mapping approaches (i.e., chromosomal location vs. eQTLs) and prioritized the tissue-specific pathways shared by SBP and DBP. Second, although we collected eQTL datasets from >40 unique tissues/cells from the GTEx database and other studies, the sample sizes for certain tissues or cells are small, resulting in small numbers of eQTLs for these cells to tissue types to be incorporated in our analysis. This limits the statistical power when used separately and most likely contributes to the lack of significant signals from many of the eQTL sets used. To alleviate this power issue, we pooled the eQTLs from related tissues or cell types in pathway analysis, which may mask certain tissue-specific signals. Third, although our findings provide multiple lines of in silico evidence to support the importance of the KDs (e.g., literature mining and BP phenotypes from mouse models), experimental validation of the novel KDs for their roles in regulating the BP GWAS genes, networks, and disease development is necessary in future studies.
Overall, high BP is a major risk factor for cardiovascular disease and has a substantial genetic contribution. Our integrative approach utilizing GWAS, functional genomics, and network modeling revealed a comprehensive system view of biological processes, tissue-specific networks, and regulators that contribute to BP regulation. The common pathways between SBP and DBP are significantly enriched for hypertension/CAD GWAS signals and the top KD subnetworks contain numerous candidate CAD causal genes (Supplementary Table 7). Collectively, these findings provide molecular level mechanistic support for the tight connection between BP and CAD risk. The processes and regulators identified from our study may open new avenues for BP-lowering and CAD therapeutics.
YZ researched and analyzed all data and drafted the manuscript. MB and LS participated in data analysis and writing of the manuscript. XS, CL, and IA contributed toward the data analysis. SK, TH, and DL contributed data and assisted with results interpretation. XY supervised study design, data analysis, and manuscript writing. All authors reviewed and edited manuscript.
XY is supported by the National Institute of Health R01DK104363 and R21NS103088, American Heart Association Scientist Development Grant (13SDG17290032), and Transatlantic Networks of Excellence Award (12CVD02) from Foundation Leducq. YZ is supported by an American Heart Association Postdoctoral Fellowship (16POST31160044). DL is funded by National Institutes of Health contract N01-HC-25195 and by the Division of Intramural Research, National Heart, Lung, and Blood Institute, National Institutes of Health, Bethesda, MD.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We thank the ICBP for providing blood pressure association results.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcvm.2019.00021/full#supplementary-material
1. Salfati E, Morrison AC, Boerwinkle E, Chakravarti A. Direct estimates of the genomic contributions to blood pressure heritability within a population-based cohort (ARIC). PLoS ONE. (2015) 10:e0133031. doi: 10.1371/journal.pone.0133031
2. Evangelou E, Warren HR, Mosen-Ansorena D, Mifsud B, Pazoki R, Gao H, et al. Genetic analysis of over 1 million people identifies 535 new loci associated with blood pressure traits. Nat Genet. (2018) 50:1412–25. doi: 10.1038/s41588-018-0205-x
4. Mäkinen VP, Civelek M, Meng Q, Zhang B, Zhu J, Levian C, et al. Integrative genomics reveals novel molecular pathways and gene networks for coronary artery disease. PLoS Genet. (2014) 10:e1004502. doi: 10.1371/journal.pgen.1004502
6. Huan T, Meng Q, Saleh MA, Norlander AE, Joehanes R, Zhu J, et al. Integrative network analysis reveals molecular mechanisms of blood pressure regulation. Mol Syst Biol. (2015) 11:799. doi: 10.15252/msb.20145399
7. Greene CS, Krishnan A, Wong AK, Ricciotti E, Zelaya RA, Himmelstein DS, et al. Understanding multicellular function and disease with human tissue-specific networks. Nat Genet. (2015) 47:569. doi: 10.1038/ng.3259
8. Shu L, Chan KH, Zhang G, Huan T, Kurt Z, Zhao Y, et al. Shared genetic regulatory networks for cardiovascular disease and type 2 diabetes in multiple populations of diverse ethnicities in the United States. PLoS Genet. (2017) 13:e1007040. doi: 10.1371/journal.pgen.1007040
9. Krishnan KC, Kurt Z, Barrere-Cain R, Sabir S, Das A, Floyd R, et al. Integration of multi-omics data from mouse diversity panel highlights mitochondrial dysfunction in non-alcoholic fatty liver disease. Cell Syst. (2018) 6:103–15. doi: 10.1016/j.cels.2017.12.006
10. Zhong H, Beaulaurier J, Lum PY, Molony C, Yang X, MacNeil DJ, et al. Liver and adipose expression associated SNPs are enriched for association to type 2 diabetes. PLoS Genet. (2010) 6:e1000932. doi: 10.1371/journal.pgen.1000932
11. Boyle AP, Hong EL, Hariharan M, Cheng Y, Schaub MA, Kasowski M, et al. Annotation of functional variation in personal genomes using RegulomeDB. Genome Res. (2012) 22:1790–7. doi: 10.1101/gr.137323.112
12. Arneson D, Bhattacharya A, Shu L, Makinen VP, Yang X. Mergeomics: a web server for identifying pathological pathways, networks, and key regulators via multidimensional data integration. BMC Genomics. (2016) 17:722. doi: 10.1186/s12864-016-3057-8
13. Shu L, Zhao Y, Kurt Z, Byars SG, Tukiainen T, Kettunen J, et al. Mergeomics: multidimensional data integration to identify pathogenic perturbations to biological systems. BMC Genomics. (2016) 17:874. doi: 10.1186/s12864-016-3198-9
14. Ehret GB, Munroe PB, Rice KM, Bochud M, Johnson AD, Chasman DI, et al. Genetic variants in novel pathways influence blood pressure and cardiovascular disease risk. Nature. (2011) 478:103. doi: 10.1038/nature10405
15. Nikpay M, Goel A, Won HH, Hall LM, Willenborg C, Kanoni S, et al. A comprehensive 1000 Genomes–based genome-wide association meta-analysis of coronary artery disease. Nat Genet. (2015) 47:1121. doi: 10.1038/ng.3396
18. Greenawalt DM, Dobrin R, Chudin E, Hatoum IJ, Suver C, Beaulaurier J, et al. A survey of the genetics of stomach, liver, and adipose gene expression from a morbidly obese cohort. Genome Res. (2011) 21:1008–16. doi: 10.1101/gr.112821.110
19. Romanoski CE, Che N, Yin F, Mai N, Pouldar D, Civelek M, et al. Network for activation of human endothelial cells by oxidized phospholipids: a critical role of heme oxygenase 1. Circul Res. (2011) 109:e27–41. doi: 10.1161/CIRCRESAHA.111.241869
21. Fehrmann RS, Jansen RC, Veldink JH, Westra HJ, Arends D, Bonder MJ, et al. Trans-eqtls reveal that independent genetic variants associated with a complex phenotype converge on intermediate genes, with a major role for the hla. PLoS Genet. (2011) 7:e1002197. doi: 10.1371/journal.pgen.1002197
22. Nica AC, Parts L, Glass D, Nisbet J, Barrett A, Sekowska M, et al. The architecture of gene regulatory variation across multiple human tissues: the muther study. PLoS Genet. (2011) 7:e1002003. doi: 10.1371/journal.pgen.1002003
23. Montgomery SB, Sammeth M, Gutierrez-Arcelus M, Lach RP, Ingle C, Nisbett J, et al. Transcriptome genetics using second generation sequencing in a caucasian population. Nature. (2010) 464:773–U151. doi: 10.1038/nature08903
24. Stranger BE, Montgomery SB, Dimas AS, Parts L, Stegle O, Ingle CE, et al. Patterns of cis regulatory variation in diverse human populations. PLoS Genet. (2012) 8:272–84. doi: 10.1371/journal.pgen.1002639
26. Dimas AS, Deutsch S, Stranger BE, Montgomery SB, Borel C, Attar-Cohen H, et al. Common regulatory variation impacts gene expression in a cell type-dependent manner. Science. (2009) 325:1246–50. doi: 10.1126/science.1174148
28. Huan TX, Liu CY, Joehanes R, Zhang XL, Chen BH, Johnson AD, et al. A systematic heritability analysis of the human whole blood transcriptome. Hum Genet. (2015) 134:343–58. doi: 10.1007/s00439-014-1524-3
29. Ardlie KG, DeLuca DS, Segrè AV, Sullivan TJ, Young TR, Gelfand ET, et al. The genotype-tissue expression (gtex) pilot analysis: Multitissue gene regulation in humans. Science. (2015) 348:648–60. doi: 10.1126/science.1262110
34. Hindorff LA, Sethupathy P, Junkins HA, Ramos EM, Mehta JP, Collins FS, et al. Potential etiologic and functional implications of genome-wide association loci for human diseases and traits. Proc Natl Acad Sci USA. (2009) 106:9362–7. doi: 10.1073/pnas.0903103106
35. Derry JMJ, Zhong H, Molony C, MacNeil D, Guhathakurta D, Zhang B, et al. Identification of genes and networks driving cardiovascular and metabolic phenotypes in a mouse f2 intercross. PLoS ONE. (2010) 5:e14319. doi: 10.1371/journal.pone.0014319
36. Wang SS, Schadt EE, Wang H, Wang XP, Ingram-Drake L, Shi W, et al. Identification of pathways for atherosclerosis in mice - integration of quantitative trait locus analysis and global gene expression data. Circul Res. (2007) 101:E11–30. doi: 10.1161/CIRCRESAHA.107.152975
37. Yang X, Schadt EE, Wang S, Wang H, Arnold AP, Ingram-Drake L, et al. Tissue-specific expression and regulation of sexually dimorphic genes in mice. Genome Res. (2006) 16:995–1004. doi: 10.1101/gr.5217506
38. Tu ZD, Keller MP, Zhang CS, Rabaglia ME, Greenawalt DM, Yang X, et al. Integrative analysis of a cross-loci regulation network identifies app as a gene regulating insulin secretion from pancreatic islets. PLoS Genet. (2012) 8:e1003107. doi: 10.1371/journal.pgen.1003107
40. Wang IM, Zhang B, Yang X, Zhu J, Stepaniants S, Zhang C, et al. Systems analysis of eleven rodent disease models reveals an inflammatome signature and key drivers. Mol Syst Biol. (2012) 8:594. doi: 10.1038/msb.2012.24
44. Menni C, Mangino M, Zhang F, Clement G, Snieder H, Padmanabhan S, et al. Heritability analyses show visit-to-visit blood pressure variability reflects different pathological phenotypes in younger and older adults: evidence from uk twins. J Hypert. (2013) 31:2356–61. doi: 10.1097/HJH.0b013e32836523c1
45. Crowley SD, Gurley SB, Herrera MJ, Ruiz P, Griffiths R, Kumar AP, et al. Angiotensin II causes hypertension and cardiac hypertrophy through its receptors in the kidney. Proc Nat Acad Sci USA. (2006) 103:17985–90. doi: 10.1073/pnas.0605545103
47. Uehata M, Ishizaki T, Satoh H, Ono T, Kawahara T, Morishita T, et al. Calcium sensitization of smooth muscle mediated by a Rho-associated protein kinase in hypertension. Nature. (1997) 389:990. doi: 10.1038/40187
48. Brænne I, Civelek M, Vilne B, Di Narzo A, Johnson AD, Zhao Y, et al. Prediction of causal candidate genes in coronary artery disease loci. Arterioscler Thromb Vasc Biol. (2015) 35:2207–17. doi: 10.1161/ATVBAHA.115.306108
51. Tu W, Eckert GJ, DiMeglio LA, Yu Z, Jung J, Pratt JH. Intensified effect of adiposity on blood pressure in overweight and obese children. Hypertension. (2011) 58:818–24. doi: 10.1161/HYPERTENSIONAHA.111.175695
52. Zhang T, Zhang H, Li S, Li Y, Liu Y, Fernandez C, et al. Impact of adiposity on incident hypertension is modified by insulin resistance in adults: longitudinal observation from the Bogalusa Heart Study. Hypertension. (2016) 67:56–62. doi: 10.1161/HYPERTENSIONAHA.115.06509
54. Eisenberg T, Abdellatif M, Schroeder S, Primessnig U, Stekovic S, Pendl T, et al. Cardioprotection and lifespan extension by the natural polyamine spermidine. Nat Med. (2016) 22:1428. doi: 10.1038/nm.4222
56. Cheng PW, Chen YY, Cheng WH, Lu PJ, Chen HH, Chen BR, et al. Wnt signaling regulates blood pressure by downregulating a GSK-3β-mediated pathway to enhance insulin signaling in the central nervous system. Diabetes. (2015) 2015:db141439. doi: 10.2337/db14-1439
58. Calabrese GM, Mesner LD, Stains JP, Tommasini SM, Horowitz MC, Rosen CJ, et al. Integrating GWAS and co-expression network data identifies bone mineral density genes SPTBN1 and MARK3 and an osteoblast functional module. Cell Syst. (2017) 4:46–59. doi: 10.1016/j.cels.2016.10.014
60. Yoda A, Adelmant G, Tamburini J, Chapuy B, Shindoh N, Yoda Y, et al. Mutations in G protein β subunits promote transformation and kinase inhibitor resistance. Nat Med. (2015) 21:71. doi: 10.1038/nm.3751
61. Warren HR, Evangelou E, Cabrera CP, Gao H, Ren M, Mifsud B, et al. Genome-wide association analysis identifies novel blood pressure loci and offers biological insights into cardiovascular risk. Nat Genet. (2017) 49:403. doi: 10.1038/ng.3768
62. Jakobsson L, Domogatskaya A, Tryggvason K, Edgar D, Claesson-Welsh L. Laminin deposition is dispensable for vasculogenesis but regulates blood vessel diameter independent of flow. FASEB J. (2008) 22:1530–9. doi: 10.1096/fj.07-9617com
63. Barallobre-Barreiro J, Oklu R, Lynch M, Fava M, Baig F, Yin X, et al. Extracellular matrix remodelling in response to venous hypertension: proteomics of human varicose veins. Cardiovasc Res. (2016) 110:419–30. doi: 10.1093/cvr/cvw075
64. Inamoto S, Kwartler CS, Lafont AL, Liang YY, Fadulu VT, Duraisamy S, et al. TGFBR2 mutations alter smooth muscle cell phenotype and predispose to thoracic aortic aneurysms and dissections. Cardiovasc Res. (2010) 88:520–9. doi: 10.1093/cvr/cvq230
65. Sagnol S, Yang Y, Bessin Y, Allemand F, Hapkova I, Notarnicola C, et al. Homodimerization of RBPMS2 through a new RRM-interaction motif is necessary to control smooth muscle plasticity. Nucleic Acids Res. (2014) 42:10173–84. doi: 10.1093/nar/gku692
66. Sabbatini AR, Barbaro NR, de Faria AP, Modolo R, Ritter AM, Pinho C, et al. Increased circulating tissue inhibitor of metalloproteinase-2 is associated with resistant hypertension. J Clin Hypertension. (2016) 18:969–75. doi: 10.1111/jch.12865
67. Wei L, Warburton RR, Preston IR, Roberts KE, Comhair SA, Erzurum SC, et al. Serotonylated fibronectin is elevated in pulmonary hypertension. Am J Physiol Lung Cell Mol Physiol. (2012) 20:L1273–9. doi: 10.1152/ajplung.00082.2012
68. López B, Querejeta R, González A, Larman M, Díez J. Collagen cross-linking but not collagen amount associates with elevated filling pressures in hypertensive patients with stage C heart failure: potential role of lysyl oxidase. Hypertension. (2012) 60:677–83. doi: 10.1161/HYPERTENSIONAHA.112.196113
69. Wickramasekera NT, Gebremedhin D, Carver KA, Vakeel P, Ramchandran R, Schuett A, et al. Role of dual-specificity protein phosphatase-5 in modulating the myogenic response in rat cerebral arteries. J Appl Physiol. (2013) 114:252–61. doi: 10.1152/japplphysiol.01026.2011
70. Pramanik K, Chun CZ, Garnaas MK, Samant GV, Li K, Horswill MA, et al. Dusp-5 and Snrk-1 coordinately function during vascular development and disease. Blood. (2009) 113:1184–91. doi: 10.1182/blood-2008-06-162180
71. Gusev NB. Some properties of caldesmon and calponin and the participation of these proteins in regulation of smooth muscle contraction and cytoskeleton formation. Biochemistry. (2001) 66:1112–21. doi: 10.1023/A:1012480829618
Keywords: blood pressure, genome wide association studies, integrative genomics, regulatory networks, key drivers
Citation: Zhao Y, Blencowe M, Shi X, Shu L, Levian C, Ahn IS, Kim SK, Huan T, Levy D and Yang X (2019) Integrative Genomics Analysis Unravels Tissue-Specific Pathways, Networks, and Key Regulators of Blood Pressure Regulation. Front. Cardiovasc. Med. 6:21. doi: 10.3389/fcvm.2019.00021
Received: 07 November 2018; Accepted: 18 February 2019;
Published: 12 March 2019.
Edited by:Tanja Zeller, Universität Hamburg, Germany
Reviewed by:Hauke Busch, Universität zu Lübeck, Germany
Joylene Elisabeth Siland, University of Groningen, Netherlands
Copyright © 2019 Zhao, Blencowe, Shi, Shu, Levian, Ahn, Kim, Huan, Levy and Yang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Xia Yang, email@example.com