We searched “asthenozoospermia” or “chronic epididymitis” in the GEO database (https://www.ncbi.nlm.nih.gov/geo/). After careful screening, the AZS datasets (GSE160749 and GSE22331) and CE dataset (GSE199903) were downloaded from the GEO database. GSE160769 contained genetic information on five AZS patients and six healthy controls. The GSE22331 dataset contained genetic information from 30 AZS patients and 30 healthy controls. Transcriptome data from an independent validation set including 3 AZS patients and 9 healthy controls were obtained from GSE26881. All AZS samples did not exclude other relevant description. Additionally, the single-cell transcriptome of the testis was obtained from GSE149512. Detailed information on the above GEO datasets is shown in Supplementary Table S1.

All expression matrices were regularly normalized and log2 transformed (if needed) by using the R package “limma” (Ritchie et al., 2015). For obtain the integrated AZS dataset, two GEO datasets (GSE160749 and GSE22331) were continually merged using Strawberry-Perl software (Version 5.30). The batch effect between arrays was then eliminated using the R package “sva” (Leek et al., 2012). The differentially expressed genes (DEGs) between the AZS and control groups were calculated based on |log2-fold change| (|log2FC|) > 0.5 and p-value <0.05. The same analytical approach was applied to DEGs between the CE and control groups. Then, the intersection of AZS DEGs and CE DEGs was performed to obtain intersecting DEGs involved in the two diseases.

To further depict the biological function and pathways implicated in the intersected genes, we performed Gene Ontology (GO) analysis and REACTOME pathway analysis by using the DAVID Bioinformatics Microarray online tool (https://david.ncifcrf.gov/) (Huang et al., 2007). GO analysis includes biological process (BP), cellular component (CC), and molecular function (MF). Visualization of GO clusters and pathway clusters was performed by R package “ggplot2” (Walter et al., 2015). p-value <0.05 was set as the cutoff criterion.

The protein interaction network was analyzed via the Search Tool for the Retrieval of Interacting Genes (STRING) database (http://string-db.org). For more potential interactions, a confidence score >0.15 was set as the cutoff criterion. Visualization of the PPI network was carried out using Cytoscape software (Version 3.7). Moreover, the MCODE plugin was used to obtain the top three gene clusters, and the cytoHubba plugin was used to obtain the top 20 genes with the highest MCC score.

We first narrowed down the number of intersected genes, and defined the co-DEGs as the intersected genes with the same expression trend (upregulation and downregulation) in the AZS and CE datasets. Subsequently, to assess the efficacy of these co-DEGs in predicting AZS and CE, we created receiver operator characteristic (ROC) curves using the R package “timeROC” (https://cran.r-project.org/web/packages/timeROC/) (Blanche et al., 2013) and evaluated the area under the curve (AUC) of these co-DEGs. The co-DEGs with AUCs >0.7 in both the AZS dataset and CE dataset were selected. Among them, hub genes were further selected and validated in the independent validation set (GSE26881) based on AUC >0.7. Visualization was carried out by the R package “ggplot2”.

The expression localization of hub genes was evaluated in immunohistochemical data downloaded from the Human Protein Atlas database (HPA) (https://www.proteinatlas.org/). Resource addresses of immunohistochemical data are shown in Supplementary Table S2.

To identify the role of hub genes in testicular function, we used GSE149512 which contained the single-cell transcriptome of testes from ten males ranging from infants to adults, and GSE211115 which contained the single-cell transcriptome of testes from nine male mice at different age stages. The expression pattern of hub genes was analyzed and visualized in the Male Health Atlas database (http://www.malehealthatlas.cn/) (Zhao L. et al., 2020).

Single-sample gene set enrichment analysis GSEA (ssGSEA), an extension of GSEA, computes enrichment scores separately for each sample-gene set pairing. This method is widely used in the calculation of immune infiltration in genetic expression datasets (Finotello and Trajanoski, 2018). We performed immune infiltration analysis for the AZS dataset and CE dataset, and calculated the enrichment score of each immune cell and function for each sample. Correlations among the above immune proportions and correlations between the hub genes and immune proportions were carried out using Spearman’s method. A p-value <0.05 was regarded as statistically significant.

For identification of the interacting genes of hub genes, a gene-gene interaction network was constructed via the online tool GeneMANIA (http://genemania.org/). The network showed the association in the fields of physical interaction, coexpression, predicted, colocalization, genetic interactions, pathway, and shared protein domains.

According to the five hub genes, the online tool miRcode (http://www.mircode.org/) (Jeggari et al., 2012) was used to predict target miRNAs in humans with default criteria. Visualization of the hub gene-miRNA network was carried out by Cytoscape software. With the use of the cytoHubba plugin, the top ten miRNAs with the highest MCC score were obtained.

To identify the association between hub genes and disease, we used the Comparative Toxicogenomics Database (CTD, http://ctdbase.org/) (Davis et al., 2017) to generate their relevancy. Male reproductive system diseases were screened, and the Inference Score was limited to >10.

Ten subjects attending Tongji Hospital were recruited in the present study. Among them, five patients were diagnosed with AZS and CE, and five subjects were healthy controls. According to the WHO laboratory manual (World Health Organization, 2021), sperm motility of less than 42% or less than 30% with progressive motility in at least two semen exams was required for AZS. Patients were diagnosed with chronic epididymitis based on the criteria requirements of CE (McConaghy and Panchal, 2016). The Institutional Research Ethics Committee of Huazhong University of Science and Technology approved our study (REC No. (2017(04)). We obtained written consent from each subject, and all procedures adhered to the Helsinki Declaration.

The medical histories of patients with AZS and CE and controls were collected. We excluded subjects with a history of varicocele, diabetes mellitus, scrotal trauma, and exposure to radiation and environmental pollutants. Chemiluminescence immunoassays (Beckman Coulter, Fullerton, United States) were used to measure the levels of FSH, LH, total T, and estradiol (E₂) in the serum. Seminal elastase was detected using enzyme-linked immunosorbent assays (ELISAs) (Milenia, Bad Nauheim, Germany). We performed an ultrasound examination for male genital condition (850 system, Hitachi, Beijing, China). Semen was harvested by masturbation after 3–5 days of abstinence, and sperm processing and testing were in strict compliance with the WHO laboratory manual (World Health Organization, 2021). We used a Percoll gradient (Solarbio, Beijing, China) and centrifugation for 20 min at 1,000 g to eliminate seminal plasma and purify sperm samples. We washed them with PBS twice, and kept sperm pellets for later measurement.

We used RNA-easyTM Isolation Reagent (Vazyme, Nanjing, China) to extract RNA from sperm samples according to the manufacturer’s instructions. And cDNA was synthesized with the Hifair^® Ⅱ 1st Strand cDNA Synthesis Kit (Yeasen, Shanghai, China). The detailed procedures of quantitative real-time polymerase chain reaction (qRT-PCR) were described previously (Chen et al., 2021). Then, primers were synthesized (Tsingke, Wuhan, China), and sequence information is shown in Supplementary Table S3. The 2^−ΔΔct relative quantification algorithm was utilized to evaluate the mRNA levels of hub genes in each sample.

The difference in mRNA levels was compared using Student’s t-test. All statistical analyses we used in the present study were carried out by R software (version 4.1.3, https://www.r-project.org/) and SPSS version 23.0 (IBM Corporation, Armonk, United States). Statistical significance was defined as a p-value less than 0.05.

The research design of this study is shown in Figure 1. With the criteria of p-value <0.05 and |log2FC| > 0.5, a total of 2,256 DEGs were selected from the AZS dataset: 800 upregulated DEGs and 1,456 downregulated DEGs. In the CE dataset, a total of 1742 DEGs were selected: 1,271 upregulated DEGs and 471 downregulated DEGs. The expression patterns of both DEG sets are shown in volcano plots and heatmaps (Figures 2A–D). To obtain the DEGs that may participate in the development of these two diseases, we took the intersection of AZS DEGs and CE DEGs. A Venn diagram showed that 109 DEGs were expressed simultaneously in the two DEG groups (Figure 2E).

To investigate inherent biological information in more detail, we performed GO and pathway analysis for the intersected DEGs. As shown in Figure 3A, for biological process, GO terms were enriched in inflammatory response, immune response, and cell adhesion. For cellular components, GO terms were enriched in extracellular region, plasma membrane, and extracellular space. For molecular function, GO terms were enriched in protein binding, transmembrane signaling receptor activity, and oligosaccharide binding. In addition, REACTOME pathway analysis showed that the DEGs were mostly enriched in the immune system, neutrophil degranulation, and interleukin-4 and interleukin-13 signaling (Figure 3B). The functional enrichment results indicated that the 109 intersecting DEGs had a close relationship with immune-related biological mechanisms.

To further explain the relationships between the intersecting genes, we used the STRING database to identify protein interactions. Figure 3C displays the PPI network, which consisted of 587 interactions. Moreover, the top three DEG clusters were calculated based on the MCODE algorithm (Figure 3D). According to the MCC algorithm, the top 20 DEGs were obtained to show their interaction in Figure 3E. Herein, CMKLR1, CD300LB, IL1B, and C3 (complement C3) were found in this top 20 DEGs interaction, and these genes are strongly associated with the activation of immune responses and the release of inflammatory mediators.

Although we screened out a batch of DEGs, the number in this batch was too large. To narrow down the range of the DEGs, we further selected 47 coexpressed DEGs (co-DEGs) with the same expression trend (upregulation and downregulation) in the AZS dataset and CE dataset (Figure 3F). With the use of the R package “pROC,” we explored the above genes with high accuracy of the diagnostic characteristics. The ROC results of the AZS dataset showed that five genes (CD300LB, CMKLR1, CCR4, B3GALT5, and CTSK) had good diagnostic values (AUC >0.7) to discriminate the samples from the patients with AZS and healthy controls. Similarly, the diagnostic values of the five genes were also good in the CE dataset (AUC >0.7) (Figures 4A,B). Based on their good diagnostic performance, these five genes were regarded as hub genes in the present study.

To further evaluate the important role of the five hub genes, we downloaded another AZS dataset as the independent validation set. After ROC evaluation, Figure 4C shows that the AUC values of the five hub genes were all greater than 0.7 in the validation set. These results indicated that the five hub genes have the potential to be biomarkers for AZS and CE.

To better reveal the condition of expression and localization for hub genes, we first explored the hub gene signature in the immunohistochemical data from the HPA database. As shown in Figure 5A, weak positive expression for CMKLR1 was mainly concentrated in the interstitial region and blood vessels beside spermatogenic tubules. For CCR4, weak positive staining was observed in the cytoplasmic region of germ cells near the basement membrane and Leydig cells. For B3GALT5, moderate positive expression was observed in the nuclei of germ cells near the basement membrane.

We investigated the expression pattern of cell types in healthy male controls with normal fertility. Various cell types in the testes from humans and mice were carefully classified according to a single-cell analysis algorithm (dimensional reduction by principal component analysis, cell clustering by FindClusters processing, and cluster annotation based on marker genes) (Figures 5B, D), and we found that hub genes had unique expression patterns in cell types. As shown in Figure 5C (human data), CMKLR1 and CD300LB were mainly enriched in macrophages due to their high average expression. CCR4 and B3GALT5 were both enriched in spermatogonia, while CTSK was predominantly enriched in Leydig cells and peritubular myoid cells. In Figure 5E (mouse data), Cd300lb and B3galt5 were strongly enriched in immune cells, and Cmklr1 and Ccr4 were preferentially enriched in spermatids/sperm. Ctsk was highly enriched in peritubular myoid cells. The data demonstrated that the hub genes may be closely related to the immune environment of the interstitial region and the entire process of spermatogenesis.

MiRNAs have been revealed to play a critical role in male-related infertility disorders. Therefore, we selected the important miRNAs around the five hub genes in CE and AZS, and constructed a miRNA-hub gene regulatory network. A total of 224 miRNA-gene interactions were found from the miRcode database, and the miRNA profile is depicted in Figure 6A. Among these miRNAs, we selected the top ten important miRNAs that have more links to hub genes: hsa-miR-873, hsa-miR-149, hsa-miR-197, hsa-miR-4644, hsa-miR-4306, hsa-miR-3473, hsa-miR-882, hsa-miR-185, hsa-miR-544-3p, and hsa-miR-544ab. A total of 42 miRNA-hub gene interactions were identified in the network (Figure 6B). Among the top ten miRNAs, hsa-miR-149 and hsa-miR-873 were both linked to all five hub genes, while each of the other eight miRNAs was linked to four hub genes. These miRNAs may function as regulators in the development of AZS and CE.

Utilizing the CTD database, we demonstrated that the five hub genes (CD300LB, CMKLR1, CCR4, B3GALT5, and CTSK) significantly affected male reproductive system diseases. The male reproductive system diseases were focused on testicular diseases, infertility (male), asthenozoospermia, oligospermia, and teratozoospermia (Figure 6C).

Before analyzing immune infiltration, we analyzed the correlation between hub genes and their interacting genes. The gene‒gene interaction results indicated that their gene functions were enriched in mononuclear cell migration, leukocyte chemotaxis, interleukin-6 production, leukocyte migration, granulocyte migration, and so on (Figure 7A).

Single-sample GSEA is a very well-established computational algorithm for calculating immune infiltration for genetic matrices. After calculating the score of 29 immune cells/functions for each sample, we identified the correlation between the 29 immune cells/functions in AZS, which is in Figure 7B, while the correlation in CE is shown in Figure 7C. In the CE dataset, most immune cells/functions had a positive influence on each other. The situation (positive influence) was the same for AZS, except for macrophages, Th1 cells, and NK cells.

Next, we further explored the correlation between the hub genes and immune infiltration. As shown in Figure 7D, for AZS samples, hub genes were significantly linked to the level of CD8⁺ cell infiltration (CMKLR1, r = 0.67, p = 0.015; CCR4, r = 0.66, p = 0.017; B3GALT5, r = 0.87, p = 0.00012) and T helper cell infiltration (CMKLR1, r = 0.68, p = 0.014; CCR4, r = 0.58, p = 0.043; B3GALT5: r = 0.64, p = 0.022). For CE samples, hub genes had an evident correlation with CD8⁺ cell infiltration (CMKLR1, r = 0.7, p = 0.031; CCR4, r = 0.65, p = 0.049; CTSK: r = 0.82, p = 0.0068; CD300LB, r = 0.76, p = 0.016) and T helper cells (CMKLR1, r = 0.83, p = 0.0056; CCR4, r = 0.81, p = 0.0082; CTSK: r = 0.9, p = 0.00088; CD300LB, r = 0.72, p = 0.024). These results suggested that CD8⁺ cells and T helper cells may act as important immune cells regulating the progression of inflammation in CE and impaired sperm function.

In the next step, we performed qRT-PCR to examine the expression of hub genes between the patients with AZS and CE and healthy controls. As shown in Supplementary Table S4, patients with AZS and CE had decreased levels of sperm progressive motility and increased levels of seminal elastase compared with the control group. Quantitative RT-PCR results indicated that there were evident differences in hub genes between the patient and control groups (p < 0.05), and the mRNA levels of hub genes (CD300LB, CMKLR1, CCR4, and B3GALT5) were significantly elevated in the sperm affected by AZS and CE (Figure 8). However, the mRNA level of CTSK was not statistically significant. The data synergistically validated that these hub genes shared the same expression trend in the sequencing results and experimental results.