Whole Genome Sequence Analysis of the Plasma Proteome in Black Adults Provides Novel Insights Into Cardiovascular Disease

Editorial, see p 371

The circulating plasma proteome plays a fundamental role in human biological function and dysfunction. Circulating proteins both mediate and respond to disease, and are frequently the targets of pharmaceutical interventions. Several recent studies have coupled genotyping and proteomic profiling to understand the genetic basis for the individual differences observed in protein levels, which are known to be heritable.^1–7 Such work has led to critical advances in our understanding of the genetic architecture of the plasma proteome and its relationship to disease, including factors specifically associated with cardiovascular risk.^4,6,7 However, initial findings were derived nearly entirely in European populations such as the Framingham Heart Study using genotyping arrays. Further, individuals with increased African ancestry are known to harbor substantially more genetic diversity than those of European ancestry,^8,9 and rare mutations found specifically among persons of African ancestry have been critical in expanding our knowledge of cardiovascular biology, as is the case for PCSK9.¹⁰ We hypothesized that coupling whole genome sequence analysis with plasma proteomics in individuals of African ancestry would greatly increase the power to identify novel genetic determinants of the plasma proteome, which would inform our understanding not only of ancestry-specific genetic variation but also of human cardiovascular biology in general.

Here we use whole genome sequence data and aptamer-based proteomic profiling of 1301 proteins on the SOMAscan platform in 1852 self-identified Black individuals from the JHS (Jackson Heart Study)¹¹ to identify novel protein quantitative trait loci (pQTLs) determining protein levels. Associations were replicated in 980 participants from MESA (Multi-Ethnic Study of Atherosclerosis)¹² and 708 from the HERITAGE Family Study (Health, Risk Factors, Exercise Training and Genetics; Table S1),¹³ and further validated using an alternate proteomic profiling platform in JHS. These data serve as the basis for an enhanced understanding of proteins highly relevant to cardiovascular homeostasis across diverse human populations.

Methods

Data Availability

Whole genomes for JHS and MESA, generated as part of the National Heart, Lung, and Blood Institute (NHLBI) TOPMed (Trans-Omics for Precision Medicine) program, are available through restricted access via the NHLBI database of Genotypes and Phenotypes (dbGaP). TOPMed accession numbers for JHS and MESA are phs000964/phs002256.v1.p1 and phs001416, respectively. Full genome-wide association study (GWAS) summary statistics for JHS (the discovery cohort) generated in this study will be available for general research use through controlled access at dbGaP accession phs001974: NHLBI TOPMed: Genomic Summary Results for the Trans-Omics for Precision Medicine program. For assistance in accessing the discovery data in JHS before full availability on dbGaP, investigators should contact the authors and follow JHS data access procedures (https://www.jacksonheartstudy.org/). GWAS data for the replication studies (MESA and HERITAGE) are fully included in the article. Individual-level proteomic and genomic data in the replication datasets are available through application to the respective cohorts.

Study Approval

The JHS study was approved by Jackson State University, Tougaloo College, and the University of Mississippi Medical Center Institutional Review Boards, and all participants provided written informed consent. All MESA participants provided written informed consent, and the study was approved by the Institutional Review Boards at The Lundquist Institute (formerly Los Angeles BioMedical Research Institute) at Harbor-University of California, Los Angeles, Medical Center, University of Washington, Wake Forest School of Medicine, Northwestern University, University of Minnesota, Columbia University, Johns Hopkins University, and University of California, Los Angeles. The human study protocols were approved by the Institutional Review Boards of Beth Israel Deaconess Medical Center, University of Washington, and the 4 clinical centers of HERITAGE.

Cohorts

The JHS, MESA, and the HERITAGE Family Study have all been previously described.^11–13 In brief, JHS is a community-based longitudinal cohort study begun in 2000 of 5306 self-identified Black individuals from the Jackson, Mississippi, metropolitan statistical area.¹¹ Included in the present study are samples collected at Visit 1 between 2000 and 2004 from 1852 individuals with whole genome sequencing (WGS)¹⁴ and proteomic profiling performed in batches (see Proteomic Profiling below).

MESA began in 2000 with 6814 men and women age 45 to 84 years recruited at 6 clinical centers across the United States. Participants were identified belonging to 4 racial/ethnic groups: Black, Hispanic, Asian, or White. Included in the present study are 980 individuals selected randomly across all 4 racial/ethnic groups with proteomic profiling from Visit 1 between 2000 and 2002 and whole genome sequence analysis.¹²

HERITAGE enrolled a combination of self-identified White and Black family units, totaling 763 sedentary participants (62% White) between the ages of 17 and 65 years in a 20-week, graded endurance exercise training study across 4 clinical centers in the United States and Canada in 1994 to 1995.¹³ Included in the present study are a random subset of 708 individuals with baseline plasma samples and genotyping.

Proteomic Profiling

Proteomic profiling by SomaScan (aptamer-based affinity platform) and Olink (antibody-based affinity platform) has been described previously.^6,15 Please see the Supplemental Methods for further details.

Genotyping and Imputation

WGS in JHS and MESA has been described previously.^14,16 Included in the present study are participants included in Freeze 6 of the TOPMed project at the Northwest Genome Center at University of Washington and the Broad Institute. Samples underwent >30× WGS. Genotype calling with vt¹⁷ and quality control were performed by the Informatics Resource Center at the University of Michigan.

Genotyping in HERITAGE was performed on the Illumina Infinium Global Screening Array. Genotypes were called using Illumina’s GenCall on the basis of the top/bottom strand method. Genotype imputation was performed using the University of Michigan Imputation Server Minimac4 to reference panel TOPMed Freeze5.¹⁸ Phasing was performed with Eagle v2.4. Sites were excluded with call rate <90%, mismatched alleles, or invalid alleles (88% of sites retained).

Statistical Analysis

All statistical methods are explained throughout the sections below.

Whole Genome Sequence Association Analysis

Across all 3 cohorts, proteomic measurements were standardized to a set of control samples (pooled plasma) that were part of each plate. The resulting values were log-transformed and scaled to a mean of 0 and SD of 1. In JHS, to account for batch effects, proteins were log-transformed and scaled within batch and then combined. In all cohorts, these log-transformed values were residualized on age, sex, batch, and principal components of ancestry 1 to 10 as determined by Genetic Estimation and Inference in Structured Samples (GENESIS).^16,19,20 In HERITAGE and MESA, measurements were also residualized on race to account for nongenetic racial effects not captured by genetic ancestry. The resulting residuals were then inverse-normalized. The association between these values and genetic variants was tested using linear mixed effects models adjusted for age, sex, the genetic relationship matrix, and principal components 1 to 10 using the fastGWA model implemented in the GCTA software package (version 1.93.2beta/gcta64).²¹ Repeat adjustment was implemented to reduce type I error and improve statistical power.²² Variants with a minor allele count less than 5 in a given cohort were excluded from analysis in that cohort. A Bonferroni-adjusted significance threshold of 3.8×10^-11 (5×10^-8/1301) was used for discovery in JHS. For variants in cis (<1 Mb from the transcription start site [TSS] of the coding gene for the associated protein),¹ variants with P values of 5×10^-6 were also considered in a separate analysis, given the biological plausibility of such associations.

Variance Explained for Each Protein

Single nucleotide polymorphism (SNP)–based heritability, hSNP2, was estimated using a linkage disequilibrium (LD)– and minor allele frequency (MAF)–stratified genomic relatedness matrix restricted maximum likelihood (GREML-LDMS) model implemented in the GCTA software. This method allows for fitting multiple genomic relatedness matrices with SNPs binned according to their regional LD and MAF.²³ It is recommended for heritability estimation on WGS data.^23,24 Using this model, we first calculated the segment-based (length of segment: 200 kb) LD scores and partitioned SNPs into 4 groups on the basis of the quartiles of the regional LD score. Genomic relatedness matrices for each of the 4 groups was then computed using SNPs binned into the corresponding group, and jointly fitted into a mixed effect model for estimating the heritability and variance. In our analysis, we allowed for a maximum of 1000 iterations. For all analyses, we adjusted for age, sex, and the first 10 principal components of genetic ancestry. Variance explained by the top performing variant (as determined by lowest P value) was estimated using the equation BETA^2×(2×AF1× (1-AF1)/VAR) where BETA was the beta estimate for the effect allele, AF1 was the allele frequency of the effect allele, and VAR was the variance of the protein residual used for WGS analysis. Variance explained by clinical covariates was estimated using linear regression of log-transformed protein level regressed on age, sex, systolic blood pressure, diabetes, use of hypertensive medication, current smoking status, and a history of coronary heart disease. Proteins whose total heritability could not be estimated by this model, which often occurs when heritability estimates are low,²³ were excluded (N=185, Figure S1, Table S2).

Defining Protein-Locus Associations and Sentinel Variants

To identify the broadest genomic region associated with a protein, we applied the following previously described algorithm:¹ a 1-Mb region around each SNP associated with a given protein was defined. Beginning with the region containing the variant with the lowest P value, overlapping regions were merged together. This was repeated until no more overlapping regions existed for the given protein. The variant with the lowest P value in each region was identified as the sentinel variant. To describe regions associated with multiple proteins, regions with sentinel variants in LD with r²≥0.8 were described as the same region, exclusively for descriptive purposes. LD was determined using SNPClip, using data from individuals of African ancestry.^25,26 Any variants within 1 Mb of the TSS for the cognate gene of a protein were considered “cis.”

Replication in MESA/HERITAGE

Associations between sentinel variants and proteins from JHS were evaluated in MESA and HERITAGE separately if associated statistics were available (if minor allele count was <5 in either cohort, that variant was not considered in that cohort). Where association statistics were available in both cohorts, the 2 cohorts were meta-analyzed using the inverse-variance weighted method using fixed effects. Validation threshold was set at P<0.05 with consistent direction of effect.

Meta-Analysis

Results from JHS, MESA, and HERITAGE were meta-analyzed together using mixed effects models in the “meta” package of R4.0.5. Only variants with a P value for association with a given protein <1×10^-5 in at least 2 of the studies were included.

Comparing With Previous pQTLs

To determine whether pQTLs were novel, we used the PhenoScanner package (version 2) for R.^27,28 For each protein-locus association identified above, we divided the locus into 1-Mb or less segments (maximum permitted by PhenoScanner API) if needed. The resulting region or regions were then passed to the phenoscanner function in R, with the following arguments: build was set to “38,” P value to 1×10^-5, catalog to “pQTL,” and proxies to “None” (query date June 28, 2020). To supplement PhenoScanner, we reviewed the literature for additional studies using SomaScan or Olink to identify the genetic architecture of the plasma proteome and identified 3 not in PhenoScanner.^2,6,7 Results from these studies were considered using the same criteria as above. If the protein linked to that region in our analysis was found to be previously associated with any variants in the region, this was considered a “previous” protein-locus association. For the subset of protein-locus associations that were previously described, we secondarily looked to see whether the sentinel SNP in JHS represented a novel genetic determinant. Sentinel SNPs were queried against both PhenoScanner and the 3 other studies to look for any variants associated with the same protein and in linkage disequilibrium with the new sentinel SNP. Again the phenoscanner function in R was used with the following arguments: build was set to “38,” P value to 1×10^-5, catalog to “pQTL,” proxies to “EUR” (because these variants were discovered in European populations), and r² to “0.5” (query date October 1, 2020).

Comparing With Previous GWAS Results

To determine overlap between clinical GWAS analyses and pQTLs in this analysis, we used the PhenoScanner package for R. All 569 sentinel SNPs as identified above were passed to the phenoscanner function in R with the following arguments: build was set to “38,” P value to “1×10^-5,” catalog to “GWAS,” r² to “0.5,” and proxies to “AFR” (query date October 15, 2020).

Comparing Results With ClinVar Data

The entirety of the ClinVar database was downloaded from the National Center for Biotechnology Information FTP site (https://ftp.ncbi.nlm.nih.gov/pub/clinvar/tab_delimited/variant_summary.txt.gz, accessed September 3, 2020).²⁹ These data were merged to all variants associated with any protein in the JHS at a P value <5×10^-6.

Variant Annotations

Reference allele frequencies from gnomAD³⁰ and variant category from GENCODE³¹ were obtained from the Functional Annotation of Variants–Online Resource (available at https://favor.genohub.org, download date July 20, 2020).³²

Results

Whole Genome Association Analysis of Proteomic Profiling

We performed whole genome association analysis between 28.1 million variants with an allele count in JHS of at least 5 and 1301 plasma protein measures in 1852 self-identified Black individuals (61% women). Proteins exhibited a wide range of estimated total heritability (median heritability = 0.33, interquartile range 0.22–0.48, Figure S1, Table S2). Imputing proteins with nonconverged heritability estimates to 0 resulted in a median heritability of 0.29 (see Methods).

At a Bonferroni-adjusted significance cutoff (5×10^-8/1301=3.8×10^-11), we identified 569 associations with 479 proteins encompassing 438 unique genetic loci (Figure 1, Table S3). Each locus is a genomic region containing at least 1 variant associated with a protein but often summarizing multiple nearby variants in varying degrees of linkage disequilibrium (see Methods). The variant with the lowest P value for association in each locus is considered the sentinel variant. Using this method, we identified 114 locus-protein associations not previously described. For previously described locus-protein associations, we identified novel sentinel variants in 217 loci.

• Download figure
• Download PowerPoint

Across these 569 associations, 329 (58%) of the sentinel variants for the given locus are within 1 Mb of the TSS of the cognate gene for that protein (termed cis), and 240 are nonlocal (termed trans). We identified an additional 183 suggestive cis associations when the P value threshold was lowered to 5×10^-6 (Table S4). The majority of cis pQTLs were close to the TSS of the cognate gene, with 90% falling within 100 kb of the TSS (Figure 2A).

• Download figure
• Download PowerPoint

The majority of proteins (70%) with a significant pQTL were associated with a single locus. Three proteins were associated with 5 different loci: Ck-beta-8-1, cyclin-dependent kinase inhibitor 1B, and apolipoprotein L1 (Figure 2B).

Patterns observed in previous studies were replicated here: most loci (388, 89%) were associated with only 1 protein, although there were several pleiotropic loci including regions near the VTN, ABO, and APOE genes (Figure 2C), all of which have been implicated in cardiovascular disease.^1,3,33–36 Sentinel variants were largely proximate to coding genes, with only 20% in intergenic regions (Figure 2D and 2E). There was a strong inverse relationship between effect size and MAF, consistent with previous pQTL studies (Figure 2F).¹

In contrast with previous studies, a significant number of sentinel variants had allele frequencies that varied substantially from those observed in European populations: 166 (36%) of the 464 identified sentinel variants in JHS had MAF <1%, whereas 65 (14%) of the variants had MAF <0.0001% among non-Finnish Europeans in gnomAD.³⁰ Many of these variants were much more common in JHS: among the 166 variants with MAF <1% in non-Finnish Europeans, 71 had MAF >5% in JHS. Figure 3 illustrates the wide disparity between allele frequencies of all 569 sentinel variants in African versus non-Finnish European populations in gnomAD.³⁰

• Download figure
• Download PowerPoint

We also completed proteomic profiling in 2 smaller cohorts, MESA (N=980, 53% women, 19% Black) and the HERITAGE Family Study (N=708, 56% women, 36% Black), each containing a subset of self-identified Black individuals, which were meta-analyzed (when possible) to validate the results. Consistent associations were observed for 90% of the 569 sentinel variants at a P value <0.05 with matching direction of effect. If a significance threshold adjusted for multiple corrections is used (P<0.05/569=8.8×10^-5), 72% replicate. Variants that did not replicate in some cases had lower MAF, falling below the minor allele count threshold of 5 in 1 of the 2 replication cohorts, reducing overall replication power (Table S3, Figure S2). Results from JHS, MESA, and HERITAGE were also meta-analyzed together. This analysis yielded 13 additional pQTLs: 9 trans and 4 cis (Table S5).

In a limited subsample of JHS participants (N=488), plasma samples were also profiled using the Olink Explore platform, which uses a completely distinct, immunoassay-based approach for protein measurement, which generally relies on polyclonal antibody conjugates.³⁷ Of the 569 sentinel variant-protein associations, 318 could be compared on the Olink platform. These associations showed a consistent effect across the 2 platforms (correlation of effect, 0.82 [95% CI, 0.78–0.85]; Figure S3). Across all 318 comparisons, the median Soma-Olink correlation was 0.62 (interquartile range, 0.35–0.74). The direction of effect matched in 86%, and 51% of associations were confirmed at a Bonferroni (0.05/318) level of significance (Table S3). There were a small number of discordant associations where effects as measured by SOMA and Olink were significant but with opposing directions of effect, such as the association between rs5744204 and lipopolysaccharide-binding protein. These may indicate platform-specific binding effects, but still support a genetic effect on protein levels as the most likely explanation, save for the unlikely possibility of opposing effects on just the binding of reagents from each platform.

Although all pQTLs are listed in Table S3, a subset of the results and information discussed in the following sections is highlighted in the Table.

Novel Genetic Determinants of Plasma Proteins Related to Thrombosis, Lipid Biology, and Myocardial Disease

To determine the novelty of the wide genomic regions identified as pQTLs by our analysis, we queried pQTL data available in PhenoScanner, a database of GWAS findings.^27,28 Of the 569 protein-locus associations, 114 (20%) had not been previously identified (Figure 1, Table S3) at a P value <1×10^-5. Of these 114 novel associations, 84 (74%) were trans associations. Sixty-two (54%) of the sentinel variants for these loci were uncommon (ie, MAF <1%) in non-Finnish European populations, but had a median MAF in JHS of 5% (interquartile range, 2%–12%). Novel pQTLs provide the opportunity to better understand biological pathways. As an example, a variant in the 5-prime untranslated region of F12, the gene for clotting factor XII, is observed to be a novel pQTL for thrombin and plasma serine protease inhibitor. This variant has previously been shown to affect thrombin generation and the coagulation cascade.³⁸

Similar to previous studies, we identify multiple pleiotropic genetic loci that affect the levels of multiple proteins. The APOE locus is 1 such well-established locus, which is known to be associated with hypercholesterolemia, atherosclerotic heart disease, and Alzheimer disease.³⁹ Our analysis reveals 6 new proteins associated with this gene at 3 distinct (r²<0.1) missense variants: rs7412, rs769455, and rs42935 (Figure 4A). These 6 proteins: b-endorphin, MMP-3 (matrix metalloproteinase-3), Sonic hedgehog, ZAP70 (zeta chain of T cell receptor associated protein kinase 70), Kelch-like ECH-associated protein 1, and matrix metalloproteinase-8, implicate new targets in understanding how APOE (apolipoprotein E) may mediate its effects. Indeed, APOE knockout mice, which develop atherosclerotic lesions that mimic human plaques, have shown reduced Zap70 activation.⁴⁰ Further, MMP-3 levels have been shown to be elevated in affected areas of the brain among those with Alzheimer disease.⁴¹

• Download figure
• Download PowerPoint

The analysis also shows a new pleiotropic locus at HPX, the gene for hemopexin. The sentinel variant, rs12117, is nearly monoallelic in European populations, but has a MAF in JHS of 5%. Six proteins are shown to be affected by this variant: bone morphogenetic protein receptor type-2, natural killer cell receptor 2B4, K-Ras, glial cell line–derived neurotrophic factor, TNFα (tumor necrosis factor-α), and tumor necrosis factor ligand superfamily member 18. Another 3 proteins are associated with other variants either in or upstream of HPX (Figure 4A). It has been posited that hemopexin protects cells from oxidative stress by clearing heme, and TNFα is known to induce HPX expression in rats as an acute phase response.⁴² Further, HPX/APOE double knockout mice had accelerated atherosclerosis related to oxidative stress and changes in macrophage function. This role of hemopexin may be particularly important in Black patients with sickle-cell disease: murine models have shown the value of heme-scavenging by hemopexin in reducing inflammation in models of sickle-cell disease.^43,44 Our findings suggest specific genetic variation may have a role in the immune functions of hemopexin. Although no members of our cohort had sickle-cell disease, 24 individuals did have both the minor allele of rs12117 and sickle-cell trait. However, no definitive interaction between these 2 variants and any protein could be identified. Unfortunately, given the low frequency of the variant in European-based GWAS, no clinical implications for rs12117 have been identified, although other variants in HPX have been linked to ulcerative colitis.⁴⁵ Further data are needed; specifically, data from patients with sickle-cell disease would be of value.

Our analysis can implicate new biology related to previously described variants as well. The variant rs2066702 in ADH1B has been identified as a risk locus for alcohol dependence across multiple ancestry-specific GWAS.⁴⁶ The same variant in our analysis is associated with levels of NAMPT (nicotinamide phosphoribosyltransferase; Figure S4A), which regulates intracellular nicotinamide adenine dinucleotide (NAD⁺) and plays a role in cardiac hypertrophy and adverse remodeling.⁴⁷ It is important that the minor allele of rs2066702 is protective of alcohol dependence, and it is this allele that is associated with higher levels of NAMPT, suggesting that alcohol use may deplete NAMPT in humans. Furthermore, previous murine studies have shown that ethanol administration diminished NAMPT levels, whereas overexpression of NAMPT was found to protect against steatosis.⁴⁸

Conversely, the associations between well-described proteins and poorly understood genes can further elucidate biology. Levels of 2 proteins, plasminogen and angiostatin (itself a fragment of plasminogen), were linked to a variant upstream of GALNT7 (Figure S4B). Plasminogen and angiostatin each have a strong cis pQTL, supporting aptamer specificity for their measurement (Tables S3 and S4). Although plasminogen and angiostatin are critical factors in clot dissolution and angiogenesis inhibition,^49,50 respectively, the biological role of GALNT7, a glycosyltransferase, has been linked by more limited evidence to cancer proliferation.⁵¹ The sentinel SNPs linked to these proteins in our analysis are monoallelic in European populations, so previous GWAS data do not exist. However, other variants at the GALNT7 locus have been linked to vascular disorders in the UK Biobank, including “Cause of death: peripheral vascular disease, unspecified” (P=1×10^-23), “Cause of death: vascular dementia, unspecified” (P=8×10^-20), and “Cause of death: chronic or unspecified with hemorrhage” (P=2×10^-17), all 3 of which are plausibly mediated by plasminogen or angiostatin.^27,28

Known Ancestry-Specific Loci Highlight Ancestry-Specific Cardiovascular Disease Pathways

Analysis of samples from individuals of greater African ancestry allows for assessment of specific loci known to be of particular clinical importance in individuals of African descent. We evaluated the proteomic signatures of 4 such well-described loci.

TTR (transthyretin) amyloidosis results from the misfolding of the transthyretin tetramer, ultimately resulting in abnormal protein deposition in myocardium and nerve tissue, leading to cardiomyopathy and neuropathy. Protein misfolding is accelerated in the presence of mutations in the TTR gene; specifically, rs76992529 encodes a V122I mutation that is found in 3% to 4% of Black individuals. In our data we show this variant to be a robust pQTL for RBP4 (retinol-binding protein 4), a binding partner of TTR.⁵² In individuals with TTR amyloidosis and overt myocardial disease (typically manifested as left ventricular thickening and diastolic dysfunction), RBP4 levels are known to be diminished—the normal transthyretin tetramer protects RBP4 from renal clearance.⁵³ However, our data show that asymptomatic carriers of this mutation have diminished RBP4 levels as well, even in the absence of reported heart failure (Figure 5A). To further explore this finding, we leveraged extensive metabolite profiling in JHS.⁵⁴ We found an unknown metabolite feature highly correlated with circulating RBP4 (Pearson correlation 0.64, [95% CI, 0.61–0.66]). As expected, the association between this metabolite and the V122I mutation was also strong (β=–0.76, P=4.6×10^-14). This metabolite feature has a mass-to-charge ratio of 269.226, which strongly suggests its identity as a dehydrated form of retinol, the binding partner of RBP4, according to the Human Metabolome Database. These data further complement and validate our proteomic association of RBP4 and TTR. Larger datasets are needed to explore the functional consequences of these proteomic and metabolomic findings.

• Download figure
• Download PowerPoint

Two alleles in the APOL1 gene (rs73885319/rs60910145 or “G1” and rs71785313 or “G2”) are linked to chronic kidney disease and cardiovascular disease in JHS and are common in individuals with African ancestry.^55–57 In JHS, rs73885319 has a MAF of 23%, whereas the variant is not present persons of European ancestry in gnomAD. In addition to being associated with levels of APOL1 (apolipoprotein L1) in our analysis, it was also the sentinel SNP determining levels of CKAP2 (cytoskeleton associated protein 2, Figure 5B). CKAP2 has been linked to tumor formation because it has a role in mitosis, but has also been observed to be upregulated in renal tubular necrosis.^58,59 In models adjusted for age, sex, body mass index, systolic blood pressure, presence of hypertension, presence of diabetes, HbA1c, and proteomic batch/plate, CKAP2 levels as measured by SOMAscan were associated with increased estimated glomerular filtration rate in JHS (β=1.16, P=0.002). Because APOL1 risk variants are associated with renal disease, this could point to a protective role for CKAP2 in response to APOL1 genetic risk, requiring further investigation as a therapeutic target.

DARC (Duffy chemokine receptor) is a binding site crucial to malarial infection with Plasmodiumvivax, but has also been shown to affect risk for cardiovascular outcomes in JHS.⁶⁰ Under positive selection in sub-Saharan Africa, the FY*O allele of this gene is thus common in individuals of African descent, although it is present in only 0.4% of individuals of non-Finnish European descent in gnomAD.³⁰ Levels of CCL14 (C-C motif chemokine ligand 14) and eotaxin have previously been linked to this gene, and to this list we now add protein S100-A9, CXCL11 (C-X-C motif chemokine ligand 11), and bactericidal permeability-increasing protein. Despite being linked to neutropenia, the Duffy-null allele has not been shown to lead to an increased risk of infection.⁶¹ However, there is evidence of a slower progression of HIV infection in the Duffy-null state.⁶² These results expand the list of inflammatory mediators affected by the Duffy-null state.

Last, the variant that causes sickle cell trait, rs334, has an allele frequency of 4% in JHS. This variant was associated with fractalkine (P=2.5×10^-6). Previous work has linked fractalkine, an inflammatory cytokine, to incident heart failure, specifically in Black individuals.⁶³

Protein Associations for Clinically Relevant Variants

Among the other 435 protein-locus pairs with previously identified pQTLs in the same region, 44 of the previous pQTLs were at P values >5×10^-8, and 177 of the previous pQTLs differed from the sentinel variants identified in JHS (r²<0.5). Thus, even in genetic regions previously linked to a given protein, many sentinel variants identified in this analysis may point to novel genetic effects when combined with existing genetic databases (Tables S6 and S7). As an example, the variant rs2234355 in the CXCR6 gene is nearly monoallelic in European populations, but is common among African populations, and thus well represented in JHS (MAF 44%). The variant has been previously shown to be protective against Pneumocystis jiorvecii infection in HIV-infected individuals, and was more common in those achieving viremic control.^64,65 Interactions between CXCR6 (C-X-C chemokine receptor type 6) and its ligand CXCL16 (C-X-C motif chemokine ligand 16) have been posited as a potential mechanism; we show this variant to be a strong (P=5.7×10^-54) sentinel pQTL for CXCL16, supporting this hypothesis. The relationship may also have cardiovascular consequences, because CXCL16 levels have been associated with acute coronary syndromes.⁶⁶

Discussion

Our data represent a comprehensive effort to understand the genetic determinants of the circulating plasma proteome using whole genome sequence analysis in individuals with greater genetic diversity than those in previous analyses. We identify numerous novel genetic determinants of a wide range of circulating proteins, many of which are important in vascular and cardiac biology. Many of these genetic variants have known clinical implications, in which case our data delineate novel biology potentially linking genetic variation to disease. As an example, the genetic mutation associated with TTR amyloidosis in persons of African ancestry, rs76992529, is shown here to be associated with RBP4 levels in persons without overt cardiomyopathy. A recent study from the BioMe database found a similar difference among persons with this mutation and without cardiomyopathy.⁶⁷ Our findings extend the small case-control biobank study to a large, well-defined prospective cohort, advancing RBP4 levels as a potential preclinical biomarker. Further studies are needed to determine if there is an interaction between this mutation, RBP4 levels, and incident cardiomyopathy.

In other cases, the proteomic associations identified represent the first meaningful annotation of a given genetic variant. Such is the case for rs12117, a missense variant in the gene for hemopexin. Despite a MAF of ≈2.6% in persons of African ancestry, little is known about this variant. Here, we describe it as a pleiotropic locus, affecting the levels of multiple inflammatory proteins. Given hemopexin’s role in heme-scavenging, identifying additional carriers, particularly those with sickle-cell disease, may offer critical insights, and the proteins identified here would be useful starting points. The paucity of genome-wide association data in diverse populations limits our ability to interrogate associations, such as rs12117, with tools such as Mendelian randomization, but hopefully highlights the need for greater inclusion of diverse populations in genetic research going forward. Greater diversity in genetic association studies will not only increase our understanding of functional genomics but may also help delineate gene-environment interactions that affect individuals of diverse ancestry. Indeed, our analysis identifies novel variants that are not particularly rare in Europeans, but are only now described in a cohort of Black Americans. This finding suggests the possibility of gene-environment interactions, including the effects of social and structural differences, which have biological/health effects at multiple levels (health care access, stress response, environmental toxins, etc).⁶⁸ Such future work is important not only for the populations themselves, but also for optimum understanding of the genomic basis of biological variability and disease susceptibility.

Future work leveraging these data may also center around the intriguing finding of genetic variants that produce opposing findings on the Soma platform compared with the Olink platform. These variants, often protein-altering, likely affect binding of one platform, but the significant opposing effects suggest they are true pQTLs. Understanding the implications of such variants on a genome-wide scale may identify functionally important gene regions and inform interpretation of binding data.

Our study has several strengths: as mentioned, it is the largest analysis of its kind in a Black population, which gives it the power to detect many novel variants. The results are compared with 2 multiethnic populations and an alternate profiling platform. Our study also has several important limitations. Although this is the largest pQTL analysis in a Black population, the sample size for genome-wide association is relatively modest compared with many GWAS. This also informs a second limitation, the use of multi-ancestry cohorts for validation rather than a population of similar ancestry to JHS. This fact is related to limited availability of proteomic data in Black persons, and the desire to maintain an adequate sample size for validation of our original findings. For example, all 980 MESA participants with proteomics are included, regardless of their racial or ethnic identification, in the hope that statistical validation can be performed on as many variants as possible. Limiting MESA to only the Black participants would have left only 190 individuals. A further limitation is aptamer specificity on the SomaScan platform. Although cis pQTLs (both from this study and others) and validation on the Olink platform can confirm aptamer specificity, off target effects may be falsely attributed as transpQTLs, though we expect most cases off nonspecificity to bias toward the null. Aptamer validation efforts beyond those included here are ongoing across many groups.^1,2,69,70

Taken together, our work highlights the importance of extending proteomics, genomics, and likely other -omics studies, to diverse populations, to identify important potential biomarkers and disease pathways not only in those populations but also in the human population at large.

Article Information

Acknowledgments

Jackson Heart Study

The authors wish to thank the staff and participants of the JHS.

Heritage

The authors thank Drs Arthur S. Leon, D.C. Rao, James S. Skinner, Tuomo Rankinen, Jacques Gagnon, and the late Jack H. Wilmore for contributions to the planning, data collection, and conduct of the HERITAGE project.

Author Contributions

D.H.K., U.A.T., A.G.B., T.J.W., J.G.W., P.N., and R.E.G contributed to the original concept of the project, planned these analyses, and formulated the methods. D.H.K., U.A.T., D.N., M.D.B., J.M.R., D.E.C., SD, L.F., S.S., D.S, Y.G., M.E.H., A.C., J.G.W., and R.E.G., collected, organized, and contributed to the quality control and management of JHS proteomic data, via both the SomaLogic platform and the Olink platform. D.H.K., U.A.T., D.N., M.D.B., J.M.R., D.E.C., S.D., L.F., S.S., D.S, R.P.T., P.D., K.D.T., Y.L., W.C.J., X.G., J.Y., Y.-D.I.C, A.W.M., S.S.R., J.I.R., and R.E.G., collected, generated, organized, and contributed to the quality control and management of MESA proteomic data. D.H.K., U.A.T., D.N., M.D.B., J.M.R., D.E.C., S.D., L.F., S.S., D.S, C.B., M.A.S., and R.E.G. collected, generated, organized, and contributed to the quality control and management of HERITAGE Family Study proteomic and genomic data. D.J. and TOPMed performed W.G.S. from JHS and MESA. D.H.K., U.A.T, A.G.B, A.P., Z.Y., P.N., and R.E.G. developed the WGS analysis pipeline and statistical methods across all proteomics data. D.H.K., U.A.T., A.G.B., P.N., J.G.W., and R.E.G. analyzed the data and wrote/revised the article.

Disclaimer

The views expressed in this article are those of the authors and do not necessarily represent the views of the NHLBI, the NIH, or the US Department of Health and Human Services.

Supplemental Material

Supplemental Methods

Figures S1–S4

Tables S1–S7

References 71–75

Footnotes