Network toxicology integrated with machine learning and SHAP analysis identifies overlapping immune signatures between Di(2-ethylhexyl) phthalate (DEHP) and Sjögren's syndrome.

Network toxicology integrated with machine learning and SHAP analysis identifies overlapping immune signatures between Di(2-ethylhexyl) phthalate (DEHP) and Sjögren’s syndrome - PMC Skip to main content An official website of the United States government Here's how you know Here's how you know Official websites use .gov A .gov website belongs to an official government organization in the United States. Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites. Search Log in Dashboard Publications Account settings Log out Search… Search NCBI Primary site navigation Search Logged in as: Dashboard Publications Account settings Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide PERMALINK Copy As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice BMC Pharmacol Toxicol . 2026 Mar 7;27:65. doi: 10.1186/s40360-026-01119-x Search in PMC Search in PubMed View in NLM Catalog Add to search Network toxicology integrated with machine learning and SHAP analysis identifies overlapping immune signatures between Di(2-ethylhexyl) phthalate (DEHP) and Sjögren’s syndrome Cheng Lili Cheng Lili 1 The First Affiliated Hospital of Anhui University of Chinese Medicine, Hefei, Anhui 230012 China Find articles by Cheng Lili 1 , Tang Zhongfu Tang Zhongfu 1 The First Affiliated Hospital of Anhui University of Chinese Medicine, Hefei, Anhui 230012 China Find articles by Tang Zhongfu 1 , Li Ming Li Ming 1 The First Affiliated Hospital of Anhui University of Chinese Medicine, Hefei, Anhui 230012 China 2 Institute of Xin’an Medicine and Modernization of Traditional Chinese Medicine, Hefei, Anhui 230012 China Find articles by Li Ming 1, 2 , Chuanbing Huang Chuanbing Huang 1 The First Affiliated Hospital of Anhui University of Chinese Medicine, Hefei, Anhui 230012 China 2 Institute of Xin’an Medicine and Modernization of Traditional Chinese Medicine, Hefei, Anhui 230012 China Find articles by Chuanbing Huang 1, 2, ✉ Author information Article notes Copyright and License information 1 The First Affiliated Hospital of Anhui University of Chinese Medicine, Hefei, Anhui 230012 China 2 Institute of Xin’an Medicine and Modernization of Traditional Chinese Medicine, Hefei, Anhui 230012 China ✉ Corresponding author. Received 2025 Oct 25; Accepted 2026 Feb 27; Collection date 2026. © The Author(s) 2026 Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/ . PMC Copyright notice PMCID: PMC13081303  PMID: 41794809 Abstract Background As a prevalent plasticizer in industrial manufacturing and daily-use products, di(2-ethylhexyl) phthalate (DEHP) has raised growing concerns regarding its safety and potential role in disease pathogenesis. Although direct epidemiological evidence linking DEHP exposure to Sjögren’s syndrome (SS) remains limited, DEHP has been reported to exert immunomodulatory and endocrine-disrupting effects that are highly relevant to the core pathological mechanisms of SS. Therefore, this study aimed to systematically explore the potential links between DEHP and SS using an integrated strategy combining toxicity prediction, network toxicology, machine learning, molecular docking, single-cell atlas interrogation, and in vitro validation. Methods The toxicity profile of DEHP was predicted using two complementary online platforms (ProTox 3.0 and ADMETlab 2.0). SS-related transcriptomic datasets were obtained from GEO, and potential DEHP targets were collected from CTD, SwissTargetPrediction, and SEA. Network toxicology analyses (intersection screening, PPI construction, and GO/KEGG enrichment) were performed to identify key genes and pathways. Machine learning models were trained to classify SS versus controls using the DEHP–SS intersect genes, and SHAP was applied to interpret key predictors. Molecular docking (AutoDock Vina) was conducted to estimate the binding propensity between DEHP and prioritized proteins. A single-cell atlas of mouse submandibular gland tissue (PanglaoDB; Mus musculus; SRA693675:SRS3206192) was queried to localize gene expression. Finally, in vitro experiments were performed in human submandibular gland epithelial cells (HSG) to validate STAT1 expression changes following DEHP exposure. Results In silico toxicity predictions suggested that DEHP may exhibit immunotoxic potential among multiple predicted toxicity endpoints. Network toxicology analysis identified 27 DEHP–SS intersect targets, and the PPI network highlighted STAT1, CCL2, CXCL10, MYD88, and IRF7 as highly connected hub genes. Across 113 machine learning models, Ridge/Elastic Net–based models demonstrated stable classification performance, and SHAP analysis prioritized five core predictors (STAT1, PHGDH, ISG15, CXCL10, and CCL2). Molecular docking suggested favorable binding energies ( ≤ − 5.0 kcal/mol) between DEHP and these five proteins. Single-cell atlas interrogation indicated that STAT1 is expressed in salivary mucous cells. In vitro experiments further demonstrated that DEHP exposure increases STAT1 mRNA and protein expression in HSG cells in a dose-dependent manner. Collectively, these results provide mechanistically plausible—yet non-causal—evidence linking DEHP to SS-relevant immune signatures. Conclusion This study identifies an overlapping immune-response signature between predicted DEHP-associated targets and SS-related transcriptomic features, highlighting STAT1 as a prioritized node within this intersection. In vitro experiments demonstrate that DEHP exposure upregulates STAT1 expression in HSG cells. These findings provide hypothesis-generating evidence that DEHP may modulate interferon-related immune pathways relevant to SS, while not establishing causality. Supplementary Information The online version contains supplementary material available at 10.1186/s40360-026-01119-x. Keywords: Di(2-ethylhexyl) phthalate, Sjögren's syndrome, Network toxicology, SHAP analysis, Machine learning Introduction Sjögren’s syndrome (SS) constitutes a widespread systemic autoimmune disorder distinguished by immune-mediated lymphocytic invasion of exocrine glandular tissues, culminating in xerostomia (oral desiccation) and keratoconjunctivitis sicca (ocular surface inflammation). This condition represents a pressing public health challenge given its profound impact on patient quality of life [ 2 ]. Epidemiological data indicate an average annual incidence of 6.0 per 100,000 individuals, with a striking ten-fold female predominance and age-associated incidence escalation [ 44 ]. The etiological landscape of SS encompasses intricate interactions among genetic predisposition, immune system dysregulation, and environmental provocations, including estrogen fluctuations, viral pathogen exposure, and environmental pollutant encounters [ 6 , 50 ]. Importantly, while environmental factors are increasingly recognized as potential triggers for SS, specific links between common plasticizers and SS remain insufficiently characterized and warrant systematic investigation. Di(2-ethylhexyl) phthalate (DEHP), functioning as a predominant plasticizing agent, ranks among the highest-volume phthalate compounds manufactured globally, finding extensive incorporation across diverse plastic formulations [ 49 ]. Human DEHP exposure occurs ubiquitously through household items, medical devices, consumer products, and food packaging materials, and its environmental detection in food matrices, air, industrial effluents, soils, and aquatic systems underscores its broad footprint [ 16 ]. Experimental and biomonitoring evidence indicates that DEHP and its metabolites can disrupt endocrine signaling and immune homeostasis, thereby potentially amplifying autoimmune-prone responses [ 32 , 34 , 51 ]. Nevertheless, direct epidemiological or causal evidence connecting DEHP exposure to SS onset is currently lacking. Accordingly, the present work is designed as a hypothesis-generating study that leverages multi-omics integration and in vitro validation to evaluate whether DEHP-targeted molecular networks overlap with SS-related transcriptomic signatures and to prioritize testable mechanistic candidates. ProTox 3.0 represents an online toxicity prediction platform employing molecular similarity and machine learning models to systematically assess 61 toxicity endpoints, including acute toxicity, organ-specific toxicity, and adverse outcome pathways [ 4 ]. ADMETlab 2.0, conversely, serves as an integrated tool for predicting compound absorption, distribution, metabolism, excretion, and toxicity (ADMET) profiles [ 48 ]. To comprehensively evaluate DEHP’s potential health impacts, this study employed both platforms, leveraging their distinct algorithmic approaches to characterize toxicological manifestations and absorption characteristics, thereby enhancing prediction interpretability. Furthermore, we applied a network toxicology framework, which synthesizes systems biology, bioinformatics, and computational toxicology to systematically interpret chemical-associated biological networks [ 19 ]. By constructing multi-dimensional compound–target–pathway networks, we aimed to elucidate key targets and core biological pathways through which DEHP may plausibly modulate SS-relevant mechanisms. Machine learning (ML) refers to computational methods that automatically identify patterns from data to build predictive models. ML algorithms are particularly adept at handling high-dimensional, non-linear complex relationships, uncovering underlying patterns that may be missed by traditional statistical analyses, and have been increasingly applied to biomarker discovery and disease classification [ 1 ]. In this study, ML techniques were employed to prioritize key genes best classify SS versus controls in independent transcriptomic cohorts. SHAP analysis was subsequently used to interpret the ML model’s decision-making process by quantifying the contribution of each gene to the final prediction, thereby improving interpretability [ 35 ]. Molecular docking was then performed to estimate the likelihood of DEHP binding to the prioritized protein targets, and single-cell RNA sequencing (scRNA-seq) atlases were queried to provide cell-type localization information for key genes [ 9 ]. In summary, this study aims to clarify the predicted toxicological characteristics of DEHP and to identify SS-relevant molecular targets and pathways potentially modulated by DEHP using an integrated framework that includes network toxicology, ML/SHAP interpretability, molecular docking, single-cell atlas interrogation, and in vitro experiments. By emphasizing biological plausibility and transparent methodological limitations, our results are intended to generate testable mechanistic hypotheses rather than establish a causal relationship between DEHP exposure and SS. Materials and methods Study design Our investigative framework initiated with computational toxicological profiling of DEHP via specialized web-based platforms. Subsequently, SS-associated datasets were retrieved from the NCBI Gene Expression Omnibus (GEO) database. Following data preprocessing, differentially expressed genes (DEGs) associated with SS were identified, and potential SS-related targets were further screened using Weighted Gene Co-expression Network Analysis (WGCNA). DEHP targets were retrieved through network toxicology approaches. The intersection between these two gene sets—SS-related targets and DEHP targets—was computed in R (e.g., ggvenn) and visualized. Subsequent functional enrichment analysis of overlapping genes was conducted within the R software environment utilizing specific packages for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses. Protein-protein interaction (PPI) networks for overlapping genes were constructed using the STRING database and subsequently imported into Cytoscape software for visualization and advanced topological analysis. The top five hub genes and subnetworks were identified employing the CytoHubba plugin within Cytoscape. To refine hub gene identification and pinpoint diagnostic biomarkers, machine learning (ML) was employed, followed by SHAP analysis to interpret model predictions and assess feature importance. Molecular docking simulations were then conducted to evaluate binding stability between hub genes (as target proteins) and DEHP. Additionally, single-cell RNA sequencing (scRNA-seq) data analysis was performed to validate expression patterns of hub genes at cellular resolution. Finally, in vitro cell experiments were performed to validate the expression of core genes in human submandibular gland epithelial cells (HSG) following exposure to different concentrations of DEHP. In summary, this integrated methodology was designed to predict DEHP’s pathogenic risk and its molecular targets in SS [ 18 , 31 ]. A schematic flowchart detailing specific procedures is presented in Fig. 1 . Fig. 1. Open in a new tab Research design process Toxicity analysis and prediction To thoroughly characterize the toxicity profile of DEHP through computational simulations, we employed two online platforms in a complementary manner: ProTox 3.0 ( https://tox.charite.de ) and ADMETlab 2.0 ( https://admetmesh.scbdd.com ). The canonical SMILES and 2D structure file (SDF) of DEHP were retrieved from PubChem( https://pubchem.ncbi.nlm.nih.gov )(Kim et al., 2023). The canonical SMILES string (CCCCC(CC)COC(= O)c1ccccc1C(= O)OCC(CC)CCCC) was used as the primary input for both webservers (default settings), whereas the 2D structure was used only for figure illustration. ProTox 3.0 was applied to predict toxicity endpoints (including organ toxicity and immunotoxicity) and to output endpoint-specific probabilities; ADMETlab 2.0 was used to obtain predicted ADMET properties and toxicity-related risks, which are reported as semi-quantitative confidence levels (e.g., “---” to “+++”) in the platform output [ 4 , 48 ]. Acquisition of SS-related targets Five SS-related datasets ( GSE23117 , GSE40611 , GSE51092 , GSE66795 , and GSE84844 ) were curated from the NCBI GEO database(Clough et al., 2024). GSE23117 (minor salivary gland) and GSE40611 (parotid gland) represent exocrine-gland tissues, whereas GSE51092 , GSE66795 , and GSE84844 are blood-based cohorts. GSE23117 , GSE40611 , and GSE51092 were designated as the discovery cohort, and GSE66795 and GSE84844 served as independent external validation cohorts. For each dataset, we downloaded the (pre)processed expression matrices and performed probe-to-gene symbol mapping using the corresponding platform annotation; when multiple probes mapped to the same gene, the average expression value was used. Expression values were log2-transformed when necessary and quantile-normalized within each dataset to improve comparability across platforms. Given that the discovery cohort includes different biological matrices (salivary gland versus blood), we treated both dataset identity and tissue type as potential sources of unwanted variation. Batch effect correction was performed in a multi-stage manner: (1) surrogate variable analysis (SVA) was used to estimate latent technical factors within the discovery cohort [ 29 ]; (2) residual batch variation was further adjusted using ComBat with an empirical Bayes framework [ 21 ], with tissue type included as a covariate. Post-correction assessments (boxplots and PCA) were used to evaluate harmonization. Because SS exhibits strong sex bias and the validation cohorts are predominantly female, we explicitly acknowledge sex-related sampling limitations and discuss their potential impact in the Discussion section. Differential gene expression analysis Differential expression analysis was performed using the limma package [ 37 ]. For the discovery cohort, a linear model was fitted with SS status as the primary factor of interest while including dataset origin and tissue type (salivary gland vs. blood) as covariates to mitigate cross-matrix confounding. Genes with an adjusted p-value < 0.05 and |log2FC| > 0.585 were considered significantly differentially expressed. Visualization included volcano plots and heatmaps generated with standard R packages (ggplot2 and pheatmap). Weighted gene co-expression network analysis (WGCNA) A scale-free co-expression network was constructed using the WGCNA package [ 27 ]. To reduce potential confounding from dataset origin and tissue type, WGCNA was performed on the batch-corrected expression matrix after regressing out dataset and tissue covariates (linear model residuals), while preserving SS status-associated variation. Briefly, (1) a merged normalized expression matrix was imported, replicate probes were collapsed with avereps, and low-variance genes were filtered (SD > 0.5); (2) sample groups were parsed to build the trait matrix; (3) missing/aberrant genes and samples were removed using goodSamplesGenes, and outliers were excluded by sample clustering with a static cut (cutHeight = 20000; minSize = 10); (4) the soft-thresholding power was selected by pickSoftThreshold using the scale-free topology criterion (signed R² ≈ 0.80) and applied to construct the adjacency matrix; (5) TOM and TOM dissimilarity (1 − TOM) were computed, genes were hierarchically clustered, and modules were identified by cutreeDynamic (deepSplit = 2; pamRespectsDendro = FALSE; minModuleSize = 60) and merged by mergeCloseModules (cutHeight = 0.25); (6) module–trait associations were assessed by correlating module eigengenes with Control/Treat traits; and (7) calculated and exported to rank candidate hub genes within trait-associated modules. Acquisition of DEHP targets To systematically identify the molecular targets of DEHP, we integrated multiple curated resources, including the Comparative Toxicogenomics Database (CTD, http://ctdbase.org ) [ 13 ], SwissTargetPrediction ( http://www.swisstargetprediction.ch ) [ 12 ], and the Similarity Ensemble Approach (SEA, https://sea.bkslab.org ) [ 23 , 45 ]. All resources were accessed in our analysis window (July 2025) and queried using DEHP identifiers (name and canonical SMILES). CTD provides curated chemical–gene/protein interactions with evidence annotations; SwissTargetPrediction and SEA are ligand-similarity-based target prediction tools. For SwissTargetPrediction, the canonical SMILES of DEHP was submitted with organism set to Homo sapiens and default parameters, predicted targets were exported and retained. For SEA, the canonical SMILES was submitted under default settings with human targets retained when available. For CTD, DEHP-associated genes/proteins were extracted from curated chemical–gene interaction records restricted to Homo sapiens, and a minimum interaction evidence threshold (interaction count ≥ 5) was applied to reduce noise. Finally, they were uniformly mapped to the official gene symbols in Uniprot ( https://www.uniprot.org/ ) and merged after removing duplicates. Screening of DEHP-SS core targets and protein-protein interaction (PPI) network construction A union set of DEGs and the genes from the key WGCNA module(s) was compiled to represent SS-related targets. Core DEHP–SS targets were identified by intersecting this SS target set with the predicted DEHP target list, and overlaps were computed entirely in R using the ggvenn package. Because predicted DEHP targets are proteins whereas transcriptomic features are mRNAs, we integrated the two layers by mapping proteins to their corresponding official gene symbols. Overlapping targets were submitted to the STRING database to construct a protein–protein interaction (PPI) network [ 41 ], using Homo sapiens as the organism and a standard confidence score threshold. Network files were exported and imported into Cytoscape (v3.9.1) for visualization [ 40 ]. Hub genes were ranked using the cytoHubba plugin in Cytoscape, and the top five genes were prioritized based on degree centrality [ 10 ]. Functional and pathway enrichment analysis Gene Ontology (GO) enrichment analysis (Biological Process, Cellular Component, and Molecular Function) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway [ 22 , 42 ] enrichment were performed using clusterProfiler in R (clusterProfiler v4.16.0 [ 47 ], together with enrichplot(v1.28.4) and org.Hs.eg.db༈V3.21.0༉. Statistical significance was evaluated using Benjamini–Hochberg adjusted p-values (p-adjusted/FDR) and p-value < 0.05 considered significant. Enrichment results were visualized as bubble plots and chord/circle plots to facilitate interpretation. Machine learning-based screening of core genes To prioritize robust transcriptomic biomarkers among the DEHP–SS intersect targets, we established a supervised machine learning (ML) framework in which the outcome variable was SS status (SS vs. control). The intersect gene expression matrix was extracted from the batch-corrected discovery cohort and split into a training set and an internal validation set using stratified sampling (typical split: ~70% training/30% validation). Model fitting and internal tuning were performed strictly within the training set to avoid data leakage: Elastic Net/Lasso/Ridge used 10-fold cross-validation in cv.glmnet (nfolds = 10) to select lambda; GBM used 10-fold CV (cv.folds = 10) to choose the optimal number of trees; glmBoost used k-fold resampling via cvrisk; LDA was trained using caret with CV (default); and XGBoost used 5-fold CV to determine nrounds. Two independent GEO cohorts ( GSE66795 and GSE84844 ) were used for external validation. A total of 11 algorithm families were evaluated: Lasso, Ridge, Elastic Net, Support Vector Machine (SVM), Random Forest (RF), glmBoost, Stepglm, Gradient Boosting Machine (GBM), Linear Discriminant Analysis (LDA), XGBoost, and Naive Bayes. The total number of 113 models reflects different hyperparameter configurations within each family (e.g., multiple Elastic Net mixing parameters, alpha = 0.1–0.9) and predefined two-stage pipelines (feature selection method + classifier); the full list is provided in Supplyment. refer.methodLists. Because several algorithms are sensitive to feature scaling, gene expression features were standardized (z-score) prior to model training; as implemented in our scripts, standardization was performed within each cohort (training and each validation cohort) using cohort-specific centering and scaling parameters. Model discrimination was primarily assessed using AUC, which was computed from continuous predicted probabilities. For threshold-dependent metrics (accuracy and F1-score), class labels were obtained using the model-specific decision rule implemented in PredictClass; for methods requiring an explicit probability cutoff (e.g., logistic regression/Stepglm, GBM, XGBoost, and plsRglm), we used a fixed threshold of 0.5 (probability > 0.5 classified as SS). Predicted-probability heatmaps were visualized using the ComplexHeatmap package. Model interpretation Given the inherent “black-box” nature of complex ML models, we applied the SHapley Additive exPlanations (SHAP) algorithm to quantify the contribution of each feature to the model predictions. This method assigns a SHAP value to each feature, enabling an interpretable assessment of its impact on the model output [ 14 ]. Molecular docking analysis Molecular docking simulations were conducted to estimate the binding propensity between DEHP and the prioritized protein targets. The 3D structure of DEHP was retrieved from PubChem in SDF format [ 24 ] and converted to the required input format after geometry minimization using the MMFF94 force field, which is well suited for small organic molecules (including phthalates) and provides chemically reasonable ligand conformations for docking. Protein structures were obtained from the RCSB Protein Data Bank (PDB), with UniProt used for target annotation and PDB cross-referencing [ 43 ]. Proteins were prepared by removing crystallographic waters/heteroatoms as appropriate, adding polar hydrogens, and assigning charges. Docking was performed using AutoDock Vina (v1.2.0 [ 15 ].Docking poses and interactions were visualized using PyMOL (The PyMOL Molecular Graphics System; Schrödinger, LLC) [ 39 ]. Single-cell RNA sequencing analysis To provide cell-type localization information for prioritized genes in salivary gland tissue, we leveraged the PanglaoDB single-cell transcriptomics database [ 17 ]. Using the dataset labeled “Submandibular Gland (SRA693675:SRS3206192)” in PanglaoDB (Mus musculus), we queried target genes and extracted their expression patterns across annotated cell types. This step was used as supportive evidence to assess whether prioritized genes are expressed in relevant glandular cell populations (e.g., epithelial/salivary mucous cells), rather than to infer disease-specific differential expression. Cell experiment validation To experimentally validate the computationally prioritized target, we focused on STAT1 because it was consistently highlighted as (i) a highly connected hub gene in the PPI network and (ii) a top predictor in the ML/SHAP analyses, and it was also detectable in salivary gland single-cell atlases. therefore, STAT1 was selected as the primary in vitro validation target. In vitro experiments were performed to assess whether DEHP exposure modulates STAT1 expression in human submandibular gland epithelial cells (HSG). Cell culture and DEHP treatment Human submandibular gland epithelial cells (HSG, purchased from Shanghai Cyagen Biosciences, HUM-iCell-g023) were cultured in DMEM/F12 medium supplemented with 10% fetal bovine serum. Cells were exposed to DEHP (Sigma-Aldrich, purity > 99%) at concentrations of 0.1 µM, 1 µM, and 10 µM for 24 h. The solvent control group was treated with an equal volume of DMSO (final concentration < 0.1%). All experiments were performed with at least three independent biological replicates ( n ≥ 3) per condition. STAT1 expression under different DEHP treatments was assessed by RT-qPCR, Western blotting, and immunofluorescence staining. PCR detection Total RNA was extracted using the TRIzol method and reverse-transcribed into cDNA, which was then analyzed by RT-qPCR. Primer sequences are listed in Table 1 . Each sample was measured in technical triplicates, and gene expression levels were calculated using the 2^–ΔΔCt method, with β-actin as the internal reference gene. Table 1. Primer sequences Gene Amplicon size (bp) Forward primer (5’→3’) Reverse primer (5’→3’) β-actin 96 CCCTGGAGAAGAGCTACGAG GGAAGGAAGGCTGGAAGAGT STAT1 198 CTGTGAAGTTGAGAGATGTGA CAGTAACGATGAGAGGACCC Open in a new tab Protein extraction and western blot analysis Total proteins were extracted with RIPA lysis buffer and quantified by the BCA method. Equal amounts of proteins were separated by SDS-PAGE, transferred to PVDF membranes, and blocked. Membranes were incubated with primary antibodies against STAT1 (1:5000, Proteintech, 10144-2-AP) and GAPDH (1:5000, Zsbio, TA-08), followed by HRP-conjugated secondary antibodies. Signals were developed using an ECL chemiluminescence system, and ImageJ was used for densitometric quantification. To address concerns about image integrity, uncropped full-length blots and molecular weight markers are provided in the Supplementary Materials. Immunofluorescence staining After cells were seeded on coverslips, they were washed twice with PBS and fixed with 4% paraformaldehyde for 15 min. Permeabilization and blocking solution was added for incubation at room temperature for 30 min. Diluted STAT1 primary antibody (1:500, Bioss, bs-1317R) was added and incubated overnight at 4 °C. Fluorescent secondary antibody was incubated at room temperature for 1 h, followed by DAPI nuclear staining for 5 min. Coverslips were mounted with anti-fade medium and imaged using a fluorescence microscope. Representative images are shown from at least three independent experiments. Results Prediction of chemical toxicity using PROTox and ADMETlab ProTox 3.0 predicted that DEHP may be active in multiple organ toxicity endpoints (e.g., hepatotoxicity, neurotoxicity, and respiratory toxicity) and also suggested potential immunotoxicity and ecotoxicity (Table 2 ). ADMETlab 2.0 further predicted toxicity-related risks including skin sensitization, carcinogenicity, and eye irritation (Table 3 ). Considering that (i) ProTox models are trained primarily on rodent toxicity data and (ii) the immunotoxicity endpoint has only moderate reported predictive performance, we interpret these results as in silico hazard signals rather than definitive evidence of human immunotoxicity. Based on these preliminary signals, we proceeded to network toxicology analyses to examine whether DEHP-associated molecular targets overlap with SS-related transcriptomic signatures. Table 2. Toxicity prediction results for DEHP from PROTox3.0 Classification Target Prediction Probability Organ toxicity Hepatotoxicity Active 0.69 Organ toxicity Neurotoxicity Active 0.87 Organ toxicity Nephrotoxicity Inactive 0.90 Organ toxicity Respiratory toxicity Active 0.98 Organ toxicity Cardiotoxicity Inactive 0.77 Toxicity end points Carcinogenicity Inactive 0.62 Toxicity end points Immunotoxicity Active 0.96 Toxicity end points Mutagenicity Inactive 0.97 Toxicity end points Cytotoxicity Inactive 0.93 Toxicity end points BBB-barrier Inactive 1 Toxicity end points Ecotoxicity Active 0.73 Toxicity end points Clinical toxicity Inactive 0.56 Toxicity end points Nutritional toxicity Inactive 0.74 Open in a new tab Table 3. Toxicity prediction results of DEHP by ADMETlab 2.0 Classification Probability Labels hERG Blockers -- H-HT --- DILI --- AMES Toxicity --- Rat Oral Acute Toxicity --- FDAMDD --- Skin Sensitization +++ Carcinogencity + Eye Corrosion --- Eye Irritation +++ Respiratory Toxicity -- Open in a new tab Tip: For the classification endpoints, the prediction probability values are transformed into six symbols: 0−0.1(---), 0.1–0.3(--), 0.3–0.5(-), 0.5–0.7(+), 0.7–0.9(++), and 0.9−1.0(+++) Note: Prediction indicates whether the endpoint is predicted to be toxic/active; Probability indicates the model-estimated confidence (0–1). The red color of the label indicates the highest probability, followed by yellow, and green has the lowest probability Identification of SS-related target genes To minimize technical variation, we integrated GSE23117 , GSE40611 , and GSE51092 as the discovery cohort and performed multi-stage batch correction as described in Section “ Acquisition of SS-related targets ”. Boxplots and PCA were used to assess data harmonization before and after correction (Fig. 2 A–D). Notably, because the discovery cohort includes different biological matrices (salivary gland vs. blood), residual clustering by dataset/tissue may persist even after correction; therefore, the integrated analysis is interpreted as capturing shared, systemic immune-related signatures rather than claiming complete cross-tissue equivalence. Fig. 2. Open in a new tab Identification of SS-related genes. ( A ) Box plot of expression values before batch correction. ( B ) Box plot of expression values after batch correction. ( C ) PCA plot before batch correction (axis labels indicate the percentage of variance explained). ( D ) PCA plot after batch correction (axis labels indicate the percentage of variance explained). ( E ) Volcano plot of the differentially expressed genes (DEGs) between SS and controls. ( F ) Heatmap of DEGs. ( G ) Gene dendrogram and co-expression module colors. ( H ) Module–trait relationships. ( I ) Gene significance across modules. ( J ) Venn diagram showing overlap between DEGs and WGCNA module genes Differential expression analysis identified 95 DEGs (59 upregulated and 36 downregulated; FDR < 0.05, |log2FC| > 0.585) between SS and controls (Fig. 2 E–F). WGCNA was then performed on the corrected expression matrix; the soft-thresholding power was set to 5, and eight co-expression modules were identified (Fig. 2 G). Module–trait analysis indicated that the pink module had the strongest association with SS status (Fig. 2 H), and module selection was based on correlation strength, statistical significance, and module-level gene significance rather than module size alone (Fig. 2 I). By taking the union of DEGs and the genes from the key WGCNA module and removing duplicate gene symbols (unique operation), we obtained 381 SS-related candidate genes; among them, 62 genes were shared between the DEG list and the WGCNA module (Fig. 2 J). Network toxicological analysis of SS-DEHP targets Potential targets of DEHP were retrieved from CTD, SwissTargetPrediction, and SEA using the criteria described in Section “ Acquisition of DEHP targets ”. After consolidating targets across resources and removing duplicates, a total of 721 unique DEHP-associated candidate targets were obtained (Fig. 3 A). Intersecting the DEHP target list with the 381 SS-related candidate genes yielded 27 overlapping targets (Fig. 3 B), which were subsequently considered for downstream network analysis.A PPI network of the 27 overlapping SS–DEHP targets was constructed using STRING and visualized in Cytoscape (Fig. 3 C-D). Based on degree ranking, the top5 hub genes were further summarized (Fig. 3 E), respectively STAT1、CCL2、CXCL10、MYD88、IRF7. Fig. 3. Open in a new tab Network toxicology analysis of DEHP–SS intersect targets. ( A ) Venn diagram showing DEHP-associated candidate targets retrieved from CTD, SwissTargetPrediction, and SEA. ( B ) Venn diagram showing the overlap between DEHP-associated targets and SS-related candidate genes, yielding 27 intersecting targets. ( C ) STRING-derived PPI network of the 27 intersecting targets. (edges represent protein–protein associations curated by STRING, and the edge thickness reflects the STRING combined confidence score (higher score indicates stronger supporting evidence). ( D ) PPT network of the 27 intersecting targets imported into Cytoscape (node size and node color were mapped to degree centrality; darker and larger nodes indicate higher degree values).( E ) Hub-gene subnetwork visualization highlighting the top5 ranked GO and KEGG enrichment analysis GO and KEGG enrichment analyses were performed on the 27 overlapping genes (DEHP targets ∩ SS-related transcriptomic candidates) to provide a functional annotation of this intersect signature. Because the input is a pre-selected gene set rather than full transcriptome-wide differential expression, enrichment is interpreted as a descriptive summary of the biological themes represented by the overlap, and not as standalone evidence that DEHP induces SS. The GO analysis revealed that for biological process (BP), the genes were significantly enriched in terms such as “response to molecule of bacterial origin,” “response to lipopolysaccharide,” and “type I interferon signaling pathway.” Regarding cellular component (CC), the genes were primarily localized to the “external side of plasma membrane,” “phagocytic vesicle,” and “endocytic vesicle.” In the molecular function (MF) category, the main enrichments were “protein kinase regulator activity,” “serine/threonine kinase regulator activity,” “Toll-like receptor binding,” and “chemokine receptor binding” (Fig. 4 A).KEGG pathway analysis indicated that the targets were involved in a spectrum of pathways, primarily including viral infections (e.g., Epstein-Barr virus infection, COVID-19, Herpes simplex virus 1 infection), immune signaling pathways (e.g., Chemokine signaling pathway, Toll-like receptor signaling pathway, IL-17 signaling pathway), and metabolic pathways (e.g., Lipid and atherosclerosis, Cholesterol metabolism) (Fig. 4 B-C). Consistent with the STRING PPI topology, many enriched terms converged on interferon-related and chemokine/cytokine signaling processes, suggesting that a subset of the overlap genes participates in interacting immune-response modules. Accordingly, we use enrichment primarily to contextualize the PPI hubs and to guide hypothesis generation and downstream validation. Fig. 4. Open in a new tab GO and KEGG enrichment analysis. ( A ) Bubble plot of GO enrichment analysis across three ontologies: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). The bubble size represents the number of genes (Count) involved in the term, the color indicates the significance level (-log10(adjusted p-value)), and the x-axis (GeneRatio) represents the ratio of significant genes associated with the term. ( B ) Circular plot of the GO enrichment analysis. Terms are color-coded by ontology (BP in green, CC in yellow, MF in purple). The number of genes enriched in each term is displayed, with darker shades indicating higher enrichment significance. ( C ) Bubble plot of KEGG pathway enrichment analysis. The bubble size corresponds to the number of genes in the pathway, and the color represents the enrichment significance (-log10(adjusted p-value)), with redder hues indicating greater significance Identification of core genes associated with DEHP–SS intersection signatures Through a comprehensive machine-learning analysis of the 27 candidate targets, we built a supervised machine-learning (ML) framework to classify Sjögren’s syndrome (SS) versus controls (Type = SS vs. Control). We evaluated 113 pre-defined model configurations spanning 11 algorithm families and multiple hyperparameter settings, and several models achieved robust discrimination in both the internal validation set and independent external cohorts (Fig. 5 A). Among them, Ridge regression showed stable generalization across datasets and was selected for interpretation. SHAP analysis was used to quantify feature contributions, highlighting STAT1, PHGDH, ISG15, CXCL10, and CCL2 as the most influential predictors (Fig. 5 B–E).ROC curves in the validation cohorts showed that STAT1 (AUC = 0.823) and CXCL10 (AUC = 0.800) provided strong discrimination between SS and controls (Fig. 5 F). The differential expression patterns of these core genes in SS were visualized using a volcano plot (Fig. 5 G). Fig. 5. Open in a new tab Identification of core genes for SS classification using DEHP–SS intersect genes. ( A ) Model performance comparison heatmap showing AUC values across cohorts (generated with the R package ComplexHeatmap). ( B ) SHAP feature importance bar plot ranking genes by mean(|SHAP|). ( C ) SHAP summary (violin/beeswarm) plot showing the distribution of SHAP values per gene; point colors reflect normalized feature values. ( D ) SHAP interaction/dependence plots illustrating feature–feature interaction effects on model output. ( E ) SHAP waterfall plot explaining a representative individual prediction by decomposing feature contributions. ( F ) ROC curves for key genes (STAT1, PHGDH, ISG15, CXCL10, and CCL2). ( G ) Volcano plot of DEGs with core genes labeled Importantly, the ML task in this section is disease-status classification (SS vs. control) using transcriptomic expression data; therefore, the AUC values reflect diagnostic discrimination within SS datasets and do not, by themselves, demonstrate an effect of DEHP exposure. We also note that PPI-derived hub genes (e.g., STAT1, CCL2, CXCL10, MYD88, IRF7) and ML/SHAP-prioritized predictors (STAT1, PHGDH, ISG15, CXCL10, CCL2) are not fully overlap, because network topology and predictive importance (contribution to classification across cohorts) capture different properties. Molecular docking validation of DEHP–core gene interactions Molecular docking was used as a supportive, hypothesis-generating step to estimate whether DEHP can adopt energetically favorable poses with a focused set of proteins prioritized by the preceding analyses. Specifically, we docked DEHP to the five ML/SHAP‑prioritized core predictors (STAT1, PHGDH, ISG15, CXCL10, and CCL2).Our results demonstrated that DEHP consistently exhibited strong binding affinity with all five target proteins, with binding energies consistently below − 5 kcal/mol (Table 4 ), suggesting stable molecular interactions. Visualization of the binding conformations (Fig. 6 ) revealed stable docking poses for all DEHP-protein complexes. Table 4. Molecular docking binding energies of DEHP with the target proteins Ligand Receptor Bind. energy [kcal/mol] DEHP STAT1 −6.4 DEHP PHGDH −6.8 DEHP ISG15 −6.1 DEHP CXCL10 −6.8 DEHP CCL2 −5.5 Open in a new tab Fig. 6. Open in a new tab Molecular docking analysis of the interactions between DEHP and the core proteins. ( A ) Docking pose of DEHP with STAT1. ( B ) Docking pose of DEHP with PHGDH. ( C ) Docking pose of DEHP with ISG15. ( D ) Docking pose of DEHP with CXCL10. ( E ) Docking pose of DEHP with CCL2. Docking indicates binding propensity under the scoring function and does not constitute experimental validation Single-cell RNA sequencing elucidates STAT1 expression To provide cell-type localization information for prioritized genes in salivary gland tissue, we queried the PanglaoDB single-cell atlas (SRA693675:SRS3206192). Among the candidate genes, STAT1 showed detectable expression within salivary mucous cells in this atlas (Fig. 7 A–B). Because only one potentially relevant submandibular-gland atlas was available and the dataset is not a Sjögren’s syndrome patient cohort, this analysis is presented as supportive localization evidence rather than disease-specific validation. We also acknowledge that additional SS-relevant markers could be prioritized for future in vitro validation when suitable human salivary gland single-cell resources become available. Fig. 7. Open in a new tab Single-cell atlas localization of STAT1 in submandibular gland tissue. ( A ) t-SNE plot of the PanglaoDB submandibular gland single-cell dataset (SRA693675:SRS3206192; Mus musculus). ( B ) Feature plot showing STAT1 expression enrichment in salivary mucous cells Experimental validation of DEHP-mediated regulation of STAT1 expression Importantly, increased STAT1 expression in this epithelial cell model represents activation of an interferon-related inflammatory signature and does not constitute direct molecular validation of DEHP-mediated Sjögren’s syndrome pathogenesis. To verify whether DEHP regulates the expression of the core target STAT1 as predicted, qPCR, Western Blot, and immunofluorescence staining analyses were performed in HSG cells. As shown in Fig. 8 A and B, Western Blot analysis revealed that compared with the control group, the relative protein expression level of STAT1 in DEHP-treated groups was significantly increased in a concentration-dependent manner. At the 10 µM DEHP treatment, STAT1 mRNA expression was upregulated by approximately 3-fold. qPCR analysis further confirmed that STAT1 mRNA expression was significantly elevated with increasing DEHP concentrations (Fig. 8 C), consistent with the trend of protein level changes. Specifically, STAT1 mRNA expression in the 10 µM DEHP group was approximately 2-fold higher than that in the control group. The results of immunofluorescence staining were consistent with the above findings (Fig. 8 D). Fig. 8. Open in a new tab Regulatory effect of DEHP on STAT1 expression in HSG cells. ( A ) Western blot showing representative STAT1 and GAPDH bands from the same membrane/exposure. ( B ) Quantification of relative STAT1 protein expression (* p < 0.05, ** p < 0.01). ( C ) RT-qPCR analysis of relative STAT1 mRNA expression (* p < 0.05, ** p < 0.01). ( D ) Immunofluorescence staining of STAT1 (green) and nuclei (DAPI, blue) in HSG cells treated with different concentrations of DEHP Discussion To explore a potential environmental–molecular link between plasticizer exposure and Sjögren’s syndrome (SS), we selected DEHP as a representative high-production-volume phthalate. The rationale is based on (i) the ubiquity of DEHP exposure in the general population and (ii) published experimental evidence indicating that DEHP and its metabolites can modulate immune signaling and endocrine pathways relevant to autoimmunity [ 32 , 34 , 51 ]. However, we emphasize that direct epidemiological evidence establishing a causal relationship between DEHP exposure and SS is currently limited. Computational toxicity profiling using two major prediction platforms suggested potential hepatotoxicity and respiratory toxicity, and flagged immunotoxicity as a possible risk signal. Notably, the immunotoxicity endpoint in ProTox 3.0 shows comparatively lower predictive performance (approximately ~ 75% on the platform’s reported statistics), which increases uncertainty, therefore, we interpret this prediction conservatively and use it only as an auxiliary indicator rather than stand-alone evidence. In addition, ProTox 3.0 does not provide a dedicated model for endocrine disruption, and thus our in silico toxicity profiling cannot be used to draw conclusions regarding endocrine-disrupting potential. Accordingly, our rationale for exploring DEHP–SS links relies primarily on the integrated transcriptomic/network analyses and mechanistic consistency with prior experimental literature, while causality and endocrine-immune mechanisms require further validation. Our network toxicology analysis identified 27 overlapping targets between the predicted DEHP target space and SS-related transcriptomic candidates, and enrichment analysis showed that these genes cluster in interferon responses, chemokine signaling, Toll-like receptor pathways, and other immune-related processes. The chemokine signaling pathway, particularly the CXCL12/CXCR4 axis, has been demonstrated to be significantly upregulated in the salivary and lacrimal glands of SS patients and disease model mice, showing a positive correlation with the degree of lymphocyte infiltration in these glands [ 26 ]. The Toll-like receptor (TLR) signaling pathway serves as a critical bridge connecting environmental stimuli to autoimmune activation. It initiates innate immune responses by recognizing microbial pathogen signals, thereby breaking immune tolerance and contributing to SS pathogenesis. In genetically susceptible individuals, environmental factors like viruses can directly activate TLR family members (mainly TLR3, TLR7, and TLR9) on the surface of salivary gland epithelial cells (SGECs), inducing the release of damage-associated molecular patterns (DAMPs) and proinflammatory cytokines such as IL-6 and type I interferons, which trigger local inflammatory responses [ 8 , 36 ]. The IL-17 signaling pathway is a central mediator regulating humoral immunity and glandular tissue destruction in SS. It acts synergistically with other cytokines to enhance B-cell function and promote tissue fibrosis, thereby driving disease progression [ 38 , 46 ]. Furthermore, studies have found that DEHP can reshape immune cell migration patterns by activating chemokine signaling pathways, with its core mechanism focusing on the dysregulation of the CXCL12/CXCR4 axis [ 30 ]. Additionally, the TLR/MyD88 signaling pathway has been identified as a central hub for DEHP to initiate innate immune responses and amplify inflammatory cascades [ 11 ]. Nevertheless, enrichment results alone do not demonstrate that DEHP promotes SS, rather they indicate that the overlap genes are embedded in immune networks that are consistent with SS biology and potentially susceptible to environmental modulation. In summary, the accumulated evidence from our analyses and the existing literature suggests that DEHP-associated molecular networks overlap with immune pathways implicated in SS, but this does not demonstrate that DEHP exposure causes or increases the risk of Sjögren’s syndrome. By integrating network topology with supervised ML/SHAP interpretation, we prioritized five core genes (STAT1, PHGDH, ISG15, CXCL10, and CCL2). Importantly, the ML task in our study is SS classification (SS vs. control) using the DEHP–SS intersect genes as candidate features; thus, the ML results demonstrate predictive association within SS transcriptomic datasets rather than a direct effect of DEHP exposure. Differential expression analysis revealed significant upregulation of STAT1, ISG15, CXCL10, and CCL2, alongside marked downregulation of PHGDH in SS. Machine learning models further confirmed the high diagnostic value of STAT1 and CXCL10. The diagnostic significance of STAT1 in autoimmune diseases is well-established. Relevant studies indicate that the Janus kinase/signal transducer and activator of transcription (JAK/STAT) signaling pathway plays a central role in cellular communication, and dysregulation of STAT1 is closely associated with the pathogenesis of Sjögren’s syndrome [ 5 ]. Chemokines are pivotal in the immunopathogenesis of SS [ 3 ]. Notably, CXCL10 expression is abnormally elevated in patients with primary SS, suggesting its potential role as both a biomarker and a novel therapeutic target [ 20 ]. SHAP interpretability analysis highlighted STAT1 and PHGDH as the most influential predictors. Phosphoglycerate dehydrogenase (PHGDH), a key rate-limiting enzyme in the serine biosynthesis pathway, catalyzes the conversion of 3-phosphoglycerate to 3-phosphohydroxypyruvate. Elevated PHGDH activity in various immune and cancer cells can be mediated by genetic amplification, post-translational modifications, increased transcription, and allosteric regulation [ 28 ]. Molecular docking provided supportive computational evidence by suggesting that DEHP can adopt energetically favorable poses with the five prioritized protein. however, docking does not validate physical binding. Similarly, the single-cell atlas analysis provides cell-type STAT1 localization and does not establish mechanistic function or disease causation. Our in vitro experiments demonstrated that DEHP exposure increases STAT1 mRNA and protein levels in HSG cells, which is consistent with STAT1 being an SS-relevant immune regulator. Together, these findings support a mechanistically plausible DEHP–STAT1 axis, but they do not demonstrate that DEHP causes SS. The JAK/STAT family comprises four intracellular JAK tyrosine kinases and seven STAT transcription factors. Cytokines, including interleukin (IL)-6, interferon (IFN)-α, and IFN-γ, primarily activate the JAK/STAT pathway by inducing STAT phosphorylation via JAK activation. The JAK/STAT1 signaling pathway is aberrantly activated in patients with primary SS [ 53 ]. Emerging evidence also indicates that the JAK/STAT pathway significantly contributes to pSS pathogenesis through direct and indirect stimulation of B cells [ 33 ]. Elevated IFN-γ levels in SS trigger salivary gland epithelial cell (SGEC) death, and IFN-γ induces SGEC ferroptosis via the JAK/STAT1 pathway [ 7 ].DEHP can dysregulate the STAT1 signaling pathway. Studies have shown that DEHP induces lipid metabolism disorders and autophagy by modulating TYK2/STAT1 [ 52 ]. Furthermore, DEHP-induced inflammation in RAW264.7 macrophages can be suppressed via the STAT pathway [ 25 ]. Therefore, the predictions of our study align with existing experimental research, and the use of multiple complementary techniques and methods reinforces the accuracy and reliability of our findings. In summary, DEHP exposure in vitro induces STAT1 upregulation, a component of interferon-related immune signaling, but this finding alone is insufficient to establish a causal role of DEHP in Sjögren’s syndrome.Future work should incorporate direct binding/target-engagement assays (e.g., CETSA, SPR/ITC) and in vivo/clinical exposure-linked studies to test causality. The novelty of this study lies in the integration of network toxicology, ML/SHAP interpretability, molecular docking, and single-cell atlas interrogation with epithelial-cell validation to generate testable hypotheses linking DEHP to SS-relevant pathways. Nevertheless, several limitations merit explicit emphasis: (1) the GEO transcriptomic cohorts do not contain measured DEHP exposure information, so our analyses cannot attribute SS signatures to DEHP exposure directly; (2) the discovery cohort integrates salivary gland and blood datasets, and residual tissue/platform effects may persist despite batch correction; (3) available cohorts are predominantly female, which reflects SS epidemiology but may limit generalizability; (4) in silico toxicity and target prediction (especially immunotoxicity and endocrine disruption) have known performance constraints and may yield false positives; (5) docking provides likelihood estimates rather than validation, and negative controls/structural analog comparisons were not included; (6) in vitro validation was restricted to STAT1 in one epithelial cell line; (7) because DEHP may influence cellular metabolism, future Western blot validations will benefit from additional loading/normalization controls (e.g., β-actin or total-protein staining) alongside GAPDH. Conclusion In summary, this study presents an integrated computational–experimental workflow suggesting that DEHP-related target networks intersect with SS-relevant immune signatures and highlighting STAT1 as a prioritized candidate within this intersection. Our ML models identify genes that predict SS status among the candidate set, and docking/scRNA-seq analyses provide supportive, non-definitive evidence regarding binding propensity and cell-type localization. In vitro experiments further show that DEHP can upregulate STAT1 in salivary gland epithelial cells. Overall, these findings generate focused hypotheses for future mechanistic and epidemiological validation. Despite the multi-layered analyses presented, further work is needed to strengthen causal inference. Future studies should incorporate exposure-linked human cohorts, appropriate animal models of SS, and target-engagement assays to validate whether DEHP (and major DEHP metabolites such as mono(2-ethylhexyl) phthalate, MEHP) directly modulates STAT1 signaling in vivo. In addition, other endocrine-disrupting chemicals (e.g., BPA) and phthalate mixtures should be evaluated to assess chemical specificity and potential mixture effects. Supplementary Information Below is the link to the electronic supplementary material. Supplementary Material 1 (2.5KB, txt) Supplementary Material 2 (5.2MB, zip) Acknowledgements We would like to acknowledge the reviewers for their helpful comments on this paper. Author contributions All authors made a significant contribution to the work reported, whether that is in the conception, Lili C study design, execution, acquisition of data, analysis and interpretation, or in all these areas; Zhongfu T took part in drafting, Ming L revising or critically reviewing the article; Chuanbing Huang agree to be accountable for all aspects of the work. Funding This work was supported by the National Natural Science Foundation of China General Program (82574970); Major and Difficult Diseases Clinical Research Project on Collaboration between Traditional Chinese and Western Medicine (ZDYN-2024-A-146); Anhui Provincial Clinical Medicine Research and Transformation Special Projects (202304295107020114, 202304295107020115); Institute of Xin’an Medicine and Modernization of Traditional Chinese Medicine, Great Health Research Institute Special Projects (2023CXMMTCM015, 2023CXMMTCM004); Anhui Provincial Graduate Quality Engineering Graduate Innovation and Entrepreneurship Practice Project (2024cxcysj121); 2024 Anhui Provincial Higher Education Institutions Key Scientific Research Project (2024AH050957); 2024 Anhui Provincial Health Commission Scientific Research Project (2024Aa30428). The funding agencies had no role in study design, data collection, and analysis, decision to publish, or preparation of the manuscript. Data availability No datasets were generated or analysed during the current study. Declarations Ethical approval Not applicable. Competing interests The authors declare no competing interests. Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. References 1. An Q, et al. A comprehensive review on machine learning in healthcare industry: classification, restrictions, opportunities and challenges. Sensors (Basel). 2023;23. [ DOI ] [ PMC free article ] [ PubMed ] 2. André F, Böckle BC. Sjögren’s syndrome. J Dtsch Dermatol Ges. 2022;20:980–1002. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 3. Antonelli A, et al. Chemokine (C-X-C motif) ligand (CXCL)10 in autoimmune diseases. Autoimmun Rev. 2014;13:272–80. [ DOI ] [ PubMed ] [ Google Scholar ] 4. Banerjee P, et al. ProTox 3.0: a webserver for the prediction of toxicity of chemicals. Nucleic Acids Res. 2024;52:W513–20. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 5. Banerjee S, et al. JAK-STAT signaling as a target for inflammatory and autoimmune diseases: current and future prospects. Drugs. 2017;77:521–46. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 6. Björk A, et al. Environmental factors in the pathogenesis of primary Sjögren’s syndrome. J Intern Med. 2020;287:475–92. [ DOI ] [ PubMed ] [ Google Scholar ] 7. Cao T, et al. Interferon-γ induces salivary gland epithelial cell ferroptosis in Sjogren’s syndrome via JAK/STAT1-mediated inhibition of system Xc(). Free Radic Biol Med. 2023;205:116–28. [ DOI ] [ PubMed ] [ Google Scholar ] 8. Chen JQ, et al. Toll-like receptor pathways in autoimmune diseases. Clin Rev Allergy Immunol. 2016;50:1–17. [ DOI ] [ PubMed ] [ Google Scholar ] 9. Cheng C et al. A review of single-cell RNA-seq annotation, integration, and cell-cell communication. Cells. 2023;12. [ DOI ] [ PMC free article ] [ PubMed ] 10. Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. 2014;8(Suppl 4):S11. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 11. Dai XY, et al. Role of toll-like receptor/MyD88 signaling in lycopene alleviated Di-2-ethylhexyl Phthalate (DEHP)-induced inflammatory response. J Agric Food Chem. 2022;70:10022–30. [ DOI ] [ PubMed ] [ Google Scholar ] 12. Daina A, et al. SwissTargetPrediction: updated data and new features for efficient prediction of protein targets of small molecules. Nucleic Acids Res. 2019;47:W357–64. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 13. Davis AP, et al. Comparative Toxicogenomics Database (CTD): update 2023. Nucleic Acids Res. 2023;51:D1257–62. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 14. Dickinson Q, Meyer JG. Positional SHAP (PoSHAP) for interpretation of machine learning models trained from biological sequences. PLoS Comput Biol. 2022;18:e1009736. [ DOI ] [ PMC free article ] [ PubMed ] 15. Eberhardt J, Santos-Martins D, Tillack AF, Forli S. AutoDock Vina 1.2.0: New docking methods, expanded force field, and Python bindings. J Chem Inf Model. 2021;61:3891–8. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 16. Erythropel HC, et al. Leaching of the plasticizer di(2-ethylhexyl)phthalate (DEHP) from plastic containers and the question of human exposure. Appl Microbiol Biotechnol. 2014;98:9967–81. [ DOI ] [ PubMed ] [ Google Scholar ] 17. Franzén O, Gan LM, Björkegren JLM. PanglaoDB: a web server for exploration of mouse and human single-cell RNA sequencing data. Database (Oxford). 2019;baz046. [ DOI ] [ PMC free article ] [ PubMed ] 18. Gao J, et al. Integrating machine learning and molecular docking to decipher the molecular network of aflatoxin B1-induced hepatocellular carcinoma. Int J Surg. 2025;111:4539–49. [ DOI ] [ PubMed ] [ Google Scholar ] 19. Geenen S, et al. Systems biology tools for toxicology. Arch Toxicol. 2012;86:1251–71. [ DOI ] [ PubMed ] [ Google Scholar ] 20. Hakbilen S, et al. The role of CXCL9, CXCL10, and CXCL13 chemokines in patients with Sjögren’s syndrome. Clin Rheumatol. 2025;44:1635–42. [ DOI ] [ PubMed ] [ Google Scholar ] 21. Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. 2007;8:118–27. [ DOI ] [ PubMed ] [ Google Scholar ] 22. Kanehisa M, Sato Y, Kawashima M, Furumichi M, Tanabe M. KEGG: integrating viruses and cellular organisms. Nucleic Acids Res. 2023;51(D1):D587–92. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 23. Keiser MJ, Roth BL, Armbruster BN, Ernsberger P, Irwin JJ, Shoichet BK. Relating protein pharmacology by ligand chemistry. Nat Biotechnol. 2007;25:197–206. [ DOI ] [ PubMed ] [ Google Scholar ] 24. Kim S, Chen J, Cheng T, Gindulyte A, He J, He S, Li Q, Shoemaker BA, Thiessen PA, Yu B, Zaslavsky L, Zhang J, Bolton EE. PubChem 2023 update. Nucleic Acids Res. 2023;51(D1):D1373–80. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 25. Kim JH, et al. N-actylcysteine inhibits diethyl phthalate-induced inflammation via JNK and STAT pathway in RAW264.7 macrophages. BMC Mol Cell Biol. 2025;26:12. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 26. Kurosawa M, et al. NF-κB2 controls the migratory activity of memory T cells by regulating expression of CXCR4 in a mouse model of Sjögren’s syndrome. Arthritis Rheumatol. 2017;69:2193–202. [ DOI ] [ PubMed ] [ Google Scholar ] 27. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 28. Lee CM, et al. PHGDH: a novel therapeutic target in cancer. Exp Mol Med. 2024;56:1513–22. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 29. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012;28:882–3. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 30. Li J, et al. Combining epidemiology, network toxicology and in vivo experimental validation to uncover takeout-mediated DEHP/MEHP toxicity in inflammatory bowel disease. Free Radic Biol Med. 2025a;239:1–13. [ DOI ] [ PubMed ] [ Google Scholar ] 31. Li NR, et al. Aspartame increases the risk of liver cancer through CASP1 protein: A comprehensive network analysis insights. Ecotoxicol Environ Saf. 2025b;294:118089. [ DOI ] [ PubMed ] [ Google Scholar ] 32. Linghu D, et al. Diethylhexyl phthalate induces immune dysregulation and is an environmental immune disruptor. J Hazard Mater. 2024;480:136244. [ DOI ] [ PubMed ] [ Google Scholar ] 33. Liu Y, et al. miR-216a-3p alleviates primary Sjögren’s syndrome by regulating the STAT1/JAK signaling pathway. Biochem Biophys Res Commun. 2025;758:151647. [ DOI ] [ PubMed ] [ Google Scholar ] 34. Nygaard UC, et al. Immune cell profiles associated with measured exposure to phthalates in the Norwegian EuroMix biomonitoring study - A mass cytometry approach in toxicology. Environ Int. 2021;146:106283. [ DOI ] [ PubMed ] [ Google Scholar ] 35. Ponce-Bobadilla AV et al. Practical guide to SHAP analysis: explaining supervised machine learning model predictions in drug development. Clin Transl Sci. 2024;17:e70056. [ DOI ] [ PMC free article ] [ PubMed ] 36. Qi W, et al. Advances in cellular and molecular pathways of salivary gland damage in Sjögren’s syndrome. Front Immunol. 2024;15:1405126. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 37. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47. [ DOI ] [ PMC free article ] [ PubMed ] 38. Schinocca C, et al. Role of the IL-23/IL-17 pathway in rheumatic diseases: an overview. Front Immunol. 2021;12:637829. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 39. Schrödinger LLC. The PyMOL Molecular Graphics System, Version 2.x. 2023. 40. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 41. Szklarczyk D, Kirsch R, Koutrouli M, Nastou K, Mehryary F, Hachilif R, Gable AL, Fang T, Doncheva NT, Pyysalo S, Bork P, Jensen LJ, von Mering C. STRING v11.5: protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2023;51(D1):D638–46. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 42. The Gene Ontology Consortium. The Gene Ontology knowledgebase in 2023. Nucleic Acids Res. 2023;51(D1):D553–62. [ Google Scholar ] 43. The UniProt Consortium. UniProt: the Universal Protein Knowledgebase in 2023. Nucleic Acids Res. 2023;51(D1):D523–31. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 44. Tian Y, et al. Advances in pathogenesis of Sjögren’s syndrome. J Immunol Res. 2021:5928232. [ DOI ] [ PMC free article ] [ PubMed ] 45. Wang Z, et al. Improving chemical similarity ensemble approach in target prediction. J Cheminform. 2016;8:20. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 46. Wei L, et al. Upregulation of IL-6 expression in human salivary gland cell line by IL-17 via activation of p38 MAPK, ERK, PI3K/Akt, and NF-κB pathways. J Oral Pathol Med. 2018;47:847–55. [ DOI ] [ PubMed ] [ Google Scholar ] 47. Wu T, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innov (Camb). 2021;2:100141. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 48. Xiong G, et al. ADMETlab 2.0: an integrated online platform for accurate and comprehensive predictions of ADMET properties. Nucleic Acids Res. 2021;49:W5–14. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 49. Yang L, et al. Exposure to di-2-ethylhexyl phthalate (DEHP) increases the risk of cancer. BMC Public Health. 2024;24:430. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 50. Zhang TP, et al. Exposure to particulate pollutant increases the risk of hospitalizations for Sjögren’s syndrome. Front Immunol. 2022a;13:1059981. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 51. Zhang Y, et al. Health risks of phthalates: A review of immunotoxicity. Environ Pollut. 2022b;313:120173. [ DOI ] [ PubMed ] [ Google Scholar ] 52. Zhang YZ, et al. Di (2-ethylhexyl) phthalate Disorders Lipid Metabolism via TYK2/STAT1 and Autophagy in Rats. Biomed Environ Sci. 2019;32:406–18. [ DOI ] [ PubMed ] [ Google Scholar ] 53. Zhong Y, et al. Screening biomarkers for Sjogren’s Syndrome by computer analysis and evaluating the expression correlations with the levels of immune cells. Front Immunol. 2023;14:1023248. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Associated Data This section collects any data citations, data availability statements, or supplementary materials included in this article. Supplementary Materials Supplementary Material 1 (2.5KB, txt) Supplementary Material 2 (5.2MB, zip) Data Availability Statement No datasets were generated or analysed during the current study. Articles from BMC Pharmacology & Toxicology are provided here courtesy of BMC ACTIONS View on publisher site PDF (6.9 MB) Cite Collections Permalink PERMALINK Copy RESOURCES Similar articles Cited by other articles Links to NCBI Databases Cite Copy Download .nbib .nbib Format: AMA APA MLA NLM Add to Collections Create a new collection Add to an existing collection Name your collection * Choose a collection Unable to load your collection due to an error Please try again Add Cancel Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top