Infertility network and hub genes for nonobstructive azoospermia utilizing integrative analysis

Overview of the transcriptomes of NOA and GO enrichment analysis

In order to elucidate the molecular pathogenesis of NOA, we firstly screened the DEGs between normal testicular biopsies and NOA cases using the limma package and then built the volcano map and heat map using the ggplot2 and pheatmap packages. The heat map of all mRNAs is shown in Supplementary Figure 1A, where significant differences in gene expression between the two groups can be seen. Subsequently, the volcano map was constructed using differentially expressed mRNAs (|log2fold change|≥2, p-value < 0.05) as shown in Supplementary Figure 1B. In total, 214 mRNAs displayed the differential expressions in NOA, including two upregulated and 212 downregulated mRNAs (|log2fold change|≥2, p-value < 0.05) (Supplementary Table 1). In addition, from the heat map and the cluster dendrogram (Supplementary Figure 3A), it can be seen that the normal testicular biopsies and NOA cases were allocated into two distinct main clusters (control and NOA).

After Metascape analysis, the Top 10 clusters with their representative enriched terms are shown in Figure 1A and Supplementary Table 2. GO analysis results were enriched in Biological Process, Cellular Components, and Molecular Functions. Biological Process clusters mainly include spermatogenesis, fertilization, cilium movement, penetration of zona pellucida, and sperm chromatin condensation. Cellular Components consist mainly of sperm part, motile cilium, and sperm fibrous sheath. For Molecular Functions, the main component is lysozyme activity. To further capture the relationships between the terms, a subset of enriched terms was selected and rendered as a network plot (Supplementary Figure 1C), where terms with a similarity > 0.3 were connected by edges. We selected the terms with the best p-values from each of the 20 clusters, with the constraint that there were no more than 15 terms per cluster and no more than 250 terms in total. From the network plot, the selected DEGs were closely related to spermatogenesis, fertilization, cilium movement, sperm part, motile cilium, and sperm fibrous sheath, further indicating the important role of DEGs in spermatogenesis-related bioprocesses. Additionally, KEGG analysis was conducted (Supplementary Table 3), and was mainly enriched in Proximal tubule bicarbonate reclamation, Aldosterone synthesis and secretion, Glycolysis/Gluconeogenesis, and Phosphonate and phosphinate metabolism. When combined with the above GO analysis, no pathways overlapped with the GO results.

Figure 1. GO enrichment and Gene Set Enrichment Analysis of NOA, and validation of hub genes. (A) Top 10 clusters with their representative enriched term based on GO enrichment analysis of DEGs. (B) NOA samples were correlated positively with gene signatures related to Steroid Hormone Biosynthesis and Ubiquitin Mediated Proteolysis pathways. (C) NOA samples were correlated negatively with gene signatures related to Glycerophospholipid Metabolism and Proximal Tubule Bicarbonate Reclamation pathway. (D) The normalized expression validation of hub genes using the GSE45887 dataset. The Y-axis expressions were normalized by log₂(TPM/10+1). (E) The normalized expression validation of hub genes using the GSE106487 dataset. The Y-axis expressions were normalized by log₂(TPM/10+1).

Gene set enrichment analysis of NOA

KEGG analysis based on the GSEA analysis method showed that the upregulated DEGs were significantly enriched in Steroid Hormone Biosynthesis and Ubiquitin Mediated Proteolysis pathways (Figure 1B), and the downregulated DEGs were enriched in Glycerophospholipid Metabolism and Proximal Tubule Bicarbonate Reclamation pathway (Figure 1C). When combined with previous data (Supplementary Table 3), the Proximal Tubule Bicarbonate Reclamation pathway may be of great importance in NOA pathogenesis.

Protein-protein interaction of NOA

To further investigate all the DEGs (Supplementary Table 1) and the potential protein levels, the STRING database was applied to reveal the core PPI network. As shown in Supplementary Figure 2A, the core PPI network was constructed by molecular action (confidence threshold: 0.4), including 199 nodes and 297 edges. These selected DEGs were found to interact closely with biological processes or cellular components, including sexual reproduction (48 genes, nodes colored with green), spermatogenesis (38 genes, nodes colored with red), and sperm flagellum (14 genes, nodes colored with blue). Subsequently, we reconstructed the PPI networks of all DEGs by molecular action (confidence threshold: 0.7) for the above biological processes or cellular components (Supplementary Figure 2B) to compare the two methods. Based on a higher confidence threshold, there were only 37 edges, which indicated that an even higher confidence threshold could further narrow the positive target interaction. Of particular note, even using the higher confidence threshold, the SPTBN2 gene still existed in PPI networks, which may suggest an underlying role for the SPTBN2 gene in NOA pathogenesis.

Weighted gene correlation network analysis of NOA

To precisely clarify the key modules and hub genes in NOA, WGCNA was used to identify clusters (modules) of highly correlated genes and correlation networks of NOA. GSM1118241 (control) and GSM1118246 (NOA) were excluded from the analysis after quality assessment (Supplementary Figure 3A). The power of β was set at 24 to ensure a scale-free network and the mean connectivity remained normal (Supplementary Figure 3B), which indicated a good scale-free topology index of our dataset analysis. Gene modules were calculated, and gray modules represent genes that cannot be clustered into any other modules (Figure 2A and Supplementary Figure 3C). Furthermore, 44 gene modules were identified by the hierarchical clustering dendrogram and Module eigengene adjacency heatmap (Figure 2A, 2B). The interactions between gene modules were then analyzed and a TOM plot of the gene network was generated based on the corresponding hierarchical clustering dendrogram and modules (Figure 2C).

Figure 2. Weighted gene co-expression network analysis (WGCNA) of genes between the control and NOA groups. (A) Hierarchical cluster tree showing co-expression modules identified by WGCNA. The X-axis represents genes in different modules marked with different colors. (B) Heatmap plot of the adjacencies in the eigengene network. Each row and column in the heatmap correspond to one module eigengene (labeled by color). In the heatmap, the blue color represents low adjacency (negative correlation), while red represents high adjacency (positive correlation). Squares of red color along the diagonal are the meta-modules. (C) Heat map plot shows the topological overlap matrix (TOM) among 400 randomly selected genes. Light colors depict a small overlap, and the red color indicates a greater overlap. The left side and the top side show the gene dendrogram and module assignment. (D, E) Two modules with the highest relative rates were selected to localize the hub genes, they are separately colored with black (D) and pale turquoise (E). The Y-axis expressions were normalized by log₂(TPM/10+1).

Identification and validation of the hub genes of NOA

As WGCNA generated a huge gene network, we narrowed two modules (Figure 2D, 2E) for the network construction to localize the hub genes using the limma package. As a result, 179 nodes and 3306 edges were screened and used for the PPI network construction by molecular action (confidence threshold: 0.9) for further analysis (Supplementary Figure 1E). The MCODE plugin of Cytoscape was used to screen three subclusters (Supplementary Figure 1E and Supplementary Table 4) and calculate K-core values of each subcluster. After that, we extracted subcluster 1 and used Cytoscape to construct the core network (Supplementary Figure 1D). From the above two networks (Supplementary Figure 1D, 1E), CHD5 and SPTBN2 were preliminarily screened as the hub genes in NOA patients. Subsequently, we cross-referenced the gene set of each subcluster set separately with the previous total DEGs set (|log2fold change|≥2, p-value < 0.05) to screen for its common component. Ultimately, we screened to obtain these two hub genes (CHD5 and SPTBN2) with the highest differential expression in NOA patients.

In order to verify the accuracy of the differential expression of above two hub genes, we tested their expression using the other dataset; the results showed that the expression was significantly reduced in the GSE45887 dataset (Figure 1D). In addition, while the current study was being performed, Chen et al. used Whole-exome sequencing to compare the whole genomes of NOA patients with healthy people, and the results show that the CHD5 gene is one of the infertility causative genes of NOA [15], which also demonstrates the reliability of the method used in this study.

Integrative analysis of hub genes with single cell transcriptome

In this section, we integrated the scRNA-seq dataset (GSE106487) with the hub genes obtained from a previous screen for integrated analysis to further clarify the differential expression of these genes in specific cell types. The scRNA-seq data were first clustered and identified into different cell types (Figure 3A–3C). During spermatogenesis, all cell types (Figure 3D) and clusters (Figure 3E) from sperm stem cells to spermatocytes were identified by specific genetic markers. From Figure 3B, 3C, more than 95% of the cells in the NOA group were identified as Sertoli cells. In addition, we compared gene expression in the control and NOA groups to further validate the findings of the microarray (Figure 1E). Furthermore, the expression of CHD5 and SPTBN2 was significantly down-regulated from different cell types in the NOA group compared to the control group (Figure 3F), which further indicated the potential role of these two hub genes in NOA pathogenesis.

Figure 3. Identification and comparison of different cell types in different samples based on the single-cell RNA-seq dataset. (A) Cell cluster distribution in the TSNE plot. (B) Different types of cells and different samples identified in tSNE plots. SSC, spermatogonial stem cells; Differentiating SPG, differentiating spermatogonia; SPC7, cell mixtures including diakinesis, metaphase, anaphase, telophase, and secondary spermatocytes. (C) Comparison of different samples in tSNE plots. (D) Multiviolin plot of expressions of specific gene markers in identified cell types. X-axis represents diverse cell types. Y-axis represents the expression of specific gene markers. (E) Multiviolin plot of expressions of specific gene markers in different clusters. X-axis represents different clusters. Y-axis represents the expressions of specific gene markers. (F) The expression validation of hub genes from diverse cell types in the NOA group compared with the control group. The Y-axis expressions were normalized by log₂(TPM/10+1).

Data collection and preprocessing

Four normal spermatogenesis control and 27 NOA case microarray expression matrix files were downloaded from the National Center for Biotechnology Information (NCBI) Sequence Read Archive with the ID GSE45885. In addition, two expression matrix files with the ID GSE45887 (including four control and 16 NOA cases) and GSE106487 (consisting of two adult normal males, seven OA males, and one NOA male patient) were also prepared for subsequent validation and integrative analysis. The clinical information of testicular samples in GSE45885, GSE45887, and GSE106487 was obtained from published literature [8, 11, 16].

Differentially expressed gene screening

The limma package was used to determine differentially expressed genes (DEGs) between normal testicular biopsies and NOA cases under the threshold of p-value < 0.05: A linear model was simply fitted to the expression matrix and further analysis was then performed using empirical Bayes. Significant DEGs between different groups were identified using “|log2fold change|≥2” and “p-value < 0.05” as the judgment threshold. The ggplot2 package and pheatmap package in R were used to build the volcano map and heat map.

GO enrichment analysis and pathway enrichment analysis

Metascape analysis (http://metascape.org) was performed to accomplish Gene Ontology (GO) analysis to depict the unique biological significance based on DEGs between different groups [30]. The Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used to determine important pathways, and was performed in combination with two websites (https://david-d.ncifcrf.gov/home.jsp and http://kobas.cbi.pku.edu.cn/kobas3). The “p-value < 0.05” and the “|log2fold change|≥2” were used as the cutoff criteria for GO and KEGG enrichment analyses.

Gene Set Enrichment Analysis

For traditional analysis using DNA microarrays, the common approach involves focusing on a handful of genes at the top and bottom of L (i.e., those showing the largest difference) to discern telltale biological clues, but this approach has a few major limitations. To overcome these analytical challenges, we used a method called Gene Set Enrichment Analysis (GSEA) that evaluates microarray data (GSE45885) at the level of gene sets. Following the standard procedure for GSEA analysis [31], we first converted the expression dataset from GSE45885 into the tab-delimited GCT format as follows: the first column displays the gene symbol, while the second column is labeled “NA.” We then populated the subsequent columns with expression values from each sample. Subsequent operations were carried out in full accordance with GSEA’s standard protocol (http://www.gsea-msigdb.org/gsea/).

Protein-protein interaction network building

Differentially expressed mRNAs (|log2fold change|≥2, p-value < 0.05) were taken into the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING). The confidence scores were set at 0.4, 0.7, and 0.9 in diverse analysis. Then the Gene network files were input into Cytoscape software and the Molecular Complex Detection (MCODE) plugin of Cytoscape was used to analyze the core modules of the protein-protein interaction (PPI) network.

Co-expression network construction

Correlation networks are increasingly being used in bioinformatics applications. In this part, WGCNA was conducted to uncover the correlation among genes using the WGCNA package. Firstly, expression data of DEGs (p-value < 0.05) was input into R software to inspect good genes and samples; GSM1118241 (control) and GSM1118246 (NOA) were excluded from the analysis after quality assessment. To ensure that the network was scale-free, the power of the β was set to 24. The minimum number of modules was 10. Hierarchical clustering dendrogram summarized the gene modules of different colors. Heat map and topological overlap matrix (TOM) plots were used to visualize the module structure. The threshold for output to Cytoscape was set to 0.02.

Hub gene selection and validation

First, the gene network files exported from the WGCNA analysis were input into the Cytoscape software. The K-core values for each subcluster were then calculated using the MCODE plugin of Cytoscape. Referring to the detailed information of the GSE45887 dataset and previous studies [8, 11, 14], we confirmed that the sample of the GSE45887 dataset was different from that in the GSE45885 dataset, so the expression of hub genes was then validated using the GSE45887 dataset.

Integrative analysis of hub gene with single-cell transcriptome

After the above validation, we conducted an integrative analysis using the single-cell RNA-seq dataset (GSE106487). For this dataset, we used the Seurat package to analyze cell transcriptome data including two control samples and one NOA sample. First, cells with less than 2000 genes and less than 10,000 transcripts were filtered and the expression of other cells was normalized using log₂ (TPM/10+1). The data were then integrated, scaled, and centered using the ScaleData function. After PCA analysis, cells were clustered using the tSNE method with the resolution set to 1. Then, all clusters were identified as different cell types by specific genetic markers, respectively. After cell type identification, we chose the same cell type from testicular tissues of healthy and NOA adult males to further characterized the differential expression of the hub genes screened in the previous step using the violin plot.

Data availability statement

The authors confirm that the data supporting the findings of this study are available within the article and its supplementary materials.