Preeclampsia (PE) is a common and severe pregnancy complication, which is one of the leading causes of maternal mortality worldwide (Ives et al., 2020). It is characterized by new-onset hypertension after 20 weeks of gestation and multiple organ dysfunction (Brown et al., 2018). Women who survive PE may have a shorter life expectancy and an increased risk of stroke, cardiovascular disease, and diabetes (Poon et al., 2019; Kvalvik et al., 2020; Pittara et al., 2021). However, the mechanism of this hypertensive disorder in pregnancy remains elusive. Increasing evidence shows that the balance between pro- and anti-inflammatory factors is an intrinsic mechanism of occurrence and development of PE.

Several studies have indicated the potential involvement of systemic inflammatory regulators in the development of PE. Various investigations have delved into the pathophysiological roles of cytokines, including tumor necrosis factor alpha (TNF-α), interleukin-6 (IL-6), interleukin-8 (IL-8), interleukin-17 (IL-17), interleukin-18 (IL-18), interferon gamma (IFN-γ), interleukin-4 (IL-4), and interleukin-10 (IL-10) in the progression of PE (Lau et al., 2013; Yang et al., 2014; Guan et al., 2023). However, due to the presence of residual confounders, the conclusive establishment of a causal relationship between cytokines and PE remains elusive.

Mendelian randomization (MR) studies stand out as a distinctive and potent statistical approach to investigate causality between exposures (e.g., circulating cytokines) and outcomes (e.g., preeclampsia) of interest (Smith and Ebrahim, 2003; Hemani et al., 2018). This method uses genetic instruments as unconfounded proxies for exposures, which can avoid residual confounding and reverse causality that is commonly present in conventional observational studies (Burgess et al., 2013). In the absence of randomized clinical trials (RCTs), MR design is an important strategy for causal inference as genetic variants are randomly assorted at meiosis, in which the procedure mimics an RCT (Lawlor et al., 2008). Furthermore, a two-sample MR design using summary statistics from a genome-wide association study (GWAS) greatly increases the statistical power of causality inference (Pierce and Burgess, 2013). Two-sample MR analysis allows researchers to evaluate the relationship between instrument exposure and instrument results in two separate population samples, thereby enhancing the applicability and effectiveness (Hartwig et al., 2016).

Prior MR studies have investigated the causal influence of systemic inflammatory regulators on various diseases. Yeung et al. discovered that the considered systemic inflammatory regulators did not impact the risk of Alzheimer’s disease (AD). In contrast, specific cytokines such as interleukin-2 (IL-2), IFN-γ, TNF-α, and IL-18 might be downstream effects of AD (Yeung and Schooling, 2021). Song et al. (2022) demonstrated that elevated levels of IL-18 were correlated with a reduced risk of acute myeloid leukemia, while IL-17 was associated with the risk of stomach cancer. Wang et al. proposed that heightened levels of monocyte-specific chemokine-3 (MCP3), vascular endothelial growth factor (VEGF), IL-10, and IL-7 were linked to an increased risk of multiple myeloma (MM), whereas lower levels of tumor necrosis factor beta (TNF-β) were strongly associated with an elevated risk of MM (Wang et al., 2022).

Nevertheless, the bidirectional association between systemic inflammatory regulators and PE has not been explored to date. In this study, we employed a two-sample MR analysis to investigate the potential causal relationship between systemic inflammatory regulators and PE and vice versa.

The summary data for the GWAS on systemic inflammatory regulators and PE were obtained from previously published studies (refer to Supplementary Table S1), obviating the need for additional ethical clearance. This study adheres to the STROBE-MR guidelines, as outlined by Skrivankova et al. (2021). The fulfillment of the three principal criteria is imperative, as delineated by VanderWeele et al. (2014) (Ives et al., 2020): instrumental variables (IVs) exhibit a robust association with the exposure (Brown et al., 2018); IVs are independent of any unidentified confounders existing between the exposure and outcome (Poon et al., 2019); and IVs exclusively influence the outcome through their impact on the exposure.

Summary-level data for 41 systemic inflammatory regulators were acquired in this study from a GWAS, encompassing 8,293 individuals of Finnish descent from 3 distinct cohort studies: the Cardiovascular Risk in Young Finns Study, the FINRISK 1997 study, and the FINRISK 2002 study (Ahola-Olli et al., 2017). To ensure robust genetic associations, adjustments were made for age, sex, body mass index, and the top 10 genetic principal components.

The genetic association data pertaining to PE were sourced from a meta-analysis that included data from the UK Biobank, FinnGen, and BioBank. This meta-analysis involved 2,355 cases and 264,887 controls of European ancestry, as reported by Sakaue et al. (2021). Notably, this study represents a multi-country GWAS meta-analysis with minimal overlap with the GWAS on inflammatory regulators. Additionally, all single-nucleotide polymorphisms (SNPs) were derived from the analysis of European populations, thereby mitigating potential biases associated with population differences.

To mitigate the risk of false-positive IVs, we opted for statistically significant criteria (p < 5 × 10⁻⁶, linkage disequilibrium (LD) r ² < 0.001, and window size = 10,000 kb) in the summary data for systemic inflammatory regulators, as recommended by Rayes et al. (2023). Following the harmonization of the selected SNPs with those associated with PE, we identified a total of 41 circulating cytokines. To ensure the robustness of the instrumental variable, we verified that the F-statistics exceeded 10, adhering to the stringent mathematical formula: F = β²_exposure/SE²_exposure. The details of the 41 systemic inflammatory regulators chosen are provided in Supplementary Table S2.

For selecting IVs of PE, LD (r ² < 0.001) and proxies were tested for the 14 SNPs estimated to be correlated to PE at the genome-wide significance level (p < 5 × 10⁻⁶). Eventually, after harmonizing the exposure–outcome datasets, nine SNPs in total were included to construct the genetic IVs for PE (Supplementary Table S3).

We systematically validated potential causal relationships between systemic inflammatory regulators and PE through a robust methodology. Primary MR analyses were conducted using the inverse-variance weighted (IVW) approach with random effects, utilizing genetic data to estimate causal effects between exposures and outcomes. Odds ratios (ORs) and 95% confidence intervals (CIs) were reported, with statistical significance set at p < 0.05.

To enhance the reliability of our MR results, we implemented various heterogeneity and sensitivity assessments. Sensitivity analysis involved the utilization of the MR–Egger, weighted median, weighted mode, and simple mode techniques, as detailed by Xu et al. (2017). The MR–Egger method, employed to assess bias due to gene pleiotropic effects, demonstrated potential susceptibility and lower precision to a weak instrumental variable bias than other methods such as IVW and weighted median. The MR–Egger intercept was employed to identify horizontal pleiotropy, with statistical significance set at p < 0.05.

To further scrutinize pleiotropy and correct horizontal pleiotropy through outlier removal with potential pleiotropy, we employed MR Pleiotropy RESidual Sum and Outlier (MR-PRESSO) as proposed by Ong and MacGregor (2019) and Verbanck et al. (2018). Cochran’s Q test was utilized to evaluate heterogeneity and horizontal pleiotropy. The analysis of Cochran’s Q test results and funnel plots served as indices to estimate individual SNP heterogeneity, in accordance with the approach outlined by Bowden et al. (2018).

Additionally, we visually presented the results through scatter plots, forest plots, and leave-one-out plots, following the methodology described by Hemani et al. (2018). This comprehensive validation process strengthens the robustness and credibility of our findings regarding the causal relationships between systemic inflammatory regulators and PE. To account for multiple testing, we applied the Bonferroni method, which requires us to calculate associations with p-values below 0.0012 (0.05/41) as strong evidence of associations. The results with p-values ranging from 0.0012 to 0.05 were regarded as suggestive associations (Shi et al., 2023). Statistical analyses were performed and visualized using the R program (version 4.3.0) with the two-sample MR package (version 0.5.6) and MR-PRESSO (version 1.0) (Hemani et al., 2018).

The MR analysis uncovered a significant association between specific systemic inflammatory regulators—TNF-α, stem cell growth factor beta (SCGF-β), IL-9, and IL-5—and the risk of PE. All SNPs exhibited F-statistics greater than 10 (range 20.35–99.619; see Supplementary Table S2). Fourteen IVs were identified for SCGF-β, while IL-9, IL-5, and TNF-α had six, five, and five IVs, respectively. The IVW method demonstrated a significant difference (p < 0.05), and consistent directional effects were observed across five methods, namely, IVW, MR–Egger, weighted median, weighted mode, and simple mode (refer to Supplementary Table S4).

The IVW analysis revealed a negative association between SCGF-β (OR = 0.89, 95% CI: 0.80–0.99, p = 0.027) and IL-5 (OR = 0.80, 95% CI: 0.65–0.98, p = 0.030), indicating a decreased risk of PE. Although increased levels of interferon gamma-induced protein 10 (IP10) were associated with a reduced risk of PE, the IVW test results were inconclusive (p > 0.05), prompting the use of the MR–Egger method (OR = 0.66, 95% CI = 0.49–0.89, p = 0.026). The MR-Egger intercept did not indicate potential horizontal pleiotropy for SCGF-β and IL-5 (p = 0.968; p = 0.732, respectively), while IP10 exhibited pleiotropy (p = 0.033) (refer to Supplementary Tables S5, S6). Given that horizontal pleiotropy violates MR assumptions, the analysis involving IP10 was considered unreliable. Additionally, Cochran’s Q values based on the IVW tests showed no significant heterogeneity for SCGF-β and IL-5 (all p > 0.05; Supplementary Table S7).

Similarly, the IVW method revealed that elevated circulating levels of TNF-α and IL-9 were causally related to an increased risk of PE (OR = 1.32, 95% CI = 1.09–1.60, p = 0.004 and OR = 1.28, 95% CI: 1.02–1.62, p = 0.033, respectively). The MR–Egger intercept did not detect potential horizontal pleiotropy for TNF-α and IL-9 (p = 0.620 and p = 0.588, respectively). Cochran’s Q values based on the IVW tests indicated no significant heterogeneity for SCGF-β and IL-5 (all p > 0.05).

The MR results, along with heterogeneity analysis, pleiotropy analysis, and sensitivity analysis results, are summarized in Supplementary Tables S5–S7. No SNP significantly influenced the overall effect of cytokines on PE in the IVW leave-one-out sensitivity analysis. Supplementary Figures S1–S5 display scatter plots, funnel plots, forest plots, and leave-one-out plots, while Figure 1 presents the forest plots of the obtained results.

Furthermore, other systemic inflammatory regulators, including cutaneous T-cell-attracting chemokine, beta nerve growth factor, VEGF, macrophage migration inhibitory factor (MIF), TNF-related apoptosis-inducing ligand, TNF-β, MCP3, stromal cell-derived factor-1 alpha, stem cell factor, interleukin-12p70, interleukin-16, platelet-derived growth factor-BB, growth-regulated protein alpha, hepatocyte growth factor, interleukin-1 receptor antagonist, monocyte chemoattractant protein-1, macrophage inflammatory protein 1b, and IL-18, did not exhibit significant associations with PE in any of the analyses.

Using nine SNPs as IVs for PE, we demonstrated that genetically predicted PE is negatively associated with MIF levels through the IVW method (OR = 0.87, 95% CI = 0.79–0.96, p = 0.007) (Supplementary Tables S8). All SNPs we selected had F-statistics greater than 10 (mean, 22.83; range, 20.97–27.69). The MR–Egger intercept did not detect potential horizontal pleiotropy (p = 0.705; Supplementary Tables S9). Cochran’s Q values based on the IVW tests showed that there was no obvious heterogeneity (p = 0.434). Detailed data are given in Supplementary Tables S10, S11. Figure 2 shows the forest plots of the above results. Scatter plots, funnel plots, forest plots, and leave-one-out plots are listed in Supplementary Tables S5.