Identification of key genes related to unilateral cryptorchidism in sheep by comprehensive transcriptomics and proteomics analyses

Abstract

Background

Cryptorchidism is one of the most common reproductive abnormalities in rams, which seriously harms the reproductive capacity of rams and causes significant economic losses to the sheep industry. However, there are few studies elucidating its hereditary predisposition in sheep.

Results

In the present study, the transcriptome and proteome of the cryptic (CT) and contralateral (CLT) testis from four unilaterally cryptorchid rams, and the normal testis (NT) from four healthy rams were analyzed using RNA-seq and TMT-based proteomics, respectively. A total of 10,357, 10,175, and 132 differentially expressed genes (DEGs) were identified between CT vs. CLT, CT vs. NT, and CLT vs. NT. Venn diagram showed that 9744 DEGs (5499 up-regulated and 4245 down-regulated) shared in CT vs. CLT and CT vs. NT. Functional enrichment analysis revealed that 5499 up-regulated DEGs were mainly involved in regulation of programmed cell death and metabolic process, while 4245 down-regulated DEGs were closely related to reproductive process, such as spermatogenesis, sexual reproduction, reproduction and male gamete generation. Furthermore, 325 overlapped genes (114 up-regulated and 211 down-regulated) between DEGs and DAPs that shared the same regulatory status were identified by combining transcriptomics and proteomics. Ten genes, including AKAP4, AKAP3, FSIP2, HSPA1L, HSPA4L, TUBB, TXNRD2, CDC42, PGK1 and HSPA1A, were identified as candidate key genes related to unilateral cryptorchidism.

Conclusion

Our results revealed that both gene and protein expression in the cryptic testis of unilateral cryptorchid rams is massively altered. Bioinformatics analysis unveiled several candidate genes and signaling pathways potentially involved in unilateral cryptorchidism. These findings provide new insights into the molecular mechanism underlying spermatogenesis failure caused by cryptorchidism.

Background

Cryptorchidism is one of the most common reproductive congenital abnormalities in males, in which one or both testicles have not descended into the scrotum and remain in the abdominal cavity or the inguinal canal area at birth [1,2,3]. This disease is characterized by testicular malformation and defective spermatogenesis, which may lead to male subfertility and even infertility in humans and animals [4, 5]. Histological analyses of unilateral cryptorchid testes in pig showed that the number of germ cells was markedly decreased, and the seminiferous tubules contained only a few spermatogonia but did not contain post-meiotic germ cells [6]. Furthermore, credible evidence suggests that bilateral cryptorchidism will cause serious damage of male reproduction than unilateral cryptorchidism [7]. A previous study reported that the incidence of oligozoospermia in patients with untreated bilateral cryptorchidism was 80%, while its incidence decreased to 50% in patients with unilateral cryptorchidism [8].

Testicular descent proceeds in two gonad migration phases, transabdominal and inguinoscrotal [5, 9, 10]. The initial transabdominal phase is dependent on the development of the fetal gubernaculum. Both INSL3 and its receptor RXFP2 have been demonstrated to play an important role in regulating the growth and development of the gubernaculum, and thus polymorphisms in these genes are closely associated with cryptorchidism in various species, including mice [11, 12], human [13, 14], dog [15] and sheep [2]. The second inguinoscrotal phase is governed by testosterone and its receptor AR, and syndromes related to the dysregulation of reproductive hormone secretion often cause cryptorchidism, such as idiopathic hypogonadotropic hypogonadism (IHH) [16]. Cryptorchidism is a complex disease caused by multiple genetic and environmental risk factors. Previous studies on cryptorchidism have primarily focused on humans and animal models, and only very few genetic risk factors have been identified [6, 17,18,19,20,21]. Recent progress in high-throughput multi-omics technologies have greatly improved our ability to identify potential risk factors and to explore the pathogenic mechanisms of complex diseases. For instance, a transcriptomic study found that massive alterations in gene expression of the inguinal testicles from unilateral cryptorchid dogs, and two risk variants identified in the region of KAT6A [17]. Recently, a large-scale genome-wide association study (GWAS) has identified a variant in HMGA2 associated with cryptorchidism in dog [22]. Additionally, a splice donor variant in DMD has been discovered to be associated with cryptorchidism in human GWAS [21].

Congenital cryptorchidism is a major form of disorder of sex development (DSD) in newborn boys, with an incidence ranging from 1.6 to 9% [5]. However, previous epidemiological investigations from the United Kingdom and the United States have found that the prevalence of cryptorchidism in sheep is only 0.5% and 0.6%, respectively, and most of them (~ 86%) are unilateral cryptorchidism [23]. Because it’s anatomical structures, hormone secretion and timing of testis descent closely resemble humans situation, sheep can serve as an appropriate large animal model for investigating the molecular pathogenesis of cryptorchidism [2]. Male fertility is a key factor in ensuring successful conception, and plays an important role in breeding systems. Compared to normal individuals, animals with bilateral cryptorchidism have been found to exhibit impaired testicular, and thereby resulted in reduced testicular size and the number of germ cells, ultimately leading to male sterility [6, 18, 24]. Despite its relatively low prevalence in sheep, cryptorchidism has serious consequences for reproduction performance, resulting in huge economic losses in the global sheep industry [25]. It is therefore essential to improve the understanding of the molecular mechanisms associated with cryptorchidism.

Hu sheep, an excellent local breed in China, is famous for its precocious puberty, perennial estrus, and high productivity, with the biggest market share in the Chinese sheep industry [26]. To maximise productivity performance, it is important to identify the cryptorchid rams as early as possible and to culled these affected individuals from breeding stock. In the present study, we comparatively analyzed and compared the gene and protein expression profile between undescended and descended testes tissues from unilateral cryptorchid rams through transcriptome and proteome data. The aim of this study was to elucidate how cryptorchidism impairs testicular development and spermatogenesis, and further to mine the putative key genes associated with unilateral cryptorchidism in Hu sheep, which will provide deeper insight into the molecular mechanisms underlying the pathogenesis of cryptorchidism in sheep.

Methods

Sample collection

The testis samples used in this study were collected as part of a larger research project, which involved the slaughter of approximately 3000 Hu sheep across nine different batches from 2018 to 2022. The animal experiments, including animal husbandry, slaughter, testicular traits and sample collection have been clearly described in previous study [26]. Briefly, each lamb was housed in a single pen under the same nutritional and management conditions. At the age of 180 days, all rams were humanely slaughtered to produce premium lamb carcasses according to halal slaughter procedure. The scrotum was incised for each ram after slaughtered, the ram with one testicle in the scrotum and the other in the abdominal cavity was diagnosed with unilateral cryptorchidism. Contralateral testis (scrotal testis, CLT, n = 4) and cryptic testis (retained testis, CT, n = 4) from rams with unilateral cryptorchidism, as well as normal testis (NT, n = 4) from healthy rams were harvested immediately after slaughter. The testis from each individual was rapidly dissected and then the testicular tissue was divided into two parts, one part was snap-frozen in liquid nitrogen and stored at -80 °C for total RNA and protein extraction, and another part was fixed in Bouin’s solution for 24 h, followed by paraffin embedding for hematoxylin and eosin (H&E) staining.

Testicular histology

Hematoxylin and eosin (H&E) staining was performed the protocol outlined in a previous study [27]. Briefly, the testis samples were fixed in 10% formalin and dehydrated by a series of alcohols, clarified in xylene, and embedded in paraffin. Then, samples were sliced serially into 5 μm sections and stained with hematoxylin and eosin by routine methods. Stained sections were examined with a light microscope (Carl Zeiss, Gottingen, Germany), and the section images were captured and analyzed with Motic Images Advanced 3.2 software (Motic, XiaMen, China).

RNA extraction and sequencing

Total RNA was extracted from each sample using Trizol reagent (TransGen Biotech, Beijing, China) according to the manufacturer’s protocol. RNA concentration and purity were measured using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). RNA degradation and contamination was monitored on 1% agarose gels, and RNA integrity was measured using a 2100 Bioanalyzer (Agilent) or LabChip GX (Perkin Elmer). 1 µg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, USA) following manufacturer’s recommendations. The concentration and insert size of the constructed library were detected by Qubit2.0 and Agilent 2100, respectively, and the effective concentration was accurately quantified by quantitative polymerase chain reaction (Q-PCR). All libraries were sequenced on an Illumina NovaSeq 6000 platform and 150 bp paired-end reads were generated. Raw sequencing data in FASTQ format were processed to filtered out adaptors and low-quality reads using TrimGalore (v0.6.10, Cutadapt, v4.6) with default parameters. The clean reads were then mapped to the sheep reference genome (Oar_v4.0, GCA_000298735.2) using HISAT2 software (v2.2.1) with default parameters. Read counts was generated by featureCounts tool (v2.0.6), and read counts were normalized to fragments per kilobase of transcript per million (FPKM). Differential gene expression analysis was performed using the DESeq2 R package (version 1.42.0). Differentially expressed genes (DEGs) were determined using a threshold of |log2FC| > 1 and FDR < 0.05. The correlation coefficient was calculated using the Pearson method by the “cor” function, and plotted by using the corrplot package in R. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the DEGs were performed by the clusterProfiler R package (version 4.10.0). GO term and KEGG pathway with an FDR < 0.05 were considered as significantly enriched.

Gene expression validation by quantitative real-time PCR (qRT-PCR)

Total RNA was reverse transcribed using RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific). The reaction conditions for mRNA reverse transcription were as follows: 25 °C for 5 min, 42 °C for 45 min and 85 °C for 5 min. Quantitative real-time PCR (qRT-PCR) was performed using the SYBRgreen PCR Master Mix (ABI, USA, 4304437) on a Bio-rad CFX384 Real-Time System (Bio–Rad, Hercules, CA, USA). All primers were designed using Primer Premier 6.0 software (Premier Biosoft, San Francisco, USA) and the primer sequences are listed in the Supplementary Table 1. The 10 µL PCR reaction mixtures contained 5 µL of 2 × TB Green Premix Ex Taq II (Tli RNaseH Plus) (Takara Bio, Dalian, China), 1 µL of cDNA (50 ng/µL), and 0.4 µL each primer (0.4 µM), and 3.2 µL ddH₂O. The reaction conditions were: 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s and 60 °C for 30 s. Fluorescence measurements were performed every 0.5 °C from 65 to 95 °C to obtain the melting curve. Relative gene expression was calculated using the 2^−∆∆Ct method with GAPDH as a reference gene [28]. Every sample was run in triplicate. Differences between groups were analyzed by paired t tests using IBM^® SPSS^® Statistics (version 25), and differences were considered as statistically significant when P < 0.05.

Protein extraction and trypsin digestion

About 100 mg testis tissue was finely ground into powder using a mortar and pestle filled with liquid nitrogen. The obtained powder was lysed in an appropriate amount of SDT lysis buffer and disrupted using a homogenizer. The samples were further ultrasonicated for 2 min before centrifuging at 16,000 g for 20 min at 4 °C. The supernatant was collected and quantified with a BCA Protein Assay Kit (Bio-Rad, USA). Protein digestion (300 µg for each sample) was performed according to the FASP procedure described by Wisniewski et al. [29].

TMT labeling and high pH reverse phase fractionation (HPRP)

100 µg of peptides from each sample were labeled with TMT10plex (Thermo Fisher Scientific) according to the manufacturer’s instructions. The labeled peptides were mixed in equal amounts, and the dried peptides were fractionated using Pierce™ High pH Reversed-Phase Peptide Fractionation Kit (Thermo Fisher Scientific). Finally, fractionated peptides were collected and combined into 10 components. The peptides of each component were dried and reconstituted with 0.1% FA for liquid chromatography tandem mass spectrometry (LC-MS/MS) analysis.

LC-MS/MS analysis and protein identification

Chromatographic separation was carried out using an Easy-nLC system (Thermo Scientific). The sample was injected into a trap column and then subjected to gradient separation through a chromatography column at a flow rate of 300 nL/min. Buffer A was 0.1% formic acid, and buffer B was 0.1% formic acid in 95% acetonitrile. The gradient was as follows: 0–2 min, 2–8% buffer B; 2–71 min, 8–28% buffer B; 71–79 min, 28–40% buffer B; 79–81 min, 40–100% buffer B; and 81–90 min, 100% buffer B. The peptide fragments were separated and analyzed using a Q-Exactive HF-X mass spectrometer (Thermo Scientific, USA) for data-dependent acquisition (DDA) mass spectrometry for 90 min. Raw data files were analyzed using the Sequest HT search algorithm in Proteome Discoverer (version 1.6.0.16, Thermo Scientific). Uniprot sheep database (48903 sequences, accessed on 06 July 2023) was used for database search. The following search parameters were used: main search mass tolerance 10 ppm and a MS/MS tolerance of 0.02 Da, trypsin as proteolytic enzyme and allowing two missed cleavages. Oxidation (M) and acetylation (Protein N-term) were specified as variable modifications, and carbamidomethylation (C) as a fixed modification. The FDR was set to 0.01 for both peptide and protein identifications.

Bioinformatics analysis

Relative protein expression levels were normalized to the median average peptide ratio. Differentially abundant proteins (DAPs) were identified as p value < 0.05 and fold changes (FC) > 1.5 or FC < 1/1.5 by using Student’s t test. Functional annotation was performed by blasting against the public databases using BLASTP (v.2.2.26, E-value ≤ 1e-5), including UniProtKB/Swiss-Prot database, Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Ontology (GO). GO and KEGG enrichment analyses were carried out with the Fisher’s exact test, and FDR correction for multiple testing was also performed. GO terms and KEGG pathways with a P < 0.05 were considered statistically significant. Protein-protein interaction (PPI) network was constructed using STRING database version 11.5 (https://cn.string-db.org/) [30], which was then visualized using Cytoscape 3.9.0 software [31]. The top five proteins with the highest degree of connectivity were selected and considered as the hub genes by using the “Network Analyzer” plugin in Cytoscape.

Correlation analysis of the transcriptomes and proteomes

To evaluate the concordance between the transcriptome and proteome in the CT vs. CLT comparison group, both transcriptome and proteome data were median normalized using log2 transformation. The Pearson’s correlation coefficient was used to calculate the correlation between transcript and protein levels, and R > 0.80 was considered a significant positive correlation [32]. Subsequently, the DEGs with the same expression trend in transcriptome and proteomics were selected for the PPI network construction.

Results

Histological assays in ram testes

Histological assays were performed to evaluate the morphological changes of the CLT and CT groups. As shown in Fig. 1a, the CLT group exhibited normal histology. The seminiferous tubules were tightly and orderly distributed, and were lined with Sertoil cells and all phases of germ cells, including spermatogonia, spermatocytes, round spermatids, and elongated spermatids (Fig. 1b). However, the seminiferous tubules of the CT group comprised only Sertoli cells and several spermatogonia but no post-meiotic germ cells, such as round spermatids and elongated spermatids (Fig. 1c).

Fig. 1

HE staining of the (a) contralateral testis (CLT), (b) normal testis (NT) and (c) cryptic testis (CT) of Hu sheep at 6 months of age. Sg, spermatogonia; Sp, spermatocyte; RS, round spermatid; ES, elongated spermatid; Sc, Sertoli cell; Lc, Leydig cell. Bar = 20 μm

Overview of transcriptomic analysis

After a stringent quality control, a total of 519.70 million clean reads from 12 samples were obtained, with an average of 43.31 million clean reads per sample. On average, the percentage with a sequencing quality score > Q30 was 93.55, and the GC content were 50.74%. The average mapping rate of clean reads to the ovine reference genome was 89.10% (range 88.11–92.58%) (Supplementary Table 2). Theses high-quality reads ensure the reliability of subsequent analyses and results. Pearson correlation coefficient among four biological replicates were all above 0.9 in each group (Supplementary Fig. 1), indicating good reproducibility among biological replicates. A gene was defined as expressed if its standardized FPKM > 0 in any one of the samples. In total, 19,651 genes were considered expressed and used for further analysis. Hierarchical cluster analysis revealed that the genes expression patterns of these 19,651 genes were highly similar between CLT and NT groups, while there was a clear difference in CT group compared to CLT and NT groups (Fig. 2a). In total, 10,357 (4422 up-regulated and 5935 down-regulated), 10,175 (5647 up-regulated and 4528 down-regulated), and 132 (58 up-regulated and 74 down-regulated) DEGs were identified between CT vs. CLT, CT vs. NT, and CLT vs. NT, respectively (Supplementary Table 3). Venn diagram showed that 9744 DEGs shared in CT vs. CLT and CT vs. NT groups, and only 40 DEGs were shared among the three comparison groups (Fig. 2b).

Fig. 2

Expression analysis and quantitative comparison of the recognized DEGs in CT vs. CLT, CT vs. NT, and CLT vs. NT. (a) Hierarchical clustering heatmap based on standardized expression values across the three comparison groups. (b) Venn diagram indicated the shared and uniquely expressed genes in three comparison groups

Functional enrichment analyses of the DEGs

GO and KEGG pathway enrichment analyses were carried out to explore the biological functions of the identified DEGs. Notably, more than 94% of the DEGs in CT vs. CLT overlapped with the DEGs in CT vs. NT. Thus, the 9744 common DEGs (5499 up-regulated and 4245 down-regulated) in both comparison groups were subjected to functional enrichment analyses, respectively. A total of 22 significantly enriched GO terms, including 10 BP, 10 CC and 2 MF, were found in up-regulated DEGs (P-adjust < 0.05, Supplementary Table 4). In the BP category, enriched GO terms were mainly associated with organonitrogen compound biosynthetic process, translation, positive regulation of programmed cell death, and positive regulation of apoptotic process. Meanwhile, 80 significantly enriched GO terms, including 53 BP, 23 CC and 14 MF, were identified when restricting the analysis to the down-regulated DEGs (P-adjust < 0.05, Supplementary Table 4). Among them, the most significant BP terms were mainly involved in reproductive process (e.g., gamete generation, spermatogenesis, and sexual reproduction), cellular component organization or biogenesis (e.g., organelle organization, nuclear division, and organelle fission). The top 10 enriched GO terms for each category of the up-regulated and down-regulated DEGs were visualized as a bubble chart, respectively (Fig. 3a, b). KEGG pathway analysis identified 138 and 6 significantly enriched pathways in the up-regulated and down-regulated DEGs, respectively (Supplementary Table 4). Notably, the up-regulated DEGs were enriched in pathways involved in cancer biology, signal transduction, and cellular functions, whereas the down-regulated DEGs were enriched in pathways involved in cell motility, cell growth and death, and carbohydrate metabolism (Fig. 3c, d).

Fig. 3

Functional enrichment analysis of DEGs that shared in the CT vs. CLT and CT vs. NT. (a) GO enrichment analysis of up-regulated DEGs. (b) GO enrichment analysis of down-regulated DEGs. (c) KEGG pathway enrichment of up-regulated DEGs. (d) KEGG pathway enrichment of down-regulated DEGs

Overview of proteomic analysis

For the proteomic analysis, a total of 62,022 unique peptides were acquired, and 7137 proteins were confidently identified and quantified in the CT and NT samples. Quality assessing of the proteomic data showed that the median length of peptides was 10 amino acids, ranging from 6 to 49 (Supplemental Fig. 2a). Proteins with a single peptide, 2–20 peptides, and > 20 peptides account for 22.57%, 71.29% and 6.14%, respectively (Supplemental Fig. 2b). The molecular weight ranged from 10 to 200 kDa (Supplemental Fig. 2c). The protein sequence coverage of 0–10%, 10–30%, 30–50% and 50–100% accounted for 43.76%, 39.32%, 13.79%, and 3.14%, respectively (Supplemental Fig. 2d). PCA analysis of the proteomic data revealed clear separation between CT and NT groups (Fig. 4a). Hierarchical clustering heatmap found that the expression profiles of identified proteins were clearly distinct between CT and NT groups (Fig. 4b). Based on the cutoff threshold for screening DAPs, 1148 DAPs were detected, of which 758 were down-regulated and 390 were up-regulated in protein abundance in the CT vs. NT (Supplemental Table 5).

Fig. 4

Comparison of protein abundance profiles in CT vs. NT. (a) Principal component analysis (PCA) plot. (b) Hierarchical cluster heatmap of the identified proteins

Functional enrichment analyses of the DAPs

The biological functions of the identified DAPs in CT vs. NT were subsequently investigated using GO and KEGG pathway enrichment analyses. GO analysis showed that the up-regulated DAPs were significantly enriched in 5 CC and 17 MF terms, while no significantly enriched BP terms were identified. The top three CC terms were chromatin, nucleosome and protein-DNA complex, and the top three MF terms were chromatin DNA binding, nucleosomal DNA binding and structural constituent of chromatin (Fig. 5a and Supplemental Table 6). Furthermore, there were 65 BP, 14 CC and 17 MF were significantly enriched in down-regulated DAPs. The most highly enriched BP, CC and MF terms were the response to microtubule-based movement, chromosome, and microtubule motor activity, respectively (Fig. 5b and Supplemental Table 6). For KEGG pathway analysis, 14 and 12 significantly enriched pathways were identified in up- and down-regulated DAPs, respectively. The up-regulated DAPs were mainly involved in metabolic pathways, carbon metabolism, carbohydrate metabolism pathways (e.g., glycolysis / gluconeogenesis, fructose and mannose metabolism, pentose phosphate pathway, and pyruvate metabolism), and amino acid metabolism pathways (cysteine and methionine metabolism, tyrosine metabolism, and glycine, serine and threonine metabolism) (Fig. 5c and Supplemental Table 6). While the down-regulated DAPs were mainly enriched in metabolic pathways, cell cycle, motor proteins and protein processing in endoplasmic reticulum pathways (Fig. 5d and Supplemental Table 6).

Fig. 5

Functional enrichment analysis of differentially abundance proteins (DAPs) in CT vs. NT. (a) GO enrichment analysis of up-regulated DAPs. (b) GO enrichment analysis of down-regulated DAPs. (c) KEGG pathway enrichment of up-regulated DAPs. (d) KEGG pathway enrichment of down-regulated DAPs

Correlation between transcriptome and proteome

To assess the consistency of gene expression at transcription and protein levels, as well as to further investigate the regulatory mechanism of DEGs on testicular development and spermatogenesis, an integrative analysis of RNA-seq and TMT-based proteomic data was conducted. As shown in Figs. 6a and 2814 of 7137 identified proteins had corresponding genes in transcriptome data, and 356 of 1148 identified DAPs had corresponding DEGs. Of which, the expression tendencies of 325 DAPs (114 up- and 211 down-regulated proteins) agreed with the transcriptomic data. Correlation analysis based on Pearson correlation showed that the fold change of the overlapped DEGs and DAPs had a high correlation (R = 0.80, Fig. 6b, Supplemental Table 7).

Fig. 6

Correlation analysis of the transcriptomes and proteomes. (a) Venn diagram of the numbers of all detected genes, DEGs, all detected proteins and DAPs in CT vs. NT. (b) Comparison of the expression between transcriptomic (y-axis) and proteomic (x-axis) profiling. Only significant changes are color coded, red and green nodes indicate the opposite trend while blue and grey nodes indicate the opposite trend of the two omics

Protein-protein interaction network analysis

To further identify the putative key genes associated with cryptorchidism, we focused on the 325 overlapped genes between DEGs and DAPs that shared the same regulatory status, and a PPI network for these genes was constructed. The resulting network consisted of 214 nodes (73 up- and 141 down-regulated genes/proteins) and 496 edges, and average number of neighbors was 4.68 (Fig. 7). The key genes were investigated separately for the up- and down-regulated gene sets using the degree algorithm in the Network Analyzer plugin in Cytoscape. Among the 141 down-regulated proteins, the top five proteins with the highest degrees of connectivity were AKAP4, AKAP3, FSIP2, HSPA1L and HSPA4L. Among the 73 up-regulated proteins, the top five proteins with the highest degrees of connectivity were TUBB, TXNRD2, CDC42, PGK1 and HSPA1A.

Fig. 7

Protein-protein interaction (PPI) network construction and analysis. The PPI network contained a total of 214 nodes (73 up- and 141 down-regulated genes/proteins) and 496 edges. The line shows the relationship between the two proteins. The size and color of the nodes is proportional to the degree of connectivity within the network

qRT-PCR validation of DEGs

To verify the reliability of the transcriptomic data, ten representative DEGs potentially related to the androgen receptor signaling pathway were selected for qRT-PCR validation, including AMH, AR, CLDN11, FSHR, GATA4, INHA, INSL3, LHCGR, SOX9 and WT1. Notably, transcriptomic profile analysis showed that these genes were significantly up-regulated in the CT group compared with the CLT and NT groups (Fig. 8a). Consistently, the RT-qPCR results also showed that the expression levels of these genes were all significantly increased in the CT group than in the NC group (Fig. 8b), indicating the accuracy of the transcriptome data.

Fig. 8

Validation of RNA-seq by qRT-PCR. (a) Differences in mRNA expression level of ten DEGs in CT vs. NT were assayed by qRT-PCR assay. (b) Comparison of Log2 fold change in ten DEGs between RNA-seq and qRT-PCR. *P < 0.05, **P < 0.01

Discussion

Cryptorchidism is one of the most common congenital abnormalities among males, and is characterized by unilateral or bilateral fetal testes descend failure. This disorder is one of the known risk factors for subfertility, infertility and testicular cancer [3]. The descent of the testes from a position near the kidneys into the constant low-temperature environment of the scrotum is an important developmental step on the path toward normal testicular development and successful spermatogenesis. It is well documented that the cryptorchid testes have a lower testicular size and the volume of the seminiferous tubules than the normal testes. In the current study, histological examination of the contralateral testis from the unilateral cryptorchid ram showed apparently normal testicular histology, the seminiferous tubules were tightly and orderly distributed, and were lined with all phases of germ cells and Sertoli cells. While, only Sertoli cells and several spermatogonia were detected, and post-meiotic germ cells was not identified in the seminiferous tubules of the cryptic testis. Similar morphological alterations have been reported in other animal models with both spontaneous and surgically induced unilateral cryptorchidism, including rats [33], boars [6, 34], horses [35] and dogs [36].

Testicular descent is a unique physiological adaptation found in placental mammals allowing optimal spermatogenesis at scrotal temperature below core body temperature [37]. Considering that the undescended testes long-term face an unusual environment, such as elevated temperature, that does not affect the descended testes. It is conceivable that gene expression profiles may vary greatly between undescended and descended testes. Consistent with this conjecture, our transcriptome results revealed that the expression patterns of most genes were significantly different in the CT group compared to the CLT and NT groups, while there was a highly similar of expression patters between the CLT and NT groups. This finding aligns with previous studies in dogs and horses, which found that the transcriptome of inguinal testicles was massively altered in comparison to contralateral descended gonads [17, 18]. Moreover, it is noteworthy that there was little change in the transcriptome of the contralateral testis (descended testicles) of the unilateral cryptorchid rams and of the scrotal testicles of the healthy rams. This observation therefore strongly suggested that the massive changes in gene expression in the cryptic testis is not governed by germinal DNA variants in the regulatory sequences, but it is related to their face an unusual environment, such as elevated temperature and intra-abdominal pressure [17].

Based on our filtering criteria (|log2FC|> 1 and adjusted P-value < 0.05), a total of 10,357, 10,175, and 132 DEGs were identified between CT vs. CLT, CT vs. NT, and CLT vs. NT, respectively. Notably, 9744 DEGs shared in CT vs. CLT and CT vs. NT groups, of which 5499 were up-regulated and 4245 were down-regulated. Among these DEGs, several have been reported to be aberrantly expressed in the cryptic testis of dogs, including INSL3, RXFP2, CYP17A1, CYP19, and AMH [15, 24, 38]. Whereas, some candidate genes associated with syndromic cryptorchidism appear to be unaltered, such as HMGA2 [22], DMD [21], and KAT6A [17]. Furthermore, many reproductive hormone-related genes were differentially expressed, such as FSHR, LHCGR, LHB, AMH, ESR1, NR3C1, INHA, INHB, INHBB, CYP17A1, and HSD17B, and all of these DEGs were up-regulated in the CT group compared to the CLT and NT groups, suggesting that hormonal imbalance may be an important factor for spermatogenesis failure in cryptorchidism. However, the cause of spermatogenesis failure in the cryptorchid testes remains controversial. Huff and his colleagues have proposed that the abnormalities observed in the cryptorchid testes may not be solely attributed to long-term thermal injury, based on biopsies analysis of both the cryptorchid and contralateral descended testes [39]. Additionally, a long-term study demonstrated that impaired transformation of spermatogonia into type A dark (Ad) spermatogonia during prepubertal period in cryptorchid boys may not be cause by elevated temperature but by severe hormonal imbalance. This result suggests that temperature stress in the cryptorchid testes may not be a major causative factor for impaired testicular development and infertility [40, 41].

Subsequently, GO and KEGG enrichment analyses were performed on 9744 DEGs shared in CT vs. CLT and CT vs. NT groups to explore the underlying biological functions and involved pathways. For up-regulated DEGs, the significantly enriched GO BP terms were mainly involved in regulation of programmed cell death and metabolic process. For down-regulated DEGs, the GO BP terms related to reproductive process were most significantly enriched, such as spermatogenesis, sexual reproduction, reproduction and male gamete generation. As mentioned above, the cryptorchid testis has been expose to complex intra-abdominal or inguinal environments lies along the high-temperature and high-pressure for a long time, which leads to spermatogenesis dysfunction due to the germ cells are extremely sensitive to heat stress [4, 42]. According to several histological studies in model animals, it worth believing that testicular development and spermatogenesis were severely impaired in the cryptorchid testis in both surgically induced cryptorchidism and congenital cryptorchidism [6, 20, 24]. Consistent with the histological findings, the up-regulated DEGs were primarily associated with the regulation of programmed cell death and metabolic process, whereas the down-regulated DEGs were mainly enriched in genes related to general reproductive processes.

The androgen receptor signaling pathway plays an important role in multiple biological processes during spermatogenesis, including the proliferation and differentiation of both germ cells and somatic cells, the formation of inter-Sertoli cell tight junctions, and the establishment of the blood-testis barrier, etc [43]. In current study, both RNA-seq and qRT-PCR results showed that all ten genes involved in androgen receptor signaling, including AMH, AR, CLDN11, FSHR, GATA4, INHA, INSL3, LHCGR, SOX9, and WT1, were significantly upregulated in the cryptorchid testes. This finding suggested that the androgen receptor signaling pathway may be activated in cryptorchidism. Claudin-11 (CLDN11) is a member of the claudin tight junction family, which is expressed during all stages of spermatogenesis and plays important roles in the formation of BTB. Previous study has shown that overexpression and relocalization of CLDN11 is closely related to impaired spermatogenesis [36].

In current study, a total of 325 overlapped genes between DEGs and DAPs that shared the same regulatory status were identified by combining transcriptomics and proteomics. These proteins were then used to construct a PPI network and to identify putative key genes associated with cryptorchidism. Of which, several gene families have been reported to be involved in spermatogenesis when we inspected this list of genes, such as DNAH, SPATA, HSP70, and AKAP gene families. In our study, four members of DNAH gene family have been identified to be differently expressed (DNAH1, DNAH2, DNAH8, DNAH12). DNAH gene family encodes for 13 axonemal dynein heavy chain proteins, which plays important roles in the structure and function of cilia and flagella that are mainly involved in cell motility and ATP hydrolysis. Genetic mutations of the DNAH family members have been linked to male infertility [44]. Spermatogenesis-associated (SPATA) gene family consists of several evolutionarily conserved genes that are specifically expressed in testis and play critical roles in regulating spermatogenesis and fertility. Several researches have found that inactivation of SPATA6 leads to acephalic spermatozoa and male sterility in humans and mice [45, 46]. Two members of the 70kD heat shock protein (HSP70) family, HSPA1L and HSPA4L have been reported to be highly expressed in post-meiosis germ cells, and appears to play a key role in spermatogenesis and sperm maturation [47, 48]. In humans, semen quality is influenced by gene expression of HSPA4L, and reduced HSPA4L expression in spermatozoa is positively related to poor sperm quality and could increase the risk for asthenozoospermia [47]. AKAP3 and AKAP4 belong to A kinase-anchoring proteins (AKAPs), which constitute one of the wider and functionally conserved family of the signal-organizing scaffolding proteins in the mammalian genome, allowing the formation of subcellular super-molecular complexes and assuring the coordination of cAMP-responsive events within a specific subcellular localization [49]. AKAP3 and AKAP4 isoforms are uniquely expressed in spermatids and spermatozoa and are involved in sperm motility and flagellar assembly. Abnormal levels of AKAP3 and AKAP4 expression affecting the flagellar structures and resulting in asthenozoospermia [49]. Our results indicate that both gene and protein expression were massively altered in the cryptic testicles of unilateral cryptorchid rams. Several candidate key genes potentially associated with unilateral cryptorchidism were identified, however, their specific functions remain unclear. Further research is needed to understand the hereditary predispositions and pathophysiological basis of cryptorchidism.

Conclusion

This study provides the first comparative analysis of transcriptome and proteome across the CT, CLT and NT groups. The results revealed that the genes expression patterns were massively altered in CT group compared to CLT and NT groups. Additionally, ten hub genes (AKAP4, AKAP3, FSIP2, HSPA1L, HSPA4L, TUBB, TXNRD2, CDC42, PGK1 and HSPA1A) and the androgen receptor signaling pathway may play an important role in unilateral cryptorchidism. These findings provide new insights into the molecular mechanism underlying spermatogenesis failure caused by cryptorchidism, and screened out several candidate key genes associated with unilateral cryptorchidism of Hu sheep.