Comparison of reproductive performance and functional analysis of spermatogenesis factors between domestic yak and semi-wild blood yak

Abstract

This study investigates differences in reproductive performance, testicular histology, and transcriptomic profiles between male Subei (SB; semi-wild) yaks and two domestic yaks, Gannan (GN) and Qinghai (QH). Key metrics including mating age, utilization time, breeding capacity, morphometric traits, and testicular indices were analyzed. SB yaks exhibited superior reproductive metrics, including earlier sexual maturity, prolonged utilization periods, and enhanced breeding capacity compared to GN and QH (P < 0.05). Morphologically, SB yaks demonstrated significantly greater body weight, and testicular dimensions. Compared with GN and QH yaks, the seminiferous tubules of SB yaks exhibited significantly larger spermatogenic cells and luminal cavities, along with a notably higher sperm density within the luminal cavity. Transcriptomic analysis identified 2,403 and 4,428 differentially expressed genes (DEGs) in GN vs. SB and QH vs. SB comparisons, respectively. Eight key genes (TPPP3, SMAD3, PAFAH1B3, BMP7, ARSA, CTNNB1, SMAD4, STAT3) and three pathways (Hippo, pluripotency regulation, TGF-β) were implicated in testicular development and spermatogenesis. These findings underscore the genetic and physiological advantages of SB yaks, offering insights for enhancing male yak reproductive performance.

The reproductive performance of male mammals is primarily reflected in testicular development, semen quality, and reproductive disorders. As a crucial indicator of fertilization capacity, semen quality significantly affects the conception rate in female animals [1]. The process of mammalian spermatogenesis is accomplished through precise coordinated regulation by numerous genes and hormones, among which genetic factors play a decisive role [2, 3]. The expression of these genes exhibits stage-specific characteristics with precise temporal and spatial regulation [4,5,6]. Numerous studies have demonstrated that spermatogenesis is a dynamic process involving the proliferation and differentiation of spermatogonial stem cells (SSCs) within the testes. The proper progression of this process directly determines semen quality and exerts a significant impact on the reproductive performance of male animals [7]. With the rapid development and widespread application of modern molecular biology technologies in the field of life sciences, investigating the regulatory mechanisms of spermatogenesis within yak testes at molecular levels, such as transcription and translation, has become one of the key research directions [8, 9]. Guowen Wang et al. revealed the key molecular mechanisms underlying the transition from gonocytes to spermatogonia and subsequent spermatogenesis by analyzing dynamic gene expression profiles in testes of domestic yaks across different age groups. They validated the functional roles of conserved genes such as CXCR4 in yaks, supporting their universal significance in mammalian sperm development [10]. Scholars have revealed key genes (Dmc1, Spo11, Tnp2, Crem) affecting male infertility in cattle-yak hybrids through comparative transcriptomic analyses of genomic differences among cattle-yak hybrids, yaks, and cattle, providing critical insights into yak reproductive physiology [11,12,13,14]. Through cross-species transcriptomics studies involving mice [15, 16], humans [17, 18], and yak [19], numerous genes influencing spermatogenesis have been identified. These discoveries not only deepen our understanding of fundamental biological mechanisms but also demonstrate significant application potential in medicine, agriculture, and ecological conservation.

Subei yaks, inhabiting the Yanchiwan National Nature Reserve, represent a unique genetic resource shaped by natural and artificial selection. While their superior growth and production traits are recognized, reproductive mechanisms remain underexplored. This study integrates histological and transcriptomic approaches to compare SB yaks with domestic GN and QH breeds, aiming to elucidate the molecular drivers of spermatogenesis in SB yaks, which holds significant research implications and value for yak breed improvement and the development of the industry.

Research methodology

Experimental animals and sample collection

The sexual maturity time, time of utilization, breeding ability, newborn weight, body weight, body height, oblique length of body, chest measurement, and body length of 4-year-old GN, QH, and SB yaks were compared by group investigation and access to records. The left testes of 4-year-old yaks were collected from the slaughterhouses of Qingyuan in Linxia, Huangyuan in Qinghai, Subei County in Jiuquan (10 yaks/copy), and the circumference of the testis, the long testis diameter, and the short testis diameter were measured and weighed. The tunica albuginea of the testes was incised with a scalpel, followed by epididymal removal. A transverse circumferential incision was made along the midsagittal plane of the testes to obtain the testicular parenchyma adjacent to the epididymal head. The parenchyma was dissected into 2 × 2 cm³ cubes and fixed in 4% paraformaldehyde for histological sectioning. Simultaneously, remaining tissue samples were flash-frozen in liquid nitrogen and transferred to a − 80 °C freezer for subsequent transcriptomic analysis and qRT-PCR validation.

H.E staining

Put the tissue sections on the dyeing rack in turn them into xylene I and xylene II in each 10 min, and in turn into different concentration gradient alcohol solutions for descending dewaxing hydration, distilled water washing 5 min, hematoxylin staining nucleus 5 ~ 6 min, distilled water rinsing 3 times (3 s/time), hydrochloric acid alcohol differentiation 2 ~ 3 s, into the tap water back to blue 10 min (running water), distilled water washing 3 min, eosin staining cytoplasm 3 min, distilled water rinsing 3 times (3 s/time), different concentration gradient alcohol solution ascending dehydration, xylene I, II in each 10 min, neutral resin glue sealing, airing at room temperature for 4 days.

RNA preparation and quality inspection

Samples of 100 mg each from the testicular tissues of GN, QH, and SB yaks were collected and finely ground in liquid nitrogen. The ground tissue was then transferred to a 1.5 mL centrifuge tube designed for non-RNA enzymes. Total RNA was extracted from the yak testicular tissue using the Trizol method. After extraction, the RNA was allowed to dry naturally at room temperature. An appropriate volume of DEPC-treated sterile water was added to dissolve the RNA completely, and the resulting solution was stored at − 80 °C in a freezer.

The purity and concentration of total RNA were detected by Nanodrop2000 ultramicro ultraviolet spectrophotometer. The samples with OD260/OD280 ≥ 1.8, OD260/230 ≥ 2.0 were qualified, and the total RNA was greater than 1 μg. The qualified total RNA can be used for subsequent experiments.

Library preparation and sequencing

RNA library construction was performed using the TruSeqTM RNA sample preparation Kit (Illumina, SanDiego, CA). Sequencing with the Illumina NovaSeq 6000 sequencer.

Raw data quality control and sequence alignment analysis

The raw data is quality controlled by fastp [20], and the low-quality data is filtered out to obtain clean reads. Then, the clean reads were compared with the yak reference genome sequence by TopHat2 and HISAT2 [21] application software to obtain mapped reads for expression calculation and subsequent transcript assembly.

Gene expression analysis

Stringtie [22] was used to reconstruct transcripts, and RSEM [23] was used to calculate the expression of all genes in each sample. Fragments Per Kilobase of transcript per Million mapped reads (FPKM) were used as a quantitative indicator.

$$\text{FPKM}={10}^{9}\text{C}/\text{NL}$$

Note: Let FPKM(A) be the expression of gene A, then C is the number of sequenced fragments compared to gene A, N is the total number of sequenced fragments compared to the reference gene, and L is the number of bases of gene A. The FPKM method eliminates the influence of the differences in gene lengths and sequencing quantities on the calculation of gene expression, and the calculated gene expression can be directly used to compare the differences in gene expression among different.

Screening of differentially expressed genes

After obtaining the number of gene read counts, the differential expression analysis of genes between different samples was performed by DESeq2 [24] application software, DEGs were identified, and then the function of DEGs was studied. In this study, |log₂ (Fold Change)|≥ 2 and FDR < 0.05 were used as the screening conditions, and the genes that met this condition could be regarded as significant DEGs.

Gene Ontoiogy (GO) enrichment analysis of DEGs

In this study, Goatools application software was used to perform GO functional enrichment analysis on DEGs. After confirming the main GO functions of DEGs, Fisher was used for accurate test. In order to control the false positive of the calculation, the P-value was corrected by Benjamini and Hochberg multiple test method. According to the GO database, according to its different functions, it is divided into three aspects: molecular function (MF), cellular component (CC) and biological process (BP). When the P-value < 0.05, it is considered that this GO function is significantly enriched.

Target gene screening

Gene cards were used to screen out the genes related to testicular development and sperm production stores greater than 5 [25] and intersected with the DEGs related to reproduction screened by GO enrichment. PCA, PLSDA, and VIP analysis of variable projection importance were performed by MetaboAnalyst (https://www.metaboanalyst.ca/). VIP > 1 and P < 0.05 were used as criteria to further screen out the main genes related to spermatogenesis in GN, QH, and SB yaks.

Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of DEGs

KEGG enrichment analysis was performed on the target gene using the KOBAS online tool. Fisher was used for an accurate test, and the P-value was corrected by BH multiple tests. When the P-value < 0.05, it is considered that there is a significant enrichment of the KEGG function.

Regulation network analysis

High-confidence protein–protein interaction (PPI) networks were constructed via the STRING database (https://cn.string-db.org/) with a confidence threshold ≥ 0.7. Network topological parameters were analyzed using Cytoscape 3.7, and genes ranking in the top 10% for both Degree and Betweenness centrality values were selected. The resulting PPI network was visualized to identify key genes highly associated with spermatogenesis regulation.

Verification of DEGs

In order to further verify the reliability of RNA-Seq data, 11 DEGs were randomly selected for quantitative verification analysis by qRT-PCR. Quantitative primers were designed by Premier5.0 application software, and the yak β-actin gene was used as the internal reference gene. The specificity of the primers was detected by NCBI Primer Blast, and then sent to Shanghai Biological Biotechnology Co., Ltd. for synthesis. The primer information is shown in Table 1.

Table 1 Primer sequences of qRT-PCR

Data analysis

The relative quantitative method was used to calculate the transcription level of the gene by 2^−△△Ct method with yak β-actin as the internal reference gene [26]. One-way ANOVA analysis of variance was performed using SPSS25.0 application software. Duncan's method was used for multiple comparisons when the difference was significant. The analysis results were expressed as Mean ± SEM. P < 0.05 indicated that the difference was significant, and P < 0.01 indicated that the difference was extremely significant.

Results and analysis

Comparative analysis of apparent performance between domestic yak and semi-wild blood yak

As shown in Table 2, the age at first mating of GN and QH yaks was significantly higher than that of SB yaks (P = 0.013). The service lifespan of SB yaks was highly significantly longer than those of GN and QH yaks (P < 0.01). Additionally, the breeding ability of SB yaks was highly significantly greater than that of QH yaks, while QH yaks exhibited a highly significant superiority in breeding capacity over GN yaks (P < 0.01).

Table 2 Comparison of reproductive performance

Table 3 showed that the birth weight, body length, and body length of SB yaks were significantly higher than those of GN yaks, and the birth weight, body length, and body length of GN yaks were significantly higher than those of QH yaks (P < 0.01). The body weight of SB yaks were significantly higher than that of QH yaks (P = 0.028). Additionally, SB yaks displayed extreme dominance over GN yaks in withers height (P = 0.006) and highly significant superiority in chest girth relative to both GN and QH yaks (P < 0.01).

Table 3 Comparison of body weight and body size indexes in 4-year-old yak

As delineated in Table 4, SB yaks and QH yaks exhibited highly significant increases in testicular circumference compared to GN yaks (P < 0.01). Specifically, SB yaks demonstrated superior testicular longitudinal diameter and testicular weight over QH yaks (P < 0.01), while QH yaks significantly surpassed GN yaks in these two metrics (P < 0.01). Furthermore, SB yaks displayed a highly significant advantage in testicular transverse diameter relative to both QH and GN yaks (P < 0.01).

Table 4 Comparison of apparent traits of yak testis

Observation on testicular tissue morphology of domestic yak and semi-wild blood yak

The testicular tissue of yaks is organized into numerous seminiferous tubules, each composed of a basement membrane formed by myoid cells and interstitial cells, along with spermatogenic cells at varying developmental stages. Sertoli cells are interspersed among the spermatogenic cells, while the luminal space of these tubules contains spermatozoa exhibiting distinct morphological forms. As shown in Fig. 1(A-C), the seminiferous tubules of SB yaks exhibit a larger diameter compared to those of GN and QH yak breeds. Figure 1 (D-F) further demonstrates that spermatogenic cells in SB yaks display a more densely arranged cellular organization and larger cell morphology relative to GN and QH counterparts, and a significantly higher sperm density in the tubular lumen compared to domestic yaks.

Fig. 1

Histomorphology observations on testis of domesticated yak and semi-wild blooded yak. Note: A, D: Histomorphology of the testis of the GN yak; B, E: Histomorphology of the testis of the QH yak; C, F: Histomorphology of the testis of the SB yak. SC: sertoli cell; LC: leydig cell; SP: spermatozoon; Sg: Spermatogenic cell. Bar = 200 μm (100x), Bar = 50 μm (400x)

Comparative analysis of testicular reference genome between domestic and semi-wild blood yak

According to Table 5, a total of 40,915,328 unique reads (89.27%) and 42,175,453 total reads (92.01%) were mapped to the reference genome in the GN₁ data. A total of 34,221,850 unique reads (89.97%) and 35,400,235 total reads (93.07%) were mapped to the reference genome in the GN₂ data. A total of 35,327,230 unique reads (89.17%) and 36,396,938 total reads (92.87%) were mapped to the reference genome in GN₃ data. A total of 35,132,313 unique reads (88.23%) and 36,203,804 total reads (90.92%) were mapped to the reference genome in QH₁ data. A total of 32,927,273 unique reads (89.76%) and 33,985,097 total reads (92.64%) were mapped to the reference genome in QH₂ data. A total of 38,050,488 unique reads (89.11%) and 39,258,579 total reads (91.94%) were mapped to the reference genome in QH₃ data. A total of 40,818,978 unique reads (89.77% %) and 42,161,107 total reads (92.72%) were mapped to the reference genome in SB₁ data. A total of 39,892,436 unique reads (89.93%) and 41,180,952 total reads (92.83%) were mapped to the reference genome in SB₂ data. A total of 31,725,343 unique reads (89.26%) and 32,784,445 total reads (92.24%) were mapped to the reference genome in SB₃ data. It could be seen that the accuracy and coverage of the sequencing results were high and could be used for subsequent data analysis.

Table 5 Comparative analysis of testicular reference genomes of yak

Distribution of gene expression between domestic and semi-wild blood yak testicular samples

The violin plot (Fig. 2) showed the distribution of gene expression between 9 testicular samples. The maximum, median, and minimum values of gene expression in 9 samples were similar, indicating that the reliability of gene expression results obtained by sequencing was high, and the gene expression levels between groups were similar.

Fig. 2

Sample expression violin plot

Correlation analysis between domestic and semi-wild blood yak testicular samples

Based on the expression information, used R (http://www.r-project.org/) to carry out principal component analysis (PCA). As illustrated in Fig. 3(A), the PCA scatter plot of samples revealed strong intra-group reproducibility, with SB testicular samples clustering tightly together, indicating high similarity within the group. Distinct inter-group differentiation was also observed. Pearson correlation coefficients were calculated to quantify expression-level relationships between any two samples, visualized as a heatmap in Fig. 3(B). All samples exhibited robust correlations (> 95%), except for GN2, which showed lower correlation with other samples (≥ 63%). It may be due to the hybridization between GN yaks and domestic yaks in different regions. This phenomenon may be attributed to the lag in local domestication processes, or it could result from varying influences of altitude, climate, temperature, and hybridization on different individuals, which has led to the formation of class group within the same breed of domestic yaks in specific regions. However, no significant differences in morphological characteristics are observed among these domestic yak populations.

Fig. 3

Correlation analysis of sample. Note: (A) The PC1 coordinate indicates the first principal component, while the percentage in parentheses indicates the contribution of the first principal component to the sample variance; the PC2 coordinate indicates the second principal component, while the percentage in parentheses indicates the contribution of the second principal component to the sample variance. B The coloured points in the graph indicate the individual samples respectively. Both the horizontal and vertical coordinates in the graph are for each sample, and the colour shade indicates the size of the correlation coefficient between the two samples. The closer the colour is to red, the higher the correlation, and the closer the colour is to blue, the lower the correlation

Basic analysis of differential genes in testis of domestic and semi-wild blood yak

In this study, raw counts were analyzed by edgeR [27] software. Based on certain standardized processing, FDR < 0.05 and |log₂FC|≥ 2 were used as criteria to screen DEGs between domestic and semi-wild blood yak testis groups. From Fig. 4(A), comparative analysis of DEGs revealed distinct transcriptional profiles among yak populations. When comparing GN yaks versus SB yaks, we identified 2,403 DEGs (687 up-regulated; 1,716 down-regulated). QH yaks versus SB yaks, with 4,428 DEGs detected (1,493 up-regulated; 2,935 down-regulated). GN yaks versus QH yaks, with 4 genes showing differential expression (2 up-regulated and 2 down-regulated).According to Fig. 4(B), 1848 common DEGs were obtained among GN, QH, and SB yaks.

Fig. 4

A Identification of DEGs. B DEGs venn diagram. Note: Horizontal coordinates: samples compared in pairs; vertical coordinates: number of DEGs; red represents up-regulated DEGs; orange represents down-regulated DEGs; blue represents no difference

GO enrichment analysis of DEGs in the testis of domestic yak and semi-wild blood yak

A total of 1848 DEGs were enriched between GN, QH, and SB yak groups. In this experiment, Goatools application software was used to perform GO enrichment analysis on these DEGs, which were significantly enriched into 1736 GO terms, including 1314 biological processes, 192 cellular components, and 230 molecular functions.

Figure 5 showed the second-level histogram of GO functional enrichment results classification. We summarized the results of GO enrichment differences and removed duplicate genes to obtain 51 DEGs (STK4, NOTCH1, LRP2, TPPP3, TBX3, NIPBL, SIN3 A, SMAD3, USP47, ROS1, SRF, CTC1, ARID4 A, INPP5B, MAP2 K1, NKX3 - 1, TGFBR1, PUM1, CLOCK, PAFAH1B3, ZBTB16, HERC2, CXADR, NOTCH2, ASH1L, BRCA2, LHFPL2, ACSL4, STAT3, BMP7, CTNNB1, ARSA, B4GALT1, TUT4, CREBRF, PSAP, BMPR2, TERF1, SMC1 A, FANCM, KIT, DLG1, PSMD13, NDC1, SMALF1, HIFFM1, NRD4) highly associated with spermatogenesis. These genes were mainly enriched in Cell, Cell part, Organelle, Organelle part, Membrane, and Macromolecular complex, and were related to biological processes such as Binding, Catalytic activity, Molecular function regulator, Transporter activity, Signal transducer activity, Molecular transducer activity, Nucleic acid binding transcription factor activity, Transcription factor activity—protein binding, and Structural molecule activity. They were involved in the physiological processes of Cellular process, Single-organism process, Metabolic process, Biological regulation, Regulation of biological process, Response to stimulus, Cellular component organization or biogenesis, Positive regulation of biological process, Localization, Multicellular organismal process, Signaling, Negative regulation of biological process, Immune system process, Locomotion, Multi—organism process, Biological adhesion, Reproductive process, Reproduction, Growth, Behavior, and Rhythmic process in GN/QH and SB yaks. Among them, TPPP3, SMAD3, BMP7, PAFAH1B3, and ARSA exhibited up-regulated expression in SB yak testes.

Fig. 5

Classification of DEGs GO enrichment analysis classification secondary histograms

Screening of target genes in domestic yaks and yaks of semi-wild origin

According to Fig. 6(A), PCA divided DEGs into two principal components, of which the PC1 contribution rate was 86.4%, the PC2 contribution rate was 9%, and the total contribution rate was 95.4%. And domestic yaks and semi-wild blood yaks have good discrimination. Figure 6(B) showed that the PLS-DA analysis results of 51 DEGs, and the contribution rate of PLS-DA model was 92.4%. Among them, the differential genes of domestic yak and semi-wild blood yak were obviously separated, indicating that the model had good classification prediction ability and stability. Figure 6(C) showed the screening of DEGs using the variable projection importance of PLS-DA, with VIP > 1 and P < 0.05 as the screening criteria. A total of 31 target genes (BMPR2, DLG1, KIT, ASH1L, INPP5B, BMPR1 A, SMC1 A, STK4, NSUN2, ACSL4, TGFBR1, LHFPL2, HERC2, SRF, STAT3, CREBRF, NIPBL, BMP7, TBX3, B4GALT1, USP47, HIF1 A, ARID4 A, PSMD13, PUM1, ARSA, CTC1, CTNNB1, FANCM, SMAD4, SMAD3) were screened out.

Fig. 6

A PCA map of spermatogenesis DEGs. B PLS-DA map of spermatogenesis DEGs. C VIP map of spermatogenesis DEGs

KEGG enrichment analysis of selected spermatogenesis-related genes

KEGG enrichment analysis was performed on the selected 31 target genes affecting yak spermatogenesis, and Fig. 7 showed the first 20 pathways with significant enrichment. The target genes were significantly enriched in the Hippo signaling pathway, Signaling pathways regulating the pluripotency of stem cells, and the TGF-beta signaling pathway.

Fig. 7

Bubble plot of significance of KEGG enrichment affecting spermatogenesis in yak

Construction and visualization analysis of DEGs PPI network in the testis of domestic and semi-wild blood yak

Using 51 screened DEGs associated with reproduction, a PPI network was constructed based on network topological parameters (Fig. 8), which identified five core regulatory genes (CTNNB1, SMAD4, STAT3, SMAD3, NOTCH1) involved in spermatogenesis. Among them, CTNNB1, SMAD4, STAT3, SMAD3 and NOTCH1 had regulatory relationships with each other. CTNNB1 had a regulatory relationship with HIF1 A, TGFBR1, BMPR2, KIT, BMPR1 A, BRCA2, NOTCH2, MAP2 K1, TBX3 and CLOCK. SMAD4 had a regulatory relationship with HIF1 A, TGFBR1, BMPR2, KIT, BMPR1 A, BRCA2 and MAP2 K1. STAT3 had a regulatory relationship with HIF1 A, TGFBR1, BMPR2, KIT, NOTCH2, MAP2 K1 and ROS1. SMAD3 had a regulatory relationship with HIF1 A, TGFBR1, BMPR2, BMPR1 A and SRF. NOTCH1 had a regulatory relationship with HIF1 A, TGFBR1, KIT, NOTCH2, and SRF.

Fig. 8

PPI network of DEGs in testis

Analysis of DEGs expression levels affecting testicular development and spermatogenesis in domestic yaks and semi-wild yak

It could be seen from Table 6 that the expression levels of TPPP3, SMAD3, BMP7, PAFAH1B3, and ARSA in SB yaks testis were significantly higher than those in GN and QH yaks (P < 0.01). The expression levels of SMAD4, CTNNB1, and STAT3 in the testis of GN and QH yaks were significantly higher than those of SB yaks (P < 0.01). However, no significant difference was no significant difference in NOTCH1 expression (P = 0.327). Consequently, eight key genes (TPPP3, SMAD3, PAFAH1B3, BMP7, ARSA, CTNNB1, SMAD4, and STAT3) influencing SB yak spermatogenesis were identified.

Table 6 Comparative analysis of DEGs expression

Verification of DEGs in testis of domestic and semi-wild blood yak by qRT-PCR

Eleven DEGs (FGFR1, MAP2 K1, STAT3, PIK3 CA, ACVR2 A, HIF1 A, BMPR1 A, ROS1, zfhx3, PAK2, ROBO2) were selected, Yak β-actin was used as the internal reference gene, the same yak testis RNA sample as RNA-Seq sequencing was used for fluorescence quantitative PCR verification. The verified trend of gene expression was consistent with the sequencing results, which proved that the results of RNA-seq transcriptome data were reliable (Fig. 9).

Fig. 9

Histogram of DEGs qRT-PCR validation

Discussions

Spermatogenesis is a complex process that involves structural and biochemical changes in the seminiferous epithelium of adult animals and promotes the spread of genetic resources [28]. Mammalian spermatogenesis is regulated by paracrine, autocrine, and endocrine pathways as well as a series of structural elements and chemical factors that regulate the activity of somatic and germ cells [29,30,31]. There are many changes in gene expression in the process of spermatocyte to spermatogonia transformation [32]. The testicular development of 3 ~ 5 month-old yak is regulated by gene dynamic expression. Comparative analysis of DEGs reveals substantial transcriptional shifts during three critical phases of spermatogenic cell maturation, particularly marking the transition from spermatogonia to spermatocytes [10]. And the communication network orchestrating cellular activities during spermatogenesis exhibits a highly intricate architecture [33, 34]. A complex regulatory interplay exists between genotypic determinants and phenotypic expression [32]. Transcriptomics is essential for understanding the regulation of genes in animal growth and development [35]. In this study, we compared the apparent characteristics and histological differences between domestic yaks and semi—wild blood yaks through measurement, follow—up investigations, data record access, H.E staining, and transcriptomics analysis. We also clarified the key genes and important pathways that affect the spermatogenesis of SB yaks. This work can assist breeders in selecting candidate genes and has significant implications for enhancing the reproductive potential of yaks.

Comparison of apparent performance and histology between domestic yak and semi-wild blood yak

The reproductive potential of bulls depends on their mating ability and fertilization ability [36]. The early development of body weight and the time to reach sexual maturity directly affect the mating ability and fertilization ability of bulls. The results of this experiment showed that the body size and body weight of SB yak are significantly or extremely significantly higher than GN and QH yak. Although there is no significant difference in sexual maturity time among the three kinds of yaks, the first mating time of SB yak is significantly lower than GN and QH yak, and the time of utilization and mating ability are significantly or extremely significantly higher than GN and QH yak. Studies have shown that larger yaks enter puberty earlier than smaller yaks [37]. Datong yak is the first artificially cultivated hybrid of wild and domestic yak, showing excellent growth characteristics and production performance. The high growth, development, and production rate of Datong yak indicate that it is feasible to improve the traits of domestic yak by using wild yak resources, and the potential importance of using wild yak genetic resources for yak breeding is clarified [38], which is consistent with the growth characteristics of semi-wild blood yak in this experiment.

The size and function of the testis were crucial indicators for evaluating and breeding male animals, and they were highly heritable [39]. The apparent traits of SB testis in this experiment are significantly higher than GN and QH yak, and the testicular tissue morphology shows that the seminiferous tubules of SB yak testis are larger than GN and QH yak. Foote et al. founded that the testicular size and body weight of Sahiwal bulls increased in a curvilinear manner from 0 to 100 months of age. They observed a significant positive correlation between testicular development and body weight. After 100 months of age, the testicular size of Sahiwal bulls stabilized, while ejaculatory volume and sperm density increased significantly. However, there was no significant effect on sperm motility, plasma membrane integrity, or normal replacement rate [40]. T.A. Bongso, M.R., et al. studied the relationship between testis size and body weight in goats aged 3 to 28 months. The results indicated a significant positive correlation between testis size and body weight. Notably, testis size increased significantly between 6 and 7 months of age, accompanied by a sudden growth of seminiferous tubules, which promoted spermatogenesis [41]. The size of the testis would increase with the size of the body [42]. Research had shown that larger mammals require proportionately larger organs to fulfill their physiological needs. The size and weight of the testis were positively correlated with both the rate of testicular sperm production and sperm reserves [43]. Testicular size was strongly linked to daily sperm production, as the amount of sperm produced per unit of testicular volume remains constant, and the overall sperm output primarily depends on testicular size [36]. A large number of studies had shown that mammalian body size and body weight were significantly positively correlated with testicular development [44, 45]. The weight and size of the testis directly affected the discharge and quality of semen [46,47,48,49,50]. Observations of body size, body weight, age of utilization, testicular characteristics, and tissue morphology indicate that SB yak exhibits strong heterosis. However, the reproductive capacity of SB yaks was significantly lower than that of GN and QH yaks. This may be due to the fact that the testicles produce sperm while simultaneously generating substantial amounts of semen and hormones, resulting in a greater dilution effect [51].

Transcriptome difference analysis between domestic yak and semi-wild blood yak

Through transcriptome screening, eight key genes impacting spermatogenesis in SB yak were identified: TPPP3, SMAD3, PAFAH1B3, BMP7, ARSA, CTNNB1, SMAD4, and STAT3. TPPP3 belongs to the tubulin polymerization-promoting protein (TPPP) family. It promotes cell proliferation by influencing the replication and separation of centrosomes during mitosis [52]. Microtubules (MTs) are polymers composed of tubulin heterodimers. These structures are distributed throughout cells and play crucial roles in diverse cellular functions. They participate in the elongation and guidance of neurons within the growth cone, the maintenance of cell shape and intracellular transport, and the formation of the mitotic spindle, which is indispensable for chromosome separation during cell division [52,53,54]. TPPP3 binds to microtubules and promotes their bundle and network formation in a concentration-dependent manner, increasing the level of microtubule acetylation [55]. Different domains of TPPP3 protein promote microtubule elongation, and TPPP3 plays a role in developing a ciliated cell basal body array, providing a loose but limited microtubule network to cope with external loads [56]. In the absence of TPPP3 protein, in addition to cell cycle arrest and down-regulation of several cyclin proteins, HeLa cells also observed abnormalities in the metaphase and/or anaphase of mitosis, forming multipolar spindles [52]. This resulted in abnormal centrosome amplification, thereby increasing the frequency of abnormal mitosis and chromosome segregation errors [57, 58]. Studies have demonstrated that TPPP3 dynamically regulates β-catenin activity in the endometrial microenvironment during embryo implantation, coordinating cellular adhesion and signaling pathways critical for successful implantation [59], and the Wnt/β-catenin signaling pathway is involved in regulating the occurrence of mature sperm cells [60, 61]. β-catenin is an important effector of the Wnt/β-catenin signaling pathway and a regulator of cell adhesion [62]. Its role in testicular development may be more as a cytoplasmic adaptor protein [63]. It is also an important molecular link that integrates the signal events necessary for the development and maturation of germ cells after sertoli cell-germ cell adhesion and meiosis. The β-catenin complex of spermatogenic cells and the apical ectoplasmic specialization of the β-catenin complex on sertoli cells on the surface of ES triggers a signal cascade to regulate the development and maturation of spermatogenic cells at the late stage of meiosis [64].

The Smad protein serves as a classic intracellular mediator in the signal—transduction process of the TGF-β superfamily. The signal transduction initiated by ligands such as activin, TGF-β, and BMP, in combination with the regulated expression and differential utilization of specific Smad proteins, plays a vital role in the accurate regulation of the responses of testicular cells to the ligands of the TGF-β superfamily [65]. The Smad protein is predominantly located in the cytoplasm of meiotic germ cells, Sertoli cells, and Leydig cells, where it serves a vital function in testicular development and spermatogenesis [66]. Smad3 and Smad4 are key signal—transduction proteins within the TGF-β signaling pathway. Signaling from the TGF-β superfamily is essential for the determination of primordial germ cells, the preservation of the integrity of fetal testicular development, and the regulation of spermatogenesis in adult testes [67, 68]. Smad3 gets activated by ligands like TGF-β, activin, nodal, myostatin, and GDF3/9. Subsequently, it forms a complex with Smad4 and translocates into the nucleus to control the expression of target genes, thereby playing a significant part in testicular development, maturation, and the maintenance of its function [69, 70]. The activation of the SMAD3 and AKT signaling pathways serves as the crucial factor in sustaining the long—term proliferation and undifferentiated state of human sperm stem cells. Through the establishment of a well-defined culture system that incorporates specific growth factors and hydrogels, human sperm stem cells have been successfully amplified and managed to remain in an undifferentiated state for as long as 2 months [71]. SMAD4 acts as a key regulator in signal transduction within the TGF-β superfamily. It experiences a cycle of ubiquitination and deubiquitination, which is regulated by the ubiquitin ligase Ecto/TIF1γ and the deubiquitinase FAM/USP9x. This dynamic regulatory process controls the transcriptional response to TGF-β by promoting the assembly and dissociation of the Smad signaling complex [72]. Double deletion of Smad4 in Sertoli and Leydig cells results in fetal testicular tubular dysplasia and mild testicular dysplasia in adulthood. These manifestations include reduced testicular volume, partial seminiferous tubule dysplasia, and decreased sperm production [73]. The experimental results demonstrated that SMAD3 expression was upregulated in SB yak testicular tissues, whereas SMAD4 expression exhibited a downregulated profile. This observation suggests that SMAD3 might compensate for the diminished SMAD4 levels at the critical regulatory element of the Fshb promoter, potentially maintaining transcriptional regulation efficiency [74]. Bone morphogenetic proteins (BMPs), which are part of the TGF-β superfamily, play a vital role in male fertility as well as the development of fetal germ cells and somatic cells post-birth [75]. Among them, BMP7 serves as a key signaling molecule in testicular development and spermatogenesis. It selectively activates particular SMAD receptor proteins, namely SMAD1 and SMAD5, and participates in the regulation of testicular development and spermatogenesis [76], Additionally, BMP7 promotes the proliferation of Sertoli cells [77], and is involved in regulating the proliferation and differentiation of spermatogonia [78]. BMP8a is capable of enhancing the proliferation of germ cells via the SMAD2/3 signaling pathway and guiding the differentiation of germ cells through the SMAD1/5/8 signaling pathway [79]. When one allele of Bmp7 is removed, it exacerbates the phenotype of Bmp8a deletion mutants in spermatogenesis and the epididymis. This finding implies that BMP7 and BMP8 transmit signals in these two systems by using the same or similar receptors [80].

Platelet-activating factor (PAF) can influence sperm motility and acrosome function and has a role in the development of the animal reproductive system [81]. The mRNA of platelet-activating factor acetylhydrolase 1b (PAFAH1B) is found in different stages of mouse testicular development and spermatogenesis. These stages involve spermatogonia, spermatocytes, and early sperm cells during meiosis as well as the early stages of haploid germ cell development [82]. Research has shown that the PAFAH1B3 protein is located on the mitotic spindle of spermatocytes undergoing meiosis, which suggests that it plays a significant part in spermatogenesis [83, 84]. The expression level of the PAFAH1B3 gene would affect the process of testicular development and spermatogenesis [85], and knockout of the PAFAH1B3 gene would lead to spermatogenesis disorder and increased apoptosis of sperm cells [86]. Arylsulfatase A (ARSA) is a lysosomal enzyme found in the male reproductive system [87]. It is present in the lysosomes of testicular Sertoli cells as well as the acrosomes of sperm cells and spermatozoa [88]. Sulfogalactosylglycerolipid (SGG) is selectively present in male germ cells. Both SGG and ARSA have a direct affinity for the zona pellucida (ZP) and play important roles in spermatogenesis and sperm-egg interaction [89]. SGG serves as the physiological substrate of lysosomal ARSA in sertoli cells. In ARSA-knockout mice, SGG accumulates in sertoli cells, leading to lysosomal accumulation. This accumulation damages the function of sertoli cells, resulting in impaired spermatogenesis and near-zero in vitro spermatogenesis [90]. Studies have shown that ARSA can degrade SGG, prevent SGG-induced cytotoxicity, and produce palmitylpalmitoylglycerol (PPG). This PPG is used to synthesize SGG in a new generation of primary spermatocytes [91]. Signal transducers and activators of transcription 3 (STAT3) is a multifunctional transcription factor. It can either activate or inhibit the transcription of target genes, thus influencing cell proliferation, survival, and apoptosis [92, 93]. STAT3 is predominantly expressed in sertoli cells [94], where it regulates the differentiation of spermatogonial stem cells through modulation of the downstream target gene Ngn3 [95]. Notably, pharmacological inhibition of STAT3 phosphorylation has been shown to downregulate Zfp637 expression, thereby impairing the normal spermatogenesis process [96]. This study revealed that in the testicular tissues of SB yaks, genes such as TPPP3, SMAD3, PAFAH1B3, BMP7, and ARSA exhibited significant upregulation, while the expression levels of CTNNB1, SMAD4, and STAT3 were markedly downregulated. These gene expression characteristics were highly consistent with the histomorphological observations of SB yak testicular tissues, specifically manifested as enlarged luminal diameters of seminiferous tubules, tightly arranged spermatogenic cell layers with increased nuclear volumes, and a significantly higher sperm density in the tubular lumen compared to domestic yaks.

In summary, this experiment compared the apparent performance and testicular tissue morphological characteristics of GN, QH, and SB yaks. During this process, we elucidated the strong advantages exhibited by male SB yaks in terms of growth performance, reproductive performance, and testicular tissue morphology. To further explore the intrinsic mechanisms of SB yaks in reproduction, we employed transcriptome analysis and successfully screened out 8 dominant genes closely related to testicular development and spermatogenesis in SB yaks. Notably, these genes showed significant enrichment in the Hippo signaling pathway, the signaling pathway regulating stem cell pluripotency, and the TGF—beta signaling pathway. This discovery not only implies that these signaling pathways may play a crucial role in the testicular development and spermatogenesis of SB yaks, but also provides a new perspective for our understanding of the male yak's reproductive physiology. It helps us to further reveal the molecular mechanisms of yak reproductive regulation and provides a solid theoretical basis for the genetic improvement and enhancement of reproductive efficiency of yaks.

Ethics statement

All animal experiments were approved by the Ethics Committee of Gansu Agricultural University, 2024 (GSAU-Eth-VMC- 2024–009).