Dissecting the dynamic cellular transcriptional atlas of adult teleost testis development throughout the annual reproductive cycle

In vertebrates, the successful proliferation and differentiation of testicular germ cells are pivotal for species preservation and survival (Shim and Anderson, 1998). The completion of the differentiation process, particularly the transition of the spermatogonia into primary spermatocyte, secondary spermatocyte, spermatid and ultimately spermatozoa, stands as a hallmark for normal spermatogenesis and the continuity of species (Nishimura and L'Hernault, 2017). This process is intricately linked to the coordination and regulation of various cell types residing within the testis microenvironment, including stromal cells, Sertoli cells, Leydig cells, macrophages and endothelial cells (Bhang et al., 2018; Crespo et al., 2016; DeFalco et al., 2015; Zhou et al., 2019). Fish, being one of the most diverse and populous vertebrates, display diverse reproductive strategies in testis development. Certain fish species, such as zebrafish and medaka (Yang and Tiersch, 2009), following sexual maturation, exhibit a continuous renewal of germ cell populations and uninterrupted production of mature spermatozoa throughout the year, resembling the pattern observed in mammals. In contrast, seasonal breeding fishes such as Nile tilapia (Kobayashi and Nagahama, 2009), Chinese tongue sole (Zhu et al., 2018) and Japanese flounder (Yang et al., 2018) maintain spermatogenic cells throughout the annual reproductive cycle in adult testes, with varying spermatogenesis activity between spawning and non-spawning seasons. Moreover, distinct morphological characteristics and cellular compositions are discernible at different development stages during the annual reproductive cycle in several season breeding fishes, such as rainbow trout (Sato et al., 2017), pike (Hoffmann et al., 1980; Idowu, 2017) and European sea bass (Valero et al., 2015). After the spawning season, the testes of these species enter the regressed (Rs) stage, wherein diploid spermatogonia become the sole germ cells capable of mitotic division, priming for the subsequent initiation of spermatogenesis and the restoration of spermatozoa production. This stringent spatiotemporal orchestration of spermatogenesis provides a framework for analyzing the physiological properties and molecular characteristics of the testis at specific spermatogenic stages, encompassing spermatogonia maintenance, proliferation, differentiation, meiosis onset and spermatogenesis.

Considerable efforts have been directed towards elucidating the developmental program of fish testes throughout the annual reproductive cycle, yielding comprehensive insights into both morphological and molecular aspects across diverse species. Despite the wealth of information derived from these investigations, our understanding of fish testicular development during the annual reproductive cycle remains limited, because these studies were confined to morphological observation, organ-level development mechanisms or specific subsets of testicular cell populations. The molecular regulation governing dynamics across diverse cell types during testis development necessitates further elucidation (Bahamonde et al., 2016; Wang et al., 2018; Dias et al., 2020; Eildermann et al., 2012; Muller et al., 2008; Ozaki et al., 2011; Yang et al., 2016). The emergence of single-cell sequencing (scRNA-seq), a pioneering technology renowned for its high throughput and precision, has opened new avenues in scientific inquiry. This transformative approach offers comprehensive biological profiles and unprecedented insights previously unattainable across various research domains, including organismal development (Samad and Wu, 2021), embryogenesis (Liu et al., 2022) and spermatogenesis (Guo et al., 2018; Shami et al., 2020). Within the realm of testicular development, scRNA-seq has unveiled the complete roadmap of testicular cells and corresponding transcriptional signatures. Additionally, unexpected cell types during testis development have been identified, and the transcriptional disparities at single-cell level across species have also been defined (Green et al., 2018). Despite the extensive application of scRNA-seq in investigating spermatogenesis in mammals (Fang et al., 2022; Deng et al., 2022), there remains a paucity of studies focusing on spermatogenesis in non-mammals, particularly in teleost (Wang et al., 2023a; Wu et al., 2021).

The black rockfish (Sebastes schlegelii), a viviparous marine teleost species, undergoes a long-term spermatozoa storage of over 4 months in the ovary. Copulation occurring in November and December, leads to fertilization taking place in the following April (He et al., 2019). In adult male black rockfish, the characteristics of testicular cell composition and quantity dynamically shift throughout the annual reproductive cycle, providing a valuable opportunity to investigate the dynamic regulation of fish spermatogenesis (Mori et al., 2003). Moreover, the genome of black rockfish has been assembled in high quality and published in previous works of our lab, which provides the efficient tools for gene annotation (He et al., 2019). In this study, scRNA-seq was employed to analyze testicular cells obtained from the Rs stage, the regenerating (Rg) stage and differentiating (Dif) stage. The comprehensive cellular atlas was constructed from a substantial cell cohort, offering insights into the proliferation and differentiation activities of spermatogonia subsets. Notably, our findings unveiled dynamic regulations and bioenergetic transitions during spermatogonia development. Re-clustering analysis further revealed a deeper level of heterogeneity and molecular events driving spermatogenesis. Additionally, the developmental trajectory of Sertoli cells and the heterogeneity of Leydig cells and macrophages during the annual reproductive cycle were investigated. Finally, an interaction network between somatic cells in the testicular micro-environment and germ cells was established. These findings contribute to a comprehensive understanding of spermatogenesis during the annual cycle development of teleost testes and represents a community resource for understanding testicular cell development in vertebrates.

Testes from adult black rockfish were collected in February (Rs stage), May (Rg stage) and July (Dif stage) (Fig. 1A) and subjected to histological analysis to evaluate changes in morphological and cellular composition (Fig. 1B; Fig. S1A). During the Rs stage, the testes exhibited sparse spermatogonia, with noticeable interstitial cell presence between seminal vesicles. Transitioning to the Rg stage, there was an increase in spermatogonia, encompassing both undifferentiated and differentiated forms, while spermatocytes remained absent. Progressing to the Dif stage, the testes displayed abundant spermatocytes, a considerable number of spermatogonia in the seminal vesicles and the occurrence of sperm. However, the area between seminal vesicles housing interstitial cells became visibly narrower with testis development (Fig. 1B).

To unravel the intricate variations in cell types and their gene expression during the entire annual developmental cycle of the testis, scRNA-seq for testes at three pivotal developmental stages was employed on a 10x Genomics platform and the comprehensive details of each dataset were summarized in Fig. S1B. Nine distinct cell types (spermatogonia, spermatocytes, spermatids, macrophages, Sertoli cells, Leydig cells, endothelial cells, stromal cells and blood cells) were identified (Fig. 1C, Fig. S1C) based on the expression patterns of representative markers (Fig. S1D,E; Table S1). The correlation analysis across each cell subset revealed a high correlation between germ cells and somatic cells (Fig. S1E).

Subsequently, temporal changes in the abundance of specific cell types were delineated (Fig. 1D), which was consistent with the histological observations (Fig. 1B), especially regarding spermatogonia, spermatocytes and spermatid. Genes enriched in each cell type were determined using the ‘Find Markers’ functionality and summarized in Table S2. The expression heatmap of maker genes vividly depicted the unique gene expression profiles characterizing each cell type, with the top 10 representative marker genes outlined in Fig. 1E and Table S2. Some of these marker genes, such as Centromere protein F (cenpf) (Handel et al., 1999), Sperm Adhesion Molecule 1 (spam1) (Chen et al., 2006) and Cytochrome P450 17A1 (cyp17a1) (Wang et al., 2019), have been previously validated for expression in respective cell type in mammals or teleost. However, the involvement of a considerable number of genes in testis development has not been validated. Lastly, Gene Ontology (GO) analysis unveiled enriched terms from each cell type, aligning with their biological characters (Guo et al., 2018; Zhang et al., 2023; Fritz and Gommerman, 2011) (Fig. 1F). This alignment underscores the accuracy and reliability of the cell type identification.

In this investigation, a series of marker genes for germ cells was identified. The validation of some germ cell marker genes was performed through in situ hybridization (ISH), including several potentially novel marker genes, such as Tudor domain containing 15 (tdrd15) and Replication protein A 14 kDa subunit (rfa3) for spermatogonia, Adenylate kinase 9 (ak9) for spermatocytes, and High choriolytic enzyme 2 (hce2) for spermatids, which belongs to the astacin metalloproteinase subfamilies. These were expanded in rockfish and predicted to be associated with the interaction between sperm and zona pellucida in black rockfish (He et al., 2019). Additionally known marker genes including Synaptonemal complex protein 1 (sycp1) for spermatocytes, which encodes a protein that localizes to the central element of the synaptonemal complex during meiosis (Nabi et al., 2022), and Centriole-associated 1 like (spatc1l) for spermatid, which is associated with mammalian fertility (Kim et al., 2018) (Fig. 2A-F) were also confirmed to be expressed in the corresponding cell types in black rockfish testis. Subsequently, re-clustering analysis of germ cells revealed six distinct germ cell subgroups: undifferentiated spermatogonia (Un-Diff), type I differentiated spermatogonia (Diff1), type II differentiated spermatogonia (Diff2), type I spermatocytes (Sc-1), type II spermatocytes (Sc-2) and spermatid (St) (Fig. 2G; Fig. S2A). Combining the datasets of black rockfish and zebrafish, germ cells from the two species clustered together (Fig. S2B). Notably, spermatogonia exhibited high correlation among different species, indicating the conservation of spermatogonia in teleost. However, different subsets of spermatogonia showed higher correlations within species rather than across species, as observed with spermatocytes and spermatid (Fig. S2C), suggesting the divergence of molecular characteristics in specific germ cell subsets among different teleost species. Throughout the testis development in black rockfish, the number of Un-Diff increased from the Rs stage to the Rg stage, but sharply decreased at the Dif stage. Conversely, a considerable number of Diff spermatogonia were identified at the Dif stage, indicating the importance of Diff in the spermatogenesis process during the Dif stage (Fig. S2D). The known markers for each cell type were presented in Fig. 2H and Table S1. Numerous markers associated with germ cell development in mammals have been demonstrated. For example, POU class 3 homeobox1 (Pou3f1), DND microRNA-mediated repression inhibitor 1 (Dnd1), and MYC proto-oncogene (c-myc; also known as Myc) have been reported as crucial for spermatogonia maintenance in mammals (Wu et al., 2010; Niimi et al., 2019; Koji et al., 1988). Indeed, a series of genes were found to be specifically expressed in spermatogonia subsets of both teleosts and humans based on scRNA-seq datasets, including DND1, piwi like RNA-mediated gene silencing 1 (PIWIL1) and c-myc (Table S3). Notably, numerous spermatogonia markers identified in mammals were not found to show obvious expressions in black rockfish, such as Gfra1 (Grasso et al., 2012), Sall4 (Lovelace et al., 2016) and Eomes (Sharma et al., 2019). Among them, gfra1 has been identified to be specifically expressed in spermatotogonia stem cells in some teleost species, such as rainbow trout (Bellaïche et al., 2014). This suggested a divergence in molecular characteristics of spermatogonia between different species. Of course, the sequencing depth and potential dropout in single-cell analysis could contribute to the low expressions or absence of certain genes in black rockfish, which remains to be resolved with additional datasets or validation experiments. Additionally, the dynamical molecular changes among different cell subsets were illustrated and the top ten markers were listed (Fig. S2E; Table S4). In mammals, the ribosome has been reported to guide piRNA formation (Sun et al., 2021), which is crucial for germ cell development. In black rockfish, the ribosome proteins, such as ribosome protein 5 (rpl5a) and ribosome protein 7 (rpl7a), were enriched in Diff spermatogonia. In addition, piRNA related genes, such as piwil1 and tdrd15, were also enriched in spermatogonia. These suggest that ribosome proteins might participate in spermatogonia development by regulating ribosome-guided piRNA pathways.

Functional enrichment analysis of marker genes in each cell type was performed to investigate their distinct functions. Representative GO terms and KEGG pathways enriched in each cell type were illustrated (Fig. S2F,G). Comparative analysis between Un-Diff and Diff spermatogonia revealed that biological processes and pathways associated with cell proliferation and cell cycle were significantly enriched in Diff, suggesting a proliferation dynamic from Un-Diff to Dif.

To gain a comprehensive insight into spermatogenesis in black rockfish, an unbiased dynamic cell trajectory analysis was conducted, ordering cells from different germ cell clusters based on pseudo-time. The analysis revealed a developmental continuity without major branching points. The precedence relationship identified based on gene expression patterns was consistent with Monocle analysis, positioning Un-Diff and spermatid at two terminals (Fig. 2I). Additionally, the distributions of germ cells from the Rs stage, Rg stage and Dif stage along pseudo-time trajectory were consistent with their developmental time points (Fig. 2J). Notably, the Diff2 from the Dif stage were positioned to connect the preceding and following stages on the trajectory. Moreover, the expression levels of representative genes specific to different cell types exhibited dynamic changes along pseudo-time trajectory (Fig. 2K).

To delve deeper into molecular events involved in germ cell development, genes enriched in cells along pseudo-time were clustered into four distinct gene expression patterns (Fig. 2L). GO enrichment analysis was performed to investigate the detailed functions of genes from different clusters. At the early stage (Group 1), germ cell-expressed genes were significantly enriched in GO categories such as ‘gene expression’ and ‘histone modification’, whereas genes in Group 2 and Group 3 were associated with GO categories such as ‘cellular metabolic process’, ‘chromosome organization’, ‘microtubule-based process’ and ‘cell cycle’. Group 4 genes were mainly expressed in spermatid, with enriched GO categories including ‘cell redox homeostasis’ and ‘glyceropholipid biosynthetic process’. Some of these molecular events, such as histone modification, have been identified to participate in germline development in mammals and be crucial for spermatogonia development and meiosis initiation, as well as being associated with cAMP pathways (Liu et al., 2010). Furthermore, statistically enriched signaling pathways corresponding to the four gene clusters were identified (Fig. 2L). Among them, PI3K-AKT pathways and the calcium signaling pathway were associated with spermatogonia proliferation and differentiation (Zhao et al., 2019a; Golpour et al., 2016). In addition, for in vitro spermatogonia culture, various growth factors, such as FGF2 (Zhang et al., 2012), LIF and GDNF (Wang et al., 2014), have been used. Here, the VEGF signaling pathway, which has been reported to be involved in spermatogonia maintenance in mammals (Sargent et al., 2016), was found to be significantly enriched in Un_Diff spermatogonia of black rockfish. So it was speculated that VEGF-related factors might be vital for stemness maintenance of fish spermatogonia.

The proliferation and differentiation of spermatogonia are pivotal in driving the spermatogenesis process, particularly amid in the dynamic changes observed during the annual cycle of teleost testis development. In rainbow trout, spermatogonia from testes at different stages exhibited a significant difference in proliferation rate (Sato et al., 2017). In black rockfish, the number of distinct spermatogonia subsets was found to significantly differ with the stage of testis development (Fig. S2B), which was further validated by the dynamical expressions of VASA protein (Fig. 3A). Subsequently, the expression analysis of spermatogonia proliferation genes in different spermatogonia subsets revealed higher expression in Diff compared with Un-Diff spermatogonia (Fig. 3B). Specifically, these genes were most abundant in Diff2 at Rg- and Dif-stage testes (Fig. 3C), indicating a proliferation expansion in Diff along with testis development. Furthermore, the higher expression levels of several genes related to cell proliferation, including Histidine triad nucleotide-binding protein 1 (hint1) (Zhang et al., 2021b), PCNA-associated factor (paf15; pclaf) (Karg et al., 2017), Marker of proliferation Ki-67 (mki67) (Zhao et al., 2018) and Chromatin assembly factor 1 subunit A (chaf1a) (Hoek and Stillman, 2003), and Proliferating cell nuclear antigen (PCNA) (Takagi et al., 2001), were validated by ISH and immunohistochemistry (IHC) (Fig. 3D,E).

Energy produced by ATP metabolism is crucial for proliferation of various cell types. In teleost, exogenous ATP has been identified to promote the proliferation of spermatogonia in rainbow trout, highlighting its importance for spermatogonia development in teleosts (Loir, 1999). Differentially expressed genes (DEGs) in black rockfish Diff2 at different stages of testis development were identified to be significantly involved in ATP synthesis and ATPase activity (Fig. S3A,B). Indeed, the genes involved these processes exhibited prominent expressions in Diff2 at the Rg and Rs stages (Fig. 3F). Additionally, the activity of ATPase in the isolated spermatogonia from testes gradually increased from the Rs stage to the Dif stage (Fig. 3G).

Concerning spermatogonia differentiation, meiosis-related genes, such as sycp1, Synaptonemal complex protein 2 (sycp2) and Meiotic recombination protein REC8 homolog (rec8a) (Watanabe and Nurse, 1999), were highly expressed in Diff2, particularly at the Dif stage (Fig. 3H,I). This indicated that Diff2 spermatogonia at the Dif stage were crucial for initiating meiosis, along with their position on the spermatogenesis trajectory (Fig. 2J). Indeed, ISH signals of cenpf and sycp2 were most prominent at the Dif stage, as was the expression of Disruption meiotic cDNA 1 (DMC1) (Kagawa and Kurumizaka, 2010) protein (Fig. 3J). Combining these findings with the highest expression of cell cycle genes in Diff spermatogonia at Dif stage (Fig. S3C), it could be speculated that the proliferation and differentiation dynamics of spermatogonia during testis development were influenced by the divergence in energy produced by ATP metabolism, resulting in changes in the mitosis and meiosis cycle.

The number of Un-Diff cells increased explosively from the Rs stage to the Rg stage, followed by a sharp decrease at the Dif stage, indicating that Un-Diff proliferation primarily contributed the abundance of germ cells at the Rg stage. To explore the heterogeneity of Un-Diff cells in detail, the expressions of Acetylcholine receptor subunit delta (achd) and pou3f1, two marker genes in Un-Diff (Fig. 2H; Fig. S2C), were examined. Interestingly, pou3f1 signals were less distributed in spermatogonia at the Rg stage compared with achd signals (Fig. 4A), suggesting the presence of distinct Un-Diff subsets. Re-clustering analysis of spermatogonia further separated Un-Diff into two cell clusters, named Un-Diff1 and Un-Diff2 (Fig. 4B). The gene expression patterns were distinctly different between Un-Diff1 and Un-Diff2 (Fig. S4A). The representative maker genes for two Un-Diff subsets, such as c-myc and Forkhead transcription factor (foxm1) (Un-Diff1), and Bone morphogenetic protein 2 (bmp2b) and START domain-containing protein 10 (stard10) (Un-Diff2), were displayed in Fig. S4B. The increase in cell number was higher in Un-Diff2 than in Un-Diff1 from the Rs stage to the Rg stage (Fig. S4C). Meanwhile, the spermatogonia proliferation genes also exhibited higher expression in Un_Diff2 (Fig. 4C), indicating a higher proliferation activity in Un-Diff2. Through DEG analysis, a total of 161 and 509 DEGs were identified in Un_Diff1 and Un_Diff2, respectively (Fig. S4D). DEGs upregulated in Un-Diff2 were mainly involved in ATP biosynthesis and oxidative phosphorylation, which was related to the TCA cycle (Fig. 4D). The expressions of representative genes are listed in Fig. 4E. DEGs upregulated in Un-Diff1 were associated with glutamine metabolism [Gamma-glutamyl hydrolase (ggh) (Chen et al., 2022) and glutathione S-transferase 1 (gst1) (Kanwar et al., 2020)] and lipoic acid metabolism [lipoyltransferase 2 (lip2) (Spalding and Prigge, 2010)], which were related to the glycolysis process (Fig. 4F,G). These findings suggested a potential variation in the regulation of the bioenergetic source from Un-Diff1 to Un-Diff2, which might involve a shift from glycolysis to the TCA cycle.

Weighted correlation network analysis (WGCNA) was performed to identify key genes and networks linked to two Un-Diff clusters. A total of ten modules were obtained based on all expressed genes. Among them, the brown module and turquoise module were significantly associated with Un-Diff1 and Un-Diff2, respectively (Fig. S4E). Co-expression networks for Un-Diff1 (pax3b and notch2) and Un-Diff2 were constructed with the top 5% weight correlation, and the top12 hub genes were identified. For example, phosphoglucomutase-1 (pgm1), a key enzyme for glycolysis (Cao et al., 2021), was identified as a hub gene in Un-Diff1. Paired box protein 3 (pax3b), Neurogenic locus notch homolog protein 2 (notch2) and Nanos C2HC-type zinc finger 3 (nanos3), the stem cell related genes (de Morree et al., 2019; Fujimaki et al., 2018; Beer and Draper, 2013), were identified as hub genes in Un-Diff1 (pax3b and notch2) and Un-Diff2 (nanos3) (Fig. 4H). Furthermore, the developmental trajectory of spermatogonia was constructed (Fig. S5A) and gene clustering analysis was performed (Fig. S5B), showing differential involvement of ATP biosynthesis and metabolism in spermatogonia development (Fig. S5C). Specific key transcription factors (TFs) were identified in different spermatogonia clusters, exerting distinct regulatory effects in each cell subset (Fig. S6A,B). The c-myc gene, a key regulator of glycolysis (Goetzman and Prochownik, 2018), was identified as a key regulator in Un-Diff1 according to the TF networks in two Un-Diff clusters analyzed by WGCNA (Fig. S6C), indicating the crucial roles of glycolysis in Un_Diff1.

Recently, the research in mouse has demonstrated a bioenergetic transition during spermatogonia development, characterized by an upregulation of mitochondria respiration-related genes and a downregulation of glycolysis-related genes. In our study, it was also found that the expression of genes related to glycolysis, such as Hexokinase-2 (hk2), L-lactate dehydrogenase A (ldha), Phosphoglycerate kinase 1 (pgk1) and Phosphoglycerate dehydrogenase (phgdh1a), were predominant in Un_Diff. Conversely, representative genes related to the mitochondrial respiration, such as Pyruvate dehydrogenase E1 component subunit alpha (pdha1a), Pyruvate dehydrogenase E1 component subunit beta (phdb) (Xu et al., 2018), ATP synthase-coupling factor 6 (atp5j; atp5pf) (Zhao et al., 2019b) and Cytochrome c oxidase subunit 8A (cox8a) (Rotko et al., 2021), were upregulated in Diff (Fig. 4I). This suggested a potential transition in energy source from glycolysis to mitochondrial respiration during spermatogonia development. Furthermore, it was found that the glycolysis-related genes were biasedly expressed in spermatogonia from testes at the Rs stage, whereas mitochondrial respiration genes showed higher expressions in spermatogonia at the Dif stage (Fig. 4J), further suggesting a potential variation in bioenergetic source transition from glycolysis to the TCA cycle in spermatogonia development during the annual reproductive cycle.

To elucidate the precise roles of mitochondrial respiration in spermatogonia development in black rockfish, spermatogonia were isolated from testes at the Rg stage and the purity was examined by VASA expression (Fig. S7A), along with genes marking germ cells and somatic cells (Fig. S7B). Subsequently, the spermatogonia were subjected to treatment with rotenone (inhibitor of mitochondria complex I) (Fig. S7C, left) and oligomycin (ATPase inhibitor, Fig. 7D, left), both resulting in the significant downregulation of the expression levels of genes related to proliferation [pcna, mki67, paf15, chaf1a and Cell division cycle-associated protein 3 (cdca3) (Uchida et al., 2012)] and differentiation [rec8, sycp1, Meiotic recombination protein SPO11 (spo11), cenpf and dmc1] (Fig. S7C,D, right). It was also found that there were fewer Diff cells in testes in the oligomycin injection group (Fig. S7E,F), suggesting that the inhibition of mitochondrial respiration could delay the spermatogonia differentiation in vivo. Of course, deeper insight into bioenergetic transitions during the development of teleost spermatogonia needs to be further investigated.

Cell type analysis of black rockfish testes unveiled a diverse array of somatic niche cells, comprising Sertoli cells, Leydig cells, stromal cells, macrophages and endothelial cells (Fig. 1). The primary somatic cells in the testes, Sertoli cells, Leydig cells and macrophages, were selected for further investigation to explore their development during the annual reproductive cycle of black rockfish testes.

Sertoli cells within the testes of black rockfish were found to be in close proximity to the spermatogonia, indicating a potential interplay with germ cells (Fig. 5A). To investigate the heterogeneity of Sertoli cells throughout the annual reproductive cycle, re-clustering analysis delineated four distinct subpopulations (Fig. 5B). Notably, Sertoli cells from the testes at the Rs stage were divided into two categories, with one primarily falling into cluster 1 and the other converging with Sertoli cells from the Rg and Dif stages (Fig. 5B), suggesting stage-specific characteristics within Sertoli cells. Examination of DEGs in Sertoli cells unveiled varied expression profiles of representative Sertoli cell-specific genes across different developmental stages (Fig. 5C). In mammals, anti-Mullerian hormone (Amh) and SRY-box transcription factor 9 (Sox9) are known marker genes highly expressed in fetal and neonatal testes, playing a crucial role in Sertoli cell maintenance (Ren et al., 2022). In zebrafish, amh has been found to maintain the homeostasis of spermatogonia, and its absence could lead to disorders of spermatogonia proliferation and differentiation (Lin et al., 2017). Similarly, in black rockfish, amh and sox9 were enriched in Sertoli cells at the Rs stage, suggesting their involvement in maintaining Sertoli cells and spermatogonia homeostasis. Moreover, Plasminogen activator inhibitor 1 (pai1; serpine1), a key regulator of the fibrinolytic system inhibiting tissue proteolysis (Stefansson et al., 2003), was expressed highly in Sertoli cells at the Rs stage. The high expression of pai1 in Sertoli cells at the Rs stage might serve to mitigate excessive tissue proteolysis and uphold the structural integrity of black rockfish testes during this stage. The expressions of several Sertoli cell markers, including pai1, Gonadal soma derived factor (gsdf) and Cholecystokinin (cckn), were validated by ISH, and the signals of these genes were specifically localized to Sertoli cells, exhibiting dynamic changes during testicular development (Fig. 5D).

To elucidate the underlying mechanisms governing the developmental dynamics of Sertoli cells, monocle analysis was performed on Sertoli cells derived from three distinct testicular developmental stages, resulting in a developmental trajectory delineating the lineage of Sertoli cells. The trajectory revealed a continuous progression from the Rs stage to the Dif stage, without major branching points (Fig. 5E), suggesting a linear developmental path of Sertoli cells during testicular development. Gene expression dynamics along pseudo-time were analyzed, leading to the identification of three distinct gene sets, each representing genes with specific expression patterns during Sertoli cell development (Fig. 5F).

Additionally, GO and KEGG enrichment analyses were performed to investigate the molecular events underlying different stages of Sertoli cell development. During the initial phase of Sertoli cell development, GO categories including ‘protein catabolic process’ and ‘protein folding’ were significantly enriched, indicating the importance of protein turnover and quality control in Sertoli cell functions at the Rs stage. As Sertoli cells progressed in their development, enriched GO terms in the Group 2 genes included ‘cell communication’ and ‘signal transduction’, suggesting the involvement of these processes in the development and functional specialization of Sertoli cells. Group 3 genes were enriched in GO terms related to ‘cellular biosynthetic process’ and ‘gene expression’, indicating the activation of biosynthetic pathways and gene regulatory mechanisms during Sertoli cell development at the Dif stage. The GnRH signaling pathway, known to promote gametogenesis and maturation in mammals and teleost (Fallah et al., 2020; Rastrelli et al., 2014), was also found to be involved in the Sertoli cells during Dif stage in black rockfish indicating its importance in maintaining spermatogenesis. (Fig. 5G). Overall, these findings provide valuable insights into the molecular events and biological processes underlying the development and maturation of Sertoli cells during the annual reproductive cycle of black rockfish testes.

Leydig cells, positioned interstitially between the seminiferous tubules in the testis, serve crucial functions including provision of nutrients to Sertoli cells, maintenance of the vascular barrier and secretion of hormones vital for supporting germ cell development (Crespo et al., 2016). In this study, Leydig cells were subjected to re-clustering analysis, yielding two discernible cell populations designated as sub1 and sub2 (Fig. 6A). DEGs analysis revealed that 83 genes were upregulated in sub1, whereas 185 genes were upregulated in sub2 (Fig. 6B). GO enrichment analysis showed that the genes upregulated in sub1 were associated with Apoptotic process and Vasculature development, whereas the enriched GO categories for sub2 upregulated genes primarily involved Cellular amide metabolic process and Peptide metabolic process (Fig. 6C). These findings indicated a molecular divergence between the sub1 and sub2 of Leydig cells. Additionally, DEGs related to steroid hormone synthesis, such as cyp17a1, Cytochrome P450 17A2 (cyp17a2), Cytochrome P450 11B (cyp11b; cyp11c1) and 3 Beta- And Steroid Delta-Isomerase 1 (hsd3b1) (Yang et al., 2022a,b), were observed among the different Leydig cell subsets (Fig. 6D). To validate the expression profiles of hsd3b1and cyp17a2 in the testis, dual fluorescence ISH was performed. The results demonstrated that the signals of hsd3b1 and cyp17a2 were not completely overlapping (Fig. 6E), indicating potential functional differences in steroid hormone synthesis between the two subsets of Leydig cells.

Testis-resident macrophages play crucial roles in tissue homeostasis, remodeling and immune response. In this study, the macrophage markers in the testes of black rockfish were compared with those in human, identifying 22 shared marker genes (Fig. S8A). These genes exhibited specific expressions in testis macrophages across both black rockfish and human, indicating their possible homologous functions in vertebrate testes (Fig. S8B). Notably, many of these genes have known roles in macrophage functions, such as Arachidonate 5-lipoxygenase-activating protein (alox5ap) (Ji et al., 2022), Granulocyte colony-stimulating factor receptor (csf3r) (Saunders et al., 2021) and C-X-C chemokine receptor type 4 (cxcr4a) (Tian et al., 2019). However, the specific functions of these genes in testis macrophages require further investigation. Subsequently, DEG analysis was performed on macrophages from testes at the Rs, Rg and Dif stages, revealing 48, 152 and 99 DEGs with distinct expression patterns, respectively (Fig. S8C). This indicated the diverse molecular characteristics of macrophages throughout testis development. Further GO analysis of DEGs aimed to uncover the potential functional divergence in macrophages with testis development. As shown in Fig. S8D, the DEGs of macrophages in Rs testes were significantly associated with proteolysis and hydrolysis, indicating potential roles in regulating testis regression. In the context of local trauma, macrophages have been linked to injury-induced muscle proteolysis (Watford, 2003). In mammals, macrophages protect sperm from immune response (Mossadegh-Keller and Sieweke, 2018) and participate in spermatogenesis by regulating sterol biosynthesis (Nes et al., 2000). Similarly, during the Rg and Dif stages, macrophages were significantly involved in sterol homeostasis and the immune system, implying their participation in spermatogenesis in black rockfish.

Somatic cells play a crucial role in the tissue microenvironment and engage in intricate interactions with spermatogenetic cells. To further investigate the interactions between germ cells and somatic cells in the testicular micro-environment, cell communication analysis was conducted using CellphoneDB, identifying ligand-receptor pairs between various somatic cells and germ cells. Among all cell types, Leydig cells exhibited robust interaction with four subsets of spermatogonia (Un-Diff1, Un-Diff2, Diff1 and Diff2), in particular with Un-Diff1 (Fig. 7A). Notably, Insulin-like peptide 3 (Insl3) in Leydig cells tightly bound to Relaxin receptor 1 (Rxfp1) in spermatogonia subsets (Fig. 7B), emphasizing the pivotal roles of Insl3 in spermatogonia development. In zebrafish, the secretion of Insl3 from Leydig cells has been demonstrated to promote the proliferation and differentiation of spermatogonia (Assis et al., 2016), as well as inhibit the apoptosis of male germ cells (Kawamura et al., 2004). Furthermore, insl3 and rxpf1 were validated to be specifically expressed in Leydig cells and germ cells in the testes of black rockfish, particularly in differentiated spermatogonia and spermatocytes at the Rg and Dif stage, respectively. (Fig. 7C). The direct interaction between Insl3 and Rxfp1 was further validated by in vitro co-immunoprecipitation (Co-IP) and bimolecular fluorescence complementation (BiFC) experiments (Fig. 7D; Fig. S9). Interestingly, although the Insl3-RXFP2 system has been extensively studied in mammalian testes (Feng et al., 2009; Ivell et al., 2020), Insl3-Rxfp1, instead of Insl3-Rxfp2, was identified in the black rockfish testes based on cell communication analysis. The potential lower affinity between Insl3 and RXFP2 in black rockfish, along with the expression dynamics of rxfp1 and rxfp2 in the testis (Fig. S10A,B), might contribute to the functional differentiation of Rxfp1 and Rxfp2 in black rockfish. Additional ligand-receptor pairs, such as CD74-MIF and CD74-APP in stromal cells, EGFR-MIF in Leydig cells and PDGFB-LRP1 in endothelial cells, were also identified and found to be associated with subsets of spermatogonia. The ligand-receptor pairs of CD74-MIF and CD74-APP have also been identified to be functional in mammalian immune response (Tanese et al., 2015; Li et al., 2022). Elucidating the specific roles of these ligand-receptor pairs in the context of spermatogonia development and overall reproductive biology in black rockfish awaits further investigation.

In conclusion, our study, using single-cell transcriptome sequencing, has provided a comprehensive roadmap of testicular development in black rockfish during the annual reproductive cycle (Fig. 8). We have identified candidate markers for germ cells and revealed the continuous nature of spermatogenesis. The dynamic regulation of proliferation and differentiation activity in spermatogonia subsets has been observed, indicating bioenergetic transitions during spermatogonia development. Furthermore, we have highlighted the heterogeneity and molecular characteristics of spermatogonia, Sertoli cells, Leydig cells, as well as macrophages, shedding light on their functional roles in testicular development. Moreover, our analysis has elucidated the crosstalk between germ cells and somatic cells, enhancing our understanding of the intricate cellular interactions within the testes. Our study contributes to a theoretical foundation of germ cell development and provides new insights into the regulation of testis development annual cycle in teleosts.

This study was approved by the College of Marine Life Sciences, Ocean University of China Institutional Animal Care and Use Committee on 10 October 2018 (project identification code: 20181010). All samples of black rockfish (21-month-old to 26-month-old) were obtained from the Qingdao Yushun Marine Technology Co., Shandong, China. At February, May and July in 2019, fish were transported to the laboratory and maintained in filtered seawater with continuous aeration and 12h light/12 h dark conditions. Individual fish were euthanized using 250 mg/L tricaine methane sulfonate (Sigma-Aldrich, MS-222). Then the testicular tissue from six individuals was dissected for the following experiments.

Each sample subjected to scRNA-seq consisted of testicular cells pooled from six adult black rockfish. The preparation of a single testicular cell suspension followed established protocols with minor adaptations outlined previously (Xu et al., 2019; Octavera and Goro, 2019). Briefly, the tunica albuginea of the testis was removed, and the testis was incubated in 10 ml Leibovitz-15 medium (Gibco Invitrogen) containing 20 mg/ml collagenase H (Roche Diagnostics), 0.25% trypsin (Worthington Biochemical), 5% fetal bovine serum (FBS; Gibco) and 0.05% DNase I (Roche Diagnostics) for 4 h at 28°C. Cell suspensions were filtered through 40 μm strainers. Following three rounds of centrifugation at 200 g and discarding the supernatant after each centrifugation, the testicular cells were suspended in 1 ml PBS. The cell number and motility were calculated using a blood counting chamber with Trypan Blue dying.

The scRNA-seq library preparation and sequencing were carried out at BGI Gene Technology Co. (https://www.bgi.com/). In summary, cells were diluted following the manufacturer's instructions, and 33.8 μl of total mixture buffer along with cells were loaded into the 10x Chromium Controller using Chromium Single Cell 3′ v2 reagents. After cDNA synthesis, 13 cycles were employed for amplification. The resulting libraries were then sequenced on a 2×150 cycle paired-end run using a BGISEQ-500 instrument. The scRNA-seq data generated in this study have been deposited in GEO under the accession number GSE252438. The scRNA-seq data of 12-month-old zebrafish testes used in this study was downloaded from https://github.com/asposato/zebrafish_testis_fertility. The gene expression table derived from scRNA-seq data and R-objects generated in this study are available from the Broad Institute Single-cell Portal (https://singlecell.broadinstitute.org/single_cell/study/SCP2464). The scRNA-seq data from humans used in this study was downloaded from GEO (GSE120508).

The raw sequencing data were processed using CellRanger (v2.20). Fastq files were run through the count application with default settings, conducting alignment (using STAR aligner) based on black rockfish reference genome (CNGB, CNP0000222) (He et al., 2019), as well as filtering and UMI counting. The resulting UMI count matrices were employed for further analysis.

The Seurat program (https://satijalab.org/seurat/, R package, v2.3.4) was employed for cell clustering analysis. Initially, UMI count matrices from each sample were loaded into R using Read 10x function and Seurat objects for each sample were constructed. The data were filtered and normalized according to the default settings. Then the data matrices from distinct samples were integrated following guidelines outlined in the tutorial (https://satijalab.org/seurat/). Then uniform manifold approximation and projection (UMAP) (Becht et al., 2019) was employed to identify the cell clusters using the combined data. The ‘FindMarkers’ function (a Wilcoxon rank sum test) was used to determine DEGs between clusters, with criteria set at a minimum expression in 25% of cells and a log fold change >0.25. The known markers used for cell type identification and references are listed in Table S1. The ‘DoHeatmap’ function was used to generate the expression heatmap for target genes and features. Pearson correlation analysis was performed by R software (R v4.0.3, https://cran.rproject.org/src/base/R-4/). The Omicshare platform (https://www.omicshare.com/tools/) was used to perform GO and KEGG analysis using DEGs. The ‘AverageExpression’ function was carried out to calculate the average expression of marker genes of different clusters, and TBtools (Chen et al., 2020) was used to generate the heatmap. Single-cell pseudo-time trajectories were constructed using the Monocle 2 package (v2.10.1) (Qiu et al., 2017) following default settings (http://cole-trapnelllab.github.io/monocle-release/). The UMI count matrix was loaded into R, then the genes (expressed in at least 1% cells with expression >0.5) were used to define the cell trajectory by the ‘reduceDimension’ and ‘orderCells’ function. DEGs (q-value<1e-200) were used for trajectory analysis. The ‘plot_pseudotime_heatmap’ function was used to generate the heatmap.

Hub genes in Un-Diff1 and Un-Diff2 were identified using WGCNA (https://github.com/Lindseynicer/WGCNA_tutorial, R package, v1.68) performed on normalized gene expression data derived from DEGs of different types of cells, as described previously (Li et al., 2016). In brief, the topological overlap matrix (TOM) was constructed with a softPower value of 12. Then, the genes were hierarchically clustered using 1-TOM as the distance measure, and the modules were determined for the resulting dendrogram. Finally, the networks (hub genes with top 5% weight values) of Un-Diff1 and Un-Diff2 were visualized by Cytoscape Software 3.7.2. TFs specifically expressed in different spermatogonia subsets were identified based on the conserved DNA-binding domains (DBDs), as described previously (Zhang et al., 2021a). Protein sequences of marker genes from different spermatogonia clusters were obtained and random sequences were removed through CD-HIT (Huang et al., 2010). Protein domains were annotated using Hidden Markov Model (HMM) profiles sourced from the Pfam database (v31.0), employing the hmmscan program within the HMMER package (v3.1b2) (Finn et al., 2016). Genes containing DBDs were screened with the e-value threshold of 10⁻⁴ and regarded as TFs. Additionally, some genes containing DBDs were screened with the threshold in Animal TFDB 3.0, such as bHLH domain (10⁻²), HMG, Homeodomain, zf-BED, zf-C2H2 domains (10⁻³) and zf-CCCH domain (10⁻²⁰). The network of TFs was constructed based on the results of WGCNA and visualized by Cytoscape Software 3.7.2.

CellphoneDB (Efremova et al., 2020) (v2.032) was used with its statistical methodology to ascertain and characterize ligand-receptor interactions among distinct cell type pairings. The heatmap of ligand-receptor pair numbers and dotplots of expression for ligand-receptor pairs between specific cell types were created by R software.

The testis samples for Hematoxylin and Eosin (H&E) staining were fixed in Bouin solution for 24 h and stored in 70% ethanol at 4°C. Subsequently, it was dehydrated with a serial ethanol (80% ethanol, 2 h; 90% ethanol, 2 h; 2× 100% ethanol, 30 min), cleared in xylene, embedded in paraffin and sectioned at 5 μm thickness. The tissue sections were stained with H&E (Solarbio). The results of H&E staining were captured using Nikon Eclipse Ti-U microscope. Each experiment was performed in triplicate. The germ cell type identification followed the germ cell characteristics of teleosts that have been demonstrated in detail in previous studies (Schulz et al., 2010).

The digoxigenin-labeled RNA probes for detected genes in this study were synthesized using a DIG RNA labelling kit (Roche). Testis samples for ISH were fixed in 4% paraformaldehyde for 24 h. Subsequently, the tissues were dehydrated with a series of methanol concentrations (30%, 50%, 70%, 80%, 90% and 100%), cleared in xylene, embedded in paraffin and sectioned at 5 μm thickness. ISH was conducted following the detailed protocols of our lab, as described previously (Wang et al., 2021). The images were viewed using Nikon Eclipse Ti-U microscope with 40× optics. The area in blue-purple was regarded as positive signal. Each experiment was performed in triplicate.

The tissue sections were deparaffinized and hydrated, treated with sodium citrate solution (Solarbio) for antigen retrieval and 3% hydrogen peroxide for blocking endogenous peroxidase. The sections were then blocked with 5% bovine serum albumin diluted by 1×PBS for 1 h and incubated with PCNA antibody (Ribeiro et al., 2017) (sc-56, Santa Cruz Biotechnology, 1:500), DMC1 antibody (Zhang et al., 2022) (ab11054, Abcam, 1:400) and VASA antibody (Zhai et al., 2022) (GTX128306, GeneTex, 1:500) at 4°C overnight. Subsequently, the HRP-conjugated rabbit IgG (CW0103, CWBIO, 1:200) or mouse IgG (CW0102, CWBIO, 1:200) were used as secondary antibodies and diluted 1:200 for incubation for 30 min at room temperature. Finally, the positive signals were visualized using the DAB Chromogenic Reagent Kit (Solarbio) and sections were counterstained with Hematoxylin. The results of IHC were observed and photographed under a Nikon Eclipse Ti-U microscope.

For in vitro assays, the spermatogonia were isolated from testes at Rg stage by density gradient centrifugation using percoll reagent, as described previously (Sato et al., 2017). Then, the purity of spermatogonia was identified through immunostaining for VASA and quantitative reverse transcription-polymerase chain reaction (qRT-PCR) for several marker genes of germ cells and somatic cells. The isolated spermatogonia were cultured in DMEM/F12 medium supplemented with 10% FBS (Gibco) and 1% penicillin and streptomycin (Solarbio). After 48 h, the cells were treated with oligomycin (an ATPase inhibitor) and rotenone (an inhibitor of mitochondria complex I) in different concentrations. Following an additional 48 h, the cells were collected for further RNA extraction and gene expression analysis. For in vivo assays, oligomycin was injected into the testes of black rockfish (554±27 g) through cloacal orifice at a dosage of 200 μM/kg. The control group was injected with an equal volume of PBS. After two injections at 15 day intervals, testes samples were collected for further analysis.

Total RNA was extracted using TRIzol reagent (Invitrogen) following the standard protocol. Genomic DNA removal and cDNA synthesis were performed with All-In-One 5× RT MasterMix (Applied Biological Materials). The subsequent qRT-PCR analysis followed established procedures as described previously (Wang et al., 2023b). Each experiment was performed in triplicate.

The Co-IP assay was performed as described in our previous study (Qu et al., 2022). Briefly, the coding sequences of insl3 and rxfp1 were cloned from black rockfish to construct recombinant Insl3-pEGFP-N1 and pBiFC-Rxfp1-VN173 vectors, respectively. Then, Insl3-pEGFP-N1 was co-transfected with pBiFC-Rxfp1-VN173 into HEK-293 T cells cultured in 12-well plates. After 48 h, cells were washed three times with PBS and collected. Co-IP was performed using BeaverBeads^® Protein A (or A/G) Immunoprecipitation Kit (BEAVER) following standard procedures. Subsequently, the results were analyzed by western blot. Antibodies used were: anti-GFP (bs-0890R, Bioss, 1:2000), anti-FLAG (bsm-33346 M, Bioss, 1:2000), goat anti-rabbit IgG, HRP conjugated (CW0103, CWBIO, 1:5000), and goat anti-mouse IgG, HRP conjugated (CW0102, CWBIO, 1:5000).

The recombinant pBiFC-Insl3-VC155 was co-transfected with pBiFC-Rxfp1-VN173 into HEK-293 T cells. The cells transfected with pBiFC-Insl3-VC155 and pBiFC-Rxfp1-VN173 were regarded as negative control, whereas cells co-transfected with pBiFC-bFos-VC155 and pBiFC-bJunVN173 were considered the positive control. The results were visualized and imaged using a Nikon Eclipse Ti-U microscope at 48 h post transfection. Each experiment was performed in triplicate.

All qRT-PCR data are presented as the mean±s.e.m. with n=3. Statistical differences between groups were calculated using one-way ANOVA by Independent-Samples t-test (SPSS 20.0, IBM). P<0.05 was considered indicative of a significant difference between groups.

We thank Dr Shang Liu from University of Chinese Academy of Sciences for his help with scRNA analysis.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.202296.reviewer-comments.pdf

This article is part of the Special Issue ‘Uncovering developmental diversity’, edited by Cassandra Extavour, Liam Dolan and Karen Sears. See related articles at https://journals.biologists.com/dev/issue/151/20.