Genetically predicted triglyceride levels and endometriosis risk: Mendelian randomization with exploratory multi-omics analyses

Endometriosis (EM) is a prevalent condition among women of reproductive age, affecting approximately 190 million women worldwide [ 1 ]. It is characterized by the presence of endometrial tissue (comprising glands and stroma) outside the uterine cavity, leading to symptoms such as dysmenorrhea, abnormal menstruation, dyspareunia, infertility, and more. Affecting approximately 10% of women [ 2 ], EM significantly impacts quality of life. These symptoms disrupt daily activities and reduce work efficiency, with pelvic pain and disease severity accounting for up to 60% of productivity loss. EM is known for its estrogen dependence, but its exact underlying causes remain incompletely understood. Currently, laparoscopic surgery is considered the gold standard for diagnosing EM. However, due to its invasive nature, it may not be suitable for early diagnosis or monitoring the progression of the disease. To date, no peripheral blood or endometrial biomarkers have suggested sufficient accuracy in diagnosing EM. Given the chronic nature of EM, long-term pharmacological management is essential. While surgical intervention can provide symptomatic relief and remove existing lesions, it does not offer a permanent solution and carries significant costs. Ahmed M Soliman et al. found that the average annual cost of medical care for patients with endometriosis was more than three times that of non-endometriosis controls ($16,573 vs $4,733; p <  0.005) [ 3 ]. At the same time, the surgical recurrence rate was high. One systematic evaluation showed a persistence or recurrence rate of 22% at 2 years and 40%–50% at 5 years after surgery [ 4 ]. Therefore, the development and investigation of effective medications for long-term management and symptom control are important aspects of managing EM [ 5 ]. These medications can help alleviate pain and improve overall quality of life for individuals with EM. In summary, further exploration is needed to better understand the pathogenesis of EM, develop improved diagnostic methods, and establish effective disease control strategies. Given the unmet need for accurate diagnosis and effective treatment options for EM, exploring the potential of existing drugs to improve the condition of patients with EM may serve as a viable alternative. Considering the large population of individuals affected by EM, utilizing currently available medications may provide relief and help manage the symptoms associated with the condition. Recent clinical research has shown a significant elevation in triglyceride levels among patients with EM compared to healthy women [ 6 ]. Meanwhile, the findings of Hongling Zhang et al. provided evidence consistent with a potential causal relationship between TG and increased risk of EM (OR = 1.112, 95% CI: 1.033–1.198, p  = 5.03 × 10 ⁻ ³). They concluded that higher serum TG levels were associated with a greater risk of EM [ 7 ]. Additionally, several studies have established a positive correlation between triglyceride levels and the severity of EM [ 8 ]. Triglycerides (TG) play a crucial role in the structural composition of living organisms and actively engage in diverse signaling pathways, exerting varying degrees of influence on cellular functions. TG is implicated in processes such as inflammation, immunity, and steroid hormone metabolism, and their close association with the onset and progression of certain diseases has been well-documented. Nevertheless, the precise nature of the association between TG and EM remains elusive. Consequently, it is imperative to investigate the genetic association between TG and EM, thus fostering a deeper comprehension of the pathogenesis of this condition and subsequently formulating a theoretical basis for its diagnosis and treatment. Currently, the main triglyceride inhibitors encompass fibrates, ω-3 PUFA, and niacin medications. Fibrates primarily target peroxisome proliferator-activated receptor-α (PPARA), leading to a reduction in triglyceride (TG) levels by 20% to 50% [ 9 ]. Fibrates function as agonists of PPARA, binding to and activating this nuclear receptor. PPARA exhibits predominant expression in tissues actively engaged in lipid metabolism, including the liver, adipose tissue, and skeletal muscle. Upon binding to PPARA, fibrates modulate the expression of genes associated with lipid and lipoprotein metabolism, leading to diverse physiological effects, such as heightened fatty acid oxidation, diminished triglyceride synthesis, and enhanced high-density lipoprotein cholesterol (HDL-C) production [ 10 ]. While fibrates primarily exert their action on PPARA, it is worth noting that PPARA itself possesses the ability to regulate the expression of numerous genes and influence diverse metabolic pathways involved in lipid and glucose metabolism. Therefore, despite fibrates having a singular primary target, their effects can transcend PPARA and indirectly engage multiple targets by means of downstream effects resulting from PPARA activation. For instance, fibrate-induced activation of PPARA upregulates the expression of genes involved in fatty acid transport, oxidation, and lipoprotein metabolism [ 11 ]. Consequently, this can instigate alterations in the activity of enzymes, transporters, and receptors implicated in these processes. Moreover, PPARA activation can modulate inflammatory and oxidative stress pathways, thereby additionally influencing gene expression and cellular processes [ 12 ]. In summary, while fibrates primarily target PPARA, they indirectly affect multiple targets by modulating gene expression and metabolic pathways associated with lipid and lipoprotein metabolism. Therefore, investigating the role of fibrates in EM can enhance our understanding of the relationship between PPARA and multiple metabolic targets, thereby providing a more scientific and rational basis for the use of fibrates. Omega-3 polyunsaturated fatty acids (ω-3 PUFA) can lower triglyceride levels in the body through various mechanisms and by targeting multiple pathways. They can reduce TG levels by approximately 25% to 30% [ 13 ]. Regarding triglyceride synthesis, omega-3 fatty acids inhibit diacylglycerol acyltransferase (DGAT), an enzyme involved in the terminal stage of triglyceride synthesis. Through the reduction of DGAT activity, omega-3 fatty acids decrease triglyceride production in the liver [ 14 ]. Simultaneously, omega-3 fatty acids upregulate the activity of enzymes implicated in fatty acid oxidation, such as acetyl-CoA carboxylase 1 (ACC1) and fatty acid synthase (FAS), facilitating the breakdown (oxidation) of fatty acids in the liver. This reduces the availability of glycerol and fatty acids for triglyceride synthesis, consequently lowering triglyceride levels [ 15 ]. Additionally, omega-3 fatty acids can enhance the activity of lipoprotein lipase (LPL), an enzyme that facilitates the breakdown of TG into fatty acids and glycerol for utilization by tissues. This, in turn, increases the rate at which TG are cleared from the bloodstream [ 16 ]. Furthermore, omega-3 fatty acids can promote a favorable distribution of lipids, which includes reducing triglyceride levels, by modulating the expression of genes involved in lipid metabolism. Additionally, omega-3 fatty acids possess anti-inflammatory properties and contribute to cardiovascular health and neuroprotection [ 17 ]. Exploring the multi-target and multi-mechanism characteristics of omega-3 fatty acids may yield novel insights and treatment strategies for addressing the chronic inflammatory state, pain, and hormonal imbalance associated with EM. Niacin is currently utilized in cases where fibrates and prescription-grade omega-3 fatty acids have been ineffective in managing TG, albeit its notable adverse effects. While the mechanism of action of niacin remains unclear, niacin appears to mitigate triglyceride and low-density lipoprotein cholesterol levels by modulating multiple targets involved in triglyceride metabolism within adipose tissue and the liver. Research indicates that niacin activates GPR109A, leading to the inhibition of lipolysis (the breakdown of stored fat) in adipose tissue. This, in turn, decreases the release of free fatty acids into the bloodstream, thereby aiding in the reduction of triglyceride synthesis in the liver [ 18 ]. Furthermore, niacin exerts a direct and non-competitive inhibitory effect on diacylglycerol acyltransferase 2 (DGAT2), a crucial enzyme involved in hepatic triglyceride synthesis. Consequently, this inhibition leads to decreased triglyceride synthesis and reduced secretion of atherosclerotic lipoproteins from the liver [ 19 ]. Moreover, niacin medications can activate the PI3K/PKB pathway, leading to the inhibition of FAS activity in the liver. As a result, triglyceride synthesis is reduced [ 20 ]. Simultaneously, niacin plays a role in energy metabolism within the body, influencing cellular redox status and ATP production [ 21 ]. Abnormalities in cell metabolism are associated with the development of EM [ 22 ]. Investigating the mechanism of action of niacin in this context can offer valuable insights into the progression of the disease. Furthermore, niacin regulates immune cell function and modulates inflammatory responses [ 23 ], while promoting tissue repair [ 24 ] and exerting antioxidant effects [ 25 ]. EM frequently involves inflammatory responses to tissue damage and heightened oxidative stress. Investigating the mechanisms by which niacin promotes tissue repair and mitigates oxidative stress may contribute to the development of novel treatment strategies for EM. In conclusion, it is evident that a variety of triglyceride inhibitor drugs are currently available. However, selecting a suitable triglyceride inhibitor for EM patients remains a challenge. To provide genetic insights into the potential relevance of triglyceride-related drug targets in EM, our study employed the Mendelian randomization (MR) method to analyze the genetic associations between triglyceride inhibitor-related targets and EM risk. Our analysis aimed to evaluate the genetic evidence supporting the potential relevance of triglyceride-related targets in EM. Fig 1 provides a detailed overview of the study rationale.

Based on the analysis conducted, a total of 350 genome-wide significant SNPs associated with TG levels and 26 SNPs associated with EM were included as instrumental variables (IVs). The inverse-variance weighted (IVW) estimation was employed to assess the association between genetically predicted TG levels and EM risk. The results of the IVW estimation indicated that genetically predicted TG levels were significantly associated with EM risk (OR = 1.1853, 95% CI: 1.0704–1.3125, p  = 0.0011). This suggests that higher genetically predicted TG levels are associated with an increased risk of EM ( Fig 3A ). (A) Genetically predicted associations of TG on EM; (B) Genetically predicted associations of EM on TG. The Fig shows the results of two data sets, A and B, involving different exposures and outcomes, demonstrating a comparison of several statistical methods. For group A, the exposure was TG and the outcome was EM. OR and p values show the effect estimates and statistical significance, respectively, for the different methods. ORs and their 95% confidence intervals (CIs) are represented as dots and lines in the Fig. For group B, the exposure was EM and the outcome was TG. Beta and p values show effect estimates and statistical significance, respectively. Nsnps represents the number of SNPs used for each method. On the other hand, the analysis did not find a significant association between EM and TG levels (β = 0.0145, 95% CI: −0.0070 to 0.0361, p  = 0.1862). This indicates that EM was not significantly associated with TG levels ( Fig 3B ). Since the MR-Egger analysis yielded significant results ( p <  0.05) in the TG-EM analysis, the leave-one-out method was utilized to test the significance of individual SNPs on the outcome. The results suggested that by eliminating SNPs related to TG-associated SNPs one by one and recalculating the effects of the remaining SNPs, the overall results remained consistent, suggesting that the pooled effects were not affected ( S1 Fig ). These findings support the hypothesis that genetically predicted higher TG levels are associated with an increased risk of EM, while EM itself does not significantly impact TG levels. We analyzed the target genes of 35 fibrates, 12 niacins, and 10 ω-3 PUFA drugs included in DGIDB and ChEMBL. Out of these, 21 genes showed genetic associations with triglyceride concentration, and among these TG-associated genes, 7 genes showed genetic associations with EM in the present analyses. Choline Fenofibrate and Ciprofibrate target TNF and SLC9A1 significantly among fibrates, while Fenofibrate targets APOA1 , VEGFA , and SCN3A . Among niacin drugs, Probucol significantly targets ITGAV and SCN3A . SCN3A is a significant target for all three ω-3 PUFA drugs. Table 3 displays the information about the included drug targets. Detailed explanations for gene exclusion or inclusion are provided in S1 Table . All significant targets associated with TG showed positive associations with CHD risk (OR > 1, p  < 0.05) ( Fig 4 ). Exposure: Indicates genes; nsnp: Number of single nucleotide polymorphisms (SNPs) associated with the gene in the sample; p value: Significance level, a value less than 0.05 is considered significant; F statistic: The F statistic, used to assess instrument strength. OR (95% CI): Odds Ratio and its 95% confidence interval. An OR greater than 1 indicates increased risk, while less than 1 indicates reduced risk. In the TG-based drug-target MR analysis, genetically proxied effects related to SCN3A, VEGFA, APOA1, TNF, and BACE1 were associated with increased EM risk (OR > 1, p < 0.05 ), whereas genetically proxied effects related to SLC9A1 and ITGAV were associated with reduced EM risk (OR < 1 , p < 0.01). In sensitivity analysis, SCN3A and BACE1 , which utilized very low density lipoprotein (VLDL) as downstream markers, yielded similar results in TG (OR < 1, p   1, p  < 0.01) ( Fig 5 ). Target genes: The genes being analyzed for their association with triglyceride levels. Exposure: The type of lipid measurements includes TG for triglycerides and HDL-C as high-density lipoprotein cholesterol. nsnp: Number of single nucleotide polymorphisms (SNPs) associated with the respective gene. p value: Significance level; values below 0.05 were considered statistically significant. F statistic: The F statistic, used to assess instrument strength. OR (95% CI): Odds Ratio and its 95% confidence interval. An OR greater than 1 indicates increased risk, while less than 1 indicates decreased risk. If the OR is less than 1 and the confidence interval does not include 1, this suggests an association with decreased risk. Narrower confidence interval indicates a more precise estimate. The PP.H4 values for all seven genes were below 0.8, indicating limited evidence that TG and EM share a common causal variant within the examined target-gene regions. In contrast, the PP.H3 value for TNF exceeded 0.8, suggesting that TNF may be associated with both traits, but through independent genetic variants. For APOA1 , VEGFA , SLC9A1 , and BACE1 , PP.H1 values were higher than the other hypotheses, suggesting stronger evidence for association with TG than with EM in these regions. Overall, the colocalization results did not provide strong support for shared causal variants between TG and EM at the examined drug-target loci ( Table 4 ). Given that fibrates serve as ligands for PPARA, we also examined the association between PPARA and the targets under investigation. The batch effects were removed from the data. For details, see S2 Fig A and B. Compared with the healthy control group, the EM group exhibited upregulated expression of SCN3A ( p  < 0.001), BACE1 , and APOA1 ( p   0.05) ( Fig 6C ). (A) Gene-expression correlation matrix of the eight selected genes in EM patients. (B) Gene-expression correlation matrix of the eight selected genes in healthy controls. (C) Box plots comparing the expression levels of the eight selected genes between EM patients and healthy controls. Red indicates a positive correlation, whereas blue indicates a negative correlation. “ns” indicates no statistical significance. Significance: p  < 0.05 (*), p <  0.01 (**), p <  0.001 (***), and p  < 0.0001 (****). Gene-expression correlation analysis revealed significant negative correlations between SCN3A and PPARA (correlation coefficient = −0.31, p  < 0.05), as well as between SCN3A and VEGFA (correlation coefficient = −0.46, p  < 0.01), within the healthy control group ( Fig 6B ). In the EM patient group, positive correlations were observed between SLC9A1 and VEGFA (correlation coefficient = 0.42, p  < 0.01), as well as between SCN3A and ITGAV (correlation coefficient = 0.27, p  < 0.05). Conversely, negative correlations were found between SLC9A1 and ITGAV (correlation coefficient = −0.36, p  < 0.01), BACE1 and TNF (correlation coefficient = −0.33, p  < 0.05), APOA1 and TNF (correlation coefficient = −0.27, p  < 0.05), APOA1 and VEGFA (correlation coefficient = −0.43, p  < 0.01), and PPARA and VEGFA (correlation coefficient = −0.34, p  < 0.05) ( Fig 6A ). The chromosomal locations of the eight genes are shown in S2 Fig C. After including PPARA , interactions were observed among seven out of the eight genes. VEGFA, TNF , and PPARA were co-mentioned textual associations (i.e., they were jointly mentioned in the PubMed abstract) as well as co-expression interactions between pairs. PPARA and APOA1 were textually associated, co-expressed, and experimentally confirmed to interact. A textual interaction relationship was indicated between APOA1 and TN F. While TNF and BACE1 exhibited both textual and co-expression interaction relationships. BACE1 and SLC9A1 showed both textual and co-expression interaction relationships, as shown in Fig 7A . (A) PPI network created with STRING. (B) MCC node-importance analysis. In Fig 7A , background colors denote specific biological processes, and different edge styles represent different types of interactions. In Fig 7B , the color coding indicates node importance. The GO enrichment analysis results for the eight genes identified the top five enriched pathways with the smallest FDR values, listed in ascending order, as follows: Regulation of localization, Regulation of transport, Regulation of vesicle-mediated transport, Regulation of biological quality, and Response to external stimulus. All genes were associated with the Regulation of localization pathway. While seven genes ( Tnf, Slc9a1, Apoa1, Ppara, Vegfa, Itgav, Scn3a ) were enriched in the regulation of transport pathway, as depicted in Fig 7A . Gene symbols in the STRING-derived network are shown in the original output format from the STRING database. Since STRING-based PPI network analysis showed no direct interaction between SCN3A and the other six genes, we utilized the GeneMANIA website to explore indirect relationships between SCN3A and the other six genes. Additionally, we employed GeneMANIA to predict potential associated genes that exhibit significant associations with these seven genes. In the regional network, PPARA , TNF , and APOA1 exhibited clear physical interaction and pathway relationships. These three genes regulated lipid transport, while PPARA, APOA1 , and ITGAV were involved in lipid transport. APOA1, VEGFA, ITGAV , and SLC9A1 collectively mediated cell-substrate adhesion, TNF and VEGFA collectively participated in the regulation of vasculature development. BACE1 exhibited Pathway connections with KLC1 to KLC4 , whereas SCN3A suggested Genetic Interactions connections with KLC4 and SLC9A1 ( Fig 8 ). Each node represents a gene. Lines connecting the nodes signify different types of interactions among the genes. The colors and styles of the edges indicate the nature of these interactions. Larger nodes indicate greater importance and more interactions. The main regulatory factors for the analyzed genes were identified using the Cytoscape plug-in iRegulon, showing that PPARA and estrogen receptor 1 (ESR1) exhibited the highest normalized enrichment score (NES) values among the transcription factors (TFs). PPARA functions as a transcription factor for 5 genes, including itself, while ESR1 acted as a transcription factor for the 7 genes (excluding TNF) . The number of motifs associated with PPARA and ESR1 significantly exceeded that of the transcription factor ranked third in terms of NES ( Fig 9 ). (A) Transcription factor regulatory network. Purple circles indicate target genes regulated by TFs, and green octagons indicate TFs. (B) Transcription factor normalized enrichment score (NES). Fig 9A shows the transcription factor regulatory network. Purple circles show target genes regulated by transcription factors. Green octagons indicate transcription factors. Connecting lines indicate regulatory relationships between transcription factors and target genes. The direction of the line indicates the direction of regulation. Fig 9B shows the normalized enrichment score (NES) of transcription factors. The table lists the NES values, the number of target genes and the number of motifs associated with different transcription factors.

This study presented MR evidence suggesting that TG-related pathways and triglyceride-lowering drug targets may be relevant to EM biology. Three classes of compounds were identified: fibrates, nicotinic acid, and ω-3 PUFA. Seven candidate drug-target-related genes were implicated: TNF , SLC9A1, APOA1, VEGFA, ITGAV, SCN3A, and BACE1 . These selected targets formed an interaction network associated with biological processes potentially relevant to EM, including lipid transport, localization, cell adhesion, and blood vessel development. PPARA and ESR1 were identified as the two TFs with the highest enrichment scores in the regulatory analysis. Considering the genetic evidence and the observed involvement of fibrate-related targets, fenofibrate-related pathways may warrant further mechanistic and translational investigation in EM.

Firstly, we utilized a two-sample univariate MR analysis to estimate the genetically predicted association of TG on EM. This analysis helped evaluate the genetic association between TG and EM and assess whether the findings were compatible with a potential causal effect. Secondly, we conducted a search in the drug-gene interaction database to identify target genes associated with triglyceride drugs. By leveraging this database and referring to previous studies on druggable genes, we aimed to characterize specific genes that are targeted by triglyceride drugs and potentially play a role in the development or treatment of EM. Thirdly, we employed the genetic variants that are mediated by the target genes and associated with TG to serve as proxies for TG-lowering drug targets. Using a two-sample MR approach, we investigated the association between genetic proxies for TG-lowering drug targets and EM risk by analyzing data from two Genome-Wide Association Studies (GWAS). Fourthly, we conducted colocalization analysis to examine whether the association between TG and EM was driven by the same genomic locus. This analysis could help determine if there was shared genetic regulation underlying both traits. Fifthly, we explored the protein interaction relationships among the selected genes and examined their enrichment functions. Additionally, we evaluated the significance and importance of the identified genes in the context of EM. Finally, we searched for transcription factors (TFs) that regulated these genes, aiming to identify the regulatory mechanisms that influenced the expression and function of the identified genes in relation to EM. The overall study design and analysis procedure are displayed in Fig 2 . The graphical summary is created by BioGDP [ 26 ]. Based on the information provided, the data used in the study were obtained from public databases and specific GWAS consortia. The sources of the GWAS data used in the analysis are as follows: TG and high-density lipoprotein (HDL) data: The data were obtained from The Global Lipids Genetics Consortium (GLGC). This consortium conducts large-scale genetic studies to unravel the genetic architecture of lipid traits. The GLGC study included 188,577 participants, predominantly of European ancestry, consistent with the majority of its constituent studies. Plasma lipid concentrations were measured after a minimum fasting period of 8 hours. Estimates were adjusted for age, age squared, sex, and population stratification, with exclusion of participants known to be taking lipid-lowering medications. The selected genetic variants collectively explained 10–14% of the total trait variance. For individual variant association estimates, additive genetic models were fitted using inverse-normal transformed trait linear regression. For pooled effect estimates, weighted meta-analysis was performed using Stouffer’s method. In the European-specific analysis (approximately 1.3 million individuals, accounting for 80% of the total sample), the identified lipid-associated signals constituted 76% of all signals. The large sample size ensures that the effect estimates in the European population are relatively stable, and the overall results exhibit high consistency. HDL cholesterol was standardized, with the genetically predicted association reflecting the effect of a 1 standard deviation increase in HDL on the outcome. Triglycerides are log-transformed due to skewed distribution, with the genetically predicted association interpreted as the effect of doubling the triglyceride level [ 27 ]. Very low-density lipoprotein (VLDL) data: The data were sourced from the IEU OpenGWAS project. The IEU OpenGWAS project is a platform that provides access to a diverse array of GWAS datasets for different traits. VLDL levels were analyzed as a continuous variable, typically reported in mmol/L (standard unit in European GWAS consortia). Sample size: 115,082 participants. EM and coronary heart disease (CHD) data: The data for these conditions were obtained from the FinnGen study. FinnGen is a research project in genomics and personalized medicine that aims to understand the genetic basis of diseases in the Finnish population. The FinnGen database includes 377,277 participants. In the FinnGen database, disease diagnoses are primarily derived from the Finnish National Health Registries and electronic health records (EHRs), encompassing inpatient and outpatient diagnoses, medication prescriptions, laboratory test results, and medical procedure documentation. All disease information is standardized using International Classification of Diseases (ICD) codes and cross-validated across multiple data sources to ensure accuracy. For specific conditions such as cancer, diagnoses are further validated through disease-specific registries (e.g., cancer registries) and pathological confirmation. Rare diseases are additionally corroborated using genomic data to enhance diagnostic precision [ 28 ]. It is important to note that the use of these data was approved by the respective institutional review boards, ensuring ethical considerations and participant privacy. Informed consent was obtained from all participants involved in the original GWAS studies. See Table 1 for detailed information on data sources. Single nucleotide polymorphisms (SNPs) that meet the genome-wide significance level (p  < 5 × 10 −8 ) were chosen as instrumental variables (IVs). Independent genetic variants were identified using a cutoff of linkage disequilibrium (LD) values (threshold set to r2 < 0.001, kb = 10,000) to ensure the independence of the IVs [ 29 ]. Proxy SNPs were not utilized in this study. SNPs lacking necessary statistical information were directly removed. The inverse variance weighting (IVW) method was utilized as the primary statistical approach. In the case of a single SNP, the genetically predicted association was estimated using the Wald ratio test. Weighted median, Simple mode, and Weighted mode methods were employed to assess model robustness. MR-Egger can identify horizontal pleiotropy in models, although it has limited statistical power [ 30 ]. Therefore, in cases where MR-Egger yielded significant results, we investigated the included SNPs using PhenoScanner, a comprehensive genotype-to-phenotype database [ 31 ]. Additionally, to exclude SNPs that may influence EM through other phenotypes, we employed the leave-one-out method to assess whether a single SNP exerted a disproportionate impact on the overall MR estimate [ 32 ]. To mitigate the possibility of reverse causation, we conducted an additional analysis to assess whether EM might influence triglyceride levels. The primary triglyceride-lowering drugs comprise fibrates, ω-3 PUFA, and niacin. Genes encoding the target proteins of these triglyceride drugs were acquired from DGIDB and ChEMBL. All targets gathered from these two sources were verified using NCBI Gene to determine if they were human genes. Non-human genes were excluded. We excluded genes that were not predicted with 90% confidence to be active according to ChEMBL’s Target Predictions. We selected SNPs that exhibited a significant association with TG levels ( p  < 5 × 10 −8 ) within a 100 Mb window surrounding the genomic region of the drug target [ 33 ]. To enhance the instrumental variable’s strength for each drug target gene, a less stringent threshold for independent clustering SNPs was employed (r 2  < 0.3, kb = 100). Restricting instrument selection to variants located within or near pharmacologically annotated target-gene regions was intended as a biologically informed strategy to strengthen the link between the selected variants and the predefined target-related pathways. Consistent with current methodological guidance, these regional instruments were treated as proxies for target-related genetic variation rather than as direct equivalents of pharmacological target inhibition or activation [ 34 , 35 ]. Positive-control MR analyses were used to assess the biological plausibility of the pharmacogenetic instrument framework in a lipid-related disease context. Because triglyceride-lowering medications are used in the management of CHD [ 36 ], and genetic studies have provided evidence consistent with a potential causal relationship between TG and CHD risk [ 37 ], CHD was selected as a positive-control outcome [ 38 ]. Results from this analysis were interpreted only as a supportive consistency check for the target-region instrument framework, rather than as validation of gene-specific pharmacological effects or evidence of therapeutic efficacy for EM [ 34 , 35 ]. In clinical studies, fibrates and niacins have been shown to reduce the risk of atherosclerosis by increasing concentrations of HDL-C and decreasing concentrations of TG in plasma [ 39 , 40 ]. Consequently, HDL-C is utilized as a downstream biomarker of the effects of fibrates and niacins. ω-3 PUFA reduces the synthesis and secretion of VLDL particles and enhances the removal of TG from VLDL and chylomicron particles through the upregulation of LPL [ 41 ]. As a result, VLDL is employed as a downstream biomarker of the effects of ω-3 PUFA. We employed colocalization analysis to investigate whether TG and EM share the same causal variant within the drug target gene region. This analysis aimed to provide additional validation regarding the robustness of including IV. In this study, we assessed the possible probabilities using the following hypotheses: H0 represents the scenario where the gene is not associated with TG or EM; H1 suggests that the gene is associated with TG but not with EM; H2 indicates that the gene is associated with EM but not with TG; H3 suggests that the gene is associated with both TG and EM, but the SNPs are independent of each other; and H4 proposes that the gene is associated with both TG and EM and that they share common SNPs. All analyses were conducted using R 4.2.2 statistical software. The ‘TwoSampleMR’ package was utilized for MR and drug target studies, while the ‘coloc’ package was employed for colocalization analysis. The Gene Expression Omnibus (GEO) database is a publicly accessible repository that collects high-throughput gene expression data from research institutions worldwide. It is maintained by the National Center for Biotechnology Information (NCBI) [ 42 ]. The extensive collection of experimental data in the GEO database is widely utilized in various research fields, including gene regulation, disease mechanisms, and drug discovery. Our objective was to investigate the expression patterns and potential regulatory relationships among genes involved in EM using the GEO database. To accomplish this, we downloaded six transcriptome datasets related to EM from the GEO database. The specific datasets we obtained for analysis are GSE23339 [ 43 ], GSE25628 [ 44 ], GSE58178 [ 45 ], GSE11691 [ 46 ], GSE7846 [ 47 ], and GSE7305 [ 48 ]. The details of each data set are shown in Table 2 . In the statistical analysis conducted using R software, we performed preprocessing steps on the gene expression profiles obtained from the downloaded datasets. The following steps were carried out: ID conversion: We performed ID conversion on the probes present in the expression profiles, converting them into corresponding gene names. This conversion enabled a more interpretable representation of the data. Log2 transformation: To standardize part of the unstandardized dataset, we applied the log2 function. This transformation helped normalize the expression values, ensuring a more symmetric distribution and facilitating subsequent statistical analyses. Extraction and merging of expression profiles: Next, we extracted the expression profiles of genes from each dataset separately and merged them. This consolidation allowed for a comprehensive analysis using a combined dataset. Batch effect removal: Batch effects were defined as technical differences resulting from the processing and measurement of samples in different batches, unrelated to any biological variation in the experiment. To enhance the accuracy of statistical inference, we used the ‘sva’ package to remove batch effects and other unwanted noise from the merged data. This step helped eliminate potential confounding factors and reduced the impact of batch effects on the results [ 49 ]. Principal Component Analysis (PCA): Prior to and after removing the batch effect, we generated a PCA graph. This analysis helped assess the reduction of batch effects by visually comparing the distribution of samples. The aim was to minimize the influence of batch effects on the results and improve the robustness of subsequent analyses. After completing these preprocessing steps, we obtained a total gene expression profile comprising 101 samples, including data from 56 EM patients and 45 healthy controls. This processed dataset serves as the basis for further research and analysis. To explore the expression patterns and mutual regulatory relationships of genes in the EM and control groups, we employed various R packages for visualization and analysis. The specific steps are as follows: Histograms of gene expression: We utilized the ‘reshape2’ package and the ‘ggpubr’ package to create histograms comparing the gene expression levels between EM patients and healthy controls. This visualization helped identify any differences in gene expression patterns between the two groups. Spearman correlation test: To assess the correlation of gene expression levels within the EM patient and healthy control groups separately, we employed the ‘corrplot’ package and performed Spearman correlation tests [ 50 ]. This analysis allowed us to evaluate the strength and direction of the relationships between gene expression levels. The results of these tests were visualized using the ‘ggplot2’ package. Visualization of gene positions and genome structure: We utilized the ‘RCircos’ package to visualize the relationship between the positions of the eight genes on the chromosome and the genome structure. This visualization approach provides insights into the spatial organization and arrangement of these genes within the genome [ 51 ]. By implementing these steps, we can gain a better understanding of the expression patterns, correlations, and genomic relationships of the genes involved in EM and control groups. The visualizations generated through these analyses facilitate the interpretation and communication of the results. To explore the direct correlation between genes, the drug targets identified through MR analysis were input into the STRING database to establish a protein-protein interaction (PPI) network. This network helped uncover the interactions among the genes of interest. In addition, a PPI network was also constructed in GeneMANIA to further examine the genes with significant indirect correlations and the functional roles they play. GeneMANIA provided insights into gene functions and interactions. The complementary nature of STRING and GeneMANIA has been suggested in several studies, with the former providing direct interactions and the latter extending indirect associations to improve the comprehensiveness and reliability of network-based functional annotation. STRING was used to initially construct the PPI network and identify core genes. GeneMANIA was further validated through enrichment analysis of the functional modules, which was cross-checked with the STRING results to reduce false positives [ 52 ]. To gain more insights from the constructed PPI networks, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment pathway analyses were performed. The thresholds for significance were set as false discovery rate (FDR)  0.01. These analyses helped identify the biological processes and pathways associated with the genes in the network. The threshold of > 0.01 for interaction strength in STRING was selected based on the database’s recommendations and established practices in network analysis [ 53 ]. In addition, a threshold of strength >0.01 preserves weak but potentially biologically significant associations while avoiding excessive stringency leading to the omission of potential pathways. In our study, this threshold allowed us to broadly capture potential interactions among lipid metabolism and inflammation-related genes, which were further validated through functional enrichment and sensitivity analyses [ 54 ]. To evaluate the importance of the target genes in the network, maximal clique centrality (MCC) analysis was conducted using the CytoHubba plugin in Cytoscape software. This analysis assessed the centrality of the target genes, indicating their significance within the network. Furthermore, the iRegulon plugin v1.3 [ 55 ] in Cytoscape software was employed to predict potential TFs associated with the key genes in the network. The enrichment score threshold was set to 3.0, area under the curve (AUC) threshold for the receiver operating characteristic (ROC) analysis was set to 0.03, ranking threshold was set to 5000, and the maximum FDR on motif similarity was set to 0.001. This analysis helped identify TFs that may regulate the expression of the key genes. By performing these analyses and utilizing the mentioned tools and thresholds, this study aimed to gain insights into the gene interactions, functions, regulatory factors, and pathways associated with the identified drug targets.

(TIF) (A) Analysis of the PCA plot prior to batch effect removal revealed distinct clustering of samples from the same dataset and clear stratification between batches. These findings indicate the presence of substantial batch effects across datasets, which could potentially confound genuine biological variations. (B) Subsequent application of the sva package for batch effect correction resulted in a more uniform distribution of samples across the two-dimensional plane, thereby substantially mitigating the batch effect and enhancing result accuracy. (C) The chromosomal locations of the eight genes. (TIF) (DOCX) The compressed archive contains Fibrate target gene information.xlsx, Niacin target gene information.xlsx, and ω-3 PUFA target gene information. (ZIP)