Unearthing of key genes driving the infertility risk of cryptorchidism via weighted gene co- expression network

doi:10.21203/rs.3.rs-1840553/v1

Download PDF

Article

Unearthing of key genes driving the infertility risk of cryptorchidism via weighted gene co- expression network

https://doi.org/10.21203/rs.3.rs-1840553/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: Boys with cryptorchidism are susceptible to infertility. However, there are few known molecular markers of cryptorchidism-related infertility. Potential genes and molecular pathways to drive cryptorchidism-related infertility were investigated in this study via bioinformatics methods.

Methods: Gene Expression Omnibus (GEO) database was utilized to download raw GSE25518 (training set) and GSE16191 (validation set) datasets. GSE25518 was used for differentially expressed genes (DEGs) identification and weighted gene co-expression network analysis (WGCNA). Hub genes were identified from the intersection of key genes from DEGs and the key module in WGCNA. Hub genes’ expression differences and diagnostic value were subsequently verified via GSE16191. Gene Ontology (GO) analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, and Gene set enrichment analysis (GSEA) were deployed for functional annotation.

Results: 761 definitive DEGs were identified while 12 hub genes (CHRM4, CCR4, LMNB1, MTNR1B, TAC3, TRHR, RAD54L, CCNA1, SMC1B, GPRC6A, SYCE2, TAS2R19) were identified finally. The hub genes were also differentially expressed and showed excellent diagnostic value to distinguish cryptorchidism from normal testes in GSE16191. Additionally, they were well correlated with genes from the hallmark_spermatogenesis gene set. And they might be involved in oxidative phosphorylation, MYC targets V1, cholesterol homeostasis, KRAS activation, etc. Hallmark_spermatogenesis gene set was significantly enriched by five hub genes (RAD54L, CCNA1, SMC1B, TAC3, MTNR1B).

Conclusion: We have identified 12 genes that might become new biomarkers for cryptorchidism-related infertility while five genes (RAD54L, CCNA1, SMC1B, TAC3, MTNR1B) deserve more attention.

cryptorchidism

infertility

biomarker

WGCNA

bioinformatics

Cryptorchidism is one of the most common connatural malformations in men. It is featured by the failure to palpate the testis in the respective hemi-scrotum[1]. Up to 30% of preterm and 1-4% of full-term male infants were estimated to suffer from cryptorchidism[2]. The conventional treatment options for cryptorchidism encompass hormonal therapy and surgery[3]. Nevertheless, the therapeutics did not appear to renovate prior pathological lesions[3].

Boys with cryptorchidism are more prone to infertility[1, 4]. The incidence of azoospermia was 13% and 89% from unilateral to untreated bilateral cryptorchidism whether or not the statistic was adjusted[5]. A history of cryptorchidism and orchidopexy was spotted in 10% of infertile men[5]. The pathological mechanism of cryptorchidism-related infertility was convoluted. Anatomical, genetic, hormonal, and environmental factors might be involved[6, 7]. ‘Mini-puberty’ refers to the period in early infancy with a brief surge in the secretion of some crucial hormones. It is vital for the maturation of hypothalamic-pituitary-gonadal axis and transformation from gonocytes into spermatogonial stem cells[8-10]. However, in undescended testes staying with a higher temperature, mini-puberty might be significantly defective, affecting postpubertal spermatogenesis and semen quality[11]. Thus, defective mini-puberty was considered to be involved in cryptorchidism-related infertility[7, 9, 12]. A unique gene expression pattern was found in the neonatal and early postnatal mouse testes, as well as dynamic regulation of DNA methylation, suggetting the complexity in the formation and differentiation of spermatogonial stem cell[13].

Many studies tried to find molecular markers that drive cryptorchidism-related infertility. For instance, PLAP-positive germs cells were found to correlate with fertility potential[14, 15]. And PLAP, c-Kit, and OCT-4 were noticed with persistent and abnormal expressions in cryptorchidism patients, which might inhibit the transformation of gonocytes[16]. The severity of infertility in cryptorchidism was found to correlate with the single nucleotide polymorphism of NR5A1, suggesting a unique genetic phenotype of cryptorchidism-related infertility characterized by NR5A1 mutation[17]. The accelerated development of high-throughput technologies such as next-generation sequencing and microarrays facilitates the massive production of biological data[18]. Hadziselimovic’s team utilized the microarray technique to perform a whole-genome analysis of biopsies of cryptorchidism and contralateral descended testicular[19, 20]. They initially found that EGR4 deficiency in cryptorchidism impacted fertility potential by impeding the formation of a spermatogonial stem cells pool[19]. Furthermore, in another microarray with a larger sample[20], they divided the cryptorchidism into two groups (low azoospermia risk and high azoospermia risk) according to the count of Ad spermatogonia and presented a transcript signature that was abnormal in cryptorchidism patients with azoospermia risk. However, the initial grouping information could not be matched on Gene Expression Omnibus (GEO) database when downloading the microarray data. Thus, a later preliminary bioinformatics mining performed by Zhou et al. only explored the gene expression differences between normal testis and undescended testis utilizing the same microarray data[21].

Unlike seeking individual genes deemed significant in isolation between groups via traditional expression differences analysis, weighted gene co-expression network analysis (WGCNA) screened gene modules significantly associated with sample traits[22]. Hundreds of disease-related genes might be identified simultaneously and act as therapeutic targets based on the global approach of WGCNA[23, 24]. To our knowledge, no study utilized WGCNA to review existing microarray data in a networked way to unearth potentially key genes in cryptorchidism and cryptorchidism-related infertility. Thus, this study aimed to screen potentially key genes as well as molecular pathways promoting cryptorchidism-related infertility based on WGCNA and other bioinformatics methods. It might contribute to clarifying the mechanism of cryptorchidism-related infertility and ultimately provide a new perspective on the prevention of cryptorchidism-related infertility.

Data collection and preprocessing

GSE25518 (training set) and GSE16191 (validation set) on the GEO database were selected and downloaded. GSE25518 dataset was composed of 19 testicular biopsies of 18 cryptorchidism boys and 4 contralateral descended testicular biopsies from patients with testicular agenesis[20]. 16 testicular biopsies of 15 cryptorchidism boys and 4 biopsies from 4 contralateral testes in patients with unilateral testicular agenesis were comprised in GSE16191 dataset[19]. The expression matrix of probes was obtained by reading the CEL files via R “Affy” package[25]. The matrix was normalized utilizing a robust multi-array average (RMA) algorithm. The annotation information on the platform was adopted to obtain the gene symbol corresponding to each probe. Then, the matrix data were log2 transformed, and missing values were filled in with the R “impute” packages.

Differentially expressed genes (DEGs) identification

DEGs between undescended testes and contralateral descended testes were screened via the R “limma” package[26]. The raw P-value was corrected to gain the false discovery rate (FDR)[27]. DEGs were identified at the threshold of adjusted P<0.01 initially. Then a minimum 2.0-fold change with FDR<0.05 was utilized for identifying definitive DEGs[28].

Co-expression network analysis

The WGCNA was conducted based on the top 50% variant genes in GSE25518 via the R “WGCNA” package[29]. To achieve a relatively high scale-free network with an independence degree exceeding 0.8, an optimal soft threshold was obtained. Then, to evaluate the connectivity of genes in the network, the adjacency matrix was converted to a topological overlap measure (TOM) matrix. And one linkage hierarchical clustering dendrogram was constructed by modules consisting of genes with analogous expression profiles. 30 was set as the least gene number per module while the threshold to merge similar modules was 0.25. To measure module-trait relationships, the gene significance (GS) was calculated as well as module significance (MS) which was determined as the average GS of all genes in a module. Generally, the module with a higher MS score was considered to be more correlated with the clinical trait[29]. Modules with |R²|>0.9 as well as P<0.01 were considered to be significantly correlated with the cryptorchidism trait and selected for subsequent analysis.

Protein-protein interaction (PPI) network construction

The STRING database (http://t.stringdb.org) and GeneMANIA database (http://genemania.org/) were used for PPI network construction while the latter was also utilized for gene function prediction[30]. Cytoscape software 3.8.2 was utilized to determine the genes’ degree and better visualize the PPI network[31]. Cytoscape plugin “cytoHubba” was adopted to identify the top 30 hub genes in the PPI network.

Identification of hub genes

Module membership (MM)>0.8 and GS>0.8 were used to determine the hub genes in the target WGCNA module. Venn diagram analysis of hub genes in the target WGCNA module and the top 30 hub genes in the PPI network of definitive DEGs was conducted to identify candidate hub genes via the public website: https://bioinfogp.cnb.csic.es/tools/venny/index.html.

Functional enrichment analyses

Gene Ontology (GO) analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, and gene set enrichment analysis (GSEA) were conducted and visualized utilizing the R “clusterProfile” package[32]. Single gene GSEA was utilized to determine the pathways associated with hub genes based on hallmark gene sets of the Molecular Signatures Database (MSigDB) which was also achieved through the R “clusterProfile” package. A P-value<0.05 was regarded as statistically significant.

External validation of hub genes

Hub genes identified from GSE25518 between the cryptorchidism group and the normal group were validated utilizing GSE16191 for expression differences. Furthermore, receiver operating characteristic (ROC) curve analysis of hub genes was performed utilizing the R “pROC” package to determine their sensitivity and specificity[33]. The likelihood that hub genes might be a biomarker was represented by the area under the curve (AUC).

DEGs identification

6480 DEGs with an adjusted P-value<0.01 were found initially. 761 genes with a minimum 2.0-fold change with FDR<0.05 were identified as definitive DEGs finally. The distribution of the gene expression profile was presented in Fig. 1A. The definitive DEGs were also visualized by a clustered heatmap (Fig. 1B).

Functional enrichment analyses of definitive DEGs

GO and KEGG analyses were performed for the definitive DEGs, and the results were presented in Fig. 2. Biological Process (BP) was mainly enriched by GO analysis of definitive DEGs. The most significant five were cellular process involved in reproduction in multicellular organism, meiotic nuclear division, meiotic cell cycle, meiotic cell cycle process, and meiosis I (Fig. 2A). Five GO terms of Cellular Component (CC) were enriched including synaptonemal complex, synaptonemal structure, condensed nuclear chromosome, lateral element, and condensed chromosome (Fig. 2B) while no GO term of Molecular Function (MF) was enriched. The KEGG analysis showed that only Non-alcoholic fatty liver disease was enriched (Fig. 2C). All the detailed results of GO and KEGG analyses were listed in Supplementary Table 1.

Co-expression network construction

After data preprocessing, a weighted gene co-expression network was produced utilizing 10230 genes of GSE25518. All 23 samples were well clustered into two groups without outliers as presented in the hierarchical clustering dendrogram (Fig. 3A). Then a scale-free network was generated after determining the soft-threshold (power of β=10) (Fig. 3B). 15 co-expression modules were exported after hierarchical clustering (Fig. 3C). There were 3508 genes in the turquoise module, 1523 in the blue module, 1308 in the brown module, 1139 in the yellow module, 660 in the green module, 522 in the red module, 364 in the black module, 315 in the pink module, 286 in the magenta module, 166 in the purple module, 95 in the greenyellow module, 92 in the tan module, 66 in the salmon module, and 60 in the cyan module. The remaining 126 genes in the grey module were not clustered into any significant module actually. Then, module eigengene (ME) was calculated for each module. Dendrogram clusters and heatmaps of ME were utilized to assess the co-expression similarity for each module (Fig. 3D). Finally, the correlation between each module and the cryptorchidism trait was calculated (Fig. 3E). Noteworthily, the turquoise module held the strongest relevancy with the cryptorchidism trait (R²=-0.96, P=9e-13).

Screening of the key module

Firstly, we performed GSEA of the whole microarray based on the hallmark_spermatogenesis gene set. This set contains genes up-regulated during the production of male gametes (sperm). It was significantly enriched, suggesting that biological processes involved in spermatogenesis might be impaired in patients with cryptorchidism (Fig. 4A). Then the overlap of genes in each module with the hallmark_spermatogenesis gene set was explored. These overlapped genes were examined again for intersection with DEGs with P-value<0.01 and definitive DEGs with at least a 2.0-fold change (Fig. 4B). The blue and turquoise modules seemed to contain most of the spermatogenesis-related genes. Among 38 spermatogenesis-related genes in the blue module, 12 were DEGs with P-value<0.01 while 6 were definitive DEGs. Among 20 spermatogenesis-related genes in the turquoise module, 18 were DEGs with P-value<0.01, and 4 were definitive DEGs. Combined with the correlation of each module and the cryptorchidism trait, we chose the turquoise module as the key module for further analysis. Finally, GO and KEGG analyses on the turquoise module were conducted (Fig. 4C). The enriched GO terms involved in BP included regulation of ion transmembrane transport, antimicrobial humoral response, antimicrobial humoral immune response mediated by antimicrobial peptide, and antibacterial humoral response. Three GO terms involved in CC were enriched: collagen-containing extracellular matrix, apical plasma membrane, and 9+0 non-motile cilium. The enriched GO terms involved in MF included peptidase regulator activity, peptidase inhibitor activity, endopeptidase regulator activity, endopeptidase inhibitor activity, and serine−type endopeptidase inhibitor activity. No KEGG pathway was significantly enriched. The GO analysis results were all provided in Supplementary Table 2.

Candidate hub genes

1277 hub genes with MM>0.8 and GS>0.8 were obtained from the turquoise module for subsequent analysis (Fig. 5A). Then, we generated a PPI network of the definitive DEGs (Fig. 5B), from which the top 30 hub genes were identified by the Cytoscape plugin “cytoHubba”. The Venn diagram was drawn to capture the common hub genes obtained by the above two approaches (Fig. 5C). Finally, 12 candidate hub genes (CHRM4, CCR4, LMNB1, MTNR1B, TAC3, TRHR, RAD54L, CCNA1, SMC1B, GPRC6A, SYCE2, TAS2R19) were identified. We also mapped the PPI network and predicted the interaction genes of these hub genes utilizing the GeneMANIA database (Fig. 5D).

Validation of candidate hub genes

GSE16191 was used to validate the expression differences of candidate hub genes. They were all differentially expressed with P<0.001 (Fig. 6A). Only LMNB1 was upregulated in the cryptorchidism group which the remaining genes were all downregulated. Additionally, we evaluated the candidate hub genes’ diagnostic value in cryptorchidism group and normal group by ROC curve. It suggested that these hub genes showed excellent diagnostic value as a biomarker between the cryptorchidism group and normal group (AUC=1) (Fig. 6B).

Correlation of the hub genes with spermatogenesis-related genes

Firstly, correlation coefficient between each hub gene was calculated, and good correlations were found (Fig. 7A). Then PPI network between hub genes and genes in the hallmark_spermatogenesis gene set was mapped (Fig. 7B). All the hub genes and 106 genes from the hallmark_spermatogenesis gene set were involved in the network. Since 33 genes in the hallmark_spermatogenesis gene set were expressed differentially with P<0.01, we draw a correlation overview map of hub genes and these genes (Fig. 7C). Pairs of genes suggested good correlations. Additionally, as two hub genes (CHRM4, CCNA1) were contained in the hallmark_spermatogenesis gene set while three genes of the gene set were not contained in the GSE25518 dataset, we calculated the correlation coefficient between hub genes and the rest 130 genes of the hallmark_spermatogenesis gene set finally. We divided the absolute value of correlation coefficient into four categories (|r|>0.8, 0.8>|r|>0.5, 0.5>|r|>0.3, 0.3>|r|). The distribution of correlation of hub genes with the spermatogenesis-related genes was presented in Fig. 7D.

Function of hub genes in cryptorchidism-related infertility

The single gene GSEA based on hallmark gene sets was utilized to analyze the function of hub genes in cryptorchidism-related infertility. Hallmark_spermatogenesis gene set was significantly enriched in five genes (RAD54L, CCNA1, SMC1B, SYCE2, and TAC3) indicating spermatogenesis-related pathways were impaired when the genes were downregulated (Fig. 8A-E). As LMNB1 seemed to have an opposite regulatory function to the other hub genes which were consistent with the expression dissonance pattern of the hub genes, we only presented the top five enriched pathways correlated with the above six genes positively or negatively (Fig. 8A-F). And we counted and visualized the occurrence frequency of the top five enriched pathways correlated with the hub genes positively or negatively (Fig. 8G-H). Hallmark_oxidative_phosphorylation, hallmark_MYC_targets_V1, hallmark_cholesterol_homeostasis, and hallmark_KRAS_signaline_DN were significantly correlated with each hub gene while LMNB1 showed opposite correlations when compared with others (Fig. 8G-H). Hallmark_reactive_oxygen_species_pathway and hallmark_DNA_repair were positively correlated with LMNB1 while they were negatively correlated with other 10 hub genes (RAD54L, CCNA1, SMC1B, SYCE2, TAC3, TRHR, CHRM4, MTNR1B, TAS2R19, CCR4) and 6 hub genes (SYCE2, TRHR, CHRM4, GPRC6A, TAS2R19, CCR4) respectively. Additionally, hallmark_unfold_protein_response and hallmark_protein_secretion were negatively correlated with GPRC6A and another 5 genes (RAD54L, CCNA1, SMC1B, TAC3, MTNR1B) respectively while hallmark_mitotic_spindle and hallmark_G2M_checkpoint were positively correlated with CCNA1. All the single gene GSEA results of each hub gene were provided in Supplementary Table 3.

Literature review on the relationship between hub genes and cryptorchidism-related infertility

To seek the association of hub genes with cryptorchidism-related infertility, we retrieved the published literature for possible proof. Though no hub gene was found exactly involved in cryptorchidism-related infertility, ten genes (CCR4, LMNB1, LMNB1, MTNR1B, TAC3, TRHR, TRHR, RAD54L, CCNA1, SMC1B, GPRC6A, SYCE2) played a key role in infertility or spermatogenesis[34-50], indicating that these hub genes might be involved in cryptorchidism-related infertility. Additionally, there was no direct evidence of the function of CHRM4 and TAS2R19 in infertility. All the results of literature proofs were presented in Table 1.

In this study, we mined two common microarrays based on WGCNA and other bioinformatics methods. 12 genes were identified that might be involved in cryptorchidism-related infertility and become new biomarkers while five genes (RAD54L, CCNA1, SMC1B, TAC3, MTNR1B) deserve more attention.

Zhou et al. had already utilized the GSE25518 for unearthing key genes in cryptorchidism-related infertility[21]. The mTOR and arachidonic acid and pathways associated with several hub genes such as GPX8 were presented to be potentially involved in the process of spermatogenesis and the appearance of cryptorchidism based on KEGG and PPT analysis of DEGs. However, whether they were involved in the pathophysiological process of adult infertility induced by early cryptorchidism was still far from clear. The DEGs between testicular biopsies from cryptorchidism and contralateral descended testicular may be a result rather than a cause because the infertility-related damage might occur in the very early stage of cryptorchidism[7, 9, 12, 13]. Hub genes of the DEGs were not always equal to hub genes in cryptorchidism-related infertility process. We have considered the association between hub genes and current known spermatogenesis-related genes obtained from the hallmark_spermatogenesis gene set in this study which strengthened the credibility that the hub genes were involved in cryptorchidism-related fertility to some extent. Remarkably, though the target turquoise module we chose to analyze was most associated with the cryptorchidism trait and might be significantly associated with spermatogenesis-related genes, no spermatogenesis-related pathway was enriched by function enrichment analysis of different WGCNA modules. It might be due to that current known spermatogenesis-related genes were dispersed in different modules and did not reach the enrichment threshold in every single module.

When considering the 12 hub genes identified in our study, though the literature review indicated that most hub genes might be involved in infertility, direct evidence about the function of these genes in cryptorchidism-related infertility is virtually lacking. It also seemed impossible to prove whether these genes were results or causes of cryptorchidism-related infertility. Actually, unraveling the role of a gene on adult infertility from the onset of cryptorchidism was not an easy task. High-throughput microarray studies which directly seek gene expression differences between infertile patients with cryptorchidism history and normal people will make the search for key genes in cryptorchidism-related infertility more accurate. What’s more, the expression pattern of target genes in early cryptorchidism and subsequent changes as well as whether the early or postoperative intervention of the gene function could affect the risk of infertility in adulthood deserve further study in the future. Among the 12 hub genes, CHRM4 and CCNA1 are contained in the hallmark_spermatogenesis gene set. But there is no direct research on the role of CHRM4 in infertility, let alone cryptorchidism-related infertility. It may be since several genes in the hallmark_spermatogenesis gene set were computationally derived[51]. Another hub gene (LMNB1) without direct literature proof of the function in infertility is a requisite component of the nuclear lamina structure. It acts pivotally in chromosome condensation, heterochromatin distribution, as well as gene splicing and expression regulation[52]. It was the only upregulated gene in the hub genes and held a reverse regulatory relationship with other hub genes. Actually, no literature about the interaction between LMNB1 and other hub genes was found after retrieval. The PPI network of hub genes might provide a perspective to explore the interaction between them. Five hub genes (RAD54L, CCNA1, SMC1B, SYCE2, and TAC3) were significantly enriched by the hallmark_spermatogenesis gene set. Interestingly, they were in the top five holding the largest number of correlation coefficients exceeding 0.5 with genes in the hallmark_spermatogenesis gene set. In the PPI network of hub genes with spermatogenesis-related genes, their degrees appeared to be higher than many genes. However, RAD54L and TAC3 had the highest and the least degree respectively of all the hub genes. It indicated that despite the existence of predictions based on prior knowledge, the hub genes, especially the above five genes, might have more unknown interactions with known genes.

The function predictions of hub genes provided a perspective to explore their role in the process of cryptorchidism-related infertility. All the hub genes were significantly correlated with the hallmark_KRAS_SIGNALING_DN and hallmark_MYC_targets_V1. Metabolic programmes were found to be activated by oncogenic mutations of KRAS promoting the biological behavior of tumors[53]. And the MYC oncoprotein could disturb cellular activities widely such as apoptosis, cell cycle, and DNA damage response, thus driving many human cancers[54, 55]. It indicated that these hub genes might be also involved in tumorigenesis. Cryptorchidism was the most common risk factor for testicular cancer. Males with cryptorchidism were three to four times more likely to develop testicular cancer than those without cryptorchidism[56, 57]. Although we intended to find key genes participating in cryptorchidism-related infertility, the 12 hub gene just presented key differences between the cryptorchidism and normal testis before exploring their function in cryptorchidism-related infertility. Actually, several genes have been proved to be closely linked to cancer. For instance, human CCR4 deadenylase, the catalytic component of the Ccr4-Not complex, was reported to impact human gastric cancer cells via regulating proliferation and tumorigenicity[58]. And RAD54L was found to be a potential marker for the progression of bladder cancer progression[59]. Furthermore, several genes were found to be associated with testicular cancer[60, 61]. Carsten et al. found that CCNA1 was highly expressed in aggressive testicular germ cell tumors[60]. Shoko’s team reported that CCR4 was positive in a primary testicular diffuse large B-cell lymphoma case with skin involvement[61]. However, though the hub genes might cause the malignant transformation of cryptorchidism, it was not in the exploration scope of this research. Single gene GSEA suggested that infertility due to cryptorchidism and subsequent malignant transformation might not be completely separate outcomes. There might be some common pathologic basis. Several forms of male-factor infertility were reported to be correlated with elevated testicular cancer risk[62], confirming the above conjecture from the side.

Other several common pathways enriched from single gene GSEA such as hallmark_oxidative_phosphorylation, hallmark_cholesterol_homeostasis, and hallmark_oxidative_phosphorylation also provided an important research perspective to examine hub genes’ function in cryptorchidism-related infertility. Oxidative stress is the consequence of the dysfunction between antioxidants and reactive oxygen species of the organism. It might be a critical cause of male infertility[63]. Some researchers found oxidative stress by cryptorchidism was involved in infertility[64, 65]. Jung et al. found that phospholipid hydroperoxidase glutathione peroxidase (PHGPx) was fundamental to maintaining male fertility under cryptorchidism[64]. And Bae’s team reported that decursin might improve infertility in experiment-induced unilateral cryptorchidism rats by reducing oxidative stress[65]. Additionally, the regulator of cellular cholesterol homeostasis Abca1 was up-regulated in Sertoli's cells of mice undergoing cryptorchidism after a 1,25-dihydroxyvitamin D supplement, contributing to spermatogenesis[66]. Despite the lack of experimental validation, the consistency of single gene GSEA results with prior literature knowledge suggested the unexplored role of these genes in the network of oxidative stress as well as cholesterol homeostasis-related signaling pathways. Actually, several genes had already been demonstrated to regulate oxidative stress and cholesterol homeostasis. For example, CCNA1 protein level was found to correlate with alleviation of hepatocyte death and oxidative stress induced by Acetaminophen after Sirt6 knockdown[67]. Additionally, the CCR4-NOT multifunctional complex was found to be involved in the resistance to the oxidative stress of cells by regulating stress-responsive genes[68]. It was also requisite for steroid hormone synthesis and cholesterol homeostasis in Drosophila[69]. On the whole, the function exploration of the discovered hub genes is far from enough, and the single gene GSEA results especially the predicted common pathways provided an important perspective for further research in the future.

An interesting thing was that several hub genes were involved in the neurological or psychiatric function or disorders such as CHRM4[70, 71], MTNR1B[72], and TAS2R19[73, 74]. Just as reported, plentiful common proteins and biochemical characteristics were shared between the human brain and testis[75]. In fact, research on several hub genes might only stay at the level of establishing correlations with human diseases while their mechanisms and regulatory networks in human diseases are far from well understood. CHRM4, cholinergic receptor muscarinic 4, was proved to have extensive links with many psychiatric disorders in the past decades[76, 77], especially schizophrenia[70, 78]. In recent years, anemias, aging, and hemolysis of myelodysplastic syndrome were noticed to be corrected by the pharmacological inhibition of CHRM4[79], indicating that many unknown functions of these hub genes should be discovered. Actually, existing function research of most hub genes did not relate to the function predicted by our study. Anyway, the exact role of these genes in the clinical prognosis of cryptorchidism still needs to be explored in-depth in the future.

The primary strength of this study included mining publicly available datasets with high-throughput gene expression profiles via WGCNA. Additionally, to our knowledge, this is the first research utilizing WGCNA to explore key genes in cryptorchidism-related infertility. And not only found the DEGs between cryptorchidism and normal testis which was also performed by the previous study[21], but we also introduced the hallmark_spermatogenesis gene set and considered the association between the hub genes and spermatogenesis-related genes which increased the reliability that the genes we discovered were involved in cryptorchidism-related infertility. Despite the above important findings, we acknowledge that there are also several limitations. First, the sample sizes of the raw datasets were relatively small. Thus, research with larger sample sizes is necessary for validating these potential biomarkers in the future. Second, experiments are also essential for results verification because only bioinformatic methods were adopted in this study. However, an external dataset was utilized for results validation in our study which ensured the accuracy of our results to a certain extent. Third, we could not determine the causal relationship between the discovered biomarkers and cryptorchidism-related infertility. Functional studies are warranted to further clarify their exact role, identify candidate drug targets, and develop potential therapies for the prevention and treatment of cryptorchidism-related infertility.

We have identified 12 genes that might be involved in cryptorchidism-related infertility and become new biomarkers while five genes (RAD54L, CCNA1, SMC1B, TAC3, MTNR1B) deserve more attention. Further studies are required to clarify the exact function of these genes in cryptorchidism-related infertility, and whether they could become candidate drug targets to prevent and treat cryptorchidism-related infertility.

Acknowledgments

This work was supported by the National Innovation and Entrepreneurship Training Program for Undergraduates (Grant No. 202110730201).

Conflict of interests

The authors declare no competing interests associated with the manuscript.

Supplementary materials

Supplementary data to this study can be found in the submitted supplementary document.

Data availability statement

Data utilized in this work were all loaded from the public database. They are publicly available.

The STRING database (http://t.stringdb.org)

GeneMANIA database (http://genemania.org/)

https://bioinfogp.cnb.csic.es/tools/venny/index.html

Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/)

Braga, L.H., A.J. Lorenzo, and R.L.P. Romao, Canadian Urological Association-Pediatric Urologists of Canada (CUA-PUC) guideline for the diagnosis, management, and followup of cryptorchidism. Can Urol Assoc J, 2017. 11(7): p. E251-e260.
Goel, P., et al., Undescended testicle: An update on fertility in cryptorchid men. Indian J Med Res, 2015. 141(2): p. 163-71.
Braga, L.H. and A.J. Lorenzo, Cryptorchidism: A practical review for all community healthcare providers. Can Urol Assoc J, 2017. 11(1-2Suppl1): p. S26-s32.
Trump, T., et al., The Impact of Targeted Education of American Urological Association Cryptorchidism Guidelines in a Rural State: Improvement Still Needed. Res Rep Urol, 2021. 13: p. 437-443.
Chung, E. and G.B. Brock, Cryptorchidism and its impact on male fertility: a state of art review of current literature. Can Urol Assoc J, 2011. 5(3): p. 210-4.
Ferguson, L. and A.I. Agoulnik, Testicular cancer and cryptorchidism. Front Endocrinol (Lausanne), 2013. 4: p. 32.
Docampo, M.J. and F. Hadziselimovic, Molecular Pathology of Cryptorchidism-Induced Infertility. Sex Dev, 2015. 9(5): p. 269-78.
Forest, M.G., et al., Hypophyso-gonadal function in humans during the first year of life. 1. Evidence for testicular activity in early infancy. J Clin Invest, 1974. 53(3): p. 819-28.
Loebenstein, M., et al., Cryptorchidism, gonocyte development, and the risks of germ cell malignancy and infertility: A systematic review. J Pediatr Surg, 2020. 55(7): p. 1201-1210.
Hadziselimovic, F., et al., Piwi-pathway alteration induces LINE-1 transposon derepression and infertility development in cryptorchidism. Sex Dev, 2015. 9(2): p. 98-104.
Hutson, J.M., et al., The regulation of testicular descent and the effects of cryptorchidism. Endocr Rev, 2013. 34(5): p. 725-52.
Kurtoğlu, S. and O. Baştuğ, Mini puberty and its interpretation. Turk Pediatri Ars, 2014. 49(3): p. 186-91.
Kubo, N., et al., DNA methylation and gene expression dynamics during spermatogonial stem cell differentiation in the early postnatal mouse testis. BMC Genomics, 2015. 16(1): p. 624.
Thorup, J., et al., The relation between adult dark spermatogonia and other parameters of fertility potential in cryptorchid testes. J Urol, 2013. 190(4 Suppl): p. 1566-71.
Thorup, J., et al., Postnatal Germ Cell Development in the Cryptorchid Testis: The Key to Explain Why Early Surgery Decreases the Risk of Malignancy. Eur J Pediatr Surg, 2018. 28(6): p. 469-476.
Vigueras-Villaseñor, R.M., et al., Analysis of POU5F1, c-Kit, PLAP, AP2γ and SALL4 in gonocytes of patients with cryptorchidism. Acta Histochem, 2015. 117(8): p. 752-61.
Ferlin, A., et al., Mutational screening of NR5A1 gene encoding steroidogenic factor 1 in cryptorchidism and male factor infertility and functional analysis of seven undescribed mutations. Fertil Steril, 2015. 104(1): p. 163-9.e1.
Raza, K., Fuzzy logic based approaches for gene regulatory network inference. Artif Intell Med, 2019. 97: p. 189-203.
Hadziselimovic, F., et al., EGR4 is a master gene responsible for fertility in cryptorchidism. Sex Dev, 2009. 3(5): p. 253-63.
Hadziselimovic, F., et al., Testicular gene expression in cryptorchid boys at risk of azoospermia. Sex Dev, 2011. 5(2): p. 49-59.
Zhou, Y., et al., Bioinformatic identification of key genes and molecular pathways in the spermatogenic process of cryptorchidism. Genes Dis, 2019. 6(4): p. 431-440.
Zhao, W., et al., Weighted gene coexpression network analysis: state of the art. J Biopharm Stat, 2010. 20(2): p. 281-300.
Chen, Y., et al., Variations in DNA elucidate molecular networks that cause disease. Nature, 2008. 452(7186): p. 429-35.
Schadt, E.E., et al., An integrative genomics approach to infer causal associations between gene expression and disease. Nat Genet, 2005. 37(7): p. 710-7.
Gautier, L., et al., affy--analysis of Affymetrix GeneChip data at the probe level. Bioinformatics, 2004. 20(3): p. 307-15.
Ritchie, M.E., et al., limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res, 2015. 43(7): p. e47.
Reiner-Benaim, A., FDR control by the BH procedure for two-sided correlated tests with implications to gene expression data analysis. Biom J, 2007. 49(1): p. 107-26.
Rao, Z.T., S.Q. Wang, and J.Q. Wang, Exploring the osteoarthritis-related genes by gene expression analysis. Eur Rev Med Pharmacol Sci, 2014. 18(20): p. 3056-62.
Langfelder, P. and S. Horvath, WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics, 2008. 9: p. 559.
Szklarczyk, D., et al., The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res, 2017. 45(D1): p. D362-d368.
Su, G., et al., Biological network exploration with Cytoscape 3. Curr Protoc Bioinformatics, 2014. 47: p. 8.13.1-24.
Yu, G., et al., clusterProfiler: an R package for comparing biological themes among gene clusters. Omics, 2012. 16(5): p. 284-7.
Robin, X., et al., pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics, 2011. 12: p. 77.
Dai, X.X., et al., The CNOT4 Subunit of the CCR4-NOT Complex is Involved in mRNA Degradation, Efficient DNA Damage Repair, and XY Chromosome Crossover during Male Germ Cell Meiosis. Adv Sci (Weinh), 2021. 8(10): p. 2003636.
Berthet, C., et al., CCR4-associated factor CAF1 is an essential factor for spermatogenesis. Mol Cell Biol, 2004. 24(13): p. 5808-20.
Gao, Y., et al., Melatonin receptor depletion suppressed hCG-induced testosterone expression in mouse Leydig cells. Cell Mol Biol Lett, 2019. 24: p. 21.
Chen, J., et al., Expression of MTNR1a and MTNR1b in different tissues and localization in testes of yak (Bos grunniens). Acta Theriologica Sinica, 2020. 40(01): p. 64-70.
Precone, V., et al., A simultaneous next-generation sequencing approach to the diagnosis of couple infertility. Minerva Endocrinol (Torino), 2021.
Yang, J.J., et al., Uncovering novel reproductive defects in neurokinin B receptor null mice: closing the gap between mice and men. Endocrinology, 2012. 153(3): p. 1498-508.
Liu, X.P., et al., [Expression of thyrotropin-releasing hormone receptor mRNA in the rat testis development]. Shi Yan Sheng Wu Xue Bao, 2001. 34(2): p. 109-14.
Messiaen, S., et al., Rad54 is required for the normal development of male and female germ cells and contributes to the maintainance of their genome integrity after genotoxic stress. Cell Death Dis, 2013. 4(8): p. e774.
Russo, A., et al., Rad54/Rad54B deficiency is associated to increased chromosome breakage in mouse spermatocytes. Mutagenesis, 2018. 33(4): p. 323-332.
Liu, D., et al., Cyclin A1 is required for meiosis in the male mouse. Nat Genet, 1998. 20(4): p. 377-80.
Haraguchi, T., et al., Cyclin and protamine as prognostic molecular marker for testicular sperm extraction in patients with azoospermia. Fertil Steril, 2009. 91(4 Suppl): p. 1424-6.
Takabayashi, S., et al., A spontaneous smc1b mutation causes cohesin protein dysfunction and sterility in mice. Exp Biol Med (Maywood), 2009. 234(8): p. 994-1001.
Greither, T., et al., Fertility Relevance Probability Analysis Shortlists Genetic Markers for Male Fertility Impairment. Cytogenet Genome Res, 2020. 160(9): p. 506-522.
De Toni, L., et al., Polymorphism rs2274911 of GPRC6A as a Novel Risk Factor for Testis Failure. J Clin Endocrinol Metab, 2016. 101(3): p. 953-61.
Oury, F., et al., Osteocalcin regulates murine and human fertility through a pancreas-bone-testis axis. J Clin Invest, 2013. 123(6): p. 2421-33.
Bolcun-Filas, E., et al., Mutation of the mouse Syce1 gene disrupts synapsis and suggests a link between synaptonemal complex structural components and DNA repair. PLoS Genet, 2009. 5(2): p. e1000393.
Bolcun-Filas, E., et al., SYCE2 is required for synaptonemal complex assembly, double strand break repair, and homologous recombination. J Cell Biol, 2007. 176(6): p. 741-7.
Liberzon, A., et al., The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst, 2015. 1(6): p. 417-425.
Camps, J., M.R. Erdos, and T. Ried, The role of lamin B1 for the maintenance of nuclear structure and function. Nucleus, 2015. 6(1): p. 8-14.
Kerk, S.A., et al., Metabolic networks in mutant KRAS-driven tumours: tissue specificities and the microenvironment. Nat Rev Cancer, 2021. 21(8): p. 510-525.
Burchett, J.B., A.M. Knudsen-Clark, and B.J. Altman, MYC Ran Up the Clock: The Complex Interplay between MYC and the Molecular Circadian Clock in Cancer. Int J Mol Sci, 2021. 22(14).
Ahmadi, S.E., et al., MYC: a multipurpose oncogene with prognostic and therapeutic implications in blood malignancies. J Hematol Oncol, 2021. 14(1): p. 121.
Gurney, J.K., et al., Risk factors for cryptorchidism. Nat Rev Urol, 2017. 14(9): p. 534-548.
Cheng, L., et al., Testicular cancer. Nat Rev Dis Primers, 2018. 4(1): p. 29.
Song, X.H., et al., Human Ccr4 and Caf1 Deadenylases Regulate Proliferation and Tumorigenicity of Human Gastric Cancer Cells via Modulating Cell Cycle Progression. Cancers (Basel), 2021. 13(4).
Mun, J.Y., et al., E2F1 Promotes Progression of Bladder Cancer by Modulating RAD54L Involved in Homologous Recombination Repair. Int J Mol Sci, 2020. 21(23).
Müller-Tidow, C., et al., Cyclin A1 is highly expressed in aggressive testicular germ cell tumors. Cancer Lett, 2003. 190(1): p. 89-95.
Nakayama-Ichiyama, S., et al., Expression of multiple cytokines and CCR4 in primary testicular diffuse large B-cell lymphoma, not otherwise specified, involving the skin. Br J Haematol, 2012. 156(4): p. 420.
Negri, L., et al., Cancer risk in male factor-infertility. Placenta, 2008. 29 Suppl B: p. 178-83.
Makker, K., A. Agarwal, and R. Sharma, Oxidative stress & male infertility. Indian J Med Res, 2009. 129(4): p. 357-67.
Jung, K.Y., et al., Phospholipid hydroperoxide glutathione peroxidase is involved in the maintenance of male fertility under cryptorchidism in mice. Reprod Toxicol, 2015. 57: p. 73-80.
Bae, W.J., et al., Protective Effect of Decursin Extracted from Angelica gigas in Male Infertility via Nrf2/HO-1 Signaling Pathway. Oxid Med Cell Longev, 2016. 2016: p. 5901098.
Hirai, T., et al., Effect of 1,25-dihydroxyvitamin d on testicular morphology and gene expression in experimental cryptorchid mouse: testis specific cDNA microarray analysis and potential implication in male infertility. J Urol, 2009. 181(3): p. 1487-92.
Zhou, Y., et al., SIRT6 as a key event linking P53 and NRF2 counteracts APAP-induced hepatotoxicity through inhibiting oxidative stress and promoting hepatocyte proliferation. Acta Pharm Sin B, 2021. 11(1): p. 89-99.
Palermo, V., et al., The C-terminal region of yeast ubiquitin-protein ligase Not4 mediates its cellular localization and stress response. FEMS Microbiol Lett, 2021. 368(15).
Zeng, J., et al., The Drosophila CCR4-NOT complex is required for cholesterol homeostasis and steroid hormone synthesis. Dev Biol, 2018. 443(1): p. 10-18.
Qvist, P., et al., The Schizophrenia-Associated BRD1 Gene Regulates Behavior, Neurotransmission, and Expression of Schizophrenia Risk Enriched Gene Sets in Mice. Biol Psychiatry, 2017. 82(1): p. 62-76.
Walker, L.C., et al., Muscarinic M(4) and M(5) receptors in the ventral subiculum differentially modulate alcohol seeking versus consumption in male alcohol-preferring rats. Br J Pharmacol, 2021. 178(18): p. 3730-3746.
Thomson, D.M., et al., Mice lacking melatonin MT2 receptors exhibit attentional deficits, anxiety and enhanced social interaction. J Psychopharmacol, 2021: p. 2698811211032439.
Dmitrzak-Weglarz, M., et al., Transcriptomic profiling as biological markers of depression - A pilot study in unipolar and bipolar women. World J Biol Psychiatry, 2021: p. 1-13.
Diószegi, J., et al., Association of single nucleotide polymorphisms with taste and food preferences of the Hungarian general and Roma populations. Appetite, 2021. 164: p. 105270.
Matos, B., et al., Brain and testis: more alike than previously thought? Open Biol, 2021. 11(6): p. 200322.
Dean, B. and E. Scarr, Possible involvement of muscarinic receptors in psychiatric disorders: a focus on schizophrenia and mood disorders. Curr Mol Med, 2015. 15(3): p. 253-64.
Levran, O., et al., African-specific variability in the acetylcholine muscarinic receptor M4: association with cocaine and heroin addiction. Pharmacogenomics, 2016. 17(9): p. 995-1003.
Rietschel, M., et al., Association between genetic variation in a region on chromosome 11 and schizophrenia in large samples from Europe. Mol Psychiatry, 2012. 17(9): p. 906-17.
Trivedi, G., et al., Muscarinic acetylcholine receptor regulates self-renewal of early erythroid progenitors. Sci Transl Med, 2019. 11(511).

Table 1 Literature review on the relationship of hub genes with cryptorchidism-related infertility.

Hub gene	Encoding protein	Proof of potential function in cryptorchidism-related infertility	Research details of proofs
CHRM4	Cholinergic Receptor Muscarinic 4	No	None
CCR4	C-C Motif Chemokine Receptor 4	Yes[34, 35]	The CNOT4 subunit of the CCR4-NOT complex is crucial in male germ cell meiosis mediating multiple important processes such as XY chromosome crossover[34]. CAF1 (part of the CCR4-NOT complex) was essential for spermatogenesis[35].
LMNB1	Lamin B1	No	None
MTNR1B	Melatonin Receptor 1B	Yes[36, 37]	Testosterone synthesis was regulated by MTNR1B through mediating multiple cell pathways[36]; MTNR1b was widely expressed in testis at different developmental stages of yak, and the expression level gradually increased with age[37].
TAC3	Tachykinin Precursor 3	Yes[38, 39]	TAC3 was identified as a possible pathogenic variant by next-generation sequencing to screen 99 couples with idiopathic infertility[38]; TACR3^-/- mice had various reproductive deficiencies that appear to be due to a central, as opposed to a gonadal, etiology[39].
TRHR	Thyrotropin Releasing Hormone Receptor	Yes[40]	Testis of rats could specifically express TRHR, and its transcription was development-dependent[40].
RAD54L	RAD54 Like	Yes[41, 42]	RAD54 had a critical role in maintaining genomic integrity of developing germ cells[41, 42].
CCNA1	Cyclin A1	Yes[43, 44]	The spermatogenesis before the first meiotic division was blocked in mice with CCNA1^-/-, resulting in infertility[43]; The transcript ratios of CCNA1 messenger RNA were significantly reduced in patients with spermatogenic disorders[44].
SMC1B	Structural Maintenance Of Chromosomes 1B	Yes[45, 46]	A spontaneous SMC1B mutation (exon 5 is missing 16 nucleotides) leads to cohesin protein malfunction and sterility in mice[45]; SMC1B was marked as a prime candidate for male fertility impairment in a predictive analysis[46].
GPRC6A	G Protein-Coupled Receptor Class C Group 6 Member A	Yes[47, 48]	Cryptorchidism during fetal life and/or low sperm production in adulthood might be attributed to exposure reduction to androgens mediated by GPRC6A inactivation or sub-function[47]. GPCR6A might become a novel susceptibility locus for human primary testicular failure[48].
SYCE2	Synaptonemal Complex Central Element Protein 2	Yes[49, 50]	Homologous synapsis from sites of recombination in normal meiosis might be promoted by the biochemical interaction between SYCE2 and RAD51[49]. SYCE2 Knockout led to the failure in synapsis and germ cell loss consequently, contributing to infertility of both males and females[50].
TAS2R19	Taste 2 Receptor Member 19	No	None

No competing interests reported.

SupplementaryMaterial.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

Unearthing of key genes driving the infertility risk of cryptorchidism via weighted gene co- expression network

Status:

Version 1

Abstract

Figures

Introduction

Materials And Methods

Results

Discussion

Conclusion

Declarations

References

Tables

Additional Declarations

Supplementary Files

Status:

Version 1