We searched gene expression profiling datasets for asthenozoospermia using the term “asthenozoospermia” from the Gene Expression Omnibus (GEO) publicly available database (https://www.ncbi.nlm.nih.gov/gds/). The inclusion criteria for the dataset are as follows: (1) the dataset includes whole genome-wide expression mRNA microarray data; (2) the samples were sperm from normozoospermic and asthenozoospermic; (3) the organism was restricted to Homo sapiens. Finally, multiple gene expression datasets, including GSE22331 (1 normozoospermia and 1 asthenozoospermia), GSE34514 (4 normozoospermias and 4 asthenozoospermias), and GSE160749 (6 normozoospermias and 5 asthenozoospermias), were downloaded and integrated for further analysis, which included a total of 11 normozoospermias and 10 asthenozoospermias ( Supplementary Table 1 ). Firstly, we conducted independent data cleaning and preprocessing on each dataset, including outlier detection, missing value filling, etc., to ensure data consistency and integrity. Additionally, the removeBatchEffect function of the R software (version 4.2.1) “limma” package was utilized for batch correction to minimize the batch effect. By reducing the dimensionality of the data, Principal Component Analysis (PCA) allowed us to visualize the overall patterns and relationships within the dataset ( Supplementary Figure S1 ). This analysis helped confirm that the integrated data remained reliable and informative, despite the potential confounding effects of batch variations.

DEGs were identified using the “limma” packages of the R language, with a significance threshold of p-value < 0.05 and |log2FC| > 0.5. Volcano plots were generated using the SangerBox online software (http://sangerbox.com/home.html) based on the “ggplot” package of R language. Venn diagrams (https://jvenn.toulouse.inrae.fr/app/example.html) were used to identify the key DEGs by intersecting the DEGs obtained from GEO databases with our previous experimental DEGs (1,128 up-regulated and 1,210 down-regulated) (13).

To explore the biological functions and pathways associated with the key DEGs, we conducted GO and KEGG analyses using the “AnnotationDbi”, “clusterProfiler”, and “ggplot2” R packages. A significance threshold of p.adjust < 0.05 was applied to determine significant terms. The results were visualized to provide insights into the enriched biological processes and pathways associated with the DEGs.

A total of 5850 genes with a variance greater than 0.5 in the integrated dataset (11 normozoospermic and 10 asthenozoospermic individuals) were selected for WGCNA. First, we calculated the similarity between samples using horizontal hierarchical clustering. Then, we constructed an adjacency matrix and determined the optimal soft threshold power (β = 14) based on scale independence and mean connectivity (scale free R^{2 =} 0.80). Next, we transformed the adjacency matrix into a topological overlap matrix (TOM) to assess network interconnectedness. The dynamic tree cutting method was applied to identify co-expressed gene modules, and a cluster dendrogram was generated to visualize the modules represented by different colors using hierarchical clustering. To ensure reliable results, we set the minimum number of genes per module to 50 and the module cutting height to 0.1. Pearson correlation analysis was performed to calculate the correlations between module eigengenes (MEs) and phenotypes, identifying the module with the highest positive correlation with asthenozoospermia. The genes most associated with the darkolivegreen and skyblue modules and asthenozoospermic traits were selected for further analysis. Additionally, gene significance (GS) and module membership (MM) values were generated, and genes with |MM| > 0.8 and |GS| > 0.4 were considered hub genes in the key module.

The gene that exhibited the most significant difference, which was found at the intersection of DEGs and the hub genes identified by WGCNA, was selected as the core gene for further investigation. To visualize the expression patterns of the intersecting genes between DEGs and the hub genes identified by WGCNA, we generated violin plots using SangerBox software (http://sangerbox.com/home.html) based on datasets and our previous experimental data. These plots provided a graphical representation of the expression levels of the selected genes, allowing for a better understanding of their differential expression in asthenozoospermic and normozoospermic groups.

To gain further insights into the functional implications of hub gene, we performed single-gene GSEA and visualization using the R packages “clusterProfiler”, “enrichplot”, and “stringr”. Firstly, we divided all samples into low-expression and high-expression groups based on the expression levels of the key gene. Next, we conducted GSEA to identify significantly different KEGG pathways between the two groups. Pearson correlation analysis was performed to assess the correlation between the hub gene and all genes in the integrated dataset, and the genes were ranked in descending order based on the correlation scores. The ranked genes were then used as the gene set to be tested, while the KEGG signaling pathway served as the predefined gene set. Finally, we examined the enrichment of the KEGG signaling pathway in the corresponding set of genes to be tested, considering a p.adjust value less than 0.05 as the cutoff criterion for statistical significance. This analysis provided valuable information about the potential functional associations and pathways related to hub gene.

To explore the potential association between the hub gene and mitochondrial autophagy, we obtained a set of 29 mitochondrial autophagy-related genes (including ATG12, FUNDC1, SQSTM1, TOMM5, UBB, ATG5, MAP1LC3A, PGAM5, SRC, UBC, MAP1LC3B, PINK1, TOMM20, ULK1, CSNK2A2, MFN1, PARK2, TOMM22, TOMM70A, VDAC1, CSNK2B, MFN2, RPS27A, TOMM40, and UBA52) from the Reactome Pathway Database (https://reactome.org). We then performed Pearson correlation analysis to assess the correlation between the key gene and these mitochondrial autophagy-related genes. The resulting correlation heatmap was visualized using the SangerBox software. p less than 0.05 was considered as cut-off criterion for statistically significant.

To investigate the potential regulatory interactions between the key gene and non-coding RNAs, we utilized online databases such as miRcode (http://www.mircode.org/), TargetScan (https://www.targetscan.org/vert_80/), and starBase (https://starbase.sysu.edu.cn/), along with our previous research findings. These resources allowed us to predict miRNAs and lncRNAs that may be associated with the key gene. We then constructed lncRNA-miRNA-mRNA networks by integrating the lncRNA-miRNA and mRNA-miRNA interactions. The resulting networks were visualized using Cytoscape software (version 3.10.0), providing a comprehensive view of the potential regulatory relationships between the key gene and non-coding RNAs.

We collected semen samples from 3 asthenozoospermia males and 3 normozoospermia males at the Hainan Women and Children’s Medical Center. The collection methods, diagnosis, and inclusion criteria followed the protocols outlined in our previous study (13). To isolate seminal plasma exosomes, the semen samples were rapidly thawed at 37°C and transferred to a new centrifuge tube. They were then centrifuged at 4°C, 2000 ×g for 30mins. The supernatant was carefully transferred to another tube and subjected to a second centrifugation at 4°C, 10,000 ×g for 45mins to remove larger vesicles. The resulting supernatant was filtered using a 0.45μm filter membrane, and the filtrate was collected. Subsequently, the filtrate was transferred to a new centrifuge tube and centrifuged at 4°C, 100,000 ×g for 70mins. After removing the supernatant, the pellet was resuspended in 10 mL of pre-cooled 1× PBS and subjected to another centrifugation at 4°C, 100,000 ×g for 70mins. Finally, the supernatant was discarded, and the exosome pellet was resuspended in 100μL of pre-cooled 1× PBS and stored at -80°C. To examine the structural characteristics and morphology of the isolated exosomes, we performed high-resolution imaging using transmission electron microscopy (TEM, HT-7700, Hitachi) ( Supplementary Figure S2A ). Additionally, the size distribution and concentration of the exosomes were determined using a NanoFCM N30E particle size analyzer ( Supplementary Figures S2B, C ).

Total RNA was extracted from exosomes of asthenozoospermia and normozoospermia using TRIZOL reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s instructions. The concentration and purity of RNA were assessed using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) by measuring the optical density (OD) at 280 and 260 nm. The integrity of the RNA was evaluated using agarose gel electrophoresis. Subsequently, the RNA was reverse transcribed into complementary DNA (cDNA) using a reverse transcription kit (Nocoprotein Co., Ltd., Suzhou, China). For qRT-PCR, the SYBR Green Master Mix (Applied Biosystems, MA, USA) and the Applied Biosystems 7500 sequence detection system were used. The cycling conditions consisted of an initial denaturation step at 95°C for 1 minute, followed by 45 cycles of denaturation at 95°C for 20 seconds, and annealing/extension at 60°C for 45 seconds. U6 was employed as an endogenous control for miRNA, while GAPDH served as the internal control for mRNA and lncRNA. The relative gene expression levels were calculated using the 2^−ΔΔCt method. The primer sequences used for qRT-PCR are provided in Table 1 .

The relative quantitative detection of m6A methylation involved the following steps: 1) Total RNA extraction and quality control: Total RNA was extracted as previously described, and quality control measures were performed to ensure RNA integrity. 2) RNA purification and fragmentation: The extracted total RNA was treated with DNase I to remove any trace genomic DNA (gDNA). Exogenous positive and negative m6A-modified fragments were introduced as controls. The RNA was then fragmented into fragments ranging from 60 to 200 base pairs (bp) by fragmenting reagents in the presence of metal ions. 3) Antibody immunoprecipitation: Dynabeads Protein A was mixed with the anti-m6A antibody to form complexes. The fragmented RNA was incubated with the bead-antibody complex, allowing for the binding of m6A-modified RNA fragments. The beads were separated from the solution using a magnet, and then washed to remove non-specifically bound RNA fragments. The m6A-modified RNA fragments were extracted using a lysis buffer. 4) Reverse transcription of RNA: A reverse transcription reaction was performed using random primers to synthesize single-stranded cDNA. This cDNA served as a quantitative template for subsequent analysis. 5) qRT-PCR: qRT-PCR was carried out using m6A-modification site-specific primers. SYBR fluorescent dye was used for detection, and ROX dye was added for inter-well signal correction. The reaction was amplified on a Real-Time PCR instrument, and the amplification and melt curves were analyzed. The primer sequences are provided in Table 2 .

Normalization and integration of the mRNA expression datasets GSE22331, GSE34514, and GSE160749 resulted in a total of 17,415 genes. A total of 2,070 DEGs were identified between asthenozoospermic and normozoospermic groups, including 159 upregulated genes and 1,911 downregulated genes. The volcano plot visualized these results, with each dot representing a gene ( Figure 2A ). Additionally, key DEGs (containing 7 upregulated genes and 55 downregulated genes) were identified by finding the intersection between the DEGs from the datasets and previous experiments ( Figure 2B ).

GO and KEGG analyses were performed to determine the potential biological functions and pathways associated with the key DEGs ( Figure 3 ). The GO analysis revealed that most of the DEGs were located in the ribosome (organellar, mitochondrial), ribosomal subunit, outer membrane (mitochondrial and organelle), and calpain complex ( Figure 3A ), which were involved in processes such as rRNA modification and methylation, membrane raft assembly, regulation of flagellated sperm motility, sperm individualization, cellularization ( Figure 3B ). Furthermore, many DEGs were found to be involved in the regulation of enzyme activity, including rRNA methyltransferase, lysophosphatidic acid acyltransferase, and glucosamine kinase ( Figure 3C ). The KEGG analysis showed that the key DEGs were mainly associated with glycerophospholipid metabolism, ubiquitin-mediated proteolysis, and protein processing in the endoplasmic reticulum ( Figure 3D ).

A total of 5,850 genes with larger variance were selected for co-expression analysis based on the variance results of gene expression. Using the Pearson’s correlation coefficient, all samples from GSE22331, GSE34514, and GSE160749 were clustered, resulting in a sample clustering tree ( Figure 4A ). With a threshold correlation coefficient of 0.80 and a soft threshold of 14, a scale-free topology module was constructed ( Figure 4B ). Based on average hierarchical clustering and dynamic tree clipping, thirty co-expression modules were created ( Figure 4C ). Correlation analysis between the modules and clinical traits revealed that the darkolivegreen (r = 0.481, p = 0.027) and skyblue (r = 0.442, p =  0.045) modules showed modest significant correlations with asthenozoospermia ( Figure 4D ). GS and MM analysis identified the darkolivegreen (correlation coefficient = 0.440, p  = 9.8 × 10⁻¹²; Figure 4E ) and skyblue (correlation coefficient = 0.460, p  = 1.2 × 10⁻¹¹; Figure 4F ) modules as important modules strongly correlated with asthenozoospermia for further analysis. According to the cutoff criteria (|MM|> 0.8 and |GS|> 0.4), the skyblue module contained 196 genes, of which 93 genes were identified as hub genes, and 141 hub genes were identified in the darkolivegreen module (218 genes) ( Supplementary Table 2 ). A Venn diagram was then developed to take the intersection between the above genes and key DEGs to identify hub genes in asthenozoospermia, which were ZNF764, COQ9 and CA5A, respectively ( Figure 4G ).

The expression levels of hub genes in asthenozoospermia and normozoospermia were validated using the integrated GO datasets (GSE22331, GSE34514, and GSE160749) and previous experimental data ( Figure 5A ). The expression level of ZNF764 was found to be higher in asthenozoospermic patients compared to normozoospermic individuals (p = 4.9e-3), while the expression of COQ9 was significantly decreased in asthenozoospermic patients (p = 3.4e-4). However, no significant differences were observed in the expression level of CA5A between the two groups (p = 0.16). Furthermore, the expression of these two genes (ZNF764 and COQ9) in seminal exosomes of asthenozoospermic and normozoospermic individuals was validated through RT-PCR ( Figure 5B ).

Furthermore, we performed GSEA to gain a more accurate understanding of the functional roles of the hub genes. Based on the designated cut-off criteria, we identified the top twenty KEGG pathways enriched in the analysis. For ZNF764, the enriched pathways included nucleocytoplasmic transport, protein processing in endoplasmic reticulum, ubiquitin mediated proteolysis, protein export, mRNA surveillance pathway, autophagy-animal, spliceosome, polycomb repressive complex, and cell cycle, etc. ( Figure 5C ). For COQ9, the enriched pathways included autophagy, nucleocytoplasmic transport, protein processing in endoplasmic reticulum, mRNA surveillance pathway, ubiquitin mediated proteolysis, ATP-dependent chromatin remodeling, polycomb repressive complex, lysosome, endocytosis, and cell cycle, etc. ( Figure 5D ).

We analyzed the correlation between expression levels of ZNF764 and COQ9 and the expression levels of 22 mitochondrial autophagy-related genes ( Table 3 ; Supplementary Figure S2 ). Significant negative correlations were found between the expression level of ZNF764 and three mitochondrial autophagy-related genes, including CSNK2A2 (cor = -0.792, p = 0.006), MFN1 (cor = -0.789, p = 0.007), and PARK2 (cor = -0.775, p = 0.009). On the other hand, the expression level of COQ9 showed significant positive correlations with 16 mitochondrial autophagy-related genes, including SQSTM1 (cor = 0.955, p < 0.001), TOMM5 (cor = 0.903, p < 0.001), UBB (cor = 0.772, p = 0.009), ATG5 (cor = 0.744, p = 0.014), MAP1LC3A (cor = 0.914, p < 0.001), PGAM5 (cor = 0.810, p = 0.005), UBC (cor = 0.804, p = 0.005), MAP1LC3B (cor = 0.989, p < 0.001), PINK1 (cor = 0.942, p < 0.001), TOMM20 (cor = 0.744, p = 0.014), PARK2 (cor = 0.858, p = 0.001), TOMM22 (cor = 0.875, p = 0.001), TOMM70A (cor = 0.846, p = 0.002), VDAC1(cor = 0.750, p = 0.013), CSNK2B (cor = 0.805, p = 0.005), MFN2 (cor = 0.948, p < 0.001), and TOMM40 (cor = 0.772, p = 0.009). Additionally, the expression level of COQ9 was significantly negatively correlated with the expression level of SRC (cor = -0.733, p = 0.016). Therefore, COQ9 was selected for further analysis.

We utilized miRcode, TargetScan, and starBase databases to predict miRNAs that are potentially associated with COQ9 and took intersections, including miR-125a-5p, miR-125b-5p, miR-148a-3p, miR-148b-3p, and miR-4139 ( Figure 6A ). Furthermore, lncRNAs potentially associated with these miRNAs were predicted using the starBase database, and a ceRNA regulatory network was constructed ( Figure 6B ). The expression levels of miR-125a-5p and LINC00893 in exosomes of asthenozoospermic and normozoospermic individuals were detected, revealing higher expression of miR-125a-5p in asthenozoospermia and significantly decreased expression of LINC00893 (p < 0.001, Figures 6C, D ). Additionally, the expression of m6A methylation level of LINC00893 was found to be significantly increased in patients with asthenozoospermia (p < 0.05, Figure 6E ).