Gene expression profile datasets related to NOA were obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/), specifically accessions GSE4797, GSE45885, GSE45887, and GSE145467. A combined dataset (GSE_combined) was created by merging GSE45885 and GSE45887 using the “sva” R package function to correct for batch effects. This combined dataset was used for differential expression analysis to identify disease-related targets. GSE4797 and GSE145467 were used as independent cohorts to validate the expression of key genes.

Genes associated with 18 PCD subtypes were identified through a comprehensive search of the MSigDB, GeneCards, and KEGG databases, along with a review of relevant literature. The PCD subtypes included: apoptosis, pyroptosis, ferroptosis, autophagy, necroptosis, cuproptosis, parthanatos, endogenous cell death, eosinophil cell death, lysosome-dependent cell death, alkaliptosis, oxytosis, neutrophil extracellular trap-associated cell death, immunogenic cell death, anoikis, paraptosis, large-scale apoptotic cell death, and intracellular engulfment-related cell death. Differential expression analysis was performed on the integrated GSE_combined dataset using the “Limma” R package with thresholds of |log2FC| > 0.65 and an adjusted p-value <0.05 to identify NOA-related targets. Venn analysis was subsequently used to reveal dysregulated PCD-related genes in NOA patients, representing the imbalanced PCD pathways involved in NOA pathogenesis.

Gene Ontology (GO) categorizes gene function into Cellular Component (CC), Molecular Function (MF), and Biological Process (BP) domains. The Kyoto Encyclopedia of Genes and Genomes (KEGG) systematically link genomic information with higher-level functional pathways. GO and KEGG enrichment analyses were performed using the R package clusterProfiler.

The STRING database (confidence score threshold ≥0.4), which integrates known and predicted protein associations including physical interactions and functional links, was used to analyze interactions among the overlapping targets. After retrieving PPI data, non-essential targets were filtered out to construct a target-protein interaction network. This interaction network (in TSV format) was imported into Cytoscape 3.9.0 for visualization. Centiscape2.2, a plugin within Cytoscape, was used to calculate topological properties (Degree, Closeness, Betweenness) to assess node importance. Topological analysis using “CentiScape 2.2”identified hub targets based on three centrality measures: Closeness Centrality ≥0.00281, Betweenness Centrality ≥240.538, and Degree Centrality ≥7.415.

This study employed a two-sample Mendelian randomization (MR) analysis to infer a potential causal relationship between genetically predicted gene expression and the risk of male infertility (MI). This methodological approach relies on three core instrumental variable (IV) assumptions: (1) Relevance–the genetic IVs must be robustly associated with the exposure (gene expression level); (2) Independence–the IVs must not be associated with any known or unknown confounders of the exposure-outcome relationship; and (3) Exclusion Restriction–the IVs must influence the outcome solely through the exposure, with no alternative biological pathways. To obtain valid IVs adhering to these assumptions, stringent criteria were applied. Genetic instruments for gene expression (exposure) were selected as cis-expression quantitative trait loci (cis-eQTLs) from the eQTLGen Consortium, a resource integrating data from 31,684 individuals predominantly of European ancestry (Kurki et al., 2023). Single-nucleotide polymorphisms (SNPs) were retained at a genome-wide significance threshold (p < 5 × 10⁻⁸) and clumped (r ² < 0.001, window size = 10,000 kb) to ensure independence, with palindromic SNPs removed. The strength of each IV was quantified using the F-statistic (F > 10 indicated sufficient strength to minimize weak instrument bias). Summary statistics for the outcome (MI) were sourced from the FinnGen Release 12 genome-wide association study, comprising 1,852 cases and 153,573 controls, all of European ancestry (Võsa et al., 2021). The primary causal estimate was derived using the Inverse-Variance Weighted (IVW) method.

To critically evaluate the validity of the MR assumptions and the robustness of findings, a suite of sensitivity analyses was performed: (1) Cochran’s Q test assessed heterogeneity across IVs (with a random-effects IVW model applied if P < 0.05); (2) the MR-Egger intercept test detected directional horizontal pleiotropy, a violation of the independence and exclusion restriction assumptions; and (3) leave-one-out analysis verified that the causal estimate was not disproportionately driven by any single influential SNP. A significant causal association was concluded if the IVW P-value was <0.05, supported by consistent results from sensitivity analyses indicating no substantial pleiotropy or heterogeneity.

The GC-1 spg (ts) mouse spermatogonial cell line was obtained from the Cell Bank of the Chinese Academy of Sciences. Cells were cultured according to the manufacturer’s protocol at 37 °C in a humidified atmosphere containing 5% CO₂. The culture medium was high-glucose Dulbecco’s Modified Eagle Medium (DMEM), and cells were seeded at a density of 1 × 10⁵ cells. The complete growth medium was supplemented with 10% fetal bovine serum (FBS). To investigate the role of TLR4 in non-obstructive azoospermia (NOA), an inflammatory model was established based on previous literature. Cells were treated with 1 μg/mL lipopolysaccharide (LPS; Sigma, L2880) for 12 h to mimic testicular inflammation associated with human infertility. Subsequently, to verify the specific role of TLR4, this inflammatory model was treated with 10 μM TAK-242, a specific TLR4 inhibitor, as reported in prior studies. Similarly, the intervention doses for the selected compounds, perfluorooctanoic acid (PFOA) and baicalein, were derived from established literature. The final concentrations used were 500 μM for PFOA and 20 μM for baicalein (Biswas and Yenugu, 2013; Ye et al., 2019; Luo et al., 2025; Wan et al., 2018; Chen et al., 2021).

Total protein was extracted using RIPA lysis buffer (Solarbio, China) supplemented with a protease inhibitor cocktail (Yeasen, China). Protein concentration was quantified using the Pierce™ BCA Protein Assay Kit (Thermo Fisher Scientific). Amounts of protein (4–6 µg per lane) were separated on 8%–12% gradient SDS-PAGE gels and transferred to PVDF membranes using a semi-dry transfer system (Bio-Rad). Membranes were blocked with 5% skim milk in TBST for 1 h at room temperature, then incubated overnight at 4 °C with primary antibodies. The primary antibodies used were anti-TLR4 (1:2,000, # 19811-1-AP, Proteintech), anti-NF-κB p65 (1:500, # 10745-1-AP, Proteintech) and anti-GAPDH (1:50,000, # 60004-1-Ig, Proteintech). After washing, incubate the membrane with HRP conjugated anti rabbit IgG secondary antibody at room temperature. Finally, visualize the protein bands and use Image Lab ™ Quantitative analysis was performed using Bio Rad software.

The Comparative Toxicogenomics Database (CTD) was utilized to identify small molecule compounds targeting the key gene, with a focus on common environmental exposures and natural bioactive products. Following initial screening, molecular docking simulations were performed to validate compound-target interactions. The molecular structures of environmental exposures and natural products were retrieved from PubChem. The three-dimensional structure of the key target protein was obtained from the AlphaFold Protein Structure Database. AutoDockTools 1.5.7 was used for structure preparation and docking simulations to predict binding modes, binding affinity (ΔG), and functional implications. A binding energy <0 kcal/mol indicates spontaneous binding, while an energy <−5.0 kcal/mol suggests a stable interaction.

Database integration identified 1,964 PCD-related genes across 18 subtypes: apoptosis (580), pyroptosis (52), ferroptosis (88), autophagy (367), necroptosis (101), cuproptosis (19), parthanatos (9), netosis (15), lysosome-dependent death (220), alkaliptosis (7), oxytosis (5), NETosis (24), immunogenic cell death (34), anoikis (338), paraptosis (66), mitotic catastrophe (8), and entosis (23). Complete details are provided in Supplementary Table S1. Batch effect correction was applied to the merged dataset (Figure 1A), yielding an integrated training set comprising 8 control and 43 NOA samples. Differential expression analysis identified 2,025 significant DEGs (1,484 downregulated and 541 upregulated in NOA versus controls) (Figure 1B). Venn analysis subsequently identified 150 dysregulated PCD-related genes, representing pathogenic PCD regulators in NOA (Figure 1C).

GO enrichment analysis of the 150 NOA-PCD genes revealed their significant involvement in crucial biological processes, including positive and negative regulation of apoptotic signaling pathway, regulation of autophagy, cellular response to nutrient levels, cellular response to extracellular stimuli, and regulation of inflammatory response (Figure 1D). KEGG pathway analysis further identified pathways closely associated with cell proliferation and programmed cell death, such as the PI3K-Akt signaling pathway, mTOR signaling pathway, HIF-1 signaling pathway, ferroptosis, cellular senescence, necroptosis, apoptosis, Toll-like receptor signaling pathway, and MAPK signaling pathway (Figure 1E). These findings collectively suggest that dysregulated PCD may contribute to NOA pathogenesis through multiple distinct pathways.

The 150 overlapping NOA-PCD target genes were imported into the STRING database for PPI analysis, applying a confidence threshold of ≥0.4. After filtering disconnected nodes, 143 NOA-PCD target proteins remained in the final network. The PPI network was visualized using Cytoscape. Nodes were arranged using a degree-based layout, where darker colors and larger diameters indicate stronger protein interactions (Figure 2A). Topological analysis identified 10 hub genes within the NOA-PCD interaction network: HIF1A, TLR4, MDM2, GPX4, SNCA, MTOR, CSNK2A2, ATG5, CTSS, and PIK3CA (Figures 2B,C). This visualization provides a clear overview of the interactions between key targets and offers valuable insights for further investigation into the molecular mechanisms linking PCD and NOA.

After excluding the core gene ATG5 due to a lack of valid instrumental variables, a two-sample MR analysis was performed on the remaining nine candidate genes. The results (Table 1) demonstrated a significant causal relationship between genetically proxied TLR4 expression and an increased risk of MI (Odds Ratio [OR]: 1.16; 95% Confidence Interval [CI]: 1.01–1.33; p = 0.032). In contrast, the other eight candidate genes showed no significant association with MI susceptibility (all p > 0.05). Sensitivity analyses confirmed the robustness of the primary finding for TLR4 (Table 1). Specifically, the MR-Egger intercept test indicated no evidence of directional pleiotropy (p = 0.66), Cochran’s Q test showed no significant heterogeneity (p = 0.99), and leave-one-out analysis demonstrated that the causal estimate was consistent and not driven by any single influential genetic variant. Consistent with the MR results, TLR4 expression was found to be significantly upregulated in NOA samples across discovery queue (GSE) (Figure 3A) and independent validation cohorts (GSE4797 and GSE145467) (Figures 3B,C). Based on this converging evidence, TLR4 was selected as the primary focus for subsequent analyses.

TLR4 expression was assessed using integrated datasets, immunohistochemistry data, and single-cell RNA sequencing data from the Human Protein Atlas (HPA) database. Our analysis revealed that TLR4 is primarily expressed within the seminiferous tubules, with lower expression levels detected in other testicular cellular compartments (Figure 3D). Consistent with this spatial localization, single-cell transcriptomic data demonstrated that in healthy testicular tissue, TLR4 expression is predominantly enriched in macrophages and endothelial cells. Moderate expression was observed in Sertoli cells, while low-level expression was detected in spermatocytes and spermatogonia. TLR4 expression was scarcely detectable in Leydig cells, peritubular cells, as well as in early and late spermatids (Figures 3E,F).

To further evaluate the role of TLR4 in the context of NOA, we established an LPS-induced inflammatory injury model using GC-1 spg (ts) spermatogonial cells (designated as the Im1spd model) to simulate the inflammatory microenvironment of NOA, following established methodology. Western blot analysis confirmed significantly increased TLR4 expression in the Im1spd model compared to control cells (Figure 4A). Subsequent functional assays demonstrated that TLR4 inhibit effectively attenuated the inflammatory damage in GC-1 spg (ts) cells: EdU incorporation assays revealed significantly protected proliferative capacity against inflammation-induced injury (Figures 4B,C). Collectively, these data highlight TLR4 as a key mediator in the pathogenesis of NOA.

To identify potential natural bioactive compounds for the treatment of NOA, we focused on TLR4-associated natural products compiled in the CTD. Our screening identified seven common natural active compounds related to TLR4: curcuNOAn, ginsenoside Rg3, resveratrol, melatonin, baicalein, icariin B, and quercetin. These compounds are well-known for their notable pharmacological properties, including anti-inflammatory, antioxidant, and anti-apoptotic activities. Subsequent molecular docking analysis confirmed that all seven natural products could bind stably to TLR4, with binding energies (ΔG) below −5 kcal/mol. The specific binding energies were as follows: curcumin (−8.4 kcal/mol; Figure 6A), ginsenoside Rg3 (−8.1 kcal/mol; Figure 6B), resveratrol (−7.8 kcal/mol; Figure 6C), baicalein (−9.2 kcal/mol; Figure 6D), melatonin (−7.1 kcal/mol; Figure 6E), icariin B (−8.1 kcal/mol; Figure 6F), and quercetin (−9.0 kcal/mol; Figure 6G). Consistent with our computational predictions, Western blot analysis confirmed that baicalein treatment significantly reduced the LPS-induced expression of both TLR4 and NF-κB p65 (Figure 6H; Supplementary Figure S1C) and alleviated LPS-induced proliferation arrest (Supplementary Figures S1B,C). Together, these findings highlight the potential of natural bioactive compounds such as baicalein in developing novel therapeutic strategies for NOA.