Causal Relationship Between Methotrexate and Chronic Obstructive Pulmonary Disease: A Study Based on Network Toxicology and Mendelian Randomization

Introduction

Chronic obstructive pulmonary disease (COPD) is a heterogeneous lung disease characterized by persistent abnormalities in the airways and/or alveoli, often leading to progressive airflow limitation. The main symptoms include chronic respiratory symptoms such as dyspnea, cough, and expectoration.¹ The pathogenesis of COPD involves multiple etiologies such as genetic susceptibility, abnormal lung development, and infections. Although there is heterogeneity in its pathophysiological processes, all share the common pathological node of aplasia of distal airway progenitor cells, ultimately leading to abnormal lung function and irreversible airflow limitation through airway structural remodeling.^1,2

Compared with the past view that regarded COPD as a self-inflicted disease mainly caused by smoking and more common in elderly men, some researchers now propose a new view: COPD is a clinical comprehensive disease caused by the interaction between genes (G) and the environment (E) during an individual’s life cycle (T).³ The G-E-T (GETomics) interaction not only damages the lungs but may also affect the normal development and aging process of the lungs.³ Therefore, although smoking cessation is crucial for the prevention of COPD, identifying various exposure factors in an individual’s life is more helpful in reducing the overall risk of developing COPD. Due to its high morbidity and mortality, COPD has now become the third leading cause of death globally.⁴ This not only brings a heavy burden to society but also seriously affects the quality of life of patients and their families. In view of this, in-depth research on the prevention strategies, pathogenesis, and treatment methods of COPD has important academic and clinical value.

Since methotrexate (MTX) was used in the treatment of rheumatoid arthritis (RA) in 1951, it has been gradually and widely used in psoriasis, ectopic pregnancy, myasthenia gravis, Crohn’s disease and a variety of tumors (such as osteosarcoma, breast cancer, bladder cancer, etc).^5–12 The clinical application dosage of MTX varies according to the diseases being treated: when treating RA, low-dose methotrexate (LD-MTX) exerts an anti-inflammatory effect mediated by adenosine; when treating tumors, high-dose methotrexate (HD-MTX) plays a therapeutic role through an anti-proliferative effect.¹³ Although MTX is generally well-tolerated, it can still have toxic side effects on the hematological system, gastrointestinal tract, liver, kidneys, and lungs. A randomized trial in 2020 showed that after using LD-MTX, the incidence of adverse events in multiple systems such as the hematological, digestive, and respiratory systems increased. Among them, the risk of adverse reactions in the lungs increased by 42%.¹⁴ In addition, RA patients using MTX have a higher risk of developing COPD than the general population, and the lifetime risk of COPD in these patients increases by 3.4%.¹⁵ However, a study involving 58,580 outpatients showed that MTX could reduce the risk of COPD exacerbation within the first six months after enrollment in the study.¹⁶ Given the differences in the above research results and the fact that the existing evidence has not fully clarified the differential effects of MTX on different populations (COPD patients and healthy people), this study used the method of Mendelian randomization (MR) combined with network toxicology to explore whether there is a causal association between MTX and the risk of developing COPD.

Network toxicology integrates the principles of network pharmacology, toxicology, and systems biology, and can construct multi-level pathways of chemical substances-targets-toxicity, providing relevant evidence for our research.¹⁷ This method can transform the complex interactions among components, diseases, and targets into a more intuitive visual presentation. Among them, molecular docking is a prediction technology of key-lock recognition based on computer simulation. It can simulate the specific binding patterns of biomolecules by constructing a three-dimensional interaction model of ligand-receptor at the atomic level, and calculate the binding free energy (ΔG) to evaluate the relative stability of these patterns.¹⁸

The core logic of Mendelian randomization (MR) is based on the principle of “random allocation of parental alleles” in Mendel’s laws of inheritance. Genetic variables are screened through the MRC-IEU Open Genome-Wide Association Study (OpenGWAS), and genetic variations are used as instrumental variables to provide anti-confounding evidence for the causal association in this study.¹⁹

The two methods have an inherent logical connection in this study: Network toxicology reveals the specific binding mechanism between MTX and key targets at the microscopic level through molecular docking technology, providing potential biological pathway explanations for MR analysis; MR, on the other hand, verifies whether this molecular interaction translates into disease risk associations at the population level through genetic instrumental variables at the macroscopic level. The two methods verify the hypothesis of this study at different levels, jointly forming a complete evidence chain from molecular mechanism to population effects.

Materials and Methods

Research Design

In this study, network toxicology was first used to systematically predict the core targets and key pathways related to the toxic side effects of MTX on COPD, providing a hypothesis orientation for subsequent verification. Among them, molecular docking mainly simulates the interaction mode at the microscopic level and evaluates the binding ability, providing key computational evidence for the direct binding events between “components-targets”. Meanwhile, we also conducted a MR analysis.

To ensure the steady progress of the research, we screened out SNPs that met the following three assumptions: (1) Association: The reproducibility and strong correlation between the working variable and the exposure factor (MTX); (2) Independence: The working variable is not related to confounding factors; (3) Exclusivity: The working variable must affect the outcome factor (COPD) through the exposure factor (MTX) and cannot affect the outcome factor (COPD) through other factors.²⁰

To verify the robustness of the research results, this study strictly followed the relevant recommendations of the “Strengthening the Reporting of Observational Studies in Epidemiology-Mendelian Randomization” (STROBE-MR) to conduct sensitivity analysis, so as to ensure the clarity and reproducibility of the research conclusions.²¹ The complete experimental design is shown in Figure 1.

Figure 1 Research design flowchart.

Data Sources and SNP Selection

In this study, the data on the exposure factor (MTX) and the outcome factor (COPD) were sourced from the OpenGWAS database (https://gwas.mrcieu.ac.uk/). For MTX (Prescription information, ID: ukbb-12209), data from 462,933 individuals of European descent were selected (2454 cases and 460,479 controls, all were mixed-gender samples), along with 9,851,867 SNPs. For COPD (ID: ukb-a-543), a total of 337,199 European patients were included (874 cases and 336,325 controls, all were mixed-gender samples), along with 10,894,596 SNPs. To ensure that there was no significant linkage disequilibrium (LD) among the selected SNPs, LD clustering analysis was performed using the European reference data from the 1000 Genomes Project, with a threshold of r² < 0.001 and a window size of 10 Mb. Meanwhile, to ensure the absence of weak instrument variable bias, the F - value should be relatively large (F - statistic > 10). After excluding confounding factors (such as demographic and behavioral factors like age, gender, and smoking, or SNPs directly related to COPD) according to the core assumptions of MR using PhenoScanner (http://www.phenoscanner.medschl.cam.ac.uk/phenoscanner), Then, conduct data harmonization: unify the allele directions, align the effect alleles, eliminate palindromic sequences and SNPs with intermediate allele frequencies, and exclude SNPs that may have pleiotropy. And, four SNPs (rs1230666, rs6605556, rs2894316, and rs9276609) were finally selected as the instrumental variables for MTX (Table 1).

Table 1 Instrumental SNPs from GWAS on Methotrexate (MTX) and Chronic Obstructive Pulmonary Disease (COPD)*

Network Toxicology

Use SuperPred (http://bioinformatics.charite.de/superpred), Swiss Target Prediction (http://swisstargetprediction.ch), and STITCH (https://cn.string-db.org/) to screen the targets of the interaction between MTX and human proteins. Meanwhile, use “COPD” as the keyword to search for the disease genes of COPD in OMIM (https://omim.org/), Gene Cards (https://www.genecards.org/), and TTD (https://db.idrblab.net/ttd/). Use the JVENN platform (https://jvenn.toulouse.inrae.fr/app/example.html) to draw a Venn diagram of the intersection of the targets of MTX and COPD. Import the intersecting targets into the Metascape database (https://metascape.org/) to perform gene ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. Among them, GO functional enrichment analysis is used to evaluate the enrichment significance level of differential genes within specific GO entries, and KEGG pathway enrichment analysis is used to clarify the enrichment of genes in metabolic pathways and signal transduction pathways. Import the screened genes into the STRING database (https://cn.string-db.org/), select the species Homo sapiens. After initially constructing the PPI network diagram, import it into Cytoscape 3.10.3 for refinement to highlight the core genes involved in COPD and prepare for subsequent molecular docking.

Based on molecular docking technology, this study predicted the binding mode and binding free energy between MTX (small molecule ligand) and protein (target receptor) through conformational sampling and binding free energy scoring. Using “Methotrexate” as the keyword, we searched in CAS Common Chemistry (https://commonchemistry.cas.org/) and downloaded the 2D structure of MTX. The downloaded file was imported into Chem3D for optimization and finally exported in mol2 format. The top nine genes (HIF1A, EP300, NFKB1, PPARG, ABL1, STAT1, HDAC1, NFE2L2, and MAPK1) in the core gene network were selected for molecular docking. To ensure accuracy, we first searched for their corresponding proteins in UniProt (https://www.uniprot.org/), selected the most appropriate gene IDs with “human” and “reviewed” as screening criteria, searched for the protein structures encoded by these genes in the AlphaFold Protein Structure Database (https://alphafold.ebi.ac.uk/) using the IDs as keywords, and downloaded them in pdb format. After removing water molecules and small molecule ligands from the structures using PyMOL 2.5.7, the structures were exported. Using Autodock Vina 1.2.5 software, after importing the above pdb and mol2 files, the spatial discretization of the active site and the pre-calculation of the potential field were achieved through the setting of distance parameters, and finally the protein-ligand molecular docking simulation was completed.

MR Analysis

Mendelian randomization (MR) is a method widely used to explore the causal effects of exposure factors on outcomes. This study mainly employed the following five statistical methods: Inverse-variance weighting (IVW), MR-Egger regression, Weighted Median, Weighted Mode, and Simple Mode. Among them, IVW makes full use of the information of all instrumental variables, so it has the highest statistical power. When all assumptions are met, IVW can provide the most accurate estimation of causal effects, so IVW was used as the main method. However, IVW has the drawback of assuming that all instrumental variables are valid (ie, there is no horizontal pleiotropy). Therefore, to ensure the robustness of the causal effect estimation, we used the MR-Egger regression to detect and correct horizontal pleiotropy, the Weighted Mode to control the bias of weak instrumental variables, Cochran’s Q to conduct heterogeneity analysis, and the leave-one-out analysis and MR-PRESSO test for sensitivity testing. These methods enhance the credibility of the IVW results.

Results

Network Toxicology

A total of 106 target genes of MTX were retrieved through SuperPred, and 20 target genes were retrieved by Swiss Target Prediction and STITCH. 9982, 149, and 85 disease genes for COPD were retrieved from Gene Cards, OMIM, and TTD, respectively. Finally, a Venn diagram of the intersection between disease genes and MTX target genes was drawn (Figure 2A). According to the aforementioned method, 95 core target genes in the intersection were screened out, and a protein-protein interaction (PPI) network containing 95 nodes and 358 edges was constructed, with an average node degree of 7.54. In the network, the greater the color difference between nodes and edges and the darker the node color, the stronger the interaction between proteins (Figure 2B). Subsequently, Cytoscape 3.10.3 software was used to further screen the core genes (Figure 2C) GO enrichment analysis and KEGG enrichment analysis were performed on the core target genes obtained through screening. In the GO enrichment analysis, the top 10 entries with the highest enrichment significance in each of the three categories of biological process (BP), cellular component (CC), and molecular function (MF) were selected to draw barplot charts. In the KEGG enrichment analysis, the top 30 entries with the highest enrichment significance were selected to draw bubble charts and lollipop charts (Figure 2D–F).

Figure 2 (A) Venn diagram of potential targets of MTX and COPD. (B) The PPI network of the potential action targets of MTX-induced COPD. (C) Hub genes analysis of MTX affecting COPD. (D) GO enrichment analysis of the core genes of MTX and COPD. (E and F) KEGG enrichment analysis of the core genes of MTX and COPD.

The results of molecular docking analysis indicate that MTX has binding potential (ΔG < −6.4 kcal/mol) with the proteins encoded by the following 9 genes: HIF1A (ΔG = −7.0 kcal/mol), EP300 (ΔG = −8.3 kcal/mol), NFKB1 (ΔG = −6.4 kcal/mol), PPARG (ΔG = −7.1 kcal/mol), ABL1 (ΔG = −7.7 kcal/mol), STAT1 (ΔG = −7.3 kcal/mol), HDAC1 (ΔG = −7.0 kcal/mol), NFE2L2 (ΔG = −6.8 kcal/mol), and MAPK1 (ΔG = −7.4 kcal/mol). Among them, EP300 shows the strongest binding energy (ΔG = −8.3 kcal/mol). The visualization results of the relevant docking conformations are shown in Figure 3.

Figure 3 The results of molecular docking of HIF1A, EP300, NFKB1, PPARG, ABL1, STAT1, HDAC1, NFE2L2 and MAPK1 with MTX are mainly shown. The yellow represents the ligand, also known as the MTX, the green represents the amino acid that binds to it, and the dotted line represents the hydrogen bond.

MR Analysis

As shown in Table 2: IVW (OR = 1.10E+00, 95% CI = 1.04E+00-1.17E+00), MR-Egger (OR = 1.09E+00, 95% Cl = 9.81E-01-1.21E+00); Weight Mode (OR = 1.10E+00, 95% Cl = 1.03E+00-1.17E+00); Simple Mode (OR = 1.10E+00, 95% Cl = 9.89E-01-1.22E+00); Weighted median (OR = 1.10E+00, 95% Cl = 1.03E+00-1.17E+00). IVW, Weight Mode, and Weighted median indicate that there is a statistically significant positive causal relationship between MTX and COPD, that is, the use of MTX increases the risk of developing COPD. We noted that the result of MR-Egger (P-value = 0.2509) was not significant, but its OR value was extremely close to that of IVW, indicating that no significant directional horizontal pleiotropy was detected. Figures 4–6 present the forest plot of the effect size and 95% confidence interval obtained through MR analysis, the scatter plot of the leave-one-out analysis results, and the scatter plot of the heterogeneity analysis. The results of the leave-one-out analysis showed that after removing each SNP one by one, the results still did not change significantly, indicating that our study is universal and robust. The heterogeneity analysis result of IVW (Q = 4.54E-01, P = 9.29E-01) indicates that there is no significant difference and the result is reliable. The MR-Egger intercept analysis (Q = 3.15E-01, P = 8.54E-01) did not find genetic pleiotropy.

Table 2 Univariate MR Analysis of Methotrexate (MTX) and Chronic Obstructive Pulmonary Disease (COPD)*

Figure 4 Scatterplot (results of MR analysis). Abbreviations: MR, Mendelian randomization; SNP, single-nucleotide polymorphism.

Figure 5 Sensitivity analysis (MTX and COPD). The red line represents the results of the sensitivity analysis of the entire dataset, and the black line represents the analysis results after excluding the target SNP. Abbreviation: MR, Mendelian randomization.

Figure 6 Forest plot (MTX and COPD). The red line represents the analysis results of the inverse-variance weighted method (IVW) and MR-Egger regression for the entire dataset; the black line represents the IVW analysis results after excluding the target SNP. Abbreviation: MR, Mendelian randomization.

Discussion

To explore the association between MTX and COPD, this study first conducted computer simulation prediction based on network toxicology, and then used the MR analysis method to integrate the OpenGWAS datasets of exposure factors and outcomes. By parallelly analyzing the results obtained from the two methods, the consistency of their effect directions was mainly evaluated.

The research results show that there is a positive correlation between MTX and COPD, which means that the application of MTX may increase the risk of developing COPD. The specific analysis methods and results are as follows: IVW (OR = 1.10E+00, 95% CI = 1.04E+00-1.17E+00, P = 0.0008), MR-Egger (OR = 1.09E+00, 95%Cl = 9.81E-01-1.21E+00, P = 0.2509); Weighted median (OR = 1.10E+00, 95%Cl = 1.03E+00-1.17E+00, P = 0.0052); Weight Mode (OR = 1.10E+00, 95%Cl = 1.03E+00-1.17E+00, P = 0.0658); Simple Mode (OR = 1.10E+00, 95%Cl = 9.89E-01-1.22E+00, P = 0.1784). The point estimates obtained by five methods are in the same direction (OR ≈ 1.10). Among them, the associations of the IVW and weighted median methods reach the level of statistical significance (P = 0.0008 and P = 0.0052), while the confidence intervals of MR - Egger, Weight mode, and Simple mode cross the null line, showing no statistical significance. This inconsistency suggests the possible presence of residual pleiotropy or insufficient statistical power, but the consistency in the direction of the effect still provides preliminary genetic clues for this association. In addition, the prediction direction of network toxicology is consistent with the results of MR, which further enhances the reliability of this study.

In the enrichment analysis, the functions in which the proteins encoded by the relevant core genes are involved have been confirmed in existing studies, such as neutrophil extracellular trap formation (NETosis). A study found that an increase in the level of NETosis can lead to an increase in the soluble interleukin-6 receptor (sIL-6R), which in turn induces COPD.²²

Multiple genes in the network core genes are significantly associated with the p38 MAPK signaling pathway. Among them, the correlations of NFKB1, STAT1, HDAC1, HDAC2, NFE2L2, and MAPK1 are particularly prominent. Studies have shown that the activation of the p38 MAPK signaling pathway can significantly inhibit the activity and expression of HDAC1 and HDAC2, leading to impaired function of the glucocorticoid receptor, thereby reducing the clinical treatment effect of COPD.²³ In addition, p38 MAPK can induce the overexpression of inflammatory factors by activating the STAT1 transcription factor, ultimately aggravating lung injury.²⁴

It is worth noting that although MTX can treat RA, during this process, MTX differentially upregulates the expression levels of pro-inflammatory cytokines IL-1β, IL-2, IL-8, IL-10, IFN-γ, and RANTES in the lungs, and at the same time downregulates the expression of anti-inflammatory factors IL-4, IL-6, IL-12, MIP-1α, and MIP-1β in a dose-dependent manner.²⁵ The above cytokines are all closely related to the p38 MAPK signaling pathway, and MTX can directly activate the p38 MAPK signaling pathway. Therefore, reducing the use of MTX can not only directly inhibit the p38 MAPK signaling pathway but also indirectly reverse the corticosteroid insensitivity of IL-8 in airway smooth muscle cells, thereby enhancing the therapeutic effect of COPD.^25,26 In addition, p38 MAPK is a key mediator of TNF-induced NF-kB-dependent gene activation. Therefore, the p38 MAPK signaling pathway and the NF-kB signaling pathway can be sequentially activated, thereby amplifying the expression of pro-inflammatory cytokines (TNF-α, IL-1β, IL-6) and chemokines, and ultimately exacerbating lung inflammation.²⁴

In summary, MTX can increase the risk of developing COPD through multiple mechanisms such as activating the p38 MAPK signaling pathway, which poses new challenges for clinicians in formulating prevention and treatment strategies for COPD. Therefore, when formulating prevention and treatment strategies for COPD, the potential risks of MTX must be fully considered to provide more effective and safe treatment options for patients.

It should be clarified that different doses of MTX can trigger different types and severity of toxic reactions. High-dose MTX can cause acute kidney injury (AKI) in 2–12% of patients, further exacerbating renal function impairment and inducing non-renal adverse events such as bone marrow suppression, hepatotoxicity, and skin lesions.²⁷ Low-dose MTX is mainly associated with liver diseases such as liver cirrhosis, gastrointestinal symptoms, mucocutaneous lesions, and pulmonary fibrosis.^11,12 Among them, pulmonary fibrosis is a very common pathological change in the development process of COPD, which can further aggravate patients’ respiratory restriction.²⁸ Therefore, from the perspectives of the safety and effectiveness of MTX, factors such as dosage, administration mode, and treatment regimen are all core variables that need to be evaluated. Based on a large number of existing studies and the results of this study, we hypothesize that the effect of MTX on COPD is also closely related to the dosage, action pathway, and treatment regimen of MTX. Reasonable regulation of the use of MTX has certain clinical significance in preventing COPD, but the potential mechanism of action between MTX and COPD still needs further exploration.

This study still has certain limitations: At the statistical methodology level, the core statistical method for MR analysis is IVW, and the OR values obtained from the five analysis methods all show a high degree of consistency. However, the P - values of MR Egger, Simple mode, and Weight mode exceed the statistical significance level (P < 0.05), which still weakens the robustness of causal inference to some extent.

At the technical verification level, network toxicology conducts predictions based on computers, so there may be defects such as database dependence, data bias, and the static model ignoring dynamic conformations. Static computer models cannot capture the dynamic conformational changes of proteins and the complete binding patterns under physiological conditions. More than 50% of the static prediction models change their accuracy categories after dynamic evaluation.²⁹

At the level of population representativeness, the sample data related to COPD and MTX in GWAS are mainly from European-descent populations. This leads to differences in allele frequencies, linkage disequilibrium patterns, and gene - environment interactions among populations, limiting the generalizability of the results. Therefore, the conclusions of this study should not be directly extrapolated to other ethnic groups.

The following three points need special clarification: First, in this study, after calibration and verification, the P - value obtained by IVW was 0.0008, which was significantly lower than the threshold. Therefore, the preset global multiple - test correction was not performed. Second, we have systematically excluded SNPs related to known COPD risk factors using PhenoScanner. However, residual confounding where the instrumental variables affect COPD through non - MTX target pathways may still exist. Third, although we adopted sensitivity analysis, there is still a possibility of reverse causation in the MR analysis results. This is not a defect of the design itself but stems from the potential reverse influence of the outcome on the exposure GWAS estimation stage.

This study provides preliminary observational genetic evidence, which should be regarded as a basis for hypothesis generation. It should be emphasized that this study can only provide an estimated value of the hypothesized causal relationship, and further real-world research is needed to confirm the causal association between MTX and COPD. Future research should focus on three aspects: First, integrate dynamic structural biology methods and multi - ancestry population data to validate and expand the current results; Second, adopt a pharmacoepidemiological design based on individual - level dose data. By constructing a dose - stratified model and a time - effect analysis framework, quantify the dose - response relationship between MTX exposure dose and the risk of COPD, providing crucial guidance for clinical drug - use decisions; Third, conduct cell and animal experiments to verify the predicted binding relationships between MTX and NFKB1, HDAC1, MAPK1, etc, and evaluate the impact of different doses of MTX on the incidence of COPD in different populations through clinical studies.

Conclusion

This study adopted a framework combining network toxicology and MR, This study adopted a combined method of network toxicology and MR, and preliminarily provided evidence from a genetic perspective for the hypothesis that there may be a causal relationship between the use of MTX and the increased susceptibility to COPD. The results of the two methods were highly consistent, which not only verified the validity of the conclusion but also provided valuable insights for clinicians to guide the management and prevention strategies of potential COPD patients. However, due to the small observed effect size and the fact that this inference is still based on observational studies and lacks clinical validation, this finding cannot be regarded as conclusive. This study only laid a statistical basis for further exploring the association mechanism between MTX and COPD.