Circulating causal protein networks linked to future risk of myocardial infarction

Domain: neuroscience
Published: 2025-12-18
Editorial outcome: minor_revision

Abstract
Variations in blood protein levels have been linked to numerous complex diseases, including cardiovascular conditions. These associations highlight the intricate interplay between local and systemic factors in cardiovascular disease development, emphasizing the need for a comprehensive, systems-level understanding of its etiology. To address this, we develop a causal network inference framework using data from one of the largest serum proteomics studies to date, comprising measurements of 7523 serum proteins in the prospective, population-based Age, Gene/Environment Susceptibility-Reykjavik Study (AGES) cohort of 5376 older adults. Using cis-acting protein quantitative trait loci (pQTLs) as instrumental variables within a causal inference framework designed to mitigate hidden confounding, we identify 185 high-confidence causal serum protein subnetworks collectively interacting with 5611 targets. Several subnetworks, many forming hierarchical frameworks of directional relationships, are significantly associated with multiple cardiometabolic traits and with future risk of myocardial infarction and its long-term complication, heart failure.

Introduction
Atherosclerotic cardiovascular disease (ACVD) is the leading cause of age-standardized deaths globally1. While lipid-lowering treatments have been shown to reduce the risk of ACVD2, residual risk persists3–5, underscoring a significant and unmet medical need. Early coronary atherosclerosis, which advances to coronary artery disease (CAD), is the primary cause of ACVD. In its most advanced stage, coronary artery plaques may rupture, manifesting clinically either through myocardial infarction (MI) or stroke. The rate of coronary plaque growth is influenced by various systemic factors across multiple organs systems, such as the immune system driving systemic inflammation, the liver regulating lipoprotein metabolism, and adipose tissue and skeletal muscle contributing to the development of obesity and type 2 diabetes (T2D)6. Other contributing factors are endocrine signaling and hemodynamic processes6. The rate at which CAD progresses depends on the interplay between these systemic factors, ultimately leading to plaque rupture, thrombus formation, and MI6,7. The complex etiology of ACVD is shaped by both genetic and lifestyle factors8, which are mediated through interactions between multiple organ systems9. The molecular disruptions across organs that contribute to ACVD have primarily been studied in isolated pathways using model systems. While these studies offer valuable insights into disease etiology and treatment, they fail to capture the systemic complexity of ACVD. In other words, the rate at which ACVD progresses to become clinically significant depends on the interplay of both local (cardiovascular) and systemic (non-cardiovascular) factors, an aspect of the disease that is often overlooked. Thus, a broader systems-level approach is required to obtain a more holistic understanding of its etiology.

Tissue-specific and cross-tissue transcriptional networks have been established as both undirected networks10–14 and directed networks15,16, operating within and across various tissues, and demonstrating robust associations with complex diseases. In contrast to undirected co-regulatory networks, directed networks have the potential to differentiate between causal and reactive nodes and elucidate how these causal nodes propagate their effects17. Gene expression quantitative trait loci (eQTLs) have been employed as genetic instruments in causal inference analysis and as priors in reconstructing transcriptional Bayesian networks15. These variants offer an effective means of natural perturbation to infer causal relationships between genes and higher-order phenotypes like disease, and even between genes themselves12,14. This has been well documented by the explosion of Mendelian Randomization (MR) studies that use genetic variants as instrumental variables for both molecular and higher-order phenotypes18. Although reconstructing directed causal networks has traditionally been computationally intensive and limited to small-scale systems, recent advancements have significantly enhanced performance16,19–21. These improvements have made the process orders of magnitude faster, enabled the coverage of a larger proportion of variance in the data, and proven especially effective when both genotype and molecular node data are available for the same sample19,20.

Proteomics has recently advanced to the point where high-throughput measurements allow for the analysis of thousands of proteins from a single tissue or biofluid sample in large population studies22,23, exposing the depth and complexity of the plasma and serum proteomes. These recent advancements have largely been driven by the highly sensitive aptamer-based affinity methods. In fact, comparisons between various proteomics platforms highlight the superior performance of the aptamer-based platform, especially in terms of detection precision and sensitivity24. Using this approach led to the identification of the first human protein co-regulatory network, reconstructed from the analysis of 4137 serum proteins measured in 5457 older adults of the prospective, population-based AGES cohort23. Furthermore, we demonstrated that the network modules for circulating proteins are under strong genetic control and are linked to a broad range of past, current and future disease states23. Notably, the structure and composition of these serum protein networks are stable and have been validated in other body fluids, such as in plasma25 and cerebrospinal fluid26. Finally, unlike solid tissue networks, the serum protein network consists of modules incorporating proteins synthesized by many, if not all, tissues across the body23,27.

In this study, we describe the reconstruction of a directed circulating Causal Protein Network (CPN) within the AGES cohort using an expanded dataset comprising 7523 serum proteins encoded by 6586 genes. Here, cis-acting protein quantitative trait loci (pQTLs) were employed as instrumental variables to highlight causal interactions between protein pairs. Applying various filters, including network size and genome linkage disequilibrium (LD) between instrumental variables, we identified 185 CPN subnetworks encompassing at least 10 protein members and interacting with a total of 5611 target serum proteins. The CPN subnetworks were examined for their relationships with each other and tested for associations with various ACVD-related outcomes. We highlight the subnetworks with the strongest associations to incident MI, along with cardiometabolic traits that contribute to the risk of ACVD.

Results
Studypopulationandanalysisoverview
Thisstudybuildsonthepopulation-basedprospectiveAGEScohortofolderadults(N=5764,meanage76.6±5.6years,agerange66–98years,57%female).Thecohortisrichlyannotatedwithdataondiseaseriskfactors,endpoints,comorbidities,genotypeinformation,anddeepserumproteomics23,28.Table1displaysthedemographic,biochemical,clinical,physiological,andanthropometricdata,aswellascardiovascularimagingmeasurementsofparticipantsintheAGESstudy,analyzedfor7523serumproteinsandstratifiedbyincidentMIafterexcludingallprevalentMIcases.Thefollow-upperiodfornewlydiagnosedMIpatientsextendedupto12yearsfrombaseline,whereperson-yearsoffollow-upwerecalculatedfromthefirstAGESvisituntilthedateofdiagnosis,death,ortheendofthefollow-upperiod,withamedianof5.6[2.8,8.2]yearsfortheincidentMIgroup(Table1).Asanticipated,severalmeasuresassociatedwithanincreasedriskofanMIeventweresignificantlyalteredintheincidentMIgroupcomparedtothenon-MIgroup(Table1).Theseincludeahigherprevalenceofmetabolicsyndrome(MetS)andtype2diabetes(T2D),aswellaselevatedcoronaryarterycalcium(CAC)andcarotidarteryplaqueburden(Table1).Furthermore,25.7%oftheincidentMIgroupdevelopedheartfailure(HF),comparedto4.6%ofthenon-MIgroup(P=3×10−55)(Table1).Figure1presentsanoverviewofthestudyanditsworkflow,includingthereconstructionofthecirculatingCPNwithintheAGESstudy,aswellaskeydiseaseendpointsandcardiometabolictraitsassociatedwithACVDthatareexaminedinthepresentstudy.AdditionaldetailsonthedataandanalysesareprovidedinSupplementaryFig.S1.

Reconstructionofthecirculatingcausalproteinnetwork
We reconstructed a global network of circulating proteins using a causal inference framework, inferring edges from pairwise protein relationships, with cis-acting pQTLs at a false discovery rate (FDR) < 5% serving as instrumental variables (Methods). A total of 5662 proteins with a cis-acting pQTL instrument, referred to as A-proteins (Fig. 1), were identified. For each A-protein, we calculated the posterior probability of it having a causal effect on the serum levels of each of the remaining 7523 proteins (Methods), referred to as B-proteins or targets, yielding approximately 42.5 million potential network edges. From all highly significant edges (FDR=1%), we defined each network regulator and its targets as a subnetwork of the global CPN. At this stage, nearly half of all network regulators had only a single target and at more permissive FDR thresholds, the proportion of A-protein with just one target decreases further (Supplementary Fig. S2). Since our primary focus was on regulatory proteins that accounted for the most variation in the serum proteome, we selected A-proteins with a minimum of 10 targets, leading to a global network consisting of 43,528 edges and 234 A-proteins. For A-proteins with multiple aptamers, we selected those with the largest number of targets, and further refinement of LD among cis-acting instruments resulted in the final CPN comprising 185 A-proteins (referred to here as network regulators), their corresponding subnetworks, and a total of 31,358 edges (Supplementary Data 1 and Fig. 2a). We identified a high number of indirect regulations (Fig. 2a), where two network regulators (x and y) are responsible for the regulation of a common set of targets (z), but y is also a target of x. Such motifs are known as feed-forward loops (FFLs) and are a common feature in biological networks29.

Independent cis instruments accounted for an average of 7.4% of the variance in protein levels across networks, with some cis signals contributing up to 84% of the variance (Methods). Furthermore, by utilizing cis signals for parental nodes (Methods, Supplementary Text), we observed a correlation between the number of regulators a target protein has, and the proportion of total variance in target protein expression explained by the cis signal of the network regulator (Spearman r=0.78) (Supplementary Fig. S3). In some cases, more than 50% of the variance in the expression of the target protein was explained solely by the cis-acting pQTLs of the parent proteins, with no contribution from local cis-components (Supplementary Text).

High robustness and edge precision in the circulating causal protein network
We found the CPN to be robust in response to hub removal (Fig. 2b), with the largest connected component size still containing more than 80% of all nodes following the removal of the top 10 largest hubs, suggesting high levels of co-regulation and biological redundancy. The CPN also exhibited the typical “scale-free” property of biological networks30, where a small number of regulators have a very high number of targets (Fig. 2c, d). We also examined the distribution of incoming edges and found that many of the network regulators with a high number of outgoing edges also have many incoming edges, further highlighting the interconnectedness of the global CPN (Fig. 2e, f). Network robustness analysis and Precision-Recall assessments for the networks were conducted across various FDR thresholds and sample sizes. This involved random sampling of AGES individuals at varying sample sizes (Methods) and demonstrated a mean area under the receiver operating characteristic (ROC) curve exceeding 90%, even when the sample size was reduced to 2000 (Fig. 2g). However, when represented as a Precision-Recall curve, the precision at an FDR of 1% declined sharply to higher recall levels for sample sizes below 3000 and was even lower at more permissive FDR thresholds (Fig. 2h and Supplementary Fig. S4). Lastly, pairwise correlations between network targets were significantly stronger (Kruskal-Wallis P-value < 10−300) compared to randomly selected protein groups of the same size (Supplementary Fig. S5). In summary, these findings emphasize the strengths of the AGES study’s large sample size, highlighting the robustness and accurate edge estimation of the identified protein network structure.

Hierarchical organization of the circulating causal protein network
We used a greedy heuristic to derive a directed acyclic graph (DAG) from the 1622 interactions between network regulators (Methods, Fig. 2a), allowing for an interrogation of the hierarchical structure between said regulators. This process led to the removal of 572 edges between network regulators, resulting in an acyclic network with 1050 interactions, and the nodes then arranged at different levels according to their distance from the roots of the network (Fig. 3). We then reintroduced the removed edges to the CPN, while maintaining the ordered layout, and found much of the hierarchical structure to be conserved within the different levels (Supplementary Fig. S6). Most edges between proteins were observed within single levels, either above or below, with a few crossing multiple levels, indicating a highly structured organization within these networks. A hierarchical structure with a small number of levels has been observed in other biological networks30,31.

Our approach infers causal relationships between protein pairs (A→B), whether directly or indirectly mediated by a third variable (e.g., A→C→B). The abundance of FFLs and within-level interactions may therefore reflect the high sensitivity of the AGES proteomics data for detecting weaker indirect causal effects. To address this, we generated a version of the directed acyclic graph (DAG) with all indirect edges removed, known as the transitive reduction of the DAG (Methods). This simplified, cascade-like network, comprising 255 edges, showed that many previously highly interconnected nodes remained prominent even after indirect edges were excluded (Supplementary Data 2).

Protein-protein interaction (PPI) networks are highly modular and feature connected hub proteins, reflecting a hierarchical biological organization32. Using the human integrated protein-protein interaction reference (HIPPIE)33 database, we identified edges in the serum protein network overlapping with direct physical PPI networks, as outlined in the Supplementary Text. Specifically, across all confidence thresholds, we found that the true CPN captured significantly more PPI edges than the average observed in random networks, including the highest confidence PPIs (z-score=24.3, P-value < 0.001) (Methods, Supplementary Text, Supplementary Fig. S7). This comparison highlights that a substantial portion of the CPN may reflect direct physical interactions, most of which have been identified in vitro from solid tissues33. However, unlike the CPN, the PPI network database cannot capture interactions that span tissue boundaries within the context of the serum.

Causal protein networks linked to ACVD related outcomes
We assessed the association of each network regulator and its corresponding CPN eigenprotein with incident MI and ACVD-related cardiometabolic traits, including MetS, T2D, CAC, carotid artery plaque burden, and incident HF (Methods). For this analysis, the first principal component (PC1) of each of the 185 subnetworks was calculated, and PC1s accounting for more than 15% of the variance within their respective subnetworks were designated as eigenproteins. Three CPN subnetworks did not meet this criterion and, therefore, did not have valid eigenproteins. Overall, 50 network regulators and 36 eigenproteins were significantly associated with incident MI (Supplementary Data 3, 4). Figure 4a–f illustrates the correlation between network regulators and their corresponding eigenproteins, highlighting their associations with various ACVD-related outcomes. For example, both the network regulator inter-alpha (globulin) inhibitor 3 (ITIH3) and the eigenprotein for that subnetwork were significantly associated with carotid plaque burden (Fig. 4c). When considering the network regulators, seven CPN subnetworks showed significant links to all six ACVD-related traits (Fig. 5a), while for the eigenproteins, only one CPN subnetwork (ITIH3) was associated with all six traits (Fig. 5b).

We ranked these subnetworks based on their associations with incident MI and other ACVD-relevant outcomes, assigning equal weight to both the network regulator and its corresponding eigenprotein (Supplementary Data 5, 6). This approach identified 25 top-ranked subnetworks (arbitrary rank score ≥ 7; Table 2). Clustering the network regulators based on their associations with these outcomes revealed that the six core traits align with the three established etiological axes of MI (Fig. 5c), as depicted in Fig. 1. In contrast, the corresponding eigenprotein associations drive a distinct clustering of phenotypes, grouping CAC, carotid plaque burden, and MI together (Fig. 5d), which represent the key outcomes associated with atherosclerosis34. The top-ranked CPN subnetworks in Table 2 were largely retained when a stricter 30% eigenprotein variance threshold was applied (Supplementary Data 7, 8). Associations of all top ranked network regulators with incident MI and/or related traits are illustrated in Supplementary Fig. S8. We note that a substantial number of the 182 subnetworks, including 62 network regulators and 121 eigenproteins, showed only one or no significant associations with ACVD-related traits (Supplementary Data 5, 6). In summary, both the top ranked network regulators and their corresponding eigenproteins capture the established multi-dimensional axes of ACVD development, albeit through slightly different mechanisms.

The top ranked subnetworks exhibit a strong degree of interconnectivity
Having identified subnetworks of the global CPN that are associated with different aspects of ACVD, we were interested in comparing structural similarities between these subnetworks. Among the 25 top ranked subnetworks (Table 2), we find examples of clusters of proteins that are co-regulated by several network regulators (Fig. 6a). Furthermore, target similarity analysis identified clusters of CPN subnetworks that are both correlated through their eigenprotein and share similar targets (Fig. 6b). Interestingly, we find a large group of negatively correlated subnetworks which shares similar targets, indicating co-regulation by distinct network regulators with opposing directional effects.

Like the global CPN, we observed multiple levels of regulation among the network regulators of the top ranked networks (Fig. 6c). The connected subnetworks form distinct modules that largely preserve this structure, even when the previously removed edges are reintroduced (Supplementary Data 1). This includes, for example, a cascade involving the network regulators C2, CFB, and GAL3ST1 (Fig. 6c), as well as the CFHR1 network regulator when applying a more stringent eigenprotein variance threshold (Supplementary Fig. S9). This is noteworthy because CFHR1, C2 and CFB are all integral components of the same complement cascade, whereas GAL3ST1 is an enzyme involved in the biosynthesis of glycosaminoglycans, which play vital roles in various biological processes, including inflammation and vascular remodeling35. The role of the complement system in inflammatory processes and cardiac health is well-established, particularly in relation to myocardial tissue injury36. Another intriguing example involves a pathway of five proteins: COL28A1, NUDT21, PTPN11, HSPA1B, and HSPA1A (Fig. 6c). Notably, all but COL28A1 are enriched in the hematopoiesis pathway (g:Profiler, P-value=0.018), essential for blood cell formation and development from hematopoietic stem cells. This process is critical for maintaining homeostasis and has broad implications in atherosclerosis37. Interestingly, mutations in the mouse Col28a1 gene are linked to abnormal blood vessel morphology, according to the MGI database38. These cascade-related proteins may collectively influence inflammation, vascular remodeling, atherosclerosis, and blood vessel integrity, supporting their observed associations with incident MI and related traits in this study.

Finally, comparison of the top ranked CPN subnetworks with the previously published serum protein co-regulatory network from the AGES study23, revealed substantial overlap between protein clusters in both network types (Supplementary Data 9, Supplementary Fig. S10, Supplementary Text). This overlap is evident in two ways: first, many CPNs align with the same co-regulatory module; second, when a single CPN intersects with multiple co-regulatory modules, these modules often belong to the same super-cluster of correlated co-regulatory modules (Supplementary Text). Overall, a significant relationship exists between the circulating CPNs and the co-regulatory networks, despite fundamental differences in the methodologies used to reconstruct them.

Replication of causal protein network architecture in an independent cohort
We performed a validation analysis of the CPN structure using the UK Biobank (UKBB) study39 as an independent dataset. There are 2186 protein assays from the UKBB Olink Explore and the AGES 7K SOMAscan profiling platforms that target the same proteins according to their annotations (Methods). Therefore, we reconstructed CPNs in both AGES and UKBB using only this subset of proteins from each platform to ensure that these networks were comparable. As correlations in protein measurements between these two different proteomic platforms have been shown to vary40, we conducted instrument selection separately in the UKBB and AGES studies (Methods), removing duplicate aptamers targeting the same protein before analysis.

We first reconstructed a subsetted CPN in AGES but restricted it to the 2186 proteins also measured in UKBB. Of these protein assays, 1835 had a valid cis-acting pQTL instrument to be used as a potential network regulator. At a 1% FDR threshold, we identified 454 proteins with at least one target and 43 with at least 10 targets in the AGES study. Among the 454 proteins with at least one target in the subsetted AGES CPN (FDR=1%), 416 were also present in the full AGES CPN constructed from all measured proteins. Furthermore, 96.3% of the edges identified in the subsetted AGES CPN overlapped with those in the full CPN. We then reconstructed a comparative CPN using data from 52,543 UKBB participants with Olink assays targeting the same 2186 proteins. We observe a significant overlap between network regulators in the two cohorts, including the 454 proteins with at least one target in AGES, of which 393 were present with at least one target in UK Biobank (Fisher’s Exact Test; OR=1.9, P=4.1×10–6). Further, we identified 629 edges shared between the AGES subsetted and UK Biobank CPNs (FDR=1%). To formalize this comparison, we applied Fisher’s Exact Test, which demonstrated that the overlap was highly significant (OR=2.9, P=1.6×10–99), using all possible edges as the background while excluding self-edges. Notably, many more targets were identified for a given CPN in UKBB than in AGES (Supplementary Data 10), likely to reflect the nearly ten fold larger sample size in UKBB.

We next sought to determine the extent to which protein targets of subnetworks reconstructed in AGES and those from UKBB overlapped. Of the 393 shared network regulators from the comparative subnetworks, we found 36 that shared at least 2 targets in the AGES and UKBB comparative CPNs, for which we assessed overlap significance with a hypergeometric test, using the full set of 2186 proteins as background. Overall, 23 subnetworks had a significant overlap in targets across both platforms (q<0.05), including three of the top-ranked subnetworks previously linked to ACVD in the full AGES CPN, namely, KLKB1, CFB and CSF3 (Supplementary Data 10 and Supplementary Fig. S11).

A recent study compared correlations between Olink and SOMAscan platforms in the Icelandic population40. From this study, we obtained correlations for 1677 of the 2186 proteins measured (mean Pearson coefficient=0.39). For 19 of the 23 significantly overlapping subnetworks, correlations were available, with a higher mean of 0.48 (Supplementary Data 10). In some cases, network regulators of significantly overlapping CPNs were not correlated across platforms (Supplementary Data 10), suggesting that even when different affinity-based technologies bind distinct epitopes or proteoforms, the network regulator remains robust, continuing to capture the same underlying biology reflected in its downstream targets.

Additional links between the top ranked networks and ACVD
We investigated several complementary lines of evidence to support the role of the top-ranked network regulators in the pathophysiology of MI and related cardiometabolic traits. Network regulators and targets from the top-ranked subnetworks showed enrichment in multiple pathways previously linked to ACVD, including a shared enrichment in the cellular response to heat shock (Supplementary Data 11, Supplementary Figs. S12–S13, Supplementary Text). STRING-based analysis revealed significantly (P-value=0.047) enriched functional and physical interactions among the top network regulators (Supplementary Fig. S14), including previously found (Fig. 6c) and previously unreported connections (e.g., ITIH3 with KLKB1, APOA5, and AFM; Supplementary Fig. S14).

Causal inference analysis of top-ranked regulators. Two-sample MR and colocalization analyses highlighted several top-ranked regulators, including APOA5, DDX39B, HSPA1A, C11orf49 and LRP4, with significant support for causal associations with ACVD-related traits (FDR<0.05; Supplementary Data 12–13). Here, APOA5 demonstrated a strong protective effect against MI and MetS, consistent with findings that APOA5 knockout mice display a fourfold elevation in plasma triglyceride levels41, a significant risk factor for ACVD42. Further, genetically determined levels of HSPA1A showed an association with MetS, indicating a causal role, whereas DDX39B was linked to both MetS and T2D. Notably, the HSPA1A gene has also been causally implicated in T2D and its microvascular complications in prior MR and colocalization analyses43. Sensitivity analyses, restricted to exposures with at least three instruments (Methods), indicated that results were largely robust across approaches (Supplementary Data 14), although significant instrument heterogeneity was observed for APOA5 in relation to MetS.

Additional context from current and prior evidence. Beyond the analyses of the present study, several top regulators, including ITIH3, HSPA1B, KEAP1, HSPA1A, C2, CSF3, FABP3, KLKB1, AFM, APOA5, PTPN11, CFB, and COL28A1, have previously been linked to ACVD or related cardiometabolic traits (Supplementary Text). For example, the top-ranked regulator ITIH3 has a gain-of-function genetic association with increased risk of MI and is expressed in vascular smooth muscle cells and macrophages within atherosclerotic plaques44. Notably, ITIH3 is positioned at the base of the data-driven CPN (Fig. 3), indicating a potential role in numerous regulatory pathways. Another example is KLKB1, which we recently demonstrated through MR analysis to be causally protective against calcific aortic valve disease45, a leading cause of heart failure via aortic stenosis46. To further characterize the top-ranked protein network regulators, we assessed their effect directions on ACVD-related traits and survival outcomes (Supplementary Data 15) and evaluated if they are druggable or targeted by available compounds using public databases (Supplementary Data 16). Table 3 summarizes the accumulated findings, offering an integrated overview of directional consistency across causal, observational, and therapeutic evidence for each regulator.

Therapeutic pathways and clinical implications. Over half of the top-ranked network regulators are potentially druggable with small molecules or biologics, although only a few have compounds currently available for either related or unrelated disease indications (Supplementary Data 16, Supplementary Text). Specifically, 13 of the 25 regulators are druggable, with 6 already targeted by licensed or investigational drugs. Building on this, the CPN analysis has identified both interconnected networks and therapeutic pathways with potential for early ACVD detection and prevention: (1) Complement-mediated inflammation: The network regulators C2 and CFB are central to the alternative complement pathway, driving vascular inflammation and atherosclerosis50. CFB is targeted by the antisense oligonucleotide Sefaxersen, in development for IgA nephropathy (Supplementary Data 16). Targeting this pathway could prevent early inflammatory progression of atherosclerotic lesions. (2) Cellular stress response and cyto-protection: KEAP1 regulates the Nrf2 oxidative stress pathway47, while HSPA1A and HSPA1B control protein folding under stress51. All three have licensed drugs targeting them (Supplementary Data 16). Modulating this pathway may protect against oxidative damage and endothelial dysfunction. (3) Hemostasis: The key target KLKB1, a regulator of plasma kallikrein in the intrinsic coagulation pathway52, is targeted by the approved drug Berotralstat against angioedema (Supplementary Data 16). Our integrative analysis, however, implies that activating KLKB1 may offer a strategy to prevent thrombotic complications in ACVD (Table 3). 4) Lipid metabolism: The network regulators APOA5 and AFM control triglyceride-rich lipoprotein metabolism and oxidative stress41,53,54, are bio-druggable (Supplementary Data 16), and could be targeted to prevent dyslipidemia-driven initiation of atherosclerosis. This network framework highlights causal pathways for ACVD prevention, with concordant data strengthening target confidence and inconsistencies pointing to opportunities for novel therapeutic discovery.

Methods
Studypopulation
Cohort participants aged 66 through 96 years at the time of blood collection were from the AGES28, a single-center, prospective, population-based study of older adults (N=5764, mean age 76.6 ± 6 years). AGES was formed between 2002 and 2006, and its participants were randomly selected from the surviving members of the established 40-year-long population-based Reykjavik study60,61. The Reykjavik study, a prospective cardiovascular survey, recruited a random sample of 30,795 adults born between 1907 and 1935 who lived in the greater Reykjavik area in 1967, that were examined in six phases from 1967 to 199660,61. The AGES measurements, which include for instance, brain and vascular imaging, are designed to assess four biologic systems: vascular, neurocognitive (including sensory), musculoskeletal, and body composition/metabolism28. All AGES participants are of European ancestry. AGES was approved by the National Bioethics Committee in Iceland, that acts as the institutional review board for the Icelandic Heart Association (approval number VSN-00-063, in accordance with the Helsinki Declaration) and by the US National Institutes of Health, National Institute on Aging Intramural Institutional Review Board.

The study comprised MI patients who met the MONICA criteria for definite MI as previously described60. The criteria for HF were based on clinical symptoms and signs, chest X-rays, and, in many cases, echocardiographic findings from hospital records, which were adjudicated by examining every record for both prevalent HF, i.e., had HF at the baseline visit, and incident HF, i.e., HF diagnosed after the baseline visit. The incident HF cases were free of HF diagnosis at the baseline visit, but who were later hospitalized and diagnosed (hospital discharge ICD-10 diagnosis codes starting with I50) with HF during the follow-up period of eight years. Each patient’s thorough medical records were subsequently adjudicated by a cardiologist to confirm the diagnosis of symptomatic HF, and the date of the incident HF event documented. Among the criteria were symptoms such as shortness of breath that could be ambulatory, and signs of pulmonary edema. Type two diabetes (T2D) was determined from self-reported diabetes, diabetes medication use, or fasting plasma glucose ≥ 7 mmol/L according to American Diabetes Association guidelines62. Metabolic syndrome (MetS) was defined by three or more of the following: Fasting glucose ≥ 5.6 mmol/L, blood pressure ≥ 140/90, triglycerides ≥ 1.7 mmol/L, high-density lipoprotein (HDL) cholesterol < 0.9 mmol/L for males or < 1.0 mmol/L for females, body mass index (BMI) > 30 kg/m2. Systolic and diastolic blood pressure were measured twice with subjects in a supine position using a Mercury sphygmomanometer. Lipoproteins and plasma glucose levels were measured on fasting blood samples. Triglycerides (TG) were measured using enzymatic colorimetry (Roche Triglyceride Assay Kit), HDL cholesterol with an enzymatic in vitro assay (Roche Direct HDL Cholesterol Assay Kit), and glucose was measured using photometry (Roche Hitachi 717 Photometric Analysis System). Coronary artery calcium (CAC) was quantified using the Agatston scoring method63, which was reviewed independently by four image analysts. Phantom-adjusted CAC was expressed as a sum score for all four coronary arteries, as previously described in greater detail64. The use of ultrasound imaging to assess the presence and severity of carotid plaque in the AGES population has been detailed elsewhere65,66. Hepatic steatosis was assessed by computed tomography (CT), serving as a non-invasive proxy for non-alcoholic fatty liver disease (NAFLD), as previously described67. We assessed overall survival (2982 events), as well as survival following incident CHD (692 events) and incident MI (299 events). For overall survival, follow-up was defined as the time from study entry in AGES until death from any cause or the end of follow-up (December 2016). For survival after incident CHD or MI, follow-up was defined as the period starting 28 days after the incident event until death from any cause or the end of follow-up.

Proteomics profiling
Blood samples were collected at the AGES-Reykjavik baseline visit after an overnight fast, and serum samples prepared using a standardized protocol and stored in 0.5mL aliquots at −80°C. Serum samples collected from the inception period of AGES, i.e., from 2002 to 2006, were used to generate proteomics data used in this study. Before the protein measurements were performed, all serum samples from this period went through their first freeze-thaw cycle. Protein levels in serum from 5376 individuals of the AGES were determined using a multiplex SOMAscan proteomic profiling platform, which employs aptamers or Slow-Off rate Modified Aptamers (SOMAmers) that bind to target proteins with high affinity and specificity. Here, 7523 aptamers mapping to 6,586 UniProt IDs were measured in a total of 8592 samples (two time points) using the SomaScan_v4.1 platform68. The aptamer-based platform measures proteins with femtomole (fM) detection limits and a broad detection range (> 8-log dynamic range) of concentration69. To prevent biases related to batch or processing time, the order of sample collection and separate sample processing for protein measurements were randomized, and all samples run as a single set at SomaLogic Inc. (Boulder, CO, US). All aptamers that passed quality control had median intra-assay and inter-assay coefficient of variation, CV < 5%. Hybridization controls were used to correct systematic variability in detection and calibrator samples of three dilution sets (20% (1:5), 0.5% (1:200), and 0.005% (1:20,000)) were included so that the degree of fluorescence was a quantitative reflection of protein concentration. The adaptive normalization by maximum likelihood (ANML) method was employed to normalize Quality Control (QC) replicates and samples using point and variance estimations from a healthy external reference population (n=1000). Consistent target specificity of aptamers was indicated by direct (through mass spectrometry) and/or indirect validation23.

Some proteins were targeted by more than one aptamer. In such cases, individual aptamers had distinct binding sites (epitopes) or binding affinity23. The single gene NPPB, for example, produces three protein products: full-length BNP, NT-proBNP, and BNP32, each of which are targeted by different aptamers. Duplicate aptamers to single pass transmembrane proteins (one to the extracellular domain and another to the intracellular loop), aptamer targeting multimers (e.g., interleukins), and duplicate aptamers generated in distinct expression systems are further examples. Finally, 233 aptamers were derived from mouse-human protein chimeras (used as SELEX input) to target proteins from both species.

Before the analyses, protein data were centered, scaled, and Box-Cox transformed70, and extreme outlier values were excluded, defined as values above the 99.5th percentile of the distribution of 99th percentile cutoffs across all proteins. Prior to reconstruction of the causal protein networks, the data were adjusted using a linear model to account for age and sex.

Genotype data and the detection of cis-acting variants
The genotype data includes assayed and imputed genotype data for 5656 AGES participants71. The genotyping arrays used were Illumina Hu370CNV and Illumina GSA BeadChip, which were quality controlled by eliminating variants with call rates < 95% and Hardy Weinberg Equilibrium (HWE) P-value < 1 × 10−6. The arrays were imputed against the Haplotype Reference Consortium imputation panel r1.1, and post-imputation quality control was performed separately for each platform. Variants with imputation quality r2 < 0.7, minor allele frequency < 0.01, as well as monomorphic and multiallelic variants, were removed before merging the platforms to generate a dataset with 7,506,463 variants for 5656 AGES individuals as previously described71. These variants were associated to each of the aptamers on the v4.1-7k serum protein panel to identify proximal (cis-acting) pQTLs, in the same way as previously described71. We applied a 300kb genomic window spanning each protein-expressing gene in the v4.1-7k serum protein panel to map out cis-acting pQTLs after accounting for the number of single-nucleotide polymorphisms (SNPs) in each window. We then corrected for multiple testing using the Storey-Tibshirani procedure for q-value estimation72. The cis-acting pQTLs serve as instrumental variables for the reconstruction of the causal serum protein networks, which are described below.

Reconstruction of the circulating causal protein network
A circulating Causal Protein Network (CPN) was reconstructed from causal relationships between serum proteins derived using a Mendelian Randomization (MR) framework73. In this model, causal relationships between protein pairs are estimated between an exposure protein (A) and outcome protein (B), using a cis-acting pQTL (E) for A as a causal instrument. We selected all proteins which had at least one valid cis-acting pQTL as our A-proteins at FDR ≤ 5%, after correcting for multiple testing using the Storey-Tibshirani procedure72. In each pairwise causal test, the lead protein regulatory SNP (pSNP) for A (lowest P-value) was selected as E to be used as an instrumental variable in accordance with the core instrumental variable assumptions73,74: (1) the instrumental variable should be strongly associated with the exposure (A), (2) the instrumental variable should only be associated with the outcome (B) through the exposure (A), and (3) the instrumental variable should not be associated with any potential confounders affecting the exposure (A) and the outcome (B).

Causal estimates between proteins were measured using the tool Findr20 (version 1.0.8) in Python (version 3.10.6) using individual-level protein expression levels for A and B and genotypes for E. Before inferring casual estimates, the data were transformed using a rank-based inverse normal transformation within the Findr package. The product of the secondary linkage test (P2) and controlled test (P5) from Findr20 was used to estimate the posterior probability (PP) of P(A→B). The secondary linkage test measures the PP of association between E and B, and the controlled test assesses the PP of the dependence between A and B following adjustment with E to exclude that E has independent effects on A and B75. We estimated P(A→B) for all A-proteins with a valid instrumental variable, where B was every other protein in our dataset. We estimated a global FDR as 1 minus the mean of all PP for P(A→B) and then filtered PP P(A→B) to achieve the desired FDR threshold, as shown previously76. Networks were reconstructed from causal interactions that fell below this FDR threshold, where parent nodes were A and child nodes were B, with edges represented as P(A→B). In instances where there were multiple aptamers targeting the same protein, we selected the A proteins with the largest number of targets as the representative aptamer for this protein. When examining hierarchy between regulators, we converted edges between A proteins to a directed acyclic graph (DAG) using a greedy heuristic as implemented in Findr20, which selects the most significant edges in an iterative fashion to avoid cycles. This is done to identify any hierarchical structure between A-proteins, and once a DAG has been constructed the removed edges can be reintroduced to examine how well any hierarchical structure is preserved within the complete networks. Network visualization was performed using Cytoscape (version 3.10.1).

As pleiotropy is a concern in MR analyses, we took steps to account for instances where two or more exposure proteins either shared or had instrumental variables in high linkage disequilibrium (LD)74. We assigned proteins to the same LD block if they shared the same lead pSNP, or if they had different pSNPs in medium or high LD (r² ≥ 0.5). For all A proteins within an LD block, we calculated the intersection of targets as a ratio of the union of targets (I) in a pairwise manner between target sets. We resolved A proteins as independent networks when I < 0.6, in cases where I > 0.6 we collapsed the union of all targets within a single unresolved network. In cases where two or more A proteins were mutual targets of each other, we considered these networks as unresolved. Where there were more than two A-proteins in a single LD block, an A-protein must have I < 0.6 across all other A-proteins to be considered independent. We did not encounter any instances of A-proteins that were shared between separate LD blocks.

We defined protein subnetworks as a single A-protein and all targets of that A-protein, i.e., a single regulator and its targets. The average expression profile of a subnetwork was estimated as an eigenvector of all proteins in a subnetwork across all samples, described here as an eigenprotein77. We calculated the first principal component (PC) from the expression profiles of all proteins within a subnetwork using principal component analysis (PCA), via the PCA function from the Python library scikit-learn78 (version 1.1.2). If PC1 explained > 15% of the variance of the subnetwork, we used this principal component as an eigenprotein for this subnetwork in downstream analyses, and for subnetworks where the variance was < 15%, PC1 was not taken forward. Three CPN subnetworks did not meet this criterium and therefore did not have a valid eigenprotein.

Network robustness analyses
To examine the impact of varying sample size on network inference, we randomly sampled subsets of individuals in different sizes for the reconstruction of networks. At each sample size threshold (n=500, 1000, 2000, 3000, 4000, 5000), we randomly sampled three subsets from the pre-processed protein expression data and reconstructed a network for each subset using the previously outlined approach. We termed the primary network, which was reconstructed using all samples, as the ground-truth network. The ground-truth network was represented as a flattened binary matrix (A, B), where 1 was indicative of the presence of an edge between A→B, and 0 for the absence of an edge. Each sub-sampled network was represented as a similar flattened matrix (A, B), but where the values were populated by the PP of P(A→B). We tested how well each sub-sampled network captured the structure of the ground-truth network by calculating the Receiver Operating Characteristic Area Under the Curve (ROC AUC) for each sub-sampled network against the ground-truth network using the roc_auc_score function from scikit-learn (version 1.1.2)78. We then calculated the mean AUC of the repeated sub-sampled networks at each sample size threshold. We also calculated Precision-Recall for each sub-sampled network against the ground-truth network using the precision_recall_curve function, also from scikit-learn78.

Variance explained model
To identify independent cis-acting protein SNPs, a ±150kb window was specified around each lead variant to define the region of interest for association testing. Using the GCTA software (version v1.94.1)79, the cojo-slct parameter was applied using a forward model selection approach. The default collinearity cutoff of 0.9 was used, and the P-value threshold was set at 0.00763, which corresponds to the highest P-value maintaining an FDR below 5%. Having identified independent cis-acting protein SNPs for every network regulator, we estimated the proportion of variance in protein expression that could be explained by cis-acting pQTLs using multiple linear regression. For each protein, we fitted a linear model in R (version 4.3.2), where the genotypes for the independent cis-acting SNPs act as the explanatory variables for protein expression. The adjusted coefficient of determination (adjusted r²) from this model was used as an estimate of the variance explained by the cis component for each protein, and the mean adjusted r² was then calculated across all network regulators.

We also estimated the variance in target protein expression that could be explained by cis-acting pQTLs from the target regulators. More specifically, the influence of cis signals on network regulators (parental nodes) was assessed in relation to the expression of 5459 target proteins within the CPN, including 162 target proteins that also served as regulators for other proteins. For every target protein, we fitted a linear model, as described previously, using the genotypes of all independent cis-acting pQTLs for the regulators (parental nodes), in addition to any cis-acting pQTLs for the target itself. We then calculated the difference in variance explained by local and parental cis-acting pQTLs combined, to that explained by cis-acting genetic variation alone. If the target protein did not have any cis-acting pQTLs, then the variance explained by the cis component was set to 0. Any unresolved networks were excluded from this analysis.

Comparison with a referenced database of protein-protein interactions
We used the human integrated protein-protein interaction reference (HIPPIE) version 2.333 to identify experimentally derived protein-protein interactions (PPIs) that have been captured by the serum CPN described in this study. We accessed all 289,112 PPIs from the HIPPIE database, which have been scored as a weighted sum, based on the reliability of the evidence underpinning the interactions and the number of studies detecting each interaction. There are 262,346 PPIs scored at medium confidence (confidence score > 0.63) and 77,630 scored at high confidence (confidence score > 0.72), based on thresholds defined by HIPPIE authors. We then calculated the overlap between edges in the CPN with PPIs from HIPPIE at the different confidence thresholds. After which, we compared the number of common interactions between HIPPIE and the CPN to common interactions between HIPPIE and the edges from random networks. These networks were generated by randomly sampling proteins from the complete set of measured AGES proteins to produce the same number of edges as in the CPN. This process was performed 10 times, and the mean number of edges captured, in addition to the standard deviation, was calculated at each of the different confidence thresholds. We then calculated a z-score comparing the number of edges captured from the true vs random networks at each confidence threshold as: z=(observed overlap-mean random overlap)/stdev random overlap. This was then converted to a P-value as: P-value=1 – cdf(z-score), using the norm.cdf() function from Scipy Stats80.

Network topology analyses
We identified weakly connected components (i.e., groups of nodes in the CPN where every node can be reached from any other node, regardless of the direction of the edges) in the complete CPN and after removal of the regulators with the most targets, applying the weakly_connected_components function from the Graphs.jl (v1.12.0) package using Julia v1.11.1. We obtained a hierarchical layout of the A-protein DAG by first defining root nodes as A-proteins without incoming edges. We then defined the level of each A-protein as the shortest distance (number of edges in a shortest path) from a root node to the protein of interest. Shortest paths were computed using the dijkstra_shortest_path function from the Graphs.jl (v1.12.0) package using Julia v1.11.1. To further simplify the structure of the DAG, we performed transitive reduction using the “transitive reduction” function from the Graphs.jl (v1.12.0) package from Julia. Transitive reduction reduces the original DAG G to a new DAG G’ with the fewest possible edges such that, if there is a path from a vertex x to a vertex y in G, there must also be a path from x to y in G’, and vice versa.

Validating the network structure across proteomic platforms and cohorts
The UKBB is a prospective population study of 502,128 individuals from the UK who have been extensively characterized by genomic and phenotypic traits39. Genotype data is available for 486,593 participants, obtained using either the Applied Biosystems UK BiLEVE Axiom Array by Affymetrix or the Applied Biosystems UKBB Axiom Array. Imputation was carried out using the HRC reference panel, yielding allelic dosages for ~96 million variants, with markers annotated using the GRCh37 assembly of the human genome58.

There are 53,013 (53.9% female, mean age=56.8 years) UKBB participants who have undergone plasma protein profiling for 2923 proteins using the antibody affinity based Olink Explore 3072 protein extension assay platform59. We processed the Olink assay data using the approach outlined by Sun et al.59, from the initial study reporting plasma measurements from UKBB, whereby assays with more than 20% missingness are removed and the remaining missing values are mean imputed across samples. The overlap between processed samples and those with corresponding genotype data yielded 52,543 samples to be taken forward for network reconstruction. The remaining protein assays were then filtered for those that had a corresponding assay in the AGES SOMAscan 7K panel, yielding 2186 comparative protein assays. These assay measurements were then corrected for age and sex by fitting to a linear model as previously described for the AGES proteomic data.

Due to differences between the SOMAscan and Olink platforms, instrument selection was carried out independently for the UKBB proteomic data. For each protein, the lead cis-acting pSNP was chosen based on the P-value, using pQTL summary statistics from Sun et al.59. Reconstruction of causal protein networks was carried out with Findr and edges selected based on FDR thresholds, as has been previously described for the AGES discovery cohort. We tested the significance of overlapping targets between comparative subnetworks by hypergeometric testing using the hypergeom() function from Scipy Stats80 in Python, while using the intersection of measured proteins between platforms as a background. We then corrected for multiple testing using the Storey-Tibshirani procedure for q-value estimation72.

Colocalization and Mendelian randomization analysis
We used colocalization analysis to identify cis-acting pQTLs that shared a common signal with six different phenotypic traits using the R package coloc (version 5.2.3)81 and publicly available GWAS summary statistics from individuals with European ancestry (Supplementary Data 10). We first identified all proteins that had a cis-acting pQTL (FDR<5%) that shared at least 1 common cis-SNP with the GWAS trait (P<5×10-8) and extracted the summary statistics for all SNPs within ±150Kb of the transcription start site for the cognate encoding gene. We then estimated the probability of a shared causal variant (PP.H4) using approximate Bayes Factor colocalization82 via the coloc.abf (coloc) function using a threshold of PP > 0.9 and > 0.5. In addition, we tested colocalization between overlapping signals of cardiometabolic GWAS traits and serum protein cis-acting pQTL using the Bayes Factor colocalization with the Sum of Single Effects (SuSiE) framework for fine mapping of genomic loci83. The SNPs within the cis window were fine-mapped for both the GWAS and pQTL signal using the Susie rss() function from the SusieRR package (version 0.12.35), while limiting the number of possible independent signals (L=5) and the maximum number of iterations (n=5000). Only loci where convergence was observed for both signals were taken forward for colocalization. Colocalization was performed using the coloc.susie() function from the Coloc R package (version 5.2.3), which performs Bayes Factor colocalization. We selected colocalized signals (H4) using the posterior probability thresholds > 0.9 and > 0.5.

The causality of selected network regulators was assessed using a forward two-sample MR approach84. Genetic instruments (cis-acting pSNPs) were identified within ±150kb of the protein-coding region based on AGES data, filtered for statistical significance (P<0.05/number of cis-region SNPs), and matched to external GWAS outcomes. These SNPs were clumped for LD (r² < 0.1 and < 0.2) using PLINK v1.985, and AGES genotype data. MR estimates were calculated using the Wald ratio for single-SNP associations and generalized weighted least squares (GWLS) for multi-SNP analyses as previously described86. Significant results (FDR<0.05) were deemed causal candidates. To ensure the robustness of our findings, we performed sensitivity analyses of the MR estimates using MR-Egger87 and weighted median estimation88. For proteins with more than three instruments available, instrument heterogeneity was assessed using Cochrans´s Q89. Horizontal pleiotropy was assessed using MR-Egger87. The generalized GWLS was performed as previously described90, while other MR analyses were performed using the TwoSampleMR91 R package.

Statistical analysis and functional enrichment analysis
The relationship between serum protein (and eigenprotein) levels and quantitative phenotypes was evaluated using linear regression controlling for age and sex, whilst the relationship between serum protein (and eigenprotein) and prevalent disease was examined cross-sectionally using logistic regression with age and sex adjustments. The associations between serum proteins (and eigenprotein) and incident disease were evaluated longitudinally using the Cox proportional-hazards model92, with a median follow-up period of 5.6 [2.8, 8.2] years for incident MI. Associations with Benjamini-Hochberg FDR < 0.05 were considered statistically significant. To identify enriched gene ontology (GO) terms and pathways within network targets, we carried out formal function enrichment analysis using gprofiler2 (version 0.2.2), the R interface of g:Profiler (R version 4.3.2). We used a custom background of all measured proteins in the AGES aptamer-based assay and corrected for multiple testing using Benjamini-Hochberg correction with a cut-off of FDR < 0.05 to identify any enriched terms across the following categories; GO:MF, GO:BP, GO:CC, KEGG, REACTOME, Wikipathways, miRTarBase, Human Protein Atlas and Human Phenotype Ontology.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this manuscript, Bankier, Gudmundsdottir and colleagues present an impressive large-scale proteomic study of 5,376 older adults from the AGES-Reykjavik cohort, in which 7,523 serum proteins were quantified. The authors reconstruct a global causal network of circulating proteins and examine their relationships with cardiometabolic traits, including myocardial infarction (MI) and heart failure (HF).\n\nThe manuscript is clearly written and presents a substantial contribution to the field. The findings are likely to inform efforts in drug development and repurposing for cardiometabolic disease. That said, I believe the study could be further strengthened through several methodological refinements and additional contextual discussion, as outlined below.\n\nMajor comments\n1) Replication. The protein network reconstruction is data-driven and could be sensitive to cohort-specific features. Has the team considered validation in an independent cohort? Confirmation of the network structure or key associations in an external dataset would considerably enhance the robustness and generalisability of the findings.\n2) One-sample Mendelian Randomisation (MR) and risk of bias The use of one-sample MR introduces potential for bias and inflation of the type I error rate (see: Burgess et al., Wellcome Open Res, 2019;4:186). Where feasible, a two-sample MR approach would be preferable, drawing upon publicly available GWAS summary statistics for the cardiometabolic traits in question. This would reduce bias, increase statistical power, and facilitate a range of sensitivity analyses (e.g., MR Egger, weighted median) that are vital for assessing the robustness of causal inference.\n3) Colocalization methodology. The authors may wish to consider the application of more advanced colocalization approaches capable of modelling multiple causal variants, such as SuSiE or PwCoCo (e.g. PMID: 35452592). These methods are better suited to complex loci and may yield more reliable insights regarding shared genetic aetiology between proteins and traits.\n4) Inclusion of positive controls. It may be instructive to include a small number of well-established biological pathways or protein complexes—such as those involved in transcriptional regulation, the cell cycle, or intracellular transport—as positive controls. Demonstrating that the network approach successfully identifies such relationships would lend confidence to the method's validity.\n5) Ancestry and generalisability As the study sample is limited to individuals of European descent, the authors should more clearly acknowledge this as a limitation. While expansion to more diverse populations may not yet be feasible, explicit recognition of this issue will help orient future work towards greater inclusivity and global relevance.\n\nMinor comment\nLine 540: It would be helpful to define the symbol “E” for the benefit of readers less familiar with network analysis terminology.\n\nIn sum, this is a high-quality manuscript addressing an important topic, with clear translational relevance. The above suggestions are intended to strengthen the rigour and transparency of the work. I look forward to seeing the revised version.", "reviewer_2": "The manuscript presents a sophisticated causal network inference framework leveraging 7,523 serum proteins measured in 5,376 elderly individuals from the AGES-Reykjavik Study. The rationale for a system-level approach and the use of pQTLs for causal protein network inference is well-articulated, and the methodology is technically sound but highly complex. Below are specific concerns and suggestions for improvement.\n\n1. Abstract\nThe statement \"we used cis-acting serum protein quantitative trait loci (pQTLs) as instrumental variables, while accounting for potential unobserved confounding factors\" requires clarification. If this refers to Mendelian randomization (MR), the phrasing is misleading—MR mitigates unobserved confounding via genetic IVs but does not explicitly \"account\" for it. Please revise to reflect MR’s assumptions accurately.\n\n2. Results (Line 129–130)\n\"The follow-up period for newly diagnosed MI patients extended up to 12 years from baseline, with recurrent cases excluded from the incident group analysis. Person-years of follow-up were calculated from the first AGES visit until the date of diagnosis, death, or the end of the follow-up period...\" Incident vs. Recurrent Cases: If follow-up was censored at the first MI diagnosis (i.e., incident cases only), recurrent MI cases should not exist in this cohort by definition.\n\n3. Discussion\nThe clinical implications of prioritized proteins (Table 2) are underdeveloped. Given the focus on ASCVD early detection and prevention, readers would benefit from:\n- A brief summary of how many proteins have licensed/investigational drugs targeting them.\n- Potential therapeutic pathways suggested by the causal network.\n\n4. Methods\nClarify definition of NAFLD/MASLD criteria.\nGenomic Window Inconsistencies: Varying window sizes are used across analyses (300 kb for cis-pQTL detection, 150 kb for variance modeling/colocalization, 500 kb for MR IVs). Justify these choices.\nColocalization Assumptions: The analysis assumes a single causal variant. For robustness, consider methods accommodating multiple causal variants (e.g., PWCoCo, SharePro) for prioritized proteins.\nThe LD clumping threshold (r² ≥ 0.2) is lenient and may retain non-independent cis-pQTLs. Address whether stricter thresholds (e.g., r² < 0.1) is not feasible here.", "reviewer_3": "Summary\nThis manuscript by Emilsson and colleagues presents a comprehensive systems biology approach to understand the molecular outcomes underlying myocardial infarction (MI) and associated cardiometabolic conditions. By analyzing serum proteomic data through Somascan with samples from the AGES-Reykjavik study and employing a causal inference framework using cis-acting pQTLs, the authors reconstruct 185 high-confidence causal protein subnetworks (CPNs). These subnetworks were then associated with incident MI, metabolic syndrome, heart failure, and related traits. The authors attempted to validate their network using internal metrics and reference protein-protein interaction databases, and lastly performed colocalization and Mendelian Randomization (MR) analyses in order to strengthen their claims.\n\nI think that this study represents an important advancement in causal proteomics and has the potential to uncover novel biomarkers and therapeutic targets for cardiovascular disease. The greater than 5,000 samples afforded by the AGES cohort is compelling. Htowever, some conceptual and methodological weaknesses should be addressed to enhance the study’s impact and interpretability. importantl, the analytical framework is strong,but the study would benefit from independent validation, greater biological resolution, and experimental follow-up.\n\nMajor Strengths\nThe AGES-Reykjavik cohort provides deep phenotyping and proteomic data, with sufficient statistical power to explore complex biological relationships because of this large, richly annotated patient cohort:. This is considered a high profile strength of the project\nThe application of Mendelian Randomization and the Findr pipeline to reconstruct directed protein networks from high-dimensional serum data is well-justified and methodologically sound. And represents and innovative use of causal inference.\nThe authors conduct extensive robustness tests, including AUC/precision-recall analysis, transitive reduction, and overlap with known PPIs, lending feasibility and credibility to their findings.\nAssociations with future MI and HF events, alongside known clinical traits, support the clinical applicability of the identified network regulators (e.g., ITIH3, APOA5) showing translational potential\n\nMajor Weaknesses and Areas for Improvement\nLack of External Validation. The causal protein networks are derived from a single cohort of elderly Icelandic participants. No independent replication is attempted in other populations or platforms. It is entirely possible to validate key subnetworks or associations (e.g., ITIH3, DDX39B, APOA5) in an external proteomic cohort even if limited in number, or more relevant to an multi-ethnic cohort. It is entirely possible as well, that there are public available datasets, such as the UK biobank which could serve as a validation informatic approach.\nInformatic analysis/Indirect Edge Ambiguity and Confounding interpretation. While transitive reduction is applied to remove some indirect edges, a large number of network edges may still reflect indirect regulation or residual confounding. A recommendation might be to incorporate formal mediation models or latent confounder correction to further refine causal edge specificity.\nFunctional Validation Missing.\nThe study is computational in nature; no experimental follow-up or perturbation assays are included to confirm network-derived predictions. Here, a suggestion would be to target a subset of high-priority regulators (e.g., ITIH3, HSPA1A, LRP4) for validation using in vitro or in vivo models, or at minimum, review existing functional studies in greater detail.\nThere is some potential ambiguity in protein-aptamer mapping, particular due to the platform that the authors used. Somascan has some inherent technical limitations. In this study, several aptamers map to the same protein family, sometimes targeting different isoforms or protein domains. Such dual identification could result in artificial interpretation or duplication of network regulators. It is recommended that authors clarify how redundancy in aptamer targeting was handled, and consider collapsing overlapping aptamers to reduce spurious edges.\nPhenotypic Associations and Effect Size Interpretation. While statistically significant, the magnitude of associations between network eigenproteins and traits is not always clinically contextualized. The authors are suggested to maybe report hazard ratios or beta coefficients with interpretation to support translational relevance.\n\nMinor Point\nThe study is missing “Negative Controls”. For instance, it would strengthen the interpretation to report null associations (e.g., subnetworks NOT linked to MI or metabolic traits) to contextualize specificity, this is important context to the manuscript."}, "author_response": "Response to Review\nWe are pleased to submit the revised version of our manuscript, entitled “Circulating causal\nprotein networks linked to future risk of myocardial infarction” (NCOMMS-25-18281-A), for\nconsideration for publication in Nature Communications. We sincerely thank the reviewers for\ntheir thoughtful and constructive feedback. In the accompanying point-by-point response\ndocument, we have carefully addressed each of the comments. These revisions also include\ncorrection of typographical errors, and adjustments to ensure full compliance with Nature\nCommunications formatting guidelines. We believe that these changes have substantially\nstrengthened the manuscript.\nReviewer comments are addressed in blue font in the response document. Additions to the revised\nmanuscript are indicated in italicized text. Line numbers cited in our responses refer to the clean\nversion of the revised manuscript, submitted alongside a tracked-changes version highlighting all\nedits.\nReviewer comments\nReviewer #1:\nIn this manuscript, Bankier, Gudmundsdottir and colleagues present an impressive large-scale\nproteomic study of 5,376 older adults from the AGES-Reykjavik cohort, in which 7,523 serum\nproteins were quantified. The authors reconstruct a global causal network of circulating proteins\nand examine their relationships with cardiometabolic traits, including myocardial infarction (MI)\nand heart failure (HF).\nThe manuscript is clearly written and presents a substantial contribution to the field. The findings\nare likely to inform efforts in drug development and repurposing for cardiometabolic disease. That\nsaid, I believe the study could be further strengthened through several methodological refinements\nand additional contextual discussion, as outlined below.\nWe thank the reviewer for his/her positive remarks.\nMajor comments\n1) Replication. The protein network reconstruction is data-driven and could be sensitive to cohort-\nspecific features. Has the team considered validation in an independent cohort? Confirmation of\nthe network structure or key associations in an external dataset would considerably enhance the\nrobustness and generalisability of the findings.\nResponse: While we previously acknowledged the current scarcity of external datasets suitable for\nour network structure replication, we fully agree with the Reviewer that independent validation\nwould strengthen the conclusion of our study. Reviewer 3 mirrors Reviewer 1’s concern (see\nbelow) and specifically recommends testing replication using the UK Biobank proteomics data.\nWe therefore evaluated the causal protein network architecture in the UK Biobank Olink dataset,\ncomprising 52,543 samples and 2,186 proteins that overlap with those measured in the AGES\nstudy, despite the previously noted modest concordance with our SomaLogic aptamer-based\nmeasurements (including only moderate Spearman correlations and discrepancies in cis-acting\npQTL signals for many analytes)1. Accordingly, while performing this analysis, we have made a\npoint of emphasizing platform heterogeneity as an important caveat for interpreting cross-cohort\n-1-\ncomparability. Based on this analysis, we have added new Supplementary Table S10 and\nSupplementary Fig. S11 (displayed at the end of this document), along with corresponding text in\nthe Results, Discussion, and Methods sections as detailed below:\nLines 297-337 in the Result section:\n“Replication of causal protein network architecture in an independent cohort\n“We performed a validation analysis of the CPN structure using the UK Biobank (UKBB) study2\nas an independent dataset. There are 2,186 assays from the UKBB Olink Explore and the AGES\n7K SOMAscan profiling platforms that target the same proteins according to their annotations\n(Methods). Therefore, we reconstructed CPNs in both AGES and UKBB using only this subset of\nproteins from each platform to ensure that these networks were comparable. As correlations in\nprotein measurements between these two different proteomic platforms have been shown to vary3,\nwe conducted instrument selection separately in the UKBB and AGES studies (Methods), removing\nduplicate aptamers targeting the same protein before analysis.\nWe first reconstructed a subsetted CPN in AGES but restricted it to the 2,186 proteins also\nmeasured in UKBB. Of these protein assays, 1,835 had a valid cis-acting pQTL instrument to be\nused as a potential network regulator. At a 1% FDR threshold we identified 454 proteins with at\nleast one target and 43 with at least 10 targets in the AGES study. Among the 454 proteins with at\nleast one target in the subsetted AGES CPN (FDR = 1%), 416 were also present in the full AGES\nCPN constructed from all measured proteins. Furthermore, 96.3% of the edges identified in the\nsubsetted AGES CPN overlapped with those in the full CPN. We then reconstructed a comparative\nCPN using data from 52,543 UKBB participants with Olink assays targeting the same 2,186\nproteins. We observed a significant overlap between network regulators in the two cohorts,\nincluding the 454 proteins with at least one target in AGES, of which 393 were present with at\nleast one target in UK Biobank (Fisher’s Exact Test; OR = 1.9, P = 4.1×10-6). Further, we\nidentified 629 edges shared between the AGES subsetted and UK Biobank CPNs (FDR = 1%). To\nformalize this comparison, we applied Fisher’s Exact Test, which demonstrated that the overlap\nwas highly significant (OR = 2.9, P = 1.6×10⁻⁹⁹), using all possible edges as the background while\nexcluding self-edges. Notably, many more targets were identified for a given CPN in UKBB than\nin AGES (Supplementary Table S10), likely to reflect the nearly tenfold larger sample size in\nUKBB.\nWe next sought to determine the extent to which protein targets of subnetworks reconstructed in\nAGES and those from UKBB overlapped. Of the 393 shared network regulators from the\ncomparative subnetworks, we found 36 that shared at least 2 targets in the AGES and UKBB\ncomparative CPNs, for which we assessed overlap significance with a hypergeometric test, using\nthe full set of 2,186 proteins as background. Overall, 23 subnetworks had a significant overlap in\ntargets across both platforms (q < 0.05), including three of the top-ranked subnetworks previously\nlinked to ACVD in the full AGES CPN, namely, KLKB1, CFB and CSF3 (Supplementary Table\nS10, Supplementary Fig. S11).\nA recent study compared correlations between Olink and SOMAlogic platforms in the Icelandic\npopulation3. From this study, we obtained correlations for 1,677 of the 2,186 proteins measured\n(mean Pearson coefficient = 0.39). For 19 of the 23 significantly overlapping subnetworks,\ncorrelations were available, with a higher mean of 0.48 (Supplementary Table S10). In some cases,\nnetwork regulators of significantly overlapping CPNs were not correlated across platforms\n(Supplementary Table S10), suggesting that even when different affinity-based technologies bind\n-2-\ndistinct epitopes, the network regulator remains robust and continues to capture the same\nunderlying biology reflected in their downstream targets.”\nLines 478-493 in the Discussion section: “Validating our CPNs is inherently challenging because\ncomplementary external data are scarce, unlike tissue-specific gene regulatory networks that can\nbe supported by gene expression QTLs and external transcription factor-target datasets4. To\naddress this challenge, we drew on existing resources such as the UKBB5,6, although its protein\ncoverage is considerably smaller than in our study, and inter-platform correlations are modest\n(mean r = 0.39), with additional differences arising from the epitope targeting inherent to affinity-\nbased technologies. Nevertheless, the identification of 23 significantly overlapping subnetworks\n(out of 36 comparable networks) across AGES and UKBB highlights that the CPN architecture is\nrobust. We further observed strong cross-cohort replication of global regulators and their\nconnections, with regulators and edges significantly overlapping beyond chance expectations. The\nUKBB’s tenfold larger sample size facilitated the detection of more targets per regulator, but its\nreduced proteomic coverage limited a more comprehensive replication of networks, emphasizing\nthe dual importance of sample size and proteomic breadth for network reconstruction and\ncomparison. Importantly, several of the top-ranked subnetworks, including KLKB1, CFB, and\nCSF3, were consistently replicated across platforms and populations, underscoring their\nbiological relevance beyond platform-specific effects.”\nLines 774-799 in Method section:\n“Validating the network structure across proteomic platforms and cohorts\nThe UKBB is a prospective population study of 502,128 individuals from the UK who have been\nextensively characterized for genomic and phenotypic traits2. Genotype data is available for\n486,593 participants, obtained using either the Applied Biosystems UK BiLEVE Axiom Array by\nAffymetrix or the Applied Biosystems UKBB Axiom Array. Imputation was carried out using the\nHRC reference panel, yielding allelic dosages for ~96 million variants, with markers annotated\nusing the GRCh37 assembly of the human genome5.\nThere are 53,013 (53.9% female, mean age = 56.8) UKBB participants who have undergone\nplasma protein profiling for 2,923 proteins using the antibody affinity based Olink Explore 3072\nprotein extension assay platform6. We processed the Olink assay data using the approach outlined\nby Sun et al.6, from the initial study reporting plasma measurements from UKBB, whereby assays\nwith more than 20% missingness are removed and the remaining missing values are mean imputed\nacross samples. The overlap between processed samples and those with corresponding genotype\ndata yielded 52,543 samples to be taken forward for network reconstruction. The remaining\nprotein assays were then filtered for those that had a corresponding assay in the AGES SOMAscan\n-3-\n7K panel, yielding 2,186 comparative protein assays. These assay measurements were then\ncorrected for age and sex by fitting to a linear model as previously described for the AGES\nproteomic data.\nDue to differences between Somalogic and Olink platforms, instrument selection was carried out\nindependently for the UKBB proteomic data. For each protein, the lead cis-acting pSNP was\nchosen based on the P-value, using pQTL summary statistics from Sun et al.6. Reconstruction of\ncausal protein networks was carried out with Findr and edges selected based on FDR thresholds,\nas has been previously described for the AGES discovery cohort. We tested the significance of\noverlapping targets between comparative subnetworks by hypergeometric testing using the\nhypergeom() function from Scipy Stats7 in Python, while using the intersection of measured\nproteins between platforms as a background. We then corrected for multiple testing using Storey-\nTibshirani procedure for q-value estimation8.”\n2) One-sample Mendelian Randomisation (MR) and risk of bias The use of one-sample MR\nintroduces potential for bias and inflation of the type I error rate (see: Burgess et al., Wellcome\nOpen Res, 2019;4:186). Where feasible, a two-sample MR approach would be preferable, drawing\nupon publicly available GWAS summary statistics for the cardiometabolic traits in question. This\nwould reduce bias, increase statistical power, and facilitate a range of sensitivity analyses (e.g.,\nMR Egger, weighted median) that are vital for assessing the robustness of causal inference.\nResponse: Two-sample Mendelian randomization (MR) was conducted and clearly described in\nboth the Results and Methods sections of the current manuscript. See for instance current Result\nsection: “We performed colocalization and two-sample MR analysis on the network regulators\nlisted in Table 2 to further investigate potential causal relationships with ACVD related outcomes\n(see Methods)”, and Method section: “The causality of selected network regulators was assessed\nusing a two-sample cis-Mendelian randomization (MR) approach…”\nWhile sensitivity analysis was not included in the previous version, it has now been performed and\nis presented in a new Supplementary Table S14, with corresponding methodological details added\nto the Results and Methods section. See also Reviewer 2 comments on MR below, related to the\nuse of different LD clumping thresholds.\nLines 356-360 in Results section: “Sensitivity analyses restricted to exposures with at least three\ninstruments (Methods), indicated that results were largely robust across approaches\n(Supplementary Table S14), although significant instrument heterogeneity was observed for\nAPOA5 in relation to MetS. Leave-one-out analysis, however, suggested that the overall causal\neffect was not driven by any single variant (data not shown).”\nLines 827-832 in Method section: “To ensure the robustness of our findings, we performed\nsensitivity analyses of the MR estimates using MR-Egger9 and weighted median estimation10. For\nproteins with more than three instruments available, instrument heterogeneity was assessed using\nCochrans´s Q11. Horizontal pleiotropy was assessed using MR-Egger9. The generalized weighted\nleast squares (GWLS) was performed as previously described12, while other MR analyses were\nperformed using the TwoSampleMR13 R package.”\n3) Colocalization methodology. The authors may wish to consider the application of more\nadvanced colocalization approaches capable of modelling multiple causal variants, such as SuSiE\nor PwCoCo (e.g. PMID: 35452592). These methods are better suited to complex loci and may\nyield more reliable insights regarding shared genetic aetiology between proteins and traits.\nResponse: In response to the Reviewer’s suggestion, we applied the SuSiE approach. We note,\nhowever, that the MR results were largely consistent with those obtained using our earlier\ncolocalization method (see also Reviewer 2’s related comment below). Compared with the original\ncolocalization analysis, many loci did not pass the SuSiE fine-mapping step, nearly two-thirds in\ntotal, and thus could not be tested for colocalization. For this reason, we present both the original\n-4-\nanalysis and the SuSiE-based results, now shown in the updated Table S13 (formerly Table S12).\nThe description in the Methods section has been revised.\nLines 801-819 in Method section:\n“Colocalization and Mendelian randomization analysis\nWe used colocalization analysis to identify cis-acting pQTLs that shared a common signal with six\ndifferent phenotypic traits using the R package coloc (version 5.2.3)14 and publicly available\nGWAS summary statistics from individuals with European ancestry (Supplementary Table S10).\nWe first identified all proteins that had a cis-acting pQTL (FDR < 5%) that shared at least 1\ncommon cis-pSNP with the GWAS trait (P < 510-8) and extracted the summary statistics for all\nSNPs within ±150 Kb of the transcription start site for the cognate encoding gene. We then\nestimated the probability of a shared causal variant (PP.H4) using approximate Bayes Factor\ncolocalization5 via the coloc.abf (coloc) function using a threshold of PP > 0.9 and >0.5.\nAdditionally, we tested colocalization between overlapping signals of cardiometabolic GWAS\ntraits and serum protein cis-acting pQTL using the Bayes Factor colocalization5 with Sum of\nSingle Effects (SuSiE) framework for fine mapping of genomic loci6. The SNPs within the cis\nwindow were fine mapped for both the GWAS and pQTL signal using the susie_rss() function from\nthe SusieR R package (version 0.12.35), while limiting the number of possible independent signals\n(L = 5) and the maximum number of iterations (n = 5,000). Only loci where convergence was\nobserved for both signals were taken forward for colocalisation. Colocalization was performed\nusing the coloc.susie() function from the Coloc R package (version 5.2.3), which performs Bayes\nFactor colocalization. We selected for colocalized signals (H4) using the posterior probability\nthresholds > 0.9 and >0.5.”\n4) Inclusion of positive controls. It may be instructive to include a small number of well-\nestablished biological pathways or protein complexes—such as those involved in transcriptional\nregulation, the cell cycle, or intracellular transport—as positive controls. Demonstrating that the\nnetwork approach successfully identifies such relationships would lend confidence to the method's\nvalidity.\nResponse: We appreciate the reviewer’s suggestion to include well-established biological\npathways or protein complexes as positive controls to validate our network approach. In the current\nversion of the manuscript, we have indeed identified edges in the causal serum protein network\nthat overlap with direct physical protein–protein interactions (PPIs) across multiple confidence\nthresholds using the HIPPIE database (z-score = 24.3, P < 0.001). These overlaps serve as a form\nof external validation, demonstrating that our network captures biologically plausible interactions.\n5) Ancestry and generalizability As the study sample is limited to individuals of European descent,\nthe authors should more clearly acknowledge this as a limitation. While expansion to more diverse\npopulations may not yet be feasible, explicit recognition of this issue will help orient future work\ntowards greater inclusivity and global relevance.\nResponse: We thank the reviewer for highlighting this point. In response, we have included a new\nsentence in the Discussion section, where we had previously addressed other limitations of our\nstudy.\n-5-\nLines 497-500 in the Discussion section: “First, all AGES participants are of European descent\n(i.e., white/Caucasian), which may limit the generalizability of our findings, as protein predictors\nand clinical indicators of heart disease can differ across populations with diverse genetic and\nenvironmental backgrounds.”\nMinor comment\nLine 540: It would be helpful to define the symbol “E” for the benefit of readers less familiar with\nnetwork analysis terminology.\nResponse: We thank the reviewer and have revised the text in the Methods section (Line 648) to\nclarify the meaning of E. Specifically; E represents the cis-acting pQTL instrument in the causal\npathway E → A → B.\nIn sum, this is a high-quality manuscript addressing an important topic, with clear translational\nrelevance. The above suggestions are intended to strengthen the rigour and transparency of the\nwork. I look forward to seeing the revised version.\n-6-\nReviewer #2:\nThe manuscript presents a sophisticated causal network inference framework leveraging 7,523\nserum proteins measured in 5,376 elderly individuals from the AGES-Reykjavik Study. The\nrationale for a system-level approach and the use of pQTLs for causal protein network inference\nis well-articulated, and the methodology is technically sound but highly complex. Below are\nspecific concerns and suggestions for improvement.\nWe thank the reviewer for his/her positive remarks.\n1. Abstract\nThe statement \"we used cis-acting serum protein quantitative trait loci (pQTLs) as instrumental\nvariables, while accounting for potential unobserved confounding factors\" requires clarification.\nIf this refers to Mendelian randomization (MR), the phrasing is misleading—MR mitigates\nunobserved confounding via genetic IVs but does not explicitly \"account\" for it. Please revise\nto reflect MR’s assumptions accurately.\nResponse: We thank the Reviewer for this helpful comment. We agree that the original phrasing\ncould be misleading in the context of MR, as our statement referred specifically to the methodology\nused for network reconstruction within the Findr framework. To avoid confusion, we have revised\nthe Abstract accordingly. The sentence now reads:\n“To reconstruct a global causal network of circulating proteins, we used cis-acting serum protein\nquantitative trait loci (pQTLs) as instrumental variables within a causal inference framework,\nwhich includes tests designed to mitigate the influence of hidden confounders.”\n2. Results (Line 129–130)\n\"The follow-up period for newly diagnosed MI patients extended up to 12 years from baseline,\nwith recurrent cases excluded from the incident group analysis. Person-years of follow-up were\ncalculated from the first AGES visit until the date of diagnosis, death, or the end of the follow-up\nperiod...\" Incident vs. Recurrent Cases: If follow-up was censored at the first MI diagnosis (i.e.,\nincident cases only), recurrent MI cases should not exist in this cohort by definition.\nResponse: We thank the Reviewer for his/her observation and would like to clarify that follow-up\nwas indeed censored at the time of the first MI diagnosis, ensuring that only first (incident) MI\nevents were included. By design, recurrent MI cases were not part of the incident cohort, and we\nhave therefore removed any reference to recurrent MI in the corresponding sentence (Lines 130-\n133).\n3. Discussion\nThe clinical implications of prioritized proteins (Table 2) are underdeveloped. Given the focus on\nASCVD early detection and prevention, readers would benefit from:\n- A brief summary of how many proteins have licensed/investigational drugs targeting them.\n- Potential therapeutic pathways suggested by the causal network.\n-7-\nResponse: We fully agree with the Reviewer and, in response, have generated a new\nSupplementary Table S16 reporting whether each top-ranked network regulator is druggable or\ntargeted by existing compounds using public databases, and have added corresponding description\nto the Supplementary Text (new section called Therapeutic pathways and clinical implications)\nand Results and Discussion sections to highlight potential therapeutic pathways and clinical\nimplications. In response to related comments from Reviewer 3, we used information from\nSupplementary Table S16 and another newly added Supplementary Table S15 (see related\ncomment by Reviewer 3 below) to generate a new Table 3 in the main text. This table expands on\nadditional links between the regulators and ACVD and integrates directional consistency across\ncausal, observational, and therapeutic evidence for each regulator.\nLines 404-426 in Results: “Therapeutic pathways and clinical implications: Over half of the top-\nranked network regulators are potentially druggable with small molecules or biologics, although\nonly a few have compounds currently available for either related or unrelated disease indications\n(Supplementary Table S16, Supplementary Text). Specifically, 13 of the 25 regulators are\ndruggable, with 6 already targeted by licensed or investigational drugs. Building on this, the CPN\nanalysis has identified both interconnected networks and therapeutic pathways with potential for\nearly ACVD detection and prevention: 1) Complement-mediated inflammation: The network\nregulators C2 and CFB are central to the alternative complement pathway, driving vascular\ninflammation and atherosclerosis15. CFB is targeted by the antisense oligonucleotide Sefaxersen,\nin development for IgA nephropathy (Supplementary Table S16). Targeting this pathway could\nprevent early inflammatory progression of atherosclerotic lesions. 2) Cellular stress response and\ncyto-protection: KEAP1 regulates the Nrf2 oxidative stress pathway16, while HSPA1A and\nHSPA1B control protein folding under stress17. All three have licensed drugs targeting them\n(Supplementary Table S16). Modulating this pathway may protect against oxidative damage and\nendothelial dysfunction. 3) Hemostasis: The key target KLKB1, a regulator of plasma kallikrein\nin the intrinsic coagulation pathway18, is targeted by the approved drug Berotralstat against\nangioedema (Supplementary Table S16). Our integrative analysis, however, implies that\nactivating KLKB1 may offer a strategy to prevent thrombotic complications in ACVD (Table 3).\n4) Lipid metabolism: The network regulators APOA5 and AFM control triglyceride-rich\nlipoprotein metabolism and oxidative stress19-21, are bio-druggable (Supplementary Table S16),\nand could be targeted to prevent dyslipidemia-driven initiation of atherosclerosis. This network\nframework highlights causal pathways for ACVD prevention, with concordant data strengthening\ntarget confidence and inconsistencies pointing to opportunities for novel therapeutic discovery.”\nLines 519-525 in Discussion: “Indeed, among 25 network regulators, 13 are druggable, with 6\nalready targeted by licensed or investigational drugs. Accordingly, the CPN reveals\ninterconnected therapeutic pathways, including complement-mediated inflammation, cellular\nstress response, hemostatic balance, and lipid metabolism, supporting both single- and multi-\npathway targeting. This framework provides a roadmap for next-generation ASCVD prevention,\nwhere concordant evidence reinforces target confidence and discrepancies highlight opportunities\nfor novel therapeutic discovery.”\n4. Methods\nClarify definition of NAFLD/MASLD criteria.\n-8-\nResponse: We have revised the Methods section to clarify the diagnostic criteria for\nNAFLD/MASLD.\nLines 583-585 in Method section: “Hepatic steatosis was assessed by computed tomography (CT),\nserving as a non-invasive proxy for non-alcoholic fatty liver disease (NAFLD), as previously\ndescribed22”\nGenomic Window Inconsistencies: Varying window sizes are used across analyses (300 kb for cis-\nacting pQTL detection, 150 kb for variance modeling/colocalization, 500 kb for MR IVs). Justify\nthese choices.\nResponse: We thank the reviewer for highlighting this point. To clarify, the 300 kb window used\nfor cis-acting pQTL detection is equivalent to the ±150 kb window applied in the variance\nmodeling and colocalization analyses, i.e., both define the same 300 kb genomic span centered on\nthe gene. For the Mendelian randomization analyses, a broader window was used to define cis-\ninstruments, based on individual disease studies conducted independently of the network\nreconstruction and the cis-window applied therein. For consistency we have repeated the MR\nanalysis using a 300 kb (±150 kb) window and updated the text accordingly (see track changes\nversion of the manuscript). Please see our response to the same reviewer's comment on LD\nclumping thresholds below. Updated results are presented in Table S13 (previously Table S12).\nColocalization Assumptions: The analysis assumes a single causal variant. For robustness,\nconsider methods accommodating multiple causal variants (e.g., PWCoCo, SharePro) for\nprioritized proteins.\nResponse: Please see our response to a similar comment from Reviewer 1 (comment 3), who\nsuggested using SuSiE. Following the recommendations of both reviewers, we revisited the\ncolocalization analysis with an additional approach. The updated results, alongside the previous\nanalyses, are presented in Table S13 (formerly Table S12), and the Methods section has been\nrevised accordingly.\nLines 801-819 in Method section:\n“Colocalization and Mendelian randomization analysis\nWe used colocalization analysis to identify cis-acting pQTLs that shared a common signal with six\ndifferent phenotypic traits using the R package coloc (version 5.2.3)14 and publicly available\nGWAS summary statistics from individuals with European ancestry (Supplementary Table S10).\nWe first identified all proteins that had a cis-acting pQTL (FDR < 5%) that shared at least 1\ncommon cis-pSNP with the GWAS trait (P < 510-8) and extracted the summary statistics for all\nSNPs within ±150 Kb of the transcription start site for the cognate encoding gene. We then\nestimated the probability of a shared causal variant (PP.H4) using approximate Bayes Factor\ncolocalization5 via the coloc.abf (coloc) function using a threshold of PP > 0.9 and >0.5.\nAdditionally, we tested colocalization between overlapping signals of cardiometabolic GWAS\ntraits and serum protein cis-acting pQTL using the Bayes Factor colocalization5 with Sum of\nSingle Effects (SuSiE) framework for fine mapping of genomic loci6. The SNPs within the cis\nwindow were fine mapped for both the GWAS and pQTL signal using the susie_rss() function from\nthe SusieR R package (version 0.12.35), while limiting the number of possible independent signals\n-9-\n(L = 5) and the maximum number of iterations (n = 5,000). Only loci where convergence was\nobserved for both signals were taken forward for colocalisation. Colocalization was performed\nusing the coloc.susie() function from the Coloc R package (version 5.2.3), which performs Bayes\nFactor colocalization. We selected for colocalized signals (H4) using the posterior probability\nthresholds > 0.9 and >0.5.”\nThe LD clumping threshold (r² ≥ 0.2) is lenient and may retain non-independent cis-acting pQTLs.\nAddress whether stricter thresholds (e.g., r² < 0.1) is not feasible here.\nResponse: In addition to narrowing the genomic window in the MR to ensure consistency across\nanalyses (see response above), we applied a more stringent LD clumping thresholds of r² < 0.1 in\naddition to the original threshold r² < 0.2. We note however that the GWLS estimator used for the\nMR analysis accounts for correlation between instruments23. Overall, the main findings remain\nbroadly consistent, with some associations reaching higher levels of significance. The updated\nresults are now provided in Table S13 (previously Table S12).\n-10-\nReviewer #3:\nSummary\nThis manuscript by Emilsson and colleagues presents a comprehensive systems biology approach\nto understand the molecular outcomes underlying myocardial infarction (MI) and associated\ncardiometabolic conditions. By analyzing serum proteomic data through Somascan with samples\nfrom the AGES-Reykjavik study and employing a causal inference framework using cis-acting\npQTLs, the authors reconstruct 185 high-confidence causal protein subnetworks (CPNs). These\nsubnetworks were then associated with incident MI, metabolic syndrome, heart failure, and related\ntraits. The authors attempted to validate their network using internal metrics and reference protein-\nprotein interaction databases and lastly performed colocalization and Mendelian Randomization\n(MR) analyses in order to strengthen their claims.\nI think that this study represents an important advancement in causal proteomics and has the\npotential to uncover novel biomarkers and therapeutic targets for cardiovascular disease. The\ngreater than 5,000 samples afforded by the AGES cohort is compelling. However, some conceptual\nand methodological weaknesses should be addressed to enhance the study’s impact and\ninterpretability. Importantly, the analytical framework is strong, but the study would benefit from\nindependent validation, greater biological resolution, and experimental follow-up.\nMajor Strengths\nThe AGES-Reykjavik cohort provides deep phenotyping and proteomic data, with sufficient\nstatistical power to explore complex biological relationships because of this large, richly annotated\npatient cohort: This is considered a high profile strength of the project\nThe application of Mendelian Randomization and the Findr pipeline to reconstruct directed protein\nnetworks from high-dimensional serum data is well-justified and methodologically sound. And\nrepresents and innovative use of causal inference.\nThe authors conduct extensive robustness tests, including AUC/precision-recall analysis,\ntransitive reduction, and overlap with known PPIs, lending feasibility and credibility to their findings.\nAssociations with future MI and HF events, alongside known clinical traits, support the clinical\napplicability of the identified network regulators (e.g., ITIH3, APOA5) showing translational\npotential.\nWe appreciate the reviewer’s positive feedback.\nMajor Weaknesses and Areas for Improvement\nLack of External Validation. The causal protein networks are derived from a single cohort of\nelderly Icelandic participants. No independent replication is attempted in other populations or\nplatforms. It is entirely possible to validate key subnetworks or associations (e.g., ITIH3,\nDDX39B, APOA5) in an external proteomic cohort even if limited in number, or more relevant to\nan multi-ethnic cohort. It is entirely possible as well, that there are public available datasets, such\nas the UK biobank which could serve as a validation informatic approach.\n-11-\nResponse: While we previously acknowledged the current scarcity of external datasets suitable for\nour network -structure replication, we fully agree with the Reviewer that independent validation\nwould strengthen the conclusion of our study. We have therefore assessed the causal protein\nnetwork architecture in the UK Biobank Olink dataset, comprising 52,543 samples and 2,186\nproteins, despite previously noted modest cross platform concordance with our -SomaLogic\naptamer-based protein measurements (including only moderate Spearman correlations and\ndiscrepancies in cis-acting pQTL signals for many analytes)1. Accordingly, while performing this\nanalysis, we have made a point of emphasizing platform heterogeneity as an important caveat for\ninterpreting cross-cohort -comparability. Based on this analysis, we have added new Supplementary\nTable S10 and Supplementary Fig. S11 (displayed at the end of the document), along with\ncorresponding text in the Results, Discussion, and Methods sections as detailed below:\nLines 297-337 in the Result section:\n“Replication of causal protein network architecture in an independent cohort\n“We performed a validation analysis of the CPN structure using the UK Biobank (UKBB) study2\nas an independent dataset. There are 2,186 assays from the UKBB Olink Explore and the AGES\n7K SOMAscan profiling platforms that target the same proteins according to their annotations\n(Methods). Therefore, we reconstructed CPNs in both AGES and UKBB using only this subset of\nproteins from each platform to ensure that these networks were comparable. As correlations in\nprotein measurements between these two different proteomic platforms have been shown to vary3,\nwe conducted instrument selection separately in the UKBB and AGES studies (Methods), removing\nduplicate aptamers targeting the same protein before analysis.\nWe first reconstructed a subsetted CPN in AGES but restricted it to the 2,186 proteins also\nmeasured in UKBB. Of these protein assays, 1,835 had a valid cis-acting pQTL instrument to be\nused as a potential network regulator. At a 1% FDR threshold we identified 454 proteins with at\nleast one target and 43 with at least 10 targets in the AGES study. Among the 454 proteins with at\nleast one target in the subsetted AGES CPN (FDR = 1%), 416 were also present in the full AGES\nCPN constructed from all measured proteins. Furthermore, 96.3% of the edges identified in the\nsubsetted AGES CPN overlapped with those in the full CPN. We then reconstructed a comparative\nCPN using data from 52,543 UKBB participants with Olink assays targeting the same 2,186\nproteins. We observed a significant overlap between network regulators in the two cohorts,\nincluding the 454 proteins with at least one target in AGES, of which 393 were present with at\nleast one target in UK Biobank (Fisher’s Exact Test; OR = 1.9, P = 4.1×10-6). Further, we\nidentified 629 edges shared between the AGES subsetted and UK Biobank CPNs (FDR = 1%). To\nformalize this comparison, we applied Fisher’s Exact Test, which demonstrated that the overlap\nwas highly significant (OR = 2.9, P = 1.6×10⁻⁹⁹), using all possible edges as the background while\nexcluding self-edges. Notably, many more targets were identified for a given CPN in UKBB than\nin AGES (Supplementary Table S10), likely to reflect the nearly tenfold larger sample size in\nUKBB.\nWe next sought to determine the extent to which protein targets of subnetworks reconstructed in\nAGES and those from UKBB overlapped. Of the 393 shared network regulators from the\ncomparative subnetworks, we found 36 that shared at least 2 targets in the AGES and UKBB\ncomparative CPNs, for which we assessed overlap significance with a hypergeometric test, using\nthe full set of 2,186 proteins as background. Overall, 23 subnetworks had a significant overlap in\ntargets across both platforms (q < 0.05), including three of the top-ranked subnetworks previously\n-12-\nlinked to ACVD in the full AGES CPN, namely, KLKB1, CFB and CSF3 (Supplementary Table\nS10, Supplementary Fig. S11).\nA recent study compared correlations between Olink and SOMAlogic platforms in the Icelandic\npopulation3. From this study, we obtained correlations for 1,677 of the 2,186 proteins measured\n(mean Pearson coefficient = 0.39). For 19 of the 23 significantly overlapping subnetworks,\ncorrelations were available, with a higher mean of 0.48 (Supplementary Table S10). In some cases,\nnetwork regulators of significantly overlapping CPNs were not correlated across platforms\n(Supplementary Table S10), suggesting that even when different affinity-based technologies bind\ndistinct epitopes, the network regulator remains robust and continues to capture the same\nunderlying biology reflected in their downstream targets.”\nLines 478-493 in the Discussion section: “Validating our CPNs is inherently challenging because\ncomplementary external data are scarce, unlike tissue-specific gene regulatory networks that can\nbe supported by gene expression QTLs and external transcription factor-target datasets4. To\naddress this challenge, we drew on existing resources such as the UKBB5,6, although its protein\ncoverage is considerably smaller than in our study, and inter-platform correlations are modest\n(mean r = 0.39), with additional differences arising from the epitope targeting inherent to affinity-\nbased technologies. Nevertheless, the identification of 23 significantly overlapping subnetworks\n(out of 36 comparable networks) across AGES and UKBB highlights that the CPN architecture is\nrobust. We further observed strong cross-cohort replication of global regulators and their\nconnections, with regulators and edges significantly overlapping beyond chance expectations. The\nUKBB’s tenfold larger sample size facilitated the detection of more targets per regulator, but its\nreduced proteomic coverage limited a more comprehensive replication of networks, emphasizing\nthe dual importance of sample size and proteomic breadth for network reconstruction and\ncomparison. Importantly, several of the top-ranked subnetworks, including KLKB1, CFB, and\nCSF3, were consistently replicated across platforms and populations, underscoring their\nbiological relevance beyond platform-specific effects.”\nLines 774-799 in Method section:\n“Validating the network structure across proteomic platforms and cohorts\nThe UKBB is a prospective population study of 502,128 individuals from the UK who have been\nextensively characterized for genomic and phenotypic traits2. Genotype data is available for\n486,593 participants, obtained using either the Applied Biosystems UK BiLEVE Axiom Array by\nAffymetrix or the Applied Biosystems UKBB Axiom Array. Imputation was carried out using the\nHRC reference panel, yielding allelic dosages for ~96 million variants, with markers annotated\nusing the GRCh37 assembly of the human genome5.\nThere are 53,013 (53.9% female, mean age = 56.8) UKBB participants who have undergone\nplasma protein profiling for 2,923 proteins using the antibody affinity based Olink Explore 3072\nprotein extension assay platform6. We processed the Olink assay data using the approach outlined\nby Sun et al.6, from the initial study reporting plasma measurements from UKBB, whereby assays\nwith more than 20% missingness are removed and the remaining missing values are mean imputed\nacross samples. The overlap between processed samples and those with corresponding genotype\ndata yielded 52,543 samples to be taken forward for network reconstruction. The remaining\nprotein assays were then filtered for those that had a corresponding assay in the AGES SOMAscan\n-13-\n7K panel, yielding 2,186 comparative protein assays. These assay measurements were then\ncorrected for age and sex by fitting to a linear model as previously described for the AGES\nproteomic data.\nDue to differences between Somalogic and Olink platforms, instrument selection was carried out\nindependently for the UKBB proteomic data. For each protein, the lead cis-acting pSNP was\nchosen based on the P-value, using pQTL summary statistics from Sun et al.6. Reconstruction of\ncausal protein networks was carried out with Findr and edges selected based on FDR thresholds,\nas has been previously described for the AGES discovery cohort. We tested the significance of\noverlapping targets between comparative subnetworks by hypergeometric testing using the\nhypergeom() function from Scipy Stats7 in Python, while using the intersection of measured\nproteins between platforms as a background. We then corrected for multiple testing using Storey-\nTibshirani procedure for q-value estimation8.”\nInformatic analysis/Indirect Edge Ambiguity and Confounding interpretation. While transitive\nreduction is applied to remove some indirect edges, a large number of network edges may still\nreflect indirect regulation or residual confounding. A recommendation might be to incorporate\nformal mediation models or latent confounder correction to further refine causal edge specificity.\nResponse: Mediation analysis was already incorporated within the Findr pipeline, and we now\npresent the results at a stringent 1% FDR threshold in a new Supplementary Table S10, which also\nincludes the replication results in the UKBB dataset. However, we note that a recognized limitation\nof mediation analysis is its susceptibility to confounding, which can lead to a high rate of false\nnegatives, a challenge that is particularly pronounced in biological networks where coregulation\nby shared upstream factors frequently occurs.\nFunctional Validation Missing.\nThe study is computational in nature; no experimental follow-up or perturbation assays are\nincluded to confirm network-derived predictions. Here, a suggestion would be to target a subset of\nhigh-priority regulators (e.g., ITIH3, HSPA1A, LRP4) for validation using in vitro or in vivo\nmodels, or at minimum, review existing functional studies in greater detail.\nResponse: Please see our responses to Reviewer 1 (comment 4), Reviewer 2 (comment 3), and the\ncurrent Reviewer’s comment for related discussion below (related to effect size/direction). In brief,\nwe have expanded on the connections among the top-ranked network regulators by integrating\nboth prior and new knowledge, along with additional interpretation, including effect direction,\nlinks to survival outcomes, established biology, and potential druggability (including licensed or\ninvestigational drugs). As detailed elsewhere, this has resulted in expanded text in the main\nmanuscript and Supplementary Text, two new Supplementary tables, and a new Table 3 in the\nmain text. Please see our detailed response below. However, given the scope, objectives, and\ndesign of the present study, we believe that experimental validation in vitro or in vivo is more\nappropriately addressed in a separate follow-up investigation and is beyond the scope of this work.\nThere is some potential ambiguity in protein-aptamer mapping, particular due to the platform that\nthe authors used. Somascan has some inherent technical limitations. In this study, several aptamers\nmap to the same protein family, sometimes targeting different isoforms or protein domains. Such dual identification could result in artificial interpretation or duplication of network regulators. It is\n-14-\nrecommended that authors clarify how redundancy in aptamer targeting was handled and consider\ncollapsing overlapping aptamers to reduce spurious edges.\nResponse: This was indeed addressed. For example, on former page 7 (current Lines 156-160) in\nthe Results section: “For A-proteins with multiple aptamers, we selected those with the largest\nnumber of targets, and further refinement of linkage disequilibrium (LD) among cis-acting\ninstruments resulted in the final CPN comprising 185 A-proteins (referred to here as network\nregulators), their corresponding subnetworks, and a total of 31,358 edges (Supplementary Table\nS1, Fig. 2A)”. In current Lines 668-670 of the Methods section, we clarified this point as well: “In\ninstances where there were multiple aptamers targeting the same protein, we selected the A\nproteins with the largest number of targets as the representative aptamer for this protein.”\nPhenotypic Associations and Effect Size Interpretation. While statistically significant, the\nmagnitude of associations between network eigenproteins and traits is not always clinically\ncontextualized. The authors are suggested to maybe report hazard ratios or beta coefficients with\ninterpretation to support translational relevance.\nResponse: This is a great point which we have addressed in the context of related comments and\nrevisions described above. Specifically, please refer to our earlier response to your comment, as\nwell as Reviewer 2’s comment 3. We agree that interpreting both effect size and direction is critical\nfor translational relevance and acknowledge that this was not fully addressed in the previous\nversion of the study. It is important to note that eigenproteins represent aggregate signals derived\nfrom multiple proteins within a module. Consequently, their associations with traits often exhibit\nweak or modest effect sizes, reflecting average variation across the module rather than strong\neffects from individual proteins, such as for specific network regulators. Therefore, effect sizes\ncannot be directly compared between network regulators and their corresponding eigenproteins.\nNonetheless, we consider the effect direction of network regulators, together with their consistency\nwith other evidence including causal inference and associations with clinical outcomes and\nsurvival to be the most critical and note that this was overlooked in the previous version of our\npaper. We have now addressed this in the revised manuscript through the new Supplementary\nTable S15 and new Table 3 (supported by data from Supplementary Tables S15 and S16) in main\ntext, as well as updates to the Results section and Supplementary text. For example, the Results\nsection titled “Additional links between top-ranked network regulators and ACVD” has been\nrewritten to reflect these updates, as shown below. The new Table 3 is shown at the end of this\ndocument.\nLines 339-403 in Result section:\n“Additional links between top-ranked network regulators and ACVD\nWe investigated several complementary lines of evidence to support the role of the top-ranked\nnetwork regulators in the pathophysiology of MI and related cardiometabolic traits. Network\nregulators and targets from the top-ranked subnetworks showed enrichment in multiple pathways\npreviously linked to ACVD, including a shared enrichment in the cellular response to heat shock\n(Supplementary Table S11, Supplementary Figs. S12–S13, Supplementary Text). STRING-based\nanalysis revealed significantly (P-value = 0.047) enriched functional and physical interactions\namong the top network regulators (Supplementary Fig. S14), including previously found (Fig. 6C)\nand novel connections (e.g., ITIH3 with KLKB1, APOA5, and AFM; Supplementary Fig. S14).\n-15-\nCausal inference analysis of top-ranked regulators: Two-sample MR and colocalization analyses\nhighlighted several top-ranked regulators including APOA5, DDX39B, HSPA1A, C11orf49 and\nLRP4, with significant support for causal associations with ACVD-related traits (FDR < 0.05;\nSupplementary Tables S12–S13). Here, APOA5 demonstrated a strong protective effect against\nMI and MetS, consistent with findings that APOA5 knockout mice display a fourfold elevation in\nplasma triglyceride levels19, a significant risk factor for ACVD24. Further, genetically determined\nlevels of HSPA1A showed an association with MetS, indicating a causal role, whereas DDX39B\nwas linked to both MetS and T2D. Notably, the HSPA1A gene has also been causally implicated\nin T2D and its microvascular complications in prior MR and colocalization analyses25. Sensitivity\nanalyses restricted to exposures with at least three instruments (Methods), indicated that results\nwere largely robust across approaches (Supplementary Table S14), although significant\ninstrument heterogeneity was observed for APOA5 in relation to MetS. Leave-one-out analysis,\nhowever, suggested that the overall causal effect was not driven by any single variant (data not\nshown).\nAdditional context from current and prior evidence: Beyond the analyses of the present study,\nseveral top regulators, including ITIH3, HSPA1B, KEAP1, HSPA1A, C2, CSF3, FABP3, KLKB1,\nAFM, APOA5, PTPN11, CFB, and COL28A1, have previously been linked to ACVD or related\ncardiometabolic traits (Supplementary Text). For example, the top-ranked regulator ITIH3 has a\ngain-of-function genetic association with increased risk of MI and is expressed in vascular smooth\nmuscle cells and macrophages within atherosclerotic plaques26. Notably, ITIH3 is positioned at\nthe base of the data-driven CPN (Fig. 3), indicating a potential role in numerous regulatory\npathways. Another example is KLKB1, which we recently demonstrated through MR analysis to\nbe causally protective against calcific aortic valve disease27, a leading cause of heart failure via\naortic stenosis28. To further characterize the top-ranked protein network regulators, we assessed\ntheir effect directions on ACVD-related traits and survival outcomes (Supplementary Table S15)\nand evaluated if they are druggable or targeted by available compounds using public databases\n(Supplementary Table S16). Table 3 summarizes the accumulated findings, offering an integrated\noverview of directional consistency across causal, observational, and therapeutic evidence for\neach regulator.\nMany of the top-ranked regulators showed consistent directionality across multiple lines of\nevidence and exhibit positive and often causal effects on ACVD traits (Table 3), including ITIH3,\nHSPA1A, DDX39B, FABP3, PTPN11, CFB, and COL28A1, whereas APOA5 and KLKB1 show\nprotective effects (Table 3, Supplementary Text). In contrast, network regulators such as KEAP1,\nC2, and CSF3 show inconsistent effect directions between the present findings and previously\nreported associations (Table 3). KEAP1 regulates the NRF2 pathway, essential for oxidative\nbalance and cardiovascular protection (Supplementary Text). Cardiomyocyte-specific Keap1\nknockout upregulates NRF2 and prevents induced cardiac dysfunction16, contrasting with the\nprotective association of circulating KEAP1 in our study (Table 3). This discrepancy may reflect\ndifferences between tissue-specific and circulating biology, where serum KEAP1 could indicate\nsystemic processes independent of NRF2 inhibition or a feedback marker of activated antioxidant\nresponses rather than a direct protective factor. A human study found C2 deficiency to be\nassociated with increased atherosclerosis risk29, suggesting a protective role for C2 (see\nSupplementary Text). In contrast, our finding that higher circulating C2 predicts increased\nincident MI risk (Table 3), may reflect differences between lifelong deficiency and elevated plasma\nlevels, with the latter potentially indicating complement activation or inflammation. Alternatively,\nelevated C2 in the present study could represent a compensatory upregulation in response to\n-16-\nsubclinical vascular injury. Further, cardiac myocytes express CSF3 under both normal and stress\nconditions, particularly following ischemic injury, where it supports cardiomyocyte survival after\na heart attack30. Notably, the positive association between circulating CSF3 levels, incident MI,\nand reduced survival (Table 3), may reflect a compensatory or adaptive response to early cardiac\ninjury or stress during subclinical disease (Supplementary Text). Finally, given the conflicting\nfindings for C11orf49 and LRP in our study and the inconsistent evidence for AFM in previous\nresearch (Table 3), the direction of their effects on ACVD etiology remains inconclusive. In\nsummary, 68% (17 of 25) of the top-ranked network regulators, this study, along with others,\nreveals additional specific connections to ACVD-related traits. These findings strengthen evidence\nfor the causal role of many of the network regulators and their targets in ACVD development and\nmay guide therapeutic strategies to modulate them.”\nLines 514-519 in Discussion section: “Our integrative analyses support a key mechanistic role for\nmany top-ranked network regulators in ACVD, with mostly consistent associations observed\nacross complementary evidence. Some regulators, including KEAP1, C2, and CSF3, showed\ndirectionally inconsistent associations, likely reflecting differences between circulating and tissue-\nspecific biology or compensatory responses. Overall, these findings reinforce the causal relevance\nof key network regulators and their targets, highlighting potential avenues for therapeutic\nintervention.”\nMinor Point\nThe study is missing “Negative Controls”. For instance, it would strengthen the interpretation to\nreport null associations (e.g., subnetworks NOT linked to MI or metabolic traits) to contextualize\nspecificity, this is important context to the manuscript.\nResponse: This is an important point, particularly given that the network reconstruction process is\ndisease-agnostic. Indeed, many of the 185 subnetworks including 62 network regulators and 121\neigenproteins, exhibited only one or no significant associations with ACVD-related traits\n(Supplementary Tables S5–S6). For instance, VTN, the regulator of the largest network\ncomprising 3,632 targets, showed no significant association with any trait (Supplementary Table\nS5). Accordingly, we have added the following sentence to the Results section:\nLines 252-255 in Results section: “We note that a substantial number of the 185 subnetworks,\nincluding 62 network regulators and 121 eigenproteins, showed only one or no significant\nassociations with ACVD-related traits (Supplementary Tables S5–S6).”\n-17-\nReferences\n1. Rooney, M.R., et al. Plasma proteomic comparisons change as coverage expands for SomaLogic\nand Olink. medRxiv (2024).\n2. Sudlow, C., et al. UK biobank: an open access resource for identifying the causes of a wide range\nof complex diseases of middle and old age. PLoS Med 12, e1001779 (2015).\n3. Eldjarn, G.H., et al. Author Correction: Large-scale plasma proteomics comparisons through\ngenetics and disease associations. Nature 630, E3 (2024).\n4. Gerstein, M.B., et al. Architecture of the human regulatory network derived from ENCODE data.\nNature 489, 91-100 (2012).\n5. Bycroft, C., et al. The UK Biobank resource with deep phenotyping and genomic data. Nature 562,\n203-209 (2018).\n6. Sun, B.B., et al. Plasma proteomic associations with genetics and health in the UK Biobank. Nature\n622, 329-338 (2023).\n7. Virtanen, P., et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat\nMethods 17, 261-272 (2020).\n8. Storey, J.D. & Tibshirani, R. Statistical significance for genomewide studies. Proc Natl Acad Sci\nU S A 100, 9440-9445 (2003).\n9. Burgess, S. & Thompson, S.G. Interpreting findings from Mendelian randomization using the MR-\nEgger method. Eur J Epidemiol 32, 377-389 (2017).\n10. Bowden, J., Davey Smith, G., Haycock, P.C. & Burgess, S. Consistent Estimation in Mendelian\nRandomization with Some Invalid Instruments Using a Weighted Median Estimator. Genet\nEpidemiol 40, 304-314 (2016).\n11. Cochran, W.G. The Combination of Estimates from Different Experiments. Biometrics 10, 101-\n129 (1954).\n12. Jonmundsson, T., et al. A proteomic analysis of atrial fibrillation in a prospective longitudinal\ncohort (AGES-Reykjavik study). Europace : European pacing, arrhythmias, and cardiac\nelectrophysiology : journal of the working groups on cardiac pacing, arrhythmias, and cardiac\ncellular electrophysiology of the European Society of Cardiology 25(2023).\n13. Hemani, G., et al. The MR-Base platform supports systematic causal inference across the human\nphenome. Elife 7(2018).\n14. Giambartolomei, C., et al. Bayesian test for colocalisation between pairs of genetic association\nstudies using summary statistics. PLoS Genet 10, e1004383 (2014).\n15. Liao, C.C., Xu, J.W., Huang, W.C., Chang, H.C. & Tung, Y.T. Plasma Proteomic Changes of\nAtherosclerosis after Exercise in ApoE Knockout Mice. Biology 11(2022).\n16. Zoccarato, A., et al. NRF2 activation in the heart induces glucose metabolic reprogramming and\nreduces cardiac dysfunction via upregulation of the pentose phosphate pathway. Cardiovascular\nresearch 121, 339-352 (2025).\n17. Hess, K., et al. Concurrent action of purifying selection and gene conversion results in extreme\nconservation of the major stress-inducible Hsp70 genes in mammals. Sci Rep 8, 5082 (2018).\n18. Wan, J., et al. Kallikrein augments the anticoagulant function of the protein C system in thrombin\ngeneration. Journal of thrombosis and haemostasis : JTH 20, 48-57 (2022).\n19. Pennacchio, L.A., et al. An apolipoprotein influencing triglycerides in humans and mice revealed\nby comparative sequencing. Science 294, 169-173 (2001).\n20. Kronenberg, F., et al. Plasma concentrations of afamin are associated with the prevalence and\ndevelopment of metabolic syndrome. Circ Cardiovasc Genet 7, 822-829 (2014).\n21. Nowicki, G.J., Ślusarska, B., Polak, M., Naylor, K. & Kocki, T. Relationship between Serum\nKallistatin and Afamin and Anthropometric Factors Associated with Obesity and of Being\nOverweight in Patients after Myocardial Infarction and without Myocardial Infarction. Journal of\nclinical medicine 10(2021).\n-18-\n22. Speliotes, E.K., et al. Genome-wide association analysis identifies variants associated with\nnonalcoholic fatty liver disease that have distinct effects on metabolic traits. PLoS Genet 7,\ne1001324 (2011).\n23. Gkatzionis, A., Burgess, S. & Newcombe, P.J. Statistical methods for cis-Mendelian randomization\nwith two-sample summary-level data. Genet Epidemiol 47, 3-25 (2023).\n24. Budoff, M. Triglycerides and Triglyceride-Rich Lipoproteins in the Causal Pathway of\nCardiovascular Disease. The American journal of cardiology 118, 138-145 (2016).\n25. Yuan, S., et al. Plasma proteins and onset of type 2 diabetes and diabetic complications: Proteome-\nwide Mendelian randomization and colocalization analyses. Cell reports. Medicine 4, 101174\n(2023).\n26. Ebana, Y., et al. A functional SNP in ITIH3 is associated with susceptibility to myocardial\ninfarction. Journal of human genetics 52, 220-229 (2007).\n27. Bortnick, A.E., et al. Plasma Proteomic Assessment of Calcific Aortic Valve Disease in Older\nAdults. J Am Heart Assoc 14, e036336 (2025).\n28. Lerman, D.A., Prasad, S. & Alotti, N. Calcific Aortic Valve Disease: Molecular Mechanisms and\nTherapeutic Approaches. European cardiology 10, 108-112 (2015).\n29. Jönsson, G., et al. Hereditary C2 deficiency in Sweden: frequent occurrence of invasive infection,\natherosclerosis, and rheumatic disease. Medicine (Baltimore) 84, 23-34 (2005).\n30. Harada, M., et al. G-CSF prevents cardiac remodeling after myocardial infarction by activating the\nJak-Stat pathway in cardiomyocytes. Nat Med 11, 305-311 (2011).\n31. Kurki, M.I., et al. FinnGen provides genetic insights from a well-phenotyped isolated population.\nNature 613, 508-518 (2023).\n32. Zhuang, L., et al. Fatty acid-binding protein 3 contributes to ischemic heart injury by regulating\ncardiac myocyte apoptosis and MAPK pathways. American journal of physiology. Heart and\ncirculatory physiology 316, H971-h984 (2019).\n33. Schachtl-Riess, J.F., et al. KLKB1 and CLSTN2 are associated with HDL-mediated cholesterol\nefflux capacity in a genome-wide association study. Atherosclerosis 368, 1-11 (2023).\n34. Araki, T., et al. Noonan syndrome cardiac defects are caused by PTPN11 acting in endocardium to\nen\nto enhance endocardial-mesenchymal transformation. Proc Natl Acad Sci U S A 106, 4736-4741\n(2009).\n35. Coan, P.M., et al. Complement Factor B Is a Determinant of Both Metabolic and Cardiovascular\nFeatures of Metabolic Syndrome. Hypertension (Dallas, Tex. : 1979) 70, 624-633 (2017).\n36. Peters, A.E., et al. Proteomic Pathways across Ejection Fraction Spectrum in Heart Failure: an\nEXSCEL Substudy. medRxiv (2023).\n-19-\nTable 3. Summary of network regulators, including findings from the present study, previously known biological\nlinks, effect directions, if druggable or status as a target of an available compound\nRegulator Evidence from MR/Colocalization in AGES, Effect direction on ACVD Druggable / available\nand prior known biology related to ACVD* and/or survival (Table S13) compound (Table S15)\nGain-of-function variant increases risk of MI. Highly\nITIH3 Risk Druggable\nexpressed in atherosclerotic lesion\nCausally linked to T2D and microvascular\nHSPA1B Risk Available compound\ncomplications\nDDX39B Causally related to T2D and MetS (AGES) Risk Not druggable\nCausally related to MetS (AGES). Causally linked\nHSPA1A Risk Available compound\nto T2D and microvascular complications\nCardiomyocyte-specific knockout mice show\nKEAP1 Protective Available compound\nimproved cardiac function\nC2 Deficiency increases atherosclerosis Risk Druggable\nPCDH8 None Risk Not druggable\nGABARAP Genetic links to hypertension and CVD Risk Not druggable\nCSH1/CSH2 None Protective Druggable\nGenetic links to diabetes and venous\nRFC4 Risk Not druggable\nthromboembolism\nIZUMO1 None Risk Not druggable\nCSF3 Promotes cardiomyocyte survival after MI Risk Druggable\nFABP3 Contributes to ischemic heart injury Risk Not druggable\nCausally protective against calcific aortic valve\nKLKB1 Protective Available compound\ndisease; influences cholesterol efflux\nDCTN2 None Risk Not druggable\nPositively associated with the metabolic syndrome,\nAFM Protective Druggable\nbut lower levels in MI patients\nAIPL1 None Protective Not druggable\nMRRF None Protective Not druggable\nCausally protective against MI and MetS (AGES).\nAPOA5 Protective Druggable\nKnockout leads to elevated plasma triglycerides\nGain-of-function mutations cause heart defects in\nPTPN11 Risk Available compound\nNoonan syndrome\nKnock-out mice show improved metabolic and\nCFB Risk Available compound\ncardiovascular features\nGAL3ST1 None Risk Not druggable\nNUDT21 None Risk Not druggable\nCOL28A1 Positively associated with HF and its subtypes Risk Druggable\nC11orf49 Causally linked to MI (AGES) Protective Not druggable\nLRP4 Causally linked to MI (AGES) Inconclusive Not druggable\n*Background on previously characterized biological roles of network regulators: ITIH326; HSPA1A/B25; KEAP116; C229;\nGABARAP31; RFC431; CSF330; FABP332; KLKB127,33; AFM20,21; APOA519; PTPN1134;CFB35; COL28A136.\n-20-\nSupplementary Fig. S11. Network visualization from the comparative CPN analysis between\nAGES and UK Biobank. Edges visualized were identified in both CPNs (FDR = 1%) and filtered\nto remove network regulators with single targets. 36 network regulators are labelled and shown\nin red and there are 281 unique target proteins shown in blue (some red regulators may also be\ntargets).\n-21-\nResponse to Review\nWe are pleased to submit the revised version of our manuscript, entitled “Circulating causal\nprotein networks linked to future risk of myocardial infarction” (NCOMMS-25-18281-B), for\nconsideration for publication in Nature Communications. We have carefully reviewed the Author\nChecklist and made all necessary changes, as reflected in the files and documents uploaded with\nthe resubmission. We thank the reviewers for their positive remarks.\nReviewer comments\nReviewer #1 (Remarks to the Author):\nThe revised manuscript is much stronger, particularly due to the replication in UKBB and the\nSuSiE colocalization analysis. I have no objections to its acceptance for publication in Nature\nCommunications.\nResponse: We appreciate the positive remarks.\nReviewer #2 (Remarks to the Author):\nNone.\nResponse: We thank the reviewer\nReviewer #3 (Remarks to the Author):\nI have no additional concerns\nResponse: We thanks the reviewer\n-1-"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The revised manuscript is much stronger, particularly due to the replication in UKBB and the SuSiE colocalization analysis. I have no objections to its acceptance for publication in Nature Communications.", "reviewer_3": "I have no additional concerns"}, "author_response": "Response to Review\nWe are pleased to submit the revised version of our manuscript, entitled “Circulating causal\nprotein networks linked to future risk of myocardial infarction” (NCOMMS-25-18281-B), for\nconsideration for publication in Nature Communications. We have carefully reviewed the Author\nChecklist and made all necessary changes, as reflected in the files and documents uploaded with\nthe resubmission. We thank the reviewers for their positive remarks.\nReviewer comments\nReviewer #1 (Remarks to the Author):\nThe revised manuscript is much stronger, particularly due to the replication in UKBB and the\nSuSiE colocalization analysis. I have no objections to its acceptance for publication in Nature\nCommunications.\nResponse: We appreciate the positive remarks.\nReviewer #2 (Remarks to the Author):\nNone.\nResponse: We thank the reviewer\nReviewer #3 (Remarks to the Author):\nI have no additional concerns\nResponse: We thanks the reviewer\n-1-"}]

Decision letter