The SMILES representations for twenty principal PFAS compounds ( Table 1 ) were obtained from the PubChem repository (https://pubchem.ncbi.nlm.nih.gov/) (18). These molecular identifiers were subsequently analyzed using ProTox version 3.0 (https://tox.charite.de/protox3/) (19) and the ChEMBL database (https://www.ebi.ac.uk/chembl/) (20) to enable comprehensive computational assessment of their toxicological profiles.

The selected PFAS compounds, represented in SMILES format and sourced from the PubChem database, were analyzed using the SwissTargetPrediction platform (http://www.swisstargetprediction.ch/)(21) and the SuperPred tool (https://prediction.charite.de) (22). In both systems, Homo sapiens was designated as the target organism to ensure biological relevance to humans, with a screening threshold set at “Probability ≥ 0.” SwissTargetPrediction utilizes machine learning algorithms that analyze molecular fingerprints alongside curated ligand–protein interaction datasets to predict probable targets. Likelihood scores are assigned based on structural and chemical feature comparisons. To maximize the inclusivity of potential targets, a minimum probability cutoff of zero was applied, thereby prioritizing sensitivity in target identification (23). Subsequently, target specificity was enhanced through systematic refinement and enrichment processes. These computational tools collectively enabled the prediction, filtration, and validation of molecular targets potentially interacting with PFAS compounds.

To uncover genes implicated in DKD, searches were conducted in the GeneCards database (https://www.genecards.org/) (24) and the Online Mendelian Inheritance in Man (OMIM) resource (https://omim.org/) (25) using the keyword “diabetic kidney disease.” After removing duplicate entries, a comprehensive set of DKD-related targets was compiled by integrating data from both databases. To enhance dataset relevance, the median “Relevance score” across all identified targets was used as a threshold for filtering (26). Employing the median as a cutoff ensures uniformity and facilitates the integration of heterogeneous gene expression datasets. This strategy strikes a balance between selectivity and biological significance, while mitigating potential biases arising from variable data volumes across platforms (27). Additionally, Venn diagram analysis was performed to identify overlapping targets shared between PFAS-associated and DKD-related gene sets. These intersecting targets were designated as candidate toxicological targets through which PFAS may contribute to the pathogenesis of DKD.

The intersecting targets identified as potential mediators of PFAS-induced DKD-represented by the overlapping region of the Venn diagram-were submitted to the STRING database (https://string-db.org/) using the “Multiple Proteins” analysis function. The species parameter was set to Homo sapiens, and the minimum required interaction score was configured to “high confidence” (>0.7). Additionally, the false discovery rate (FDR) threshold was set to a stringent level to minimize false positives. Applying a high confidence cutoff (>0.7) effectively excludes low-quality or speculative interactions, retaining only experimentally validated or strongly supported associations (28, 29). Concurrently, stringent FDR control significantly reduces the likelihood of false-positive signals arising from random errors, thereby enhancing the reliability and reproducibility of the interaction network.

The results obtained from the STRING database were imported into Cytoscape software (version 3.10.3) to calculate network parameters for each node and to optimize the visualization of molecular interactions (30). A PPI network was constructed based on the topological characteristics of the nodes. Using the CytoNCA plugin “Centiscape”, key topological parameters for each gene node were calculated, including Closeness Centrality, Betweenness Centrality, and Degree Value. Hub targets were identified based on the following three criteria:

The selection of these three parameters reflects their complementary roles in network topology analysis: degree centrality quantifies the importance of a node based on the number of its direct connections, thereby identifying locally dense hub nodes (31); betweenness centrality measures the extent to which a node functions as a critical intermediary facilitating global information flow within the network (32); and closeness centrality evaluates the ability of a node to efficiently disseminate information across the entire network (33). The combined application of these metrics enables a more comprehensive and nuanced assessment of node significance. By integrating multiple topological indicators, this approach minimizes selection bias and enhances the robustness and reliability of the identified hub targets (34).

To elucidate the biological functions of potential targets involved in PFAS-induced DKD, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the DAVID database (https://david.ncifcrf.gov). The GO analysis encompassed three domains: biological processes (BP), cellular components (CC), and molecular functions (MF), providing a comprehensive overview of the primary biological roles of these targets. KEGG enrichment analysis identified key signaling pathways associated with PFAS-induced DKD targets, applying a FDR threshold of < 0.05 to ensure statistical robustness. Furthermore, KEGG pathway enrichment was specifically conducted on the subset of selected hub targets to delineate their involvement in DKD-related pathways, thereby underscoring critical signaling cascades implicated in disease pathogenesis.

To validate the expression patterns of the identified hub targets, immunohistochemical (IHC) staining data were retrieved from the Human Protein Atlas (HPA; https://www.proteinatlas.org/). Protein expression levels were evaluated across various organs by examining staining intensity, the identity of positively stained cell types, and tissue-specific distribution. Additionally, single-cell transcriptomic analyses were conducted to characterize target gene expression across distinct cell populations, including endocrine and epithelial cells. Uniform Manifold Approximation and Projection (UMAP) was employed to visualize gene expression distribution among heterogeneous cell types. Subcellular localization of the corresponding proteins was further assessed using fluorescence microscopy, providing spatial resolution of protein enrichment within key cellular compartments such as the nucleus, cytoplasm, mitochondria, and endoplasmic reticulum.

Differential expression of hub genes between normal renal tissues and those affected by chronic kidney disease (CKD) was validated using the Nephroseq database (https://www.nephroseq.org/). To corroborate these findings in the context of DKD, gene expression datasets were obtained from the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/), with particular focus on the GSE30122 dataset comprising 19 DKD samples and 50 normal controls, enabling assessment of differential gene expression between diseased and healthy renal tissues.

To explore the relationships among immune cells during infiltration, the immune microenvironment was profiled using CIBERSORT (35), enabling comparison of relative immune cell subsets and immune scores in the GSE30122 dataset. A significance threshold of P < 0.05 was applied for selecting features for further analysis. Correlation matrices of immune infiltration were then generated, categorizing correlation strength as weak (absolute coefficients between 0.10 and 0.39), moderate (0.40 to 0.69), and strong (ranging from 0.70 to 0.89).

GSEA was performed using GSEA software (version 3.0), obtained from the Broad Institute website (http://software.broadinstitute.org/gsea/index.jsp) (36). Based on the expression levels of core genes, samples were stratified into high-expression (≥50%) and low-expression (<50%) groups. Curated gene sets (c2.cp.kegg.v7.4.symbols.gmt) were sourced from the Molecular Signatures Database (MSigDB, http://www.gsea-msigdb.org/gsea/downloads.jsp) (37) for pathway and molecular mechanism analyses. The minimum and maximum gene set sizes were set to 5 and 5000, respectively, with 1,000 permutations performed to ensure robust statistical assessment. Pathways achieving a nominal P < 0.05 and a FDR < 0.25 were considered statistically significant.

Discovery Studio, a specialized software for molecular simulations in the life sciences, is commonly employed in drug development and biomacromolecular modeling, particularly for analyzing protein and antibody structures (38). In this study, PFAS compounds were employed as ligands for molecular docking simulations. Receptor proteins were selected based on hub targets exhibiting the highest degree centrality within the PPI network. Three-dimensional structures of the PFAS molecules were obtained from the PubChem database in SDF format. Corresponding 3D protein structures of the primary toxic targets were retrieved from the RCSB Protein Data Bank (http://www1.rcsb.org/), specifying Homo sapiens as the source organism and restricting resolution to 0–3.0 Å. Using Discovery Studio 2019, receptor structures were preprocessed by adding hydrogen atoms, removing crystallographic water molecules, and reconstructing missing loop regions. Potential binding pockets were subsequently predicted. Ligands were prepared through hydrogenation and conformational analysis. Molecular docking simulations were performed with the CDOCKER module to predict the optimal binding orientations between PFAS ligands and selected targets (39). CDOCKER interaction energies were calculated, with higher scores indicative of stronger ligand–receptor binding affinities.

Targets for 20 major PFAS compounds were predicted using the SwissTargetPrediction and SuperPred databases. Following data integration and removal of duplicates, a total of 424 unique PFAS-related targets were identified. Gene targets associated with DKD were obtained from the GeneCards and OMIM databases, resulting in 9,999 unique targets after filtering. Intersection analysis using Venny 2.1.0 revealed 86 overlapping targets ( Figure 1A ), which represent potential toxicological mediators through which PFAS may contribute to the pathogenesis of DKD.

To analyze these common targets, a PPI network was constructed by submitting the 86 shared targets to the STRING database ( Figure 1B ), with parameters set to “Multiple proteins,” species designated as Homo sapiens, and a high-confidence minimum interaction score. This yielded a network comprising 83 nodes and 572 edges ( Figure 1C ), illustrating the intricate interactions between PFAS-associated toxic targets and the pathophysiology of DKD. To further elucidate the mechanisms underlying PFAS-induced toxicity in DKD, the network was imported into Cytoscape 3.9.0, where the CytoHubba plugin was employed for clustering analysis. The top 10 hub targets, identified based on degree centrality, included MMP9, BCL2, CYP3A43, ACE, HNF4A, HSP90AA1, AGTR1, MMP2, AGTR2, and HMGCR ( Figure 1D ), suggesting their pivotal roles in mediating PFAS-induced DKD toxicity.

GO and KEGG pathway enrichment analyses were conducted on the 86 candidate targets associated with PFAS-induced DKD toxicity. GO analysis revealed that, in terms of BP, these targets are primarily involved in transport, regulation of biological processes, and lipid metabolic processes ( Figure 1E ). Regarding CC, the targets are predominantly associated with intrinsic components of membranes and plasma membrane parts ( Figure 1F ). For MF, the targets are mainly linked to transition metal ion binding, zinc ion binding, and lyase activity ( Figure 1G ). KEGG pathway enrichment identified the top 10 significantly enriched pathways, with the most prominent including nitrogen metabolism, PPAR signaling, endocrine resistance, insulin resistance, and AMPK signaling pathways ( Figure 1H ). Notably, the majority of these pathways are related to metabolic regulation ( Figures 1I, J ). Collectively, these findings provide valuable insights into the potential mechanisms by which PFAS contribute to DKD toxicity, emphasizing key biological processes and signaling pathways involved in disease progression.

Validation using the HPA database revealed that the hub targets implicated in PFAS-induced DKD toxicity-namely MMP9, BCL2, CYP3A43, ACE, HNF4A, and HMGCR-were predominantly localized in renal tubules, with no detectable expression in glomeruli. HSP90AA1 and AGTR1 were undetected in both tubules and glomeruli. MMP2 demonstrated moderate expression in the renal tubules but was expressed at low levels in the glomeruli. Immunohistochemical data for AGTR2 are currently unavailable in the HPA database ( Figure 2A ). Single-cell transcriptomic analysis indicated the expression of these targets across germ cells, neuronal cells, blood immune cells, and endocrine cells ( Figure 2B ). Subcellular localization analysis further showed that HSP90AA1, MMP9, and CYP3A43 were primarily localized in the cytosol; BCL2 was enriched in mitochondria; AGTR1, ACE, and MMP2 were predominantly found in vesicles; and HNF4A was localized within the nucleoplasm ( Figure 2C ).

Analysis of the Nephroseq database revealed that MMP9, HSP90AA1, BCL2, and CYP3A43 were significantly upregulated in CKD, whereas AGTR1, ACE, HNF4A, MMP2, and HMGCR exhibited significant downregulation, with all differences reaching statistical significance (P < 0.05). Expression data for AGTR2 were not available in this database ( Figure 2D ). In the GSE30122 dataset, several genes-including MMP9, HSP90AA1, BCL2, MMP2, AGTR1, and HMGCR-were upregulated in DKD, whereas ACE, HNF4A, CYP3A43, and AGTR2 were downregulated. Among these, the differential expression of CYP3A43, ACE, MMP2, and AGTR2 was statistically significant (P < 0.05) ( Figure 2E ).

The immune microenvironment was characterized using CIBERSORT to assess correlations of immune cell infiltration within the GSE30122 dataset ( Figure 3A ). Moderate correlations were observed between plasma cells and resting mast cells, as well as between regulatory T cells (Tregs) and resting dendritic cells. Weak correlations were detected between activated natural killer (NK) cells and M1 macrophages, and between monocytes and M1 macrophages. Additionally, CD8+ T cells, activated mast cells, and follicular helper T cells were significantly decreased in DKD compared to controls (P < 0.05) ( Figure 3B ).

GSEA was performed on hub targets to elucidate potential pathogenic mechanisms. MMP9 was predominantly enriched in the allograft rejection pathway, while HSP90AA1 was associated with limonene and pinene degradation. BCL2 showed significant enrichment in the TGF-β signaling pathway. High CYP3A43 expression correlated with enrichment in olfactory transduction, whereas AGTR1 was primarily enriched in nitrogen metabolism. Elevated ACE expression was linked to maturity onset diabetes of the young, and HNF4A was enriched in pyruvate metabolism. MMP2 was associated with extracellular matrix (ECM) receptor interaction, and high AGTR2 expression correlated with Type II diabetes mellitus pathways. Finally, HMGCR was enriched in alanine, aspartate, and glutamate metabolism. Collectively, these findings suggest that these hub targets may play critical roles in immune-mediated renal injury ( Figure 3C ).

Molecular docking heatmap analysis demonstrated strong binding affinities between PFAS compounds and the top 10 hub targets ( Figure 4A ), including MMP9 (PDB ID: 2OVX), HSP90AA1 (PDB ID: 2QF6), BCL2 (PDB ID: 4IEH), CYP3A43 (PDB ID: 4NY4), AGTR1 (PDB ID: 4YAY), ACE (PDB ID: 4X5K), HNF4A (PDB ID: 1PZL), MMP2 (PDB ID: 3AYU), AGTR2 (PDB ID: 5UNF), and HMGCR (PDB ID: 1HW8). The top ten docking poses with the lowest binding energies were selected for visualization using PyMOL (40). Notably, perfluorooctanoic acid (PFOA) interacted with MMP9 residues GLN-402, ALA-189, and LEU-188, while perfluorononanoic acid (PFNA) engaged with HSP90AA1 residues GLN-402, ALA-189, LEU-188, GLU-208, HIS-411, and HIS-401. FTS bound to BCL2 at ARG-66 and GLY-104 residues, and perfluoropropanoic acid (PFPrA) interacted with AGTR1 residues ASN-117, LEU-77, LYS-126, THR-90, GLY-78, CYS-79, ARG-85, GLY-87, and GLY-89. Perfluoroheptanoic acid (PFHpA) bound CYP3A43 at TYR-35 and ARG-167 residues. The interaction between NMeFOSA and ACE involved GLY-237, whereas GenX engaged HNF4A at ILE-53, ARG-52, and ILE-19 residues. Perfluoropentanoic acid (PFPeA) bound AGTR2 at GLU-559, HIS-63, ASN-755, LYS-691, ILE-746, SER-745, and TYR-749 residues. Perfluorodecanoic acid (PFDA) interacted with MMP2 residues GLU-202, ALA-205, GLN-206, and TYR-204, and perfluorohexane sulfonate (PFHxS) bound HMGCR at LEU-477 ( Figure 4B ).