Single-cell analysis identifies critical regulators of spermatogonial development and differentiation in cattle-yak bulls

Spermatogenesis is a continuous process in which functional sperm are produced through a series of mitotic and meiotic divisions and morphological changes in germ cells. The aberrant development and fate transitions of spermatogenic cells cause hybrid sterility in mammals. Cattle-yak, a hybrid animal between taurine cattle ( Bos taurus ) and yak ( Bos grunniens ), exhibits male-specific sterility due to spermatogenic failure. In the present study, we performed single-cell RNA sequencing analysis to identify differences in testicular cell composition and the developmental trajectory of spermatogenic cells between yak and cattle-yak. The composition and molecular signatures of spermatogonial subtypes were dramatically different between these 2 animals, and the expression of genes associated with stem cell maintenance, cell differentiation and meiotic entry was altered in cattle-yak, indicating the impairment of undifferentiated spermatogonial fate decisions. Cell communication analysis revealed that signaling within different spermatogenic cell subpopulations was weakened, and progenitor spermatogonia were unable or delayed receiving and sending signals for transformation to the next stage in cattle-yak. Simultaneously, the communication between niche cells and germ cells was also abnormal. Collectively, we obtained the expression profiles of transcriptome signatures of different germ cells and testicular somatic cell populations at the single-cell level and identified critical regulators of spermatogonial differentiation and meiosis in yak and sterile cattle-yak. The findings of this study shed light on the genetic mechanisms that lead to hybrid sterility and speciation in bovid species.


INTRODUCTION
Fertility relies on efficient spermatogenesis, which is a complex cellular differentiation process involving 3 major phases: mitosis of spermatogonia, meiosis of spermatocytes and spermiogenesis (de Rooij and Russell, 2000, Hess and Renato de Franca, 2008).Bulls produced by crossing taurine cattle (Bos taurus) and yak (Bos grunnines) display complete sterility due to spermatogenic arrest at the meiosis stage (Li et al., 2020b).Female cattle-yak hybrids are fertile and widely used for milk production (Barsila et al., 2015); however, the sterile phenotype of the males prevents attempts to fix desired traits such as high milk yield through hybridization in yak breeding.The genetic basis underlying hybrid male sterility in cattle-yak has been extensively studied in recent years, yet the candidate genes that have functional roles in the reproductive isolation and speciation of these 2 species remain undefined.Spermatogenic defects in interspecific hybrid animals are often associated with various problems in meiosis.Meiosis prophase I contains 5 successive events, including leptotene, zygotene, pachytene, diplotene, and diakinesis (Rivenzon, 1954).Silver and blue fox hybrids are sterile because of the developmental arrest of primary spermatocytes at the pachytene stage due to asynapsis (Yang et al., 2020, Bugno-Poniewierska et al., 2021).In an interspecific cross between the domestic cat and the Jungle cat, meiotic sex chromosome inactivation failed to occur, resulting in the death of the pachytene spermatocytes (Bredemeyer et al., 2021).Meiosis failure is also the primary cause of spermatogenic arrest in cattleyak, but the underlying cause appears to be different.For example, meiotic chromosomes formed synapsis,

Single-cell analysis identifies critical regulators of spermatogonial development and differentiation in cattle-yak bulls
and a few diplotene spermatocytes were detected in the testis of cattle-yak (Li et al., 2020b).The major defects in these hybrid animals appear to arise from problems associated with DNA double-strand break repair (Wu et al., 2023).Additionally, recent findings from transcriptome analysis of yak and cattle-yak testicular tissues indicated that genes regulating spermatogonial fate decisions are dysregulated in hybrid animals, particularly those directing spermatogonial differentiation and meiosis entry (Endo et al., 2015, Cai et al., 2017).Together, these findings suggest that the homeostasis of the spermatogonial population and meiosis progression are impaired in cattle-yak.
The spermatogonial population contains spermatogonial stem cells (SSCs) and transit-amplifying progenitors, which are capable of committing to meiosis after several rounds of division (Brinster and Zimmermann, 1994).The spermatogonial lineage is specified from a transient germ cell population termed gonocytes or prespermatogonia during neonatal and prepubertal development (Orth et al., 2000).This transition occurs at 3-5 mo of age in Holstein bulls (Curtis and Amann, 1981) and at 5-8 mo of age in yak bulls (Wang et al., 2019a).A foundational SSC pool is established to maintain lifelong spermatogenesis under the guidance of niche factors.SSCs differentiate to provide progenitors, which enter the differentiation pathway and progress through successive differentiation steps to enter meiosis (de Rooij and Russell, 2000).Interactions among different germ cell subpopulations and communications between germ cells and somatic cells are pivotal for spermatogonial differentiation and meiosis (Raverdeau et al., 2012, Endo et al., 2017).Because of differences in cell composition and differentiation trajectory, the molecular signatures of different spermatogenic cell subtypes and the cellular interactions of testicular cells in cattle-yak are not well understood.
Single-cell RNA sequencing (scRNA-seq) can effectively compensate for the heterogeneity of gene expression in single cells ignored by conventional sequencing and thus obtain the expression profile of the whole transcriptome at the single-cell level (Tang et al., 2009).Recently, a single-cell RNA-seq survey of spermatogenic cells in cattle-yak revealed that interactions between germ cells and their niche environment were abnormal (Mipam et al., 2023); however, due to limitations in cell number and in-depth analysis, the subtypes of spermatogenic cells and genes regulating cell-to-cell commutations in spermatogenesis remain unclear.In the present study, we performed scRNA-seq analyses on the testes of 12-mo-old yak and age-matched cattle-yak to obtain gene expression profiles of different spermatogenic lineages and somatic cells.Using comparative analysis, we identified abnormally expressed genes that lead to unbalanced spermatogonial fate decisions and meiosis arrest in cattle-yak.In particular, the results of these analyses revealed that intercellular signaling was dysregulated in cattle-yak, highlighting the intrinsic defects related to spermatogonial self-renewal, differentiation, and meiosis in hybrid sterility.

Animals
All animal experiments were performed in accordance with the Guide for the Care and Use of Laboratory Animals and were approved by the Animal Ethics and Welfare Committee at the Northwest Institute of Plateau Biology, Chinese Academy of Sciences.Because the timing of spermatogenic cell development was similar between yak and cattle-yak, we selected 4 12-moold yak and cattle-yak (male cattle crossed with female yak) for single-cell RNA sequencing analysis in Qinghai Province, China.The testis was obtained by standard castration and washed 3 times with PBS after removing the epididymis and tunica vaginalis.Fresh testicular samples were cut into small pieces: one part was used for cell suspension, and the other part was fixed in Bouin's fluid for histological analysis.

Histological analysis
Testes were fixed in Bouin's solution for 12 h and then dehydrated and embedded in paraffin as described previously (Zhao et al., 2023).The fixed tissues were cut into 5 μm sections by a microtome (Leica RM2235).Sections were deparaffinized and stained with hematoxylin and eosin (H&E).After drying, the slices were sealed with neutral balsam (Biotopped) and analyzed under a Leica microscope (Germany).

Preparation of single-cell suspensions
Testicular tissues were digested to obtain a singlecell suspension in type IV collagenase (Sigma, c5138-1 G), 0.25% trypsin (Gibco, 25200-072) and 1 mg/mL hyaluronidase (Sigma, H3506-1 G).During this period, pipetting continuously to completely release the germ cells in the seminiferous tubules.One equivalent of FBS (Gibco, 16000-044) was added to the cell suspension to stop digestion.The cells were then filtered and centrifuged at 1500 rpm for 8 min.The cells were resuspended in 1 mL of DPBS (Gibco, 2534480) and treated with 3 mL of red blood cell lysis buffer (Sorolebol, R1010) to remove erythrocytes.The cells were then resuspended in cold sample buffer provided by the Single-Cell Capture and cDNA Synthesis Kit (Cat.No. 633731) in the

Construction and sequencing of scRNA-seq libraries
We constructed single-cell libraries according to the single-cell 3′ whole transcriptome amplification (WTA) protocol (Doc ID: 210966).The cell suspensions were loaded into a special U-shaped groove using a BD Rhapsody Cartridge Kit (Cat.No. 633733) and naturally precipitated in the micropores.Then, an excess of magnetic beads with cell labels and unique molecular identifier (UMI) sequences were loaded to ensure that most magnetic beads fell into the micropores to capture the mRNA released after cell lysis.After the magnetic beads were extracted and washed, the standard sequencing library was constructed by reverse transcription with a BD Rhapsody cDNA Kit (Cat.No. 633773), and cDNA was amplified with a BD Rhapsody WTA Amplification Kit (Cat.No. 633801).Libraries were sequenced at 150 bp on a NovaSeq sequencer at Novogene Bioinformatics Technology Co., Ltd.

Single-cell RNA-seq analysis
We used a newly assembled domestic yak genome from our previous works as the reference genome (v2.7.1a, https: / / github .com/alexdobin/ STAR) (Gao et al., 2022).After importing the matrices into Seurat (v.4.0.4,https: / / github .com/satijalab/ seurat), we obtained distributions of gene counts, UMI counts, and mitochondrial transcript levels.Based on previous single-cell analysis methods (He et al., 2023, Jia et al., 2023, Tao et al., 2023, Li et al., 2024), the low-quality cells were removed by the following conditions: all genes were expressed in at least 3 cells, each cell contained at least 200 genes, and the percentage of mitochondria was less than 20%.The filtered data were normalized by "LogNormalize" and calculated for highly variable genes by "FindVariableFeatures" (the threshold was set at 2,000).We then performed the analysis by principal component analysis, calculated the nearest neighbor graph based on Euclidean distance in principal component analysis space and selected 15 principal components.After that, we used the "RunTSNE" function and set the resolution parameter to 0.6 to cluster the cells.The cell types of each cluster were identified by feature plots and violin plots using accepted wellknown marker genes.To further separate Leydig cells from peritubular myoid cells (PTMCs), after cell cluster analysis of Leydig cells and PTMCs, the cell type of each cluster was identified again by feature plots and violin plots using accepted well-known marker genes.The differentially expressed genes (DEGs) for each cell type were determined by log2(fold change) > 0.25 using the default Wilcoxon rank-sum test in the "FindAll-Markers" function.

Integration analysis of scRNA-Seq data
The matrices of the same yak and cattle-yak cell types were extracted and integrated, and the "RunTSNE" function was used to perform the cell cluster analysis again.The resolution parameter of the combined germ cells was 2.0, and that of the other cell types was 0.6.The first cell cluster analysis of the integrated germ cells revealed 18 clusters, 4 of which (clusters 7, 8, 13, and 15) highly expressed Sertoli cell, Leydig cell and PTMC markers (Supplementary Table 5).To remove somatic cells, we removed these 4 clusters and performed germ cell clustering analysis again, with the resolution parameter set to 2.0.DEGs between the same cell types were identified by the "FindMarker" function with the parameters "min.pct = 0.25," "logfc.threshold = 0.20" and "p_val < 0.05."DEGs were subjected to GO and KEGG analyses with g:Profiler (https: / / biit .cs.ut.ee/gprofiler/ gost), and the cutoff for significantly enriched terms was p_val < 0.05.The results of functional enrichment were visualized by the ggplot2 package (v3.3. 5, https: / / r -graph -gallery .com/ggplot2package .html).Heatmaps were generated with https: / / www .bioinformatics.com.cn(last accessed on 10 Nov 2023), an online platform for data analysis and visualization.

Cell trajectory analysis
The Monocle2 package (Qiu et al., 2017) (v2. 18. 0, http: / / cole -trapnell -lab .github.io/monocle -release/ docs/ ) was used to reconstruct the germ cell differentiation trajectory.We identified DEGs that varied in pseudotime order and set q_val < 0.01 for pseudotime analysis, applied the "DDRTree" method of the "re-duceDimension" function for dimensionality reduction, and used the "orderCells" function to sort cells and complete trajectory construction.Visual tracks were generated in pseudotime order and colored according to pseudotime values, cell states, and original identities.Pseudotime heatmaps were generated using the "plot_pseudotime_heatmap" function based on DEGs among cell clusters in pseudotime with q_val <0.1.We examined the associations of the expression of germ cell marker genes with cell states by identifying the genes involved.

Integration analysis of scRNA-Seq data from 12-moold and 24-mo-old cattle and yak
We downloaded the published scRNA-Seq data of 24-mo-old yak and cattle-yak testicular tissues (Mipam et al., 2023) and used the same methods and parameters as before to identify germ cells and somatic cells.After that, we integrated 12-mo-old and 24-mo-old germ cells and used the "RunTSNE" function with the resolution parameter set to 0.6.The first cell cluster analysis of the integrated germ cells revealed 16 clusters, 4 of which (clusters 4, 5, 13, 14 and 15) highly expressed Sertoli cell, Leydig cell, PTMC and immune cell markers (Supplementary Table 20).We removed these 5 clusters and performed germ cell clustering analysis again, with the resolution parameter set to 0.6.DEGs between the same cell types were identified by the "FindMarker" function with the parameters "min.pct = 0.25," "logfc.threshold = 0.20" and "p_val < 0.05."DEGs between 12-mo-old and 24-mo-old spermatogenic cells were identified by the "FindMarker" function with the parameters "min.pct = 0.25," "logfc.threshold = 0.20" and "p_val < 0.05."

Intercellular communication network analysis
Intercellular communication analysis was performed using the CellChat package (v1.5.0) (Jin et al., 2021).Since this package only provides databases for mouse and human ligand-receptor interactions, we used g:Profiler (https: / / biit .cs. ut.ee/gprofiler/gost) to convert yak gene IDs to human gene IDs and selected the database for human gene IDs.The conversion rate was 93.05% and 90.16% for yak and cattle-yak.We first performed intercellular communication analysis on data from 12-mo-old yak and cattle-yak.The "com-puteCommunProb" and "aggregateNet" functions were used to compute communication probabilities and infer cellular communication networks.Afterward, the analyzed data were integrated to compare the differences between 12-mo-old yak and cattle-yak.The "compare-Interactions" function was used to compare the number and strength of interactions between the cellular communication networks of 12-mo-old yak and cattle-yak.Based on the intensity of outgoing and incoming signals in 2D space, the "netAnalysis_signalingRole_scatter" function identified cell populations that significantly varied between 12-mo-old yak and cattle-yak in terms of transmitted or received signals.The "rankNet" function was used to compare the overall information flow of each signal path, and the "netAnalysis_signal-ingRole_heatmap" function was used to compare the outgoing or incoming signals.The "netAnalysis_contribution" function was used to calculate the contribution of each ligand-receptor to the overall signaling pathway (the ratio of the total communication probability of the inferred network for each ligand-receptor pair to the total communication probability of the signaling pathway).

scRNA-Seq analysis of testicular cells from 12-moold yak and cattle-yak
Previous studies have shown that the seminiferous tubules of 12-mo-old yak and cattle-yak contain enriched numbers of spermatogonia and spermatocytes (Li et al., 2020b, Wu et al., 2023); therefore, we conducted scRNA-seq analysis on these samples to eliminate the interference from advanced germ cells in yak.H&E staining confirmed that the seminiferous tubules of both yak and cattle-yak contain morphologically similar germ cells (Figure 1 A).After quality control, unbiased cluster projections onto t-distributed stochastic neighbor embedding (t-SNE) analysis identified 16 and 13 clusters of 13385 yak and 10750 cattle-yak testicular cells, respectively (Fig. S2 A).To further define the cell type of each cluster, we used marker genes and DEGs to examine the molecular identities of each testicular type.Sertoli cells (SOX9), PTMCs and Leydig cells (ACTA2 & DCN), endothelial cells (CDH5), immune cells (RGS1), and germ cells (DAZL) were detected in yak and cattle-yak testes (Figure 1 B, Supplementary Table 1 and 2).The proportions of germ cells, PTMCs, Leydig cells, endothelial cells and immune cells were lower, while the proportions of Sertoli cells were greater in the cattle-yak testis than in the yak (Table 1).Each cell type exhibited its own gene expression signatures, as revealed by the heatmap of the significant DEGs (Figure 1 C, Supplementary Table 3 and 4;).The number of DEGs in each cell type in yak and cattle-yak were as follows: Sertoli cells (1019, 908), PTMCs and Leydig cells (764, 754), endothelial cells (563, 647), immune cells (524, 424), and germ cells (3091,3241).

Differences in gene expression among different spermatogenic cell subpopulations
Next, we examined the differences in gene expression in the spermatogenic cells of cattle-yak.A total of 477 genes were identified to be differentially expressed, among which 255 were downregulated and 222 were upregulated (Figure 3 A, Supplementary Table 10).Functional enrichment analysis of the DEGs revealed that the upregulated genes were mostly associated with the cell cycle, while the downregulated genes were associated with cell differentiation (Figure 3 B, Supplementary Table 11).The upregulated genes included DMRT1 (Zhang et al., 2016), ASB9 (Li et al., 2023), and EZH2 (Cai et al., 2020), which are associated with spermatogonia self-renewal.Downregulated genes, including PHF13 (Bördlein et al., 2011), COL1A1 (Chen et al., 2012) and STRA8 (Zhou et al., 2008), were associated with spermatogonial differentiation.Several genes controlling cell proliferation and apoptosis, including SPOCD1 (Liu et al., 2018), HMGB3 (Zheng et al., 2018), KDM1B (Wang et al., 2019b) and SIK3 (Ponnusamy et al., 2020), were also found to be differentially expressed.

Comparison of testicular somatic cells in yak and cattle-yak
Gene expression in germ cells is a coordinated process partially directed by signals from somatic cells.To examine the potential influences of somatic cells on germ cell fate decisions, we extracted and integrated the 4 types of somatic cells.A total of 8498 yak and 9380 cattle-yak Sertoli cells (SOX9) were identified (Figure 4 A).Interestingly, DEG analysis revealed that only 7 genes were downregulated (Figure 4 A, Supplementary Table 13;).These data indicated that Sertoli cell gene expression in yak and cattle-yak was similar and that Sertoli cells of cattle-yak could support normal spermatogenesis.Because the PTMCs and Leydig cells were grouped together, we further separated PTMCs and Leydig cells based on the expression of IGF1 and ACTA2 (Fig. S4), which share common progenitors during development and play crucial roles in spermatogenesis (Zhao et al., 2021, Tao et al., 2023).After removing the mixed Sertoli cells and endothelial cells, a total of 240 and 48 Leydig cells were identified in yak and cattle-yak, respectively (Figure 4 B).One hundred and 20 genes (60 upregulated and 60 downregulated) were differentially expressed in the Leydig cells of cattleyak (Figure 4 B, Supplementary Table 14).Similarly, among 935 and 388 PTMCs from yak and cattle-yak, respectively, only 34 genes were differentially expressed (Figure 4 C, Supplementary Table 15).
Recent studies have provided evidence that endothelial cells of the testis also play crucial roles in spermatogenesis (Bhang et al., 2018).Next, we analyzed testicular endothelial cells and immune cells using a similar approach.There were 610 and 178 endothelial cells in the yak and cattle-yak samples, respectively (Figure 4 D).Overall, 330 genes were downregulated and 11 were upregulated in the cells from cattle-yak (Figure 4 D, Supplementary Table 16).The downregulated genes were enriched in protein metabolic process, cell differentiation, response to chemical, and regulation of cell communication (Supplementary Table 17).Similarly, 858 and 164 immune cells were detected in yak and cattle-yak, respectively (Figure 4 E).DEG analysis revealed that 34 genes were upregulated and 531 were downregulated in the cattle-yak (Figure 4 Zhang et al.: Single-cell analysis identifies… E, Supplementary Table 18;).The terms enriched for downregulated genes identified by GO analysis included "cell communication," "signal transduction," "regulation of cellular metabolic process," and "response to chemical" (Supplementary Table 19).These data suggested that the cellular composition was dramatically different between yak and cattle-yak (Figure 4 F, Table 3); however, the gene expression at the single-cell resolution did not differ significantly among Sertoli cells, Leydig cells and PTMCs.

Intercellular communication network analysis identified abnormal signaling molecules in the testis of cattle-yak
To clarify the underlying intercellular communication that drives spermatogonial fate transitions in yak and to identify candidate factors that cause the abnormal development of spermatogenic cells in cattle-yak, we analyzed the number and strength of intercellular communication networks from the scRNA-seq data using CellChat (Figure 5 A and B).In yak, all spermatogonial subtypes communicate with each other, and sper- matogonia at late developmental stages communicate with spermatocytes.In contrast, all the spermatogonia failed to interact with spermatocytes, suggesting that communication between different spermatogonial subtypes differed among cattle-yak.The overall number and strength of intercellular communication between germ cell populations were lower in cattle-yak than in yak (Figure 5 C and D).The number of communications and strength decreased by 65.00% and 53.40%, respectively.Specifically, the ability of spermatogonial cell subpopulations to communicate with Diff-SPG1 was enhanced in the cattle-yak (Figure 5 D).
Further analysis indicated that Undiff-SPG1 played a weaker role in the communication network while Undiff-SPG2 acted more as a sender than as a receiver.Diff-SPG1 and Diff-SPG2 were the main receivers in the communication network (Figure 5 E).The SPCs had a high capacity for both receiving and transmitting signals in yak.Compared with those in yak, the roles of Undiff-SPG1 and Diff-SPG1 in cattle-yak in the communication network were similar, while the ability of Undiff-SPG2 and Diff-SPG2 to receive signals was reduced.SPCs had a reduced ability to both receive and send signals.Ligand-receptor interaction analysis revealed that the epidermal growth factor (EGF) signaling pathway was enriched in yak but not detectable in cattle-yak, and the fibroblast growth factor (FGF) signaling pathway was reduced in cattle-yak (Figure 5 F).EGF signals were emitted by Undiff-SPG1, Undiff-SPG2, Diff-SPG1 and Diff-SPG2 and were received only by Undiff-SPG1 (Figure 5 G).Specifically, the HBEGF-EGFR ligand receptor pair participates this communication networks (Figure 5 I).The FGF signaling pathways from Undiff-SPG2, Diff-SPG1 and Diff-SPG2 to undifferentiated spermatogonia were only present in yak (Figure 5 H).FGF10-FGFR1 was the main contributor to this signaling pathway (Figure 5 J).However, FGF10-FGFR1 was absent, and only FGF8-FGFR1 was present in cattle-yak cells.These data suggested that the germ cell subpopulation of the cattle yak had a reduced ability to receive and send signals and that the EGF and FGF signaling pathways associated with spermatogenesis were absent or dysregulated.
Signaling communication between somatic cells and germ cells in the testicular microenvironment dictates spermatogenic cell development.Next, we utilized a similar approach to examine the number and strength of interactions between germ cells and various types of testicular somatic cells (Figure 6 A and B).Compared with those of yak, the overall number and strength of intercellular communication between cell populations were greater in cattle-yak (Figure 6C and D).In particular, the number of communications from PTMCs and Leydig cells to germ cells and reversed interactions were increased.The greatest increase in the strength of cell communication was from Leydig cells to germ cells.Among the cell populations of yak testes, Sertoli cells had the strongest ability to send and receive signals (Figure 6 E).Germ cells were better at receiving signals than at sending signals, while endothelial cells were the opposite.Immune cells primarily act as receivers of signals, while the main function of Leydig cells is to send signals (Figure 6 E).Compared with those in yak communication networks, the greatest changes in the cattle-yak communication network were in germ cells, which have an increased ability to receive signals.
The pleiotrophin (PTN) signaling pathways that regulate spermatogenesis (Tian et al., 2022, Rajachandran et al., 2023) were enriched in cattle-yak (Figure 6 F).Consistent with the finding that Sertoli cells and Leydig cells can secrete PTN proteins (Tian et al., 2022, Rajachandran et al., 2023), ligands were also sent by these 2 cell types in yak and cattle-yak (Figure 6 G).Further comparisons revealed enhanced PTN signals from Leydig cells to germ cells in cattle yaks.The NCL and SDC receptors identified in previous studies (Tian et al., 2022, Rajachandran et al., 2023)were present in both yak and cattle-yak, with PTN-NCL playing a major role (Figure 6 H).In addition, NCL was upregulated in all the cattle-yak germ cells and Undiff-SPG1 (Supplementary Tables 10 and 12).Together, these data indicated that there were abnormalities in the communication between niche cells and germ cells in the cattle-yak group, with germ cells receiving more signals from Leydig cells.

Significant differences existed in spermatogenic cells from 12-mo-old and 24-mo-old yak and cattleyak
For a more comprehensive analysis of gene expression dynamics in spermatogenic cells in adult yak and cattle-yak, we obtained published data from scRNA-seq data from 24-mo-old yak and cattle-yak testis (Mipam et al., 2023).After quality control, we identified germ cells (DDX4), Sertoli cells (INHA), PTMCs, Leydig cells (ACTA2 & DCN), endothelial cells (CDH5), and immune cells (RGS1) according to the DEGs of each cluster (Fig. S5A, Supplementary Tables 20 and 21).We extracted 1652 and 388 germ cells from 24-mo-old yak and cattle-yak, respectively.After integration, we obtained a total of 3622 cells, which contained 15 and 13 clusters of cells in yak and cattle-yak, respectively (Fig. S5 B).We classified these cells into 9 types based on the DEGs of each cluster (Figure 7 A, Supplementary Table 23).The integrated spermatogonia included 4 subpopulations (Undiff-SPG1, Undiff-SPG2, Diff- SPG1, and Diff-SPG2).Moreover, we identified 5 types of spermatocytes and spermatids in yak: early primary spermatocytes (highly expressing YTHDC2 and SYCP2), late primary spermatocytes (highly expressing INSL6 and LYAR), diakinesis to secondary spermatocytes (highly expressing TMEM89 and TEKT), early spermatids (highly expressing CHD5 and TEKT3), and late spermatids (highly expressing TNP2 and SPEM1) (Fig. S5 C).In cattle-yak, there were only 4 types of spermatogonia and early primary spermatocytes, while late primary spermatocytes to late spermatids were missing.

DISCUSSION
Deciphering the genetic causes of hybrid sterility via transcriptome analysis is important for understanding the process of reproductive isolation.The male offspring of cattle and yak hybrids have been extensively studied, and a large amount of RNA-seq data have been generated using whole testicular tissues (Cai et al., 2017).However, the data generated from these studies ignored differences in specific cell types.We performed scRNA-seq analysis to detect differences in the subtypes of spermatogenic cells and analyzed molecular interactions between germ cells and somatic cells in yaks and cattle-yaks.The information gained from these analyses not only provides a valuable data set for mining functional genes regulating spermatogenesis in yak but also reveals the developmental trajectory and impairments of spermatogonial fate decisions in cattleyak.
A key finding of the present study is that the composition and developmental trajectory of spermatogonial lineages are different between yak and cattle-yak.The germ cells of 12-mo-old yak and cattle-yak were divided into 5 different subtypes, including 2 types of undifferentiated spermatogonia (Undiff-SPG1&2), 2 types of differentiated spermatogonia (Diff-SPG1&2) and sper-     To sec.SPC, Early SPT and Late SPT denote undifferentiated spermatogonia 1, undifferentiated spermatogonia 2, differentiated spermatogonia 1, differentiated spermatogonia 2, early primary spermatocytes, latr primary spermatocytes, diakinesis to secondary spermatocytes, early spermatids and late spermatids, respectively.(B) Volcano plot showing DEGs between 12-mo-old and 24-mo-old yak and cattle-yak.Red denotes upregulated genes, blue denotes downregulated genes, and gray denotes genes for which |avg_log2FC| < 0.2.(C) GO analysis of upregulated and downregulated genes in germ cells between 12-mo-old and 24-mo-old yak and cattle-yak.Red represents upregulated transcripts, while blue represents downregulated transcripts.The dot size indicates the number of expressed cells.The vertical axis represents the terms, and the horizontal axis represents the transformed FDR (-log10P value).(D) The t-SNE plot of integrated germ cells from yak and cattle-yak.(E) Sankey dots show a subset of DEGs and their corresponding functions in 5 germ cell populations.Undiff-SPG1, Undiff-SPG2, Diff-SPG1, Diff-SPG2 and Early Pri.SPC denotes undifferentiated spermatogonia 1, undifferentiated spermatogonia 2, differentiated spermatogonia 1, differentiated spermatogonia 2 and early primary spermatocytes, respectively matocytes (SPCs).In yaks with normal spermatogenesis, Undiff-SPG1 represented the most primitive cell state that likely contained enriched SSCs.This population was reduced by more than 90% in cattle-yak.A reduction in the undifferentiated spermatogonial pool was detected in cattle-yak, and it was proposed that the gonocyte to spermatogonial transition was impaired, resulting in defects in the proliferation and expansion of the foundational undifferentiated spermatogonia (Wu et al., 2023).Another notable change was the increase in the ratio of differentiated spermatogonia.The relative proportion of differentiated spermatogonia was increased by 2-fold in cattle-yak, indicating defects in spermatogonial fate transitions, likely from premature formation of progenitor spermatogonia from the undifferentiated compartment or a block in spermatogonial differentiation.It has been reported that the proliferation and differentiation of undifferentiated spermatogonia are impaired in cattle-yak (Mipam et al., 2023).Spermatogenic cells isolated from yak responded to RA treatment with enhanced expression of genes associated with differentiation status (Mipam et al., 2023).In contrast, RA treatment of spermatogonia from cattle-yak failed to modulate the expression of these genes, further illustrating impairments in the fate transition of the spermatogonial population.RA-STRA8 signaling induces spermatogonial differentiation (Bowles et al., 2006).The expression of STRA8 begins in late undifferentiated spermatogonia, persists in differentiated spermatogonia, and decreases greatly or terminates after cells enter meiosis (Schug et al., 2007).RA concentrations in 10-mo-old yak and cattleyak testes did not differ (Wu et al., 2023); however, STRA8 expression decreased in cattle-yak germ cells.Therefore, we hypothesized that the ability of undifferentiated and differentiated spermatogonia to receive and respond to RA signals is impaired in cattle-yak.Dissecting the causes of spermatogonial differentiation will be the focus of future studies.
After integrating the germ cells of yak and cattleyak at 24 mo of age, we found that the composition of yak germ cells at different ages varied greatly, with spermatocytes and spermatozoa accounting for 96.34%.The proportion of spermatogonia in 24-mo-old yaks was 3.73%, which was much lower than that in 12-mo-old yak.However, we found that the cattle-yak at different ages were similar in composition, and all contained 4 spermatogonial subpopulations and early primary spermatocytes.We compared the gene expression levels of 5 germ cell subpopulations of integrated yak and cattleyak and detected abnormal expression of genes related to cell proliferation and apoptosis, spermatogonial cell differentiation, and meiotic entry.These data indicated that defects in spermatogonial fate decisions emerged in 12 mo-old cattle-yak.
In addition to detecting potential factors that cause defects in spermatogenic cell development, this study investigated the transcriptome landscape of 6 testicular cell types and identified genes associated with undifferentiated spermatogonial fate.Phf13 is expressed in PLZF-positive undifferentiated spermatogonia (Bördlein et al., 2011), and TSPAN33 is highly enriched in primitive undifferentiated spermatogonia (Sohni et al., 2019).LFNG-mediated Notch signaling regulates neural stem cell quiescence, division and fate determination (Semerci et al., 2017).The neurotransmitter receptor GABBR1 is a key factor regulating the proliferation of hematopoietic stem/progenitor cells (Shao et al., 2021).CNN3 can regulate the function of neural stem cells, and CNN3 deficiency increases the number and size of newly formed neutrospheres without interfering with their differentiation potential (Junghans and Herzog, 2018).A key regulator of the cell cycle, FBXW7, plays a role in maintaining somatic stem cells, embryonic stem cells and induced pluripotent stem cells (Takeishi and Nakayama, 2014).These genes may be potential SSC markers for yak and other closely related ruminant species.
Cell-to-cell communication is essential for the maintenance of homeostasis in individual cells and multicellular organs (Qiao et al., 2020).Compared with yak, the overall number and strength of cell communication between germ cell populations in cattle-yak were reduced, implying a weaker ability to signal within the germ cell subpopulations.As a result, germ cell subpopulations are delayed or unable to receive and send signals to transform into the next state, resulting in a slowed or even stalled spermatogenesis process.Niche-germline communication affects self-renewal regulation, proliferation rates, metabolism and transitions between the differentiation states of germ cells (Kanatsu-Shinohara andShinohara, 2013, Guo et al., 2018).We inferred  that abnormal communication between somatic cells and germ cells was also one of the causes of male sterility in cattle-yak.EGF promotes the self-renewal and proliferation of stem cells and regulates germ cell development (Niederberger et al., 1993, Yu et al., 2019).The FGF signaling pathway is involved in spermatogenesis (Cotton et al., 2006), and FGF is associated with the maintenance of SSCs (Kubota et al., 2004).Enhanced niche-germline communication in cattle-yak was mainly reflected by the enrichment of the MK and PTN signaling pathways.The PTN signaling pathway, which regulates spermatogonial function, has been identified in both sheep and human testis scRNA-seq data (Tian et al., 2022) and mouse testis spatial transcriptomics (Rajachandran et al., 2023).The absence and dysregulation of these signals may imply abnormalities in the regulation of SSC population self-renewal and subsequent spermatogenesis in cattle-yak.

CONCLUSIONS
Overall, this study revealed the transcriptome landscape of testicular cells of yaks and cattle-yaks at the prepubertal stage.We described the developmental trajectories and cellular compositions of spermatogenic cells, including 4 spermatogonial subpopulations and spermatocytes.The spermatogenic cell composition of cattle-yak was different from that of yak, with a decrease in the proportion of undifferentiated spermatogonia and an increase in the proportion of differentiated spermatogonia.By comparing gene expression levels in different germ cell subpopulations of yak and cattleyak, we identified genes associated with cell proliferation and apoptosis, spermatogonial cell differentiation, and meiosis entry.These analyses demonstrated that spermatogonial differentiation and meiotic entry were abnormal in cattle-yaks, mainly due to dysregulated gene expression and cellular communication, resulting in slow or even stagnant spermatogenesis.Moreover, the data reported in this study provide a reference for screening SSC markers and valuable data for understanding abnormally expressed genes associated with spermatogonial cell self-renewal, differentiation, and meiotic entry.
Zhang et al.:  Single-cell analysis identifies… BD Rhapsody Single-Cell Analysis System for singlecell library construction and sequencing.

Figure 1 .
Figure 1.Single-cell transcriptome profiling and cluster identification of 12-mo-old yak and cattle-yak testicular cells.(A) H&E staining of 12-mo-old yak and cattle-yak testes.Red arrows indicate spermatogonia, black arrows indicate spermatocytes, and blue arrows indicate Sertoli cells.(B) t-SNE plots showing the identified cell types (cells are colored according to the 5 cell types) of 12-mo-old yaks and cattle-yaks.(C) The expression patterns of the known marker genes projected on the t-SNE plot.The colors from blue to gray represent the expression levels from high to low, respectively.(D) Heatmap of the marker genes for each cell type.Different colors at the top of the heatmap represent different cell types.The right of the heatmap represents the DEGs between each cell type.
Figure 2. Single-cell transcriptome profiling and cluster identification of 12-mo-old yak and cattle-yak germ cells.(A) T-SNE plot of the 2 groups of germ cells after integration.(B) Pseudotime-ordered analysis of germ cells.The left panel is colored according to the pseudotime order.The middle panel shows cells from yak.The right panel shows cells from cattle-yak.(C) Clusters of genes that were differentially expressed across pseudotimes from integrated germ cell data sets are shown as heatmaps according to the expression color code noted at the bottom.Key genes and GO analysis results are noted on the right (see Supplementary Tables6 and 7).(D) Expression patterns of key genes that regulate spermatogonial development among integrated germ cells according to pseudotime (cell coloring id according to the state).(E) The stacked area shows the proportion of germ cells in each cell type in 12-mo-old yaks and cattle-yaks.Undiff-SPG1, Undiff-SPG2, Diff-SPG1, Diff-SPG2 and SPCs denote undifferentiated spermatogonia 1, undifferentiated spermatogonia 2, differentiated spermatogonia 1, differentiated spermatogonia 2 and spermatocytes, respectively.
Figure 3. Visualization of DEGs and functional enrichment analysis of germ cells from 12-mo-old yaks and cattle-yaks.(A) Volcano plot showing DEGs in germ cells of 12-mo-old yaks and cattle-yaks.Red denotes upregulated genes, blue denotes downregulated genes, and gray denotes genes for which |avg_log2FC| < 0.2.(B) GO analysis of upregulated and downregulated genes in germ cells between 12-moold yaks and cattle-yaks.Red represents upregulated transcripts, while blue represents downregulated transcripts.The dot size indicates the number of expressed cells.The vertical axis represents the terms, and the horizontal axis represents the transformed FDR (-log10P value).(C) The number of DEGs in 5 germ cell subtypes of 12-mo-old yaks and cattle-yaks.Red represents upregulated transcripts, while blue represents downregulated transcripts.(D) Sankey dots show a subset of DEGs and their corresponding functions in 5 germ cell populations.Undiff-SPG1, Undiff-SPG2, Diff-SPG1, Diff-SPG2 and SPCs denote undifferentiated spermatogonia 1, undifferentiated spermatogonia 2, differentiated spermatogonia 1, differentiated spermatogonia 2 and spermatocytes, respectively.
Figure 4. Single-cell transcriptome profiling and cluster identification of 12-mo-old yak and cattle-yak somatic cells.(A-E) T-SNE plots of integrated 12-mo-old yak and cattle-yak Sertoli cells, Leydig cells, PTMCs, endothelial cells and immune cells (top).Volcano plots showing DEGs in Sertoli cells, Leydig cells, PTMCs, endothelial cells and immune cells in 12-mo-old yaks and cattle-yaks.Red denotes upregulated genes, blue denotes downregulated genes, and gray denotes genes for which |log2(fold change)| < 0.2.(F) Donut chart showing the percentage of each cell type in yak and cattle-yak testes.The table shows specific percentages.
Zhang et al.: Single-cell analysis identifies…

Figure 5 .
Figure 5. Communication between germ cell populations of 12-mo-old yaks and cattle-yaks.(A-B) Circos plots showing the communication between 5 germ cell populations of 12-mo-old yaks and cattle-yaks.The arrows and the colors of the lines indicate the direction.The loops indicate autocrine circuits.Edge thickness indicates the number (A) or strength (B) of interactions between cell populations.(C) The bar graph shows the overall number of interactions between 12-mo-old yaks and cattle-yaks (left).The circos plot shows the difference in the number of interactions between 12-mo-old yaks and cattle-yaks (right).The blue line indicates a decreased number of communications in the cattle-yak compared with that in the yak, while the red line indicates an increased number of communications in the cattle-yak.(D) The bar graph shows the overall strength of interactions between 12-mo-old yaks and cattle-yaks (left).The Circos plot shows the difference in the strength of interactions between 12-mo-old yaks and cattle-yaks (right).The blue line indicates a decreased strength of communications in the cattle-yak compared with that in the yak, while the red line indicates an increased strength of communications in the cattle-yak.(E) The incoming and outgoing communication capacities of each cell population in 12-mo-old yak (left) and cattle-yak (right).(F) Overall information flow of each signaling pathway in 12-mo-old yak and cattle-yak.The red signaling pathways were more enriched in 12-mo-old yak, and the black signaling pathway was equally enriched in 12-mo-old yak and cattle-yak.(G) Heatmap of outgoing and incoming signals from EGF signaling pathways in 12-moold yak.The vertical coordinate is the sender, and the horizontal coordinate is the receiver.(H) Heatmap of outgoing and incoming signals from FGF signaling pathways in 12-mo-old yak (left) and cattle-yak (right).The vertical coordinate is the sender, and the horizontal coordinate is the receiver.(I) Relative contribution of each ligand-receptor to the overall communication network of the EGF signaling pathway in 12-mo-old yak.(J) Relative contribution of each ligand-receptor to the overall communication network of the FGF signaling pathway in 12-mo-old yak (left) and cattle-yak (right).Undiff-SPG1, Undiff-SPG2, Diff-SPG1, Diff-SPG2 and SPCs denote undifferentiated spermatogonia 1, undifferentiated spermatogonia 2, differentiated spermatogonia 1, differentiated spermatogonia 2 and spermatocytes, respectively.

Figure 6 .
Figure 6.Communication between cell populations in the testes of 12-mo-old yak and cattle-yak.(A-B) Circos plots showing the communication between 6 cell populations in 12-mo-old yak and cattle-yak testes.The arrows and the colors of the lines indicate the direction.The loops indicate self-secreting circuits.Edge thickness indicates the number (A) or strength (B) of interactions between cell populations.(C) The bar graph shows the overall number of interactions between 12-mo-old yak and cattle-yak (left).The circos plot shows the difference in the number of interactions between 12-mo-old yak and cattle-yak (right).The blue line indicates a decreased number of communications in the cattle-yak compared with that in the yak, while the red line indicates an increased number of communications in the cattle-yak.(D) The bar graph shows the overall strength of interactions between 12-mo-old yak and cattle-yak (left).The Circos plot shows the difference in the strength of interactions between 12-mo-old yak and cattle-yak (right).The blue line indicates a decreased strength of communications in the cattle-yak compared with that in the yak, while the red line indicates an increased strength of communications in the cattle-yak.(E) The incoming and outgoing communication capacities of each cell type in 12-mo-old yak (left) and cattle-yak (right).(F) Overall information flow of each signaling pathway in 12-mo-old yak and cattle-yak.The black signaling pathways are equally enriched in 12-mo-old yak and cattle-yak, and the green signaling pathways are more enriched in 12-mo-old cattle-yak.(G) Heatmap of outgoing and incoming signals from PTN signaling pathways in 12-mo-old yaks (left) and cattle-eyads (right).The vertical coordinate is the sender, and the horizontal coordinate is the receiver.(H) Relative contribution of each ligand-receptor to the overall communication network of the PTN signaling pathway in 12-mo-old yaks.

Table 1 .
Zhang et al.:Single-cell analysis identifies… The number and percentage of cells of each cell type in 12-mo-old yak and cattle-yak

Table 2 .
The number and percentage of cells in each germ cell type of combined yak and cattle-yak

Table 3 .
The number of cells in each cell type of 12-mo-old yak and cattle-yak after purification

Table 4 .
The number and percentage of cells of each germ cell type in 24-mo-old yaks and cattle yaksCell typeUndiff-SPG1 Undiff-SPG2 Diff-SPG1 Diff-SPG2 Early pri.SPC Late pri.SPC Diak.To sec.SPC Early SPT Late SPT