5 month-old adult AB line zebrafish (Danio rerio) were prepared and treated with tricaine methanesulfonate (MS-222, 0.25%) on ice for 15 min before experiment. The study was conducted in line with the Chinese law for the Protection of Animals, and the animals were treated properly. The testes samples collected from those zebrafish were first washed by PBS three times, and then were digested in 10 ml 0.25% trypsin at 37°C for 15 min, during which process the tissues were pipetted up and down every 3 min. The cell suspension was first filtered through 70 μm nylon mesh after stopping the digestion process by DMEM (with 10% FBS), and then were centrifuged at the speed of 1,000 rpm for 10 min. After centrifugation, the cells were suspended in 1 ml DMEM (with 10% FBS) and were further filtered through 40 μm nylon mesh. Before loading onto the 10x chromium chip, the cells were washed by BSA DPBS (0.04%, three times) and were resuspended to a concentration of 800∼1,000 cells/μl (viability >85%). The single cell library preparation was carried out according to the manufacturer’s protocol (Chromium Single Cell 3′ Reagent Kits V3 Chemistry). The cDNA library was further size filtered and sequenced on the Illumina NovaSeq 6,000 System following 150 bp paired-end sequencing protocol.

Cell ranger software (v5.0.0) from 10X genomics was used to preprocess the raw sequencing data, and primary UMI filtering was carried out by Cell ranger following criteria “No N contained in the UMI,” “all base quality should be no less than 10,” and “not a homopolymer sequence.” For cell ranger processing, zebrafish genome and gene annotation from Ensembl were used as reference datasets (Ensembl Gene 102). Then, R package Seurat (V4.0.0) was used for further data filtering, data normalization, cell clustering, and cluster-level marker gene discovery. The UMI matrix data was filtered at two dimensions: gene filtering and cell filtering. For gene filtering, only the genes that were detected at least in 5 cells were retained. For cell filtering, the retained cells should contain the number of expressed genes within the range between 200 and 5,000, and the proportion of UMIs from mitochondria should be less than 5%. The raw UMI matrix was normalized by sctransform, which fitted the data through regularized negative binomial regression. Three types of features were used for data regression: total number of expressed genes, total number of UMIs, and the percentage of UMIs from mitochondria. In order to cluster cells, the principal components (PC) were computed by principal components analysis (PCA), and the first fifteen PCs were used for cluster finding. K nearest neighbors (KNN) method was applied to cell clustering, where the K was set to be 20. Finally, the marker genes for each cell cluster were statistically computed by Wilcoxon Rank Sum test, and only the genes detected in at least 25% of the cells (in either of the two tested cell clusters) were included for testing. The marker genes were defined as the genes with p value less than 0.01 and the log transformed (base 2) fold change larger than 0.25.

To annotate the cell clusters identified in this study, the testicular cell-type marker genes were collected from previous studies in zebrafish (Leal et al., 2009; Chen et al., 2013; Assis et al., 2016; Safian et al., 2016; Lin et al., 2017; de Castro Assis et al., 2018; Takemoto et al., 2020) or other species (Hermann et al., 2018). In the further analyses, only the collected marker genes showing cell-type specific expression patterns in our single cell studies were used for annotation. The final list of known marker genes and their sources can be found in Supplementary Table S1.

In order to study the cell-cell communication among cell populations, the Ligand-Receptor (LR) pairs were analyzed, which followed the procedure with minor modifications described in one recently published paper (Shi et al., 2021) about LR studies in mouse. In detail, CellTalkDB (Shao et al., 2020) was selected as the database for LR analyses, which contained 3,398 manually curated human LR pairs. The orthologous genes between human and zebrafish were downloaded from Ensembl biomart (Ensembl Genes 102), and then 13,004 one-one orthologous genes were retained for LR analyses. SingleCellSignalR (Cabello-Aguilar et al., 2020) was applied to the inference of intercellular LR interaction networks, where the default LR database was replaced by CellTalkDB and only the genes found in the zebrafish-human one-one orthologous gene list were included for the analysis.

Adult zebrafish (5 months after fertilization) were anesthetized in 0.25% ms-222 (ethyl 3-aminobenzoic; Sigma, E10505) in filtered system water. Zebrafish testes were fixed with 4% paraformaldehyde (Sigma, P6148) overnight at 4°C. Then, the fixed samples were embedded in OCT (Sakura, 4583) after dehydration of sugar, sectioned at 10 μm thickness, and stained with hematoxylin and eosin (BBI, E607318).

Probe fragments were first amplified from zebrafish testes cDNA libraries (primers information in Supplementary Table S2) and then inserted into pGEM-T-easy vector (Promega, A1360). Digoxigenin-labeled antisense RNA probes were made by using DIG-RNA labeling Kit (Roche, 11175025910). Zebrafish testes were cut into 12 μm slices by cryotome under temperature between −25°C and −23°C. Slides were first fixed with 4% PFA, then digested with 0.1% Proteinase K, followed by washing in PBS. Lastly, the digested slides were first incubated with the digoxygenin-labelled RNA probes at 68°C overnight and then incubated with alkaline phosphatase-conjugated anti-DIG antibody, followed by treatment with AP-substrate NBT/BCIP solution (Roche, 11681451001) for visualization. Brightfield images of sections were obtained with Nikon Eclipse Ni-U microscope.

To visualize the cell clusters and the gene expression patterns in two-dimensional graphs, Uniform Manifold Approximation and Projection (UMAP) was applied to dimensional reduction, where the first 15 PCs from PCA were used as the input for UMAP. ClusterProfiler (Yu et al., 2012) (version 3.18.0) was used for GO enrichment analyses, and the p values were corrected by Benjamini-Hochberg (BH) method. Only the GO terms with BH corrected p value less than 0.05 were retained as significant ones.

The raw sequence data reported in this paper had been deposited in the Genome Sequence Archive (Wang et al., 2017) in National Genomics Data Center (Members, 2021), China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences, under accession number CRA003925 that are publicly accessible at https://bigd.big.ac.cn/gsa. All the codes used for this study can be found in https://github.com/gangcai/zebTestis.

Testes from five adult male zebrafish were pooled for 10X genomics single RNA sequencing, initially there were 16,032 single cells detected, and after quality filtering 14,315 single cells were remained and used for further analyses (Figure 1; Supplementary Figure S1A). 10 different cell types were identified, including spermatogonia subpopulation 1 and 2 (c3, c5), spermatocyte (c1), early round spermatid (c4), middle round spermatid (c6), late round spermatid (c0), elongated spermatid (c2), Leydig (c9), Sertoli (c7), and a small proportion of red blood cells (c8) (Figure 2A). In detail, the cell type with highest proportion was late round spermatid, and the second highest one was spermatocyte, which were 3,889 and 2,396 single cells respectively (Supplementary Figure S1B). Notably, the proportion of two somatic cells (Leydig and Sertoli) is relatively small compared to the germ cells, in total, we identified 58 Leydig cells and 133 Sertoli cells (Supplementary Figure S1B). The top three cell types with largest number of detected genes are SPG2, early round spermatid and SPG1, and the average number of genes is also positively corrected with the average number of UMIs (Figure 2B; Supplementary Figures S1C,D). Lower gene expression was detected in late-stage cell types of spermatogenesis (elongated spermatids, and late round spermatids) (Supplementary Figure S1D). All the cell types identified contain lower proportions of expressed mitochondrial genes (0.5–2% on average, Figure 2B), which indicates rare contamination of apoptotic, stressed, or lower quality cells.

The cell types were annotated based on known marker genes examined in zebrafish testis or the testis of other species (Supplementary Table S1; Supplementary Figure S2). For example, star was reported as the marker gene of Leydig cell (Lin et al., 2017), and our scRNA-seq data clearly showed that it was specifically expressed in cell cluster c9. Furthermore, based on previous studies, sycp2 had been recognized as a marker gene for spermatocyte cells (Takemoto et al., 2020), while dazl was mainly expressed in SPG and weakly detected at later spermatogenesis stages (Chen et al., 2013). In this study, sycp2 was found to be significantly enriched in c1, and dazl was discovered to be highly enriched in c3 and c5. Supplementary Table S3

Next, the cluster level marker genes were examined, in total, thousands of marker genes were identified for each cell type, including 1,449 marker genes for SPG1, 1,271 marker genes for SPG2, 418 genes for Sertoli cells, and 412 genes for Leydig cells (Figure 3A; Supplementary Table S3). The top 10 marker genes based on the rankings of gene expression foldchange can be found in Figure 3B, including star for Leydig (c9), sycp2 for Spermatocyte (c1), hbba1/hbba2 for red blood cell (c8).Based on the marker genes, GO enrichment analysis was performed: “cilium assembly,” “cilium movement,” “cilium organization” GO terms were significantly enriched in spermatids (c4, c6, c0, c2), and “DNA recombination,” “chromatin organization,” “mRNA metabolic process” GO terms were significantly enriched in early-stage cell types of spermatogenesis (c1, c3, c5) (Figure 3C). Notably, “ribosome assembly” and other ribosome related terms were enriched in non-spermatids cell types, which might indicate lack of translation events in spermatids.

The cluster-level marker genes provide a way to find novel spermatogenesis marker genes. We found edrf1 (erythroid differentiation regulatory factor 1) was specifically expressed at elongated spermatids, and larp1b (La ribonucleoprotein 1B) RNA was enriched in round spermatids (Figure 3B). We further carried out RNA in situ hybridization (ISH) to validate the novel testicular cell type specific marker genes. Firstly, for better recognition of each cell types in zebrafish testis, we annotated the cell types based on the cell shapes and relative locations in the H&E staining image (Figure 3D). As we expected, the ISH staining for larp1b and edrf1 RNA showed positive signal in round spermatids and elongated spermatids respectively, and the ISH staining for krt18a.1 and sumo3b showed positive signal in Sertoli and Spermatogonia (SPG) cells respectively (Figure 3E). Furthermore, we also provided all the novel marker genes for each zebrafish testicular cell type in Supplementary Table S3.

For zebrafish spermatogonia cells, two subpopulations were found in this study: SPG1 and SPG2 (Figure 2A). In total, 656 genes were significantly differentially expressed between SPG1 and SPG2 (Figure 4A, Supplementary Table S4). We found histone related genes involved in heterochromatin assembly were enriched in SPG2, such as hist1h2a6, hist1h4l.6, hist1h4l.11 (Figure 4B). Further GO enrichment analysis based on SPG subpopulation specifically expressed genes confirmed that the GO terms such as “condensed chromosome,” “DNA packaging complex” were significantly enriched in SPG2 (Figure 4C). In human, there were three types of spermatogonia: Type A dark, Type A pale, and Type B, where type A pale went through division every seminiferous epithelial cycle (Kolasa et al., 2012). In zebrafish, there were also three types of spermatogonia: type A undifferentiated (A_und), type A differentiated (A_diff), and type B spermatogonia (Schulz et al., 2010). The enrichment of GO term “condensed chromosome” suggests more histone proteins involved for chromosome DNA packaging, as shown in Figure 4B histone protein genes such as hist1h2a6, hist1h4l.16 were significantly enriched in SPG2. Further ISH experiment for hist1h4l.6 RNA confirmed its highly specific expression pattern in spermatogonia cells (Figure 4D). The enriched GO terms and high expression level of histone protein genes in SPG2 suggests that SPG2 was actively undergoing mitosis and differentiation, which was in line with the definition of type A differentiated spermatogonia in zebrafish (A_diff) and type B spermatogonia. In contrast to SPG2, SPG1 was not enriched with histone protein genes and DNA packaging complex (Figure 4C), which indicates that SPG1 was not actively involved in mitosis, and it might be the type A undifferentiated spermatogonia cells (A_und). New marker genes for A_und were found based on our scRNA-seq datasets, such as eno3, e2f5, and ripply2 (Figure 4B). Interestingly, ENO3 gene might be involved in muscle regeneration, and had been recognized as one of the marker genes during muscle regenerative phase of maturation (Forcina et al., 2020). Furthermore, E2F5 gene belongs to E2F family of transcription factors, and it might play a role as a tumor suppressor gene and control the cell cycle process (Chen et al., 2009).

Lastly, we studied the influence of somatic cells on zebrafish germ cells and the spermatids. As shown in Figure 5A, more LR interactions were observed between Leydig-Germ cells than Sertoli-Germ cells, which demonstrated a stronger influence of Leydig cell on germ cells than Sertoli cells. However, based on previous studies, Sertoli cell played important roles in germ cell differentiation, survival, development, and physiological functioning (Schulz et al., 2010), which made us expect more LR interactions emitted from Sertoli than Leydig. In our study (Figure 5A), both of the two somatic cells (Sertoli and Leydig, as LR interaction signal emitter) had most abundant paracrine interactions with SPG2 (as LR interaction signal receiver), while the number of LR interactions emitted from Leydig (116) is much more than the number of LR interactions emitted from Sertoli (12). Among two spermatogonia cells, both of Leydig and Sertoli cells exhibit higher influence on SPG2 than SPG1. The detail LR interactions between somatic-SPG cells were shown by chord diagrams (Figures 5B–E). Calm3a-plpp6 was predicted with highest LR interaction score (LRscore) between Sertoli (emitter) and SPG1 (receiver) (Figure 5B), and hsp90aa1.2-erbb2 showed the strongest paracrine interaction between Sertoli and SPG2 (Figure 5C). The top LR interactions between Leydig and SPG1 include calm3a-adcy8, vegfaa-grin2bb, hsp90b1-tlr7, appb-FP236542.1 (Figure 5D), and the top LR interactions between Leydig and SPG2 include appb-cav1, calm3a-mylkb, hsp90aa1.2-cftr, adma-gpr182 (Figure 5E). We also analyzed other paracrine interactions between each testicular cell, and the full list of LR interactions and their LR interaction scores can be found in Supplementary Table S5. Our result based on the comparative analyses of the paracrine interactions between somatic-germ cells suggested the importance for the future study of the paracrine influence of Leydig on Germ cells.