Transcriptomic Analysis of Testicular Gene Expression in Normal and Cryptorchid Horses

Simple Summary Cryptorchidism is a common congenital malformation that results in impaired fertility in horses. The high abdominal temperature and the effects of this disease lead to differences in gene expression between retained testes and descended testes (DTs). Here, we focus on the genetic effects of cryptorchidism. All the differentially expressed genes (DEGs) between undescended testes (UDTs) and DTs were analyzed in this study. A total of 84 DEGs were associated with functions related to sperm development and male reproductive performance. Our study has provided fundamental transcriptomic data for future studies on equine testes and cryptorchidism. Abstract Testes produce sperm, and investigations into gene expression in the testes will enhance the understanding of the roles of testicular genes in male reproduction. Cryptorchidism, the failure of one or both testes to descend into the scrotal sac, is a common congenital malformation in horses. The major clinical consequence of this abnormality is impaired fertility. The aim of this study was to analyze the expression patterns of testicular genes and to identify the differentially expressed genes (DEGs) in testes between cryptorchid and normal horses. In this study, the gene expression patterns in equine testes and the DEGs between mature descended testes (DTs) and undescended testes (UDTs) were identified by RNA-seq and validated by real-time qPCR. Our results provide comprehensive transcriptomic data on equine testes. The transcriptomic analysis revealed 11 affected genes that were downregulated in UDTs, possibly as a result of the higher temperature in the abdomen than in the scrotal sac. These 11 genes have previously been associated with male reproduction, and their downregulation might explain the impaired fertility of cryptorchid horses. Two homozygous missense mutations detected in horses with cryptorchidism were absent in normal horses and were listed as potential pathogenic mutations; these mutations should be verified in the future.


Introduction
Cryptorchidism, or impaired testicular descent, is the failure of one or both testes to descend into the scrotal sac and is a common congenital malformation in horses. The major clinical consequence of this abnormality is impaired fertility [1] and a significantly increased risk of testicular malignancy due to the negative effects of elevated temperatures on the seminiferous tubules [2,3]. Undescended testes

Gene Expression Analysis
The equine genome EquCab2.0 (GCF_000002305.2) and the Y chromosome gene sequences from [24] were utilized as the reference genome and sequences for read mapping using TopHat, and the unique mapped reads were further analyzed. The read count was estimated with htseq-count. The expression levels of the assembled transcripts were normalized, and the DEGs were identified based on their size factors (SFs), relative log expression (RLE), and trimmed means of M-values (TMMs) with the DESeq2 R package (1.8.1) (https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/ doc/DESeq2.html). The fold changes (in log2 scale) and the corresponding p-values and padj values relative to the significance thresholds were estimated for the DEGs according to the normalized gene expression levels. Based on the expression levels, genes with a padj <0.05 and a |log2-fold change| >1 were considered differentially expressed in this study. The functional associations of the DEGs were analyzed, and the protein functions were predicted with the Database for Annotation, Visualization, and Integrated Discovery (DAVID) 6.8 (https://david.ncifcrf.gov/). All searches were conducted with a minimum false discovery rate (FDR) of <0.05 as the threshold.

Real-Time qPCR Verification of the mRNA Sequencing Results
Total RNA was reverse transcribed into cDNA by PrimeScript RT Enzyme in a reaction including gDNA Eraser (RR047, TaKaRa, Dalian, China). Each 10 µL of real-time qPCR included 5 µL of SYBR Green Real-time PCR Master Mix, 0.8 µL of cDNA, and 0.4 µL of each primer. The PCR conditions consisted of 1 cycle at 94 • C for 30 s followed by 40 cycles at 95 • C for 5 s, 61-67 • C for 30 s, and 65 • C for 5 s on a Bio-Rad CFX96 instrument; fluorescence was acquired at 95 • C in the single mode. The relative expression levels were determined using the 2 −∆∆Ct method with GAPDH as the control.

Overview of the Sequenced RNAs from Horse Testes
This study used RNA-seq to compare the gene expression of normal testes and retained testes. The RNA-seq generated a total of 73.15 million raw reads, with an average of 7.31 million reads, and a total of 70.80 million clean reads, with an average of 7.08 million (range: 5.98 to 8.88 million) clean Animals 2020, 10, 102 4 of 13 reads for each sample ( Table 1). The read alignment showed that 87.53% of the reads (61,997,273 reads) mapped to the equine genome (EcuCab2) on average. Of the mapped reads, 84.45% mapped to unique positions, and 1.41% mapped to multiple positions (Table 1). Only the uniquely mapped reads were considered in subsequent analyses. We conducted a principal component analysis (PCA) on the 10 samples and found that the two UDTs (CKY2b and GU4b), the five normal DTs of Guanzhong horses (GU1, GU2, GU3, GU4a, and GU5), and the three normal DTs of Chakouyi horses (CKY1, CKY2a, and CKY3) tended to form clusters. The testes were classified as normal DTs (PC1 < 10) or UDTs (PC1 > 10) based on PC1 (Figure 1). The current study included three comparisons: A comparison of the gene expression and DEGs between Guanzhong and Chakouyi horses (group 1), a comparison of the gene expression and DEGs between the DTs and UDTs of cryptorchids (group 2), and a comparison of the gene expression and DEGs between the DTs of normal horses and the UDTs of cryptorchid horses (group 3).  (Table 1). Only the uniquely mapped reads were considered in subsequent analyses. We conducted a principal component analysis (PCA) on the 10 samples and found that the two UDTs (CKY2b and GU4b), the five normal DTs of Guanzhong horses (GU1, GU2, GU3, GU4a, and GU5), and the three normal DTs of Chakouyi horses (CKY1, CKY2a, and CKY3) tended to form clusters. The testes were classified as normal DTs (PC1 < 10) or UDTs (PC1 > 10) based on PC1 ( Figure  1). The current study included three comparisons: A comparison of the gene expression and DEGs between Guanzhong and Chakouyi horses (group 1), a comparison of the gene expression and DEGs between the DTs and UDTs of cryptorchids (group 2), and a comparison of the gene expression and DEGs between the DTs of normal horses and the UDTs of cryptorchid horses (group 3). PCA of horse testicular tissues. CKYa represents CKY2a; CKY represents CKY1 and CKY3; GU represents GU1, GU2, GU3, and GU5; GUa represents GU4a; CKY2b represents CKY2b; and GUb represents GU4b.

Top Genes Expressed in Horse Testicular Tissue
The top 30 genes expressed in the testicular tissue in the three comparisons are shown in Figure  2. Strikingly, eight genes, PRM1, UBC, TNP1, ODF2, HSP90AA1, YBX3, GSTM3, and ACTG1, were included among the top 10 genes in all three comparisons, suggesting that these eight genes are highly expressed in horse testes. As shown in Figure 2B,C, compared to normal DTs, UDTs exhibited downregulation of some genes, including TNP1, DNAAF1, CALM3, CLGN, SPA17, and RPGRIP1; these results reveal the existence of different expression patterns between normal and retained testes. PCA of horse testicular tissues. CKYa represents CKY2a; CKY represents CKY1 and CKY3; GU represents GU1, GU2, GU3, and GU5; GUa represents GU4a; CKY2b represents CKY2b; and GUb represents GU4b.

Top Genes Expressed in Horse Testicular Tissue
The top 30 genes expressed in the testicular tissue in the three comparisons are shown in Figure 2. Strikingly, eight genes, PRM1, UBC, TNP1, ODF2, HSP90AA1, YBX3, GSTM3, and ACTG1, were included among the top 10 genes in all three comparisons, suggesting that these eight genes are highly expressed in horse testes. As shown in Figure 2B,C, compared to normal DTs, UDTs exhibited downregulation of some genes, including TNP1, DNAAF1, CALM3, CLGN, SPA17, and RPGRIP1; these results reveal the existence of different expression patterns between normal and retained testes.

Differential Gene Expression Analysis
The RNA-seq technique enabled analysis of the differential expression profiles via analysis of the transcript abundance with a high sensitivity. When a p < 0.05 and a |log2 fold change| ≥1 were used as cutoffs, a total of 400, 5959, and 5324 DEGs, respectively, were identified in the three comparisons defined above. A total of 191 DEGs were downregulated in the Guanzhong horses compared to the Chakouyi horses while 219 DEGs were upregulated ( Figure 3A). Sets of 3341 and 2618 transcripts were downregulated and upregulated in the CKY2b and GU4b testes, respectively, compared to the CKY2a and GU4a testes ( Figure 3B). A total of 3049 genes were downregulated in UDTs (CKY2b and GU4b) compared to normal horse testes (CKY1, CKY3, GU1, GU2, and GU3) while 2275 genes were upregulated in these tissues ( Figure 3C). Many genes showed differential expression between DTs and UDTs, as shown in Figure 3B,C, indicating that the expression of these genes might be associated with the reproductive performance of cryptorchids. The top 10 genes that were upregulated and downregulated in the three groups are shown in Tables S3-S5.

Differential Gene Expression Analysis
The RNA-seq technique enabled analysis of the differential expression profiles via analysis of the transcript abundance with a high sensitivity. When a p < 0.05 and a |log2 fold change| ≥1 were used as cutoffs, a total of 400, 5959, and 5324 DEGs, respectively, were identified in the three comparisons defined above. A total of 191 DEGs were downregulated in the Guanzhong horses compared to the Chakouyi horses while 219 DEGs were upregulated ( Figure 3A). Sets of 3341 and 2618 transcripts were downregulated and upregulated in the CKY2b and GU4b testes, respectively, compared to the CKY2a and GU4a testes ( Figure 3B). A total of 3049 genes were downregulated in UDTs (CKY2b and GU4b) compared to normal horse testes (CKY1, CKY3, GU1, GU2, and GU3) while 2275 genes were upregulated in these tissues ( Figure 3C). Many genes showed differential expression between DTs and UDTs, as shown in Figure 3B,C, indicating that the expression of these genes might be associated with the reproductive performance of cryptorchids. The top 10 genes that were upregulated and downregulated in the three groups are shown in Tables S3-S5.

Real-Time qPCR Validation of the RNA-Seq Results
To validate the gene expression results obtained by sequencing, the expression levels of 13 DEGs, namely, ATP1A4, ROPN1, NME8, CATSPER3, CATSPER1, AKAP4, MNS1, CABS1, LRGUK, TSGA10, PRM1, CAPZA3, and ENKUR, were analyzed through qPCR. Detailed information on the primers for the 13 genes is given in Table S6. The log2 fold change values determined by RNA-seq and the log2 values determined by qPCR through the 2 −ΔΔCt method with normalization to GAPDH for the comparisons of these 13 genes are shown in the electronic Supplementary Materials in Table  S7. The expression patterns of these 13 genes were consistent with the RNA-seq results according to a Pearson correlation coefficient of 0.9272 and p < 0.0001 ( Figure 4); quantification of the differential expression levels of the genes indicated significant findings and validated the repeatability and reproducibility of the gene expression data in this study.

Real-Time qPCR Validation of the RNA-Seq Results
To validate the gene expression results obtained by sequencing, the expression levels of 13 DEGs, namely, ATP1A4, ROPN1, NME8, CATSPER3, CATSPER1, AKAP4, MNS1, CABS1, LRGUK, TSGA10, PRM1, CAPZA3, and ENKUR, were analyzed through qPCR. Detailed information on the primers for the 13 genes is given in Table S6. The log2 fold change values determined by RNA-seq and the log2 values determined by qPCR through the 2 −∆∆Ct method with normalization to GAPDH for the comparisons of these 13 genes are shown in the electronic Supplementary Materials in Table S7. The expression patterns of these 13 genes were consistent with the RNA-seq results according to a Pearson correlation coefficient of 0.9272 and p < 0.0001 ( Figure 4); quantification of the differential expression

Functional Associations of the DEGs
Various genes cooperate with each other to exercise their biological functions. Accordingly, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis and Gene Ontology (GO) analysis for terms in the biological process, cellular component, and molecular function categories to further elucidate the biological functions of the DEGs. The DEGs were significantly enriched for 5, 5, and 8 GO terms and 5, 11, and 10 KEGG pathways that met the criterion of an FDR < 0.05 in the three comparison groups ( Figure 5). The most significantly enriched GO terms were the extracellular space term (GO: 0005615) with 34 associated genes, the extracellular exosome term (GO: 0070062) with 600 associated genes, and the extracellular space term (GO: 0005615) with 218 associated genes ( Figure  5). In particular, the DEGs enriched for the spermatogenesis, sperm motility, and acrosomal vesicle terms are associated with sperm structure and function. The KEGG pathway terms showing the highest levels of significance were the herpes simplex infection term (ecb05168) with 18 enriched DEGs, the pathways in cancer term (ecb05200) with 127 enriched DEGs, and the PI3K-Akt signaling pathway term (ecb04151) with 106 enriched DEGs ( Figure 5).

Functional Associations of the DEGs
Various genes cooperate with each other to exercise their biological functions. Accordingly, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis and Gene Ontology (GO) analysis for terms in the biological process, cellular component, and molecular function categories to further elucidate the biological functions of the DEGs. The DEGs were significantly enriched for 5, 5, and 8 GO terms and 5, 11, and 10 KEGG pathways that met the criterion of an FDR < 0.05 in the three comparison groups ( Figure 5). The most significantly enriched GO terms were the extracellular space term (GO: 0005615) with 34 associated genes, the extracellular exosome term (GO: 0070062) with 600 associated genes, and the extracellular space term (GO: 0005615) with 218 associated genes ( Figure 5). In particular, the DEGs enriched for the spermatogenesis, sperm motility, and acrosomal vesicle terms are associated with sperm structure and function. The KEGG pathway terms showing the highest levels of significance were the herpes simplex infection term (ecb05168) with 18 enriched DEGs, the pathways in cancer term (ecb05200) with 127 enriched DEGs, and the PI3K-Akt signaling pathway term (ecb04151) with 106 enriched DEGs ( Figure 5).

Effects of Cryptorchidism on Gene Expression
Our investigations focused on determining the effects of cryptorchidism on testicular gene expression levels in horses. An examination of the overlap of the three groups revealed that 4162 DEGs were shared by groups 2 and 3 when the DEGs between breeds were excluded ( Figure 6). These DEGs were assigned to five GO terms and eight KEGG pathways with the criterion of an FDR < 0.05 (Figure 7). It is worth noting that within the GO cellular component category, the acrosomal vesicle (GO: 0001669) and sperm principal piece (GO: 0097228) terms were enriched while within the biological process category, the spermatogenesis (GO: 0007283) and sperm motility (GO: 0030317) terms were enriched (Figure 7). There were 84 DEGs associated with these functions that might influence sperm development and male reproductive performance. Thus, our results support the hypothesis that high abdominal temperature or other factors associated with cryptorchidism impact sperm production and relevant gene expression levels.

Effects of Cryptorchidism on Gene Expression
Our investigations focused on determining the effects of cryptorchidism on testicular gene expression levels in horses. An examination of the overlap of the three groups revealed that 4162 DEGs were shared by groups 2 and 3 when the DEGs between breeds were excluded ( Figure 6). These DEGs were assigned to five GO terms and eight KEGG pathways with the criterion of an FDR < 0.05 (Figure 7). It is worth noting that within the GO cellular component category, the acrosomal vesicle (GO: 0001669) and sperm principal piece (GO: 0097228) terms were enriched while within the biological process category, the spermatogenesis (GO: 0007283) and sperm motility (GO: 0030317) terms were enriched (Figure 7). There were 84 DEGs associated with these functions that might influence sperm development and male reproductive performance. Thus, our results support the hypothesis that high abdominal temperature or other factors associated with cryptorchidism impact sperm production and relevant gene expression levels.

Single-Nucleotide Polymorphisms (SNPs) May Cause Cryptorchidism in Horses
Our investigations identified 1,032,196 SNPs in total, of which 146,610 were located in exonic regions. There were 48,083 missense mutations and 77,143 synonymous mutations. Herein, 46 SNPs, including 2 homozygous mutations and 44 heterozygous mutations, were detected in horses with unilateral cryptorchidism but not in normal horses. These variations occurred within 29 genes and are listed in Table 2.

Single-Nucleotide Polymorphisms (SNPs) May Cause Cryptorchidism in Horses
Our investigations identified 1,032,196 SNPs in total, of which 146,610 were located in exonic regions. There were 48,083 missense mutations and 77,143 synonymous mutations. Herein, 46 SNPs, including 2 homozygous mutations and 44 heterozygous mutations, were detected in horses with unilateral cryptorchidism but not in normal horses. These variations occurred within 29 genes and are listed in Table 2.

Discussion
In this study, we compared gene expression levels between UDTs and DTs to determine the genetic effects of cryptorchidism in horses. The statistical power of the analysis of the present research was not so high, because only two samples, belonging to two individuals of Guanzhong and Chakouyi horse breeds, showing unilateral cryptorchidism, were available for the study. This is due to the rarity of the cryptorchidic condition being found in horses. A total of 84 DEGs were associated with functions that might influence sperm development and male reproductive performance; among them, 11 genes were downregulated and were associated with more than two functional terms in our study. These genes have been previously reported to be associated with the fibrous sheaths (FSs) of spermatozoa, nuclear morphology maintenance during meiotic prophase, sperm construction, and male infertility (Table S7). In particular, three of them, CATSPERD, CATSPER1, and CATSPER3, are sperm ion channel proteins involved in sperm hyperactivated motility, which facilitates sperm penetration through the zona pellucida [25][26][27][28]. The expression analysis in the current study, showed that CatSper1-4 was expressed in horse testes, consistent with previous research [25,29,30]. Furthermore, we found lower expression levels of CatSper1-4 in UDTs than in normal mature DTs. Other studies carried out on ejaculated sperm have found that CatSper transcript levels are positively correlated with sperm motility [31,32]; this association could be the result of essential roles of sperm CatSper transcripts in sperm functions, such as motility, capacitation, and the acrosomal reaction [31,33]. Therefore, cryptorchidism leads to reduced expression levels of the CatSper1-4 genes in UDTs and impaired fertility in horses.
The gene expression levels of AKAP3 and AKAP4 were significantly lower in UDTs than in DTs in this study. AKAP3 and AKAP4 are the most abundant structural protein components of the sperm FS [34]. It is generally accepted that the FS provides mechanical support and flexibility to the sperm flagellum and provides a scaffold for signaling enzymes and glycolysis that is necessary for hyperactivated motility during capacitation [34][35][36][37]. No research on the relationship between AKAP gene expression and cryptorchidism has been reported. Based on this study, AKAP3 and AKAP4 are thought to be downregulated in retained testes as a result of the abdominal environment.
ROPN1 has also been localized to the FS and binds to AKAP3 and AKAP4 [34,38]. ROPN1L and ROPN1 compensate for each other in the absence of the opposing protein, possibly to maintain AKAP3 incorporation in the FS. We found that the ROPN1 and ROPN1L genes both showed significantly lower expression in UDTs than in DTs, indicating that retained testicles have higher frequencies of sperm with dysplasia or do not produce sperm at all. SPAG6 deficiency causes sperm motility defects and morphological abnormalities [39]. The SPAG6 gene plays an important role in maintaining flagellar structure and function in mammals [39]. Therefore, cryptorchidism might significantly downregulate the expression of the ROPN1, ROPN1L, and SPAG6 genes.
SNP calling was used to clarify the genetic pathogenesis of equine cryptorchidism. A large number of SNPs were screened; however, only 46 missense mutations were detected in horses with unilateral cryptorchidism but not in normal horses. Unfortunately, none of them were consistent with the causal mutations reported previously. Among the 46 missense SNPs, two were homozygous, located within the PRICKLE3 and PPP1R42 genes, although no relationship was shown between cryptorchidism and these two genes. The PRICKLE3 gene has been found to regulate ciliogenesis and cilia growth in gastrocoel roof plates [40]. PPP1R42 contains a protein phosphatase-1-binding site and translocates from the apical nucleus to the centrosome at the base of the flagellum during spermiogenesis [41]. Therefore, the homozygous SNPs within the PRICKLE3 and PPP1R42 genes are considered to be related to horse anorchism, although further verification is needed.

Conclusions
In conclusion, DEGs between normal DTs and UDTs were identified, and many genes were upor downregulated in UDTs as a result of cryptorchidism. The decreased expression of 11 genes is considered to be related to the impaired fertility of horses with cryptorchidism. In addition, two homozygous SNPs in the PRICKLE3 and PPP1R42 genes are proposed to cause cryptorchidism in horses. A limitation of the present study is the low sample size, since only two individuals showing unilateral cryptorchidism were available for analysis, but this condition is rare in horses. An increased sample size for analysis is aimed for in further research to corroborate the results of the present study.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-2615/10/1/102/s1. Additional supporting information may be found in the online version of this article. Table S1: Information on the horse testicular samples. Table S2: Top 10 unigenes that were upregulated or downregulated in the testes in Guanzhong horses compared with Chakouyi horses. Table S3: Top 10 unigenes with higher or lower expression in UDTs than in DTs. Table S4: Top 10 unigenes that were upregulated or downregulated in horses with UDTs compared with normal horses with DTs. Table S5: Information on the primers used in qPCR. Table S6: Log2 |fold change| values determined by RNA-seq and log2 values determined by qPCR through the 2 −∆∆Ct method for 13 verified genes. Table S7: Functions of 11 genes affected by increased temperature. Figure