Insight into 2,3,7,8-tetrachlorodibenzo-p-dioxin-induced disruption of zebrafish spermatogenesis via single cell RNA-seq.

2,3,7,8-Tetrachlorodibenzo-p-dioxin (TCDD) is a potent and environmentally persistent endocrine disrupting chemical. Our previous work demonstrated the latent reproductive maladies of early-life TCDD exposure in zebrafish. Zebrafish acutely exposed to low, environmentally relevant levels of TCDD (50 pg/mL) during two windows of sexual differentiation in development (1 hour of exposure at 3 and 7 weeks postfertilization) were later infertile, showed a reduction in sperm, and exhibited gene expression consistent with an altered microenvironment, even months after exposure. Due to the highly heterogeneous cell- type and -stage landscape of the testes, we hypothesized various cell types contribute markedly different profiles toward the pathology of TCDD exposure. To investigate the contributions of the diverse cell types in the adult zebrafish testes to TCDD-induced pathology, we utilized single-cell RNA-seq and the 10x Genomics platform. The method successfully captured every stage of testicular germ cell development. Testes of adult fish exposed during sexual differentiation to TCDD contained sharply decreased populations of late spermatocytes, spermatids, and spermatozoa. Spermatogonia and early spermatocyte populations were, in contrast, enriched following exposure. Pathway analysis of differentially expressed genes supported previous findings that TCDD exposure resulted in male infertility, and suggested this outcome is due to apoptosis of spermatids and spermatozoa, even years after exposure cessation. Increased germ cell apoptosis was confirmed histologically. These results provide support for an environmental exposure explanation of idiopathic male infertility.

Analyzing the reproductive consequences of environmental exposures throughout the lifespan is a critical step in understanding adult male infertility. With most human male infertility considered idiopathic, these findings provide a toxicogenomic profile of the sequelae of peripubertal exposure to dioxin and, thus dioxin-like contaminants, wherein spermatids and spermatozoa are diminished in testes exposed to TCDD, while earlier germ cell populations are enriched. This manuscript serves as a useful reference for inevitable future scRNA-seq studies on zebrafish testes, which are comparable in function to mammalian testes, yet express distinct genes during spermatogenesis.

Introduction

Spermatogenesis is the vital process of germ cell development in the seminiferous tubules of the testes. Spermatogonial stem cells (SSCs) self-renew and generate daughter spermatogonia (SPG), which undergo meiosis to form spermatocytes (SPCs). A second meiosis followed by dramatic morphological remodeling results in the haploid sperm cells essential for sustaining human life. This is a tightly controlled, orchestrated process; thus, disruption early in life can result in major defects. Exposure to endocrine-disrupting compounds (EDCs) during critical stages of development can lead to reproductive dysfunction later in life, as the endocrine system is poised to respond to miniscule changes in hormone signaling. The psychological and economic burdens of exposure-related reproductive outcomes such as infertility are detrimental to men and couples alike (1, 2). Infertility treatment globally was a $14 billion industry in 2020, and expected to rise 10% each year (3). Male infertility accounts for about 50% of couples’ infertility (4), and is often considered idiopathic, hence an active area of research.

It is widely understood that EDCs in the environment contribute to male reproductive dysfunction (5). One environmental EDC notoriously damaging to fertility (6) is 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD): a potent, enduring byproduct of industrial manufacturing of pesticides and other products. TCDD is the most toxic of the class of “dioxins,” exerting the strongest influence on the aryl hydrocarbon (AhR) endocrine disruption pathway, earning it the status of model dioxin ligand (7), and therefore, a suitable compound to additionally represent the potential effects of the dioxin-like compounds routinely found in the general population, such as polychlorinated biphenyls (PCBs) and furans (8).

The developmental origins of health and disease (DOHaD) hypothesis posits that early-life exposure to EDCs can have long-ranging effects (9). Furthermore, the sequelae of exposure can vary based on the timing of exposure. An illustrative case is the major TCDD explosion event in Seveso, Italy, where men who were exposed to high TCDD levels during the explosion between the ages of 1 to 9 years old, but not those exposed ages 10 or older, later displayed reduced sperm and sperm motility (10). However, less is known about the long-term reproductive effects of low, environmentally relevant levels of exposure to EDCs such as dioxin during critical windows of susceptibility during development. In a rare human study measuring the influence of background levels of TCDD on sperm parameters, Mínguez-Alarcón et al. (11) took TCDD levels in peripubertal 8- to 9-year-old boys, measured semen parameters 10 years later, and found lowered sperm concentration, count, and motility associated with higher background TCDD levels. The link between early-life low-dose dioxin exposure and adult infertility has, in contrast, been well-studied in rodents and aquatic life (12–16). The NIH-approved model organism, the zebrafish (Danio rerio), has emerged as a useful tool to study these effects of exposure across the lifespan due to its quick generation time and genomic similarity to humans. The process of gonadal differentiation in zebrafish begins around 20 days postfertilization (dpf) (17) and is completed around 42 dpf (18). At the outset, zebrafish possess a juvenile ovary (19), which can differentiate to testes during this window, relying on a mixture of genetic and environmental factors for direction of the binary gonad fate (20). Due to the presence of gonocytes at 3 weeks postfertilization (wpf), zebrafish testes at this exposure timepoint resemble prenatal rather than postnatal human testes (21). By 7 wpf, gonocytes have been replaced by SPG (22), resulting in a landscape more similar to human postnatal testes. We have previously shown brief perinatal exposure to low levels of TCDD (1 hour, 50 pg/mL) at critical points of gonadal differentiation and maturation in zebrafish (21 and 49 dpf) led to decreased reproductive capacity and altered spermatogenesis as evidenced by changes in elicited fecundity (15) and altered cell populations of the seminiferous tubule (16).

The testes contain a rich diversity of cell types and cells in various stages of development, and dysfunction or arrest in one cell type alone can be associated with infertility (23–28). The single cell RNA sequencing (scRNA-seq) pipeline has emerged as a valuable tool for uncovering individual cellular functions in thousands to millions of cells, an advancement over the bulk RNA-seq method of averaging gene expression across all cells in a tissue. Cell type-specific transcriptome dynamics in heterogeneous tissue can be undetected with methods that require analyses of pooled, presumably uniform cells. Previous scRNA-seq studies have investigated the neonatal (29) and adult (30–38) testes in mammalian species, and one other adult zebrafish testes dataset exists as part of a scRNA-seq atlas (39). We examined the adult zebrafish testes in detail from an exploratory and toxicological standpoint to characterize the baseline testicular landscape in this emerging model organism, and to understand the long-term effects of early-life exposure to an environmental contaminant.

Given our histological findings that TCDD altered populations of specific testicular cell types but not others (16, 40) in infertile fish, we sought to further illuminate TCCD-induced testicular cell population alterations at the single-cell level, and indicate the cellular dynamics of the spermatozoa contributing to their loss. We employed scRNA-seq using the 10x Genomics platform to query the transcriptional landscape of the testes in adult zebrafish that previously underwent early-life exposure to TCDD (16). We found TCCD induced significant loss of testicular cells from the pachytene SPC stage through mature spermatozoa, with a concomitant enrichment of SPG and earlier stage SPCs. Despite differences between zebrafish and mammals, there was similarity among marker genes from our data and published datasets from mammals. Given the zebrafish's increasing popularity in toxicogenomics research, this work will serve as a resource for future single cell and transcriptomic research in the zebrafish testes, and contribute to understanding the roots of early developmental exposure-induced long-term reproductive outcomes such as infertility due to EDCs.

Results and Discussion

The scRNA-seq pipeline can be successfully employed for characterization of zebrafish testes

To understand later-life effects on testicular development following early-life exposure to TCDD, while retaining important cell type-specific differential information, we used the 10x Genomics platform to capture the transcriptomes of ∼7,500 single testicular cells from zebrafish exposed to either TCDD or a control vehicle solution of 0.1% dimethyl sulfoxide (DMSO) during sexual differentiation and maturation. Since we previously observed that TCDD-induced changes in certain testicular germ cell populations, but not others, may contribute to observed male infertility (16), single cell RNA-seq (scRNA-seq) was critical to understanding, from a mechanistic standpoint, how exposure disrupts the dynamic process of spermatogenesis. We sought to extend our base of knowledge of cell type-specific population alterations following exposure, and to establish a transcriptional roadmap of each cell type in the zebrafish testes prior to and following exposure. Fig. 1(a) illustrates the overall workflow. The method successfully captured an average of 31.9% of the transcriptome for each biological replicate of control (N = 3) and TCDD (N = 2) testes. Further characterization of the datasets and quality control (QC) metrics are detailed in Figures S1 and S2 (Supplementary Material) and Tables S1 and S2 (Supplementary Material).

Fig. 1.

Single cell transcriptome profiling of TCDD-exposed zebrafish testes. (a) Illustration of experimental workflow. wpf: weeks postfertilization. (b) t-distributed stochastic neighbor embedding (tSNE) plot of 10 unbiased clusters of combined control and TCDD testicular datasets (N = 7,428 cells) encapsulating every stage of germ cell development from SSC to spermatozoa. Each dot represents a single cell. Cluster cell types (legend) were discerned by marker gene expression. SSC: spermatogonial stem cell. SPG: spermatogonia. SPC: spermatocyte. RS: round spermatid. (c) Expression profiles of select marker genes in testicular cells projected onto tSNE plots.

The 10x platform recovered the full lineage of gametogenic cell types and developmental stages in the testes from SSCs to spermatozoa. T-distributed stochastic neighbor embedding (tSNE) plots reduced dimensionality and aided visualization of the dataset, which projected ten clusters of this range of cell populations (Fig. 1b). Each dot in a cluster represents a single cell. Cell types represented in each cluster were identified by cell-type-specific marker expression patterns established in the literature; these clusters are denoted by color in the tSNE. Expression of marker genes characteristic of each general cell type are projected onto the tSNE clusters in Fig. 1(c), and functional expression is highlighted in Fig. 2. The size of the dot represents the relative number of cells in the cluster expressing the gene. Predictably, clusters follow waves of expression corresponding to overarching developmental processes in spermatogenesis. SPG exhibit transcription associated with proliferation and differentiation; meiosis genes are expressed mostly strongly by SPC; sperm architecture genes come online in anticipation of remodeling the late SPC into RS. The male epigenome in zebrafish is similarly inherited by the offspring as in humans, however, in zebrafish there is an absence of some canonical mammalian proteins such as protamines (41). The dot plot suggests that zebrafish testes make use of histones, methyltransferases, transcription factors, and retrotransposon silencers to accomplish the epigenetic maintenance necessary to germ line conservation. The gene encoding DMRT1, a DNA-binding protein involved in male sex determination in zebrafish (42), is strongly present as expected in SSCs (43). Histone expression largely disappears after the ES stage, when sperm have fully condensed.

Fig. 2.

Expression of marker genes and functional suites of genes supporting spermatogenesis, by cluster. Dot plot illustrates the % of cells in each cluster expressing each transcript. Human orthologues: zte38/HORMAD2; hist1h2a2/H2AX.

Established markers of cell type and/or of cellular function facilitated identification in conjunction with spatial relations within the tSNE, as the temporal order of cell development in the testes is retained in the plot: cluster 1 contains a mix of SSCs and undifferentiated type A spermatogonia (SPG-A; dnd1) (44); cluster 2: differentiated type B spermatogonia (SPG-B; pcna) (45); cluster 3: preleptotene SPCs; cluster 4: leptotene and zygotene SPCs (L/Z SPC; sycp3) (46); cluster 5: zygotene and pachytene SPCs (Z/P SPC; HORMAD2; D. rerio zte38) (47); cluster 6: pachytene and diplotene SPCs (P/D SPC; sycp1) (48); cluster 7: round spermatid (RS) and elongating spermatids (ES; armc3) (32), cluster 8: ES (cfap52) (49); whereas clusters 9 and 10 are spermatozoa (tssk6; spag8) (50, 51). As spermatogenesis is a continuous developmental process, most clusters represent a spectrum of multiple stages within a cell type, and cannot be defined computationally by a discrete stage. The presence of two spermatozoa clusters is interpreted elsewhere by Guo et al. (31) to reflect only one relatively similar sperm population, albeit with varying levels of RNA degradation status amongst the sperm cells during their transcriptional quiescence. Cluster 3 (preleptotene SPC) lacked classical preleptotene markers of meiosis initiation, expressing mostly histone-related genes instead. Cluster 3’s most highly expressed marker was D. rerio si: ch211-113a14.24, an orthologue of H1, H2, and H5 (average log2-fold change (LFC) 2.18), and expressed the H3-like transcript (D. rerio LOC108190746) at an average LFC of 1.36. Shiraishi et al. (52) reported histone H3.5 localized to preleptotene SPCs (D. rerio hist2h3c, average LFC 1.25), and Ueda et al. (53) report a testis-specific H3 isoform (H3t) that is required for meiosis initiation.

It should be noted that zebrafish and other anamniote testes differ from mammalian testes in spermatogenic cycle length (45), Sertoli cell organization (54), and histone dynamics (55, 56); and the sperm within do not require acrosomes (57) nor do they require hypercompaction to the extent mammalian sperm do (56). However, zebrafish testes have the same overall processes and developmental stages as mammals, offering a recognizable spermatogenic marker gene landscape. Further, 71% of zebrafish genes are homologous to human and mouse genome (58). Despite the absence of many classical mammalian marker genes (e.g, evolutionarily lost meiotic gatekeeping gene STRA8 and protamine genes PRM1/PRM2; absence of definitive sex-determining chromosomes or gene SRY), the zebrafish testes cell types nevertheless exhibited heterogeneous expression programs exhibiting established markers, enabling identification of the tSNE clusters. The full marker dataset can be found in Dataset S3 (Supplementary Material).

Somatic cells such as Sertoli, Leydig, and immune cells were not present in the datatset. This absence is likely due to the cell loss inherent in the scRNA-seq protocol (59) combined with the bioinformatic requirements defining a tSNE cluster of a cell type. Our input of 5,000 cells per sample was well within the 10x Genomics guidelines of 500 to 10,000. Given the previously published proportion of overall somatic cells in zebrafish testes of ∼1.4% (39), there would likely be less than 30 somatic cells in our sample, which borders the cluster threshold of 30 cells; this number diminishes for specific somatic cell types such as immune cells. Additionally, we have shown scRNA-seq underdetects immune cell gene expression as compared to bulk RNA-seq, all other factors held constant (60). This could contribute to the observed lack of immune cell presence.

Surveying testicular cell type density shows wide variation among, between species

One question that arises when characterizing a scRNA-seq dataset is whether the makeup of the tissue, i.e. the relative number of cells from each cell type in the total sample, follows an expected pattern seen in vivo, lest the analysis be misled by an unrepresentative sample due to dropout (59) or other factors leading to false negative results. In pursuit of a scRNA-seq testes cell type population distribution reference, we supply Table S3 (Supplementary Material), summarized in Fig. 3, noting the relative contribution of SSC/undifferentiated SPG (u-SPG), differentiated SPG (d-SPG), SPC, RS, and ES/sperm to the overall sample of 12 datasets of unselected steady-state spermatogenesis for which this information was explicitly reported, could be derived from supplemental information, or was provided upon request from the authors.

Fig. 3.

Summary graph depicting the cell type contributions (%) to overall cell population from 13 datasets from 11 peer-reviewed scRNA-seq studies in testes. Average % of each cell type observed in published studies, separated by species. Mouse:  N = 5; human:  N = 5; “Other” includes one sheep and one macaque study; zebrafish:  N = 2. Table S1 (Supplementary Material) lists the breakdowns in each dataset.

The 13 datasets span five species (five human datasets (29–33); five mouse (32, 34–37); one sheep (38); one macaque (33); one zebrafish (39)), and vary widely among and between species. The relative contribution of each cell type population to the total dataset is seldom explicitly reported in the main text (of the identified studies, 23%). However, as research using scRNA-seq continues to grow, the authors encourage the scientific community to report this information going forward, particularly considering the large variation in published works on testes of the same species. Our observed distributions are presented as part of Table 1.

Table 1.

Cluster cell population contribution to total cell population in control and TCDD-treated datasets. Relative distribution of clusters in control and treated samples. The % contribution is the number of cells in a given cluster divided by total cells in the population. P-value from two-tailed t test; rounded to nearest 1/1,000th. P < 0.05*, P < 0.01**.

Cluster	Cell type	% contribution from cluster to whole population—control dataset	% contribution from cluster to whole population—TCDD dataset	% change in contribution compared to control	P-value
1	SSC/SPG-A	4.13 ± 4.25	5.12 ± 3.19	+24.11	0.692
2	SPG-B	11.80 ± 1.97	37.63 ± 2.36	+219	0.001 **
3	Preleptotene SPC	7.23 ± 2.95	23.24 ± 4.06	+221.51	0.013 *
4	Leptotene/zygotene SPC	10.60 ± 0.83	20.83 ± 2.50	+96.47	0.005 **
5	Zygotene/pachytene SPC	10.96 ± 0.81	7.42 ± 2.01	−32.28	0.062
6	Pachytene/diplotene SPC	15.78 ± 3.57	1.36 ± 0.93	−91.41	0.012 *
7	Late RS/early ES	9.72 ± 1.90	0 ± 0	−100	0.006 **
8	ES/spermatid	7.31 ± 1.15	0.11 ± 0.04	−98.45	0.003 **
9	Spermatozoa	14.78 ± 4.78	4.22 ± 0.62	−71.44	0.055
10	Spermatozoa	7.71 ± 3.18	0.08 ± 0.12	−99.02	0.041 *

Role of ncRNA and transcripts of unknown function in zebrafish testes

The testes express more noncoding RNA (ncRNA) than other tissues (61), and more ncRNAs are specific to testes (62). Various ncRNA are upregulated at different stages of development. Their main roles are maintaining stemness in SSC and recruiting chromatin-modifying proteins to condense the head during elongation (63). Zebrafish testes display ncRNA accordingly. The top 50 marker genes (FindMarkers function in Seurat) in each cluster contain “predicted” D. rerio ncRNA sequences (RefSeq mRNA accession numbers prefaced with NR or XR; Dataset S2, Supplementary Material), with the exception of cluster 5 (Z/P SPC), spanning 24 unique ncRNAs. Clusters 10 (spermatozoa) and 7 (RS/ES) express the most ncRNA diversity (six and seven transcripts, respectively). A total of one or two long noncoding RNA (lncRNA) transcripts are expressed in half of the 10 clusters’ top 50 DEGs, beginning in SSC/SPG-A through SPG-B, then resurfacing in P/D SPC and spermatozoa. The lncRNA common to the largest share of clusters (six, nine, and 10 (P/D SPC; spermatozoa) is si: ch211-155m12.5, the function of which is yet unknown. Likewise, many top markers in either dataset are either uncharacterized loci or human orthologues. Some are predicted by similarity to share function, but lack experimental evidence (typically Gene Symbols beginning with zgc, zmp, si: dkey(p), si: rp, si: zfos, si: ch, LOC; Dataset S4, Supplementary Material). This information will be valuable to future research investigating unexplored genes in the underexploited zebrafish testes.

One fascinating such transcript is the predicted vicilin-like seed storage protein At2g18540 (D. rerio LOC100332973). Also known as PAP85 in Arabidopsis thaliana, it is the most highly expressed marker gene in cluster 10 (spermatozoa; average LFC 2.51). PAP85 function in testes, or in eukaryotes at large, is unstudied. Vicilins are found in spermatophyte seeds (64), and the embryo and endosperm of developing seeds (65). The function of this specific gene and others with similar homology is unknown, but the function of seed storage proteins such as vicilin is to store amino acids in the ER, and protect them for later harvest (64). Further, an emerging function of PAP85 is RNA virus accumulation (66). PAP85’s facilitation of accumulation of viral RNA has a comparison to the RNA stored during this period in eukaryotic germ cell development as well. In seeking to validate the expression of this unexpected transcript, in situ hybridization probe design was infeasible due to cross-detection of predicted golga6l22 (D. rerio LOC101884685/si: ch73-367j5.3), the second-highest expressing marker in cluster 10 (average LFC 2.46). The mysterious function of PAP85 clearly requires more prokaryotic and eukaryotic research.

TCDD exposure differentially depletes and enriches populations of particular cell types

In keeping with previous findings (16), the TCDD-exposed testes dataset displays strikingly fewer spermatozoa populations (clusters 9 and 10) than controls. Further, every germ cell cluster from diplotene SPC onward is virtually bereft of cells (Fig. 4) In Table 1, we report TCDD-driven population shifts by the change in proportional contribution of each cluster to the entire population; absolute numbers of cells per cluster are reported in Table S2 (Supplementary Material). Interestingly, the depletion seen in late germ cells is not observed for earlier cell types; on the contrary, the SPG-B proportion increased by 219%. Early SPCs up to the leptotene/zygotene SPC stage also saw a ∼1- to 2-fold population increase in the dataset. The shift to depletion occurs around the zygotene/pachytene SPC stage, at which stage the population begins to decrease insignificantly, before the significant decrease beginning in late SPCs. The SSC population was singularly unresponsive to exposure in terms of contribution to the total population. Another way to interpret the changes in cell type contributions to the overall dataset is that early germ cell proportion is increased due to TCDD exposure simply because the sperm contribution is decreased in the sample we analyzed. Our scRNA-seq findings are largely consistent with TCDD-induced lowered sperm counts in animal models (67), and histological examination by our lab (16), wherein TCDD exposure induced a significant reduction in seminiferous tubule area coverage by spermatozoa concomitant with a significant increase in coverage by SPG.

Fig. 4.

tSNE plots of cells from control fish testes (left) and cells from TCDD-exposed testes (right). Clusters are labeled as previously shown.

TCDD exposure disrupts gene expression in multiple pathways during spermatogenesis

TCDD exposure during gonadal differentiation has previously been associated with reduced phenotypic spawning success in adult zebrafish (15). In the present study, multiple molecular pathways related to infertility were enhanced in response to exposure. Following TCDD exposure, 126 genes were significantly upregulated; 950 downregulated (≥ 1 or ≤ 1 average LFC, P < 0.01),  compared to controls. The software Ingenuity Pathway Analysis (IPA; Qiagen Bioinformatics, Redwood City, CA) categorized the biological processes disrupted by exposure (Table 2). This analysis used the “pseudo bulk-seq” dataset, where read counts from all clusters in each condition (exposed or controls) are collapsed (68), which allows a global comparison of pathway disruption following TCDD exposure. This is especially useful when individual clusters lack a sufficient number of strongly up- or downregulated DEGs to perform high-confidence pathway analysis on a cluster-by-cluster basis (cluster-level DEGs are provided in Dataset S5, Supplementary Material). Pseudo bulk-seq analysis serves as a proof-of-concept tool, where affected pathways concurred with previously observed phenotypes in exposed males (16). Exposed testes exhibit upregulated pathways involved with fertility defects. Sperm disorder was the most highly enriched pathway, with top upregulated genes: mov10l1, slc12a2, and prdx4; top downregulated: ube2w, mieg1, dnah1, and cfap91. Other specific pathways of infertility upregulated in the TCDD dataset included teratozoospermia and male germ cell apoptosis. Laterality defects at childbirth such as situs inversus totalis are associated with defective sperm quality (69); these laterality defect pathways were also increased. Likewise, pathways of gametogenesis, sperm movement, and formation of cilia were downregulated in TCDD-exposed cells. A summary of these IPA results and others can be found in Table 2. Dataset S1 (Supplementary Material) contains the full IPA report; Dataset S2 (Supplementary Material) contains the pseudo bulk-seq DEGs.

Table 2.

IPA-generated list of “Diseases and Bio Functions” enriched or suppressed in adult testes following peripubertal TCDD exposure. Differentially expressed genes were required to have an average log2-fold change of ≥ 1 or ≤ 1 and P < 0.01 to be included for analysis. The complete list can be found in Dataset S1 (Supplementary Material).

Rank	Annotation	z-score	P-value
1	Sperm disorder	3.731	8.92E-05
13	Teratozoospermia	2.646	2.37E-05
16	Oligozoospermia	2.538	3.55E-04
21	Laterality defect	2.243	5.45E-08
27	Heterotaxy or ciliopathy	2.236	5.24E-34
31	Situs inversus totalis	2	2.56E-05
50	Apoptosis of germ cells	0.762	9.84E-05
51	Apoptosis of gonadal cells	0.64	6.02E-05
67	Gametogenesis	−0.555	1.05E-05
80	DNA recombination	−1.146	5.05E-06
83	Formation of cilia	−1.465	4.90E-27
84	Homologous recombination	−1.525	1.93E-04
88	Meiosis of germ cells	−1.718	9.05E-05
91	Cell movement of sperm	−2	8.12E-08

Early TCDD exposure is associated with testicular apoptosis over a year after brief exposure

IPA also indicated downregulated pathways of meiosis, and specifically of recombination. Failure of chromosomes to successfully recombine during meiosis can contribute to apoptosis. Taken together, these pathways suggest a partial arrest of germ cells during meiosis as an explanation for previously observed infertility phenotype. IPA returned several pathways related to germ cell death, with downregulation of ubiquitin genes playing an important role. The ubiquitin-driven apoptosis pathways and the importance of ubiquitination in proper sperm maturation led us to hypothesize that TCDD was inducing sperm apoptosis. In healthy sperm, histones are widely ubiquitinated and, thus removed to begin the process of nuclear condensation (70). ube2b specifically is expressed in RSs for this purpose (71), with deficient mice displaying teratozoospermia and infertility (72). tSNEs illustrating the loss of ube2b and multiple other ubiquitin-expressing late germ cells in the TCDD dataset is shown in Figure S3 (Supplementary Material). To test the in silico apoptosis predictions, apoptosis was quantified with a cleaved caspase-3 (cC3) assay (Fig. 5).

Fig. 5.

TCDD exposure is associated with testicular apoptosis. Representative images of H&E stains (a) and (b) and cleaved caspase-3 (cC3) immunolabeling (c) and (d) of control and TCDD sections, respectively, at 40x magnification. Quantification of % apoptotic total cells (P = 2.88E-05) (e), SPG (P = 0.2338) (f), SPCs (P = 2.17e-04) (g), spermatids (P = 0.01246) (h), and spermatozoa (P = 1.85e-04) (i). N = 8 unique regions from 60x images of three fish (controls); nine unique regions from images of three fish (TCDD). * = P < 0.05; ** = P < 0.01; and *** = P < 0.001 (student's t test).

Similar to previous histological quantifications using area coverage of each cell type (μm²) (16), there was no change in the number of SPC or spermatids in TCDD-exposed testes sections. In the present histological study, the exposure also resulted in no significant change in SPG or spermatozoa number (Table S4, Supplementary Material). Interestingly, while there was no significant association between exposure and total cell population changes via H&E staining, the proportion of cC3-stained cells significantly increased following exposure in SPC, spermatids and spermatozoa. Sperm are still produced in TCDD-exposed fish, but are terminating in the seminiferous tubule. At first glance, these results seem to conflict with our in situ scRNA-seq results that clearly depict an absence of both spermatid and spermatozoa clusters following exposure. However, inherent in the scRNA-seq protocol is tissue dissociation to obtain single cells. This disruptive process may have been a sufficient insult added to the TCDD-induced precarity of the late germ cells to cause them to apoptose before reaching sequencing, which would explain their presence in histological samples and the control scRNA-seq dataset, yet absence in sequencing data and subsequent clustering. Leal et al. (45) note that apoptotic germ cells are rarely detected in adult zebrafish under normal conditions. Testicular apoptosis in TCDD-exposed organisms is understudied, but existing evidence is in conflict, with some findings of increased apoptosis following adult rodent exposure (73) and tissue incubation (74), but null results in offspring following gestational exposure (75).

Summary

Spermatogenesis is a conserved process across species, yet the zebrafish testes’ distinct morphology and unique gene expression warrant an in-depth characterization of the transcriptomic landscape in each cell type. ScRNA-seq successfully differentiated nine cell types in consecutive stages of development. This dataset serves as a testicular atlas for the rapidly growing use of this model organism. Zebrafish are renowned in particular for their contribution to developmental toxicology. We show that exposure to the EDC TCDD during sexual differentiation and maturation bestows permanent effects on the cellular makeup of the zebrafish testes. Spermatids and spermatozoa apoptose in adult male fish that were exposed for only 1 hour at 3 and 7 wpf, the windows during which zebrafish gonads commit to a binary fate. This period is heavily influenced by environmental conditions. The indifference of SPG to TCDD-induced apoptosis suggests that rather than total tissue toxicity, TCDD exposure confers a cryptozoospermia-like phenotype, where germ cells progress normally early in the cycle of spermatogenesis, but arrest or apoptose during spermatozoal maturation. DEG and pathway analysis suggest partial arrest during meiosis may contribute to the permanent developmental postponement. This work substantially builds on previous work histologically showing a decrease in spermatozoa, and phenotypic male spawning failure by revealing significant testicular disruption long after exposure cessation.

Materials and Methods

Fish husbandry

As described in Meyer et al. (40), Zebrafish (AB strain) were maintained on a 14:10 hour light/dark cycle (76) in reverse osmosis water buffered with Instant Ocean salts (60 mg/L, Aquarium Systems, Mentor, OH), with temperatures maintained at 27 to 30°C. Fish were fed twice daily. Adult fish were raised on a recirculating system at a maximum density of five fish per liter. Fish were spawned to produce embryos and following the first exposure (3 wpf) were raised in beakers with daily water changes of 40% to 60% at a density of 5 fish per 400 mL beaker between 3 and 6 weeks, and five fish per 800 mL beaker between 6 and 9 wpf. Animal use protocols were approved by the Institutional Animal Care and Use Committees at Wayne State University and the University of Wisconsin-Madison, according to the National Institutes of Health Guide to the Care and Use of Laboratory Animals.

TCDD exposure

As described in Baker et al. (15), TCDD (> 99% purity; Chemsyn) was used as a 0.4 ng/mL stock solution in DMSO. Zebrafish were exposed as previously described (15) at 3 wpf and again at 7 wpf to water-borne TCDD (50 pg/mL; 0.155 nM) or vehicle (0.1% DMSO) for 1 hour each time in small glass beakers with gentle rocking. The number of fish per volume of dosing solution was 1 fish/mL at 3 wpf, and due to growth between 3 and 7 wpf, 1 fish/2mL at 7 wpf.

Testes isolation and enzymatic dissociation of testes

Adult (1.5-year-old (+/−1 month)) male zebrafish were euthanized in tricaine methanesulfonate (1.67 mg/mL; Fisher Scientific, Waltham, MA) for 10 minutes. Testes were dissected, and excess tissue removed from testes in 4°C PBS (Gibco, Waltham, MA). Testes were minced in 4°C PBS, then centrifuged at RT for 5 minutes at 500 g. PBS was removed, and 100 uL of digestion media (100 uL Leibovitz's L-15 medium; MilliporeSigma, Burlington, MA), 1 uL bovine serum albumin (New England BioLabs, Ipswich, MA), 1 uL DNAseI (Zymo Research, Irvine, CA), and 1 mg collagenase Type II (Worthington Biochemical Corporation, Lakewood, NJ)) were added. Tissue was shaken at 280 rpm for 1.5 hours at 27°C, with manual disruption via wide-bore pipetting every 15 minutes. Cells were centrifuged at RT for 5 minutes at 500 g, digestion media aspirated, and cells resuspended in RT PBS and placed on ice. Dead cells were removed with a Dead Cell Removal Kit (Miltenyi Biotec, Bergisch Gladbach, Germany). Cell viability of 90% was determined using the BioVision Live/Dead Cell Viability Assay Kit (BioVision Inc., Milpitas, CA), according to manufacturer's instructions. Approximately 2 million cells were immediately harvested.

Library preparation and sequencing

Single cell libraries were constructed using the 10x Chromium Controller v2 chemistry following the Chromium Single Cell 3' region protocol (CG00052; 10x Genomics, Pleasanton, CA) (77), with an input of 5,000 cells per sample. Libraries were sequenced on a NovaSeq 6000 (Illumina, San Diego, CA) at a depth of 50,000 reads/cell.

Data processing

Cell Ranger v6.0.1 (10x Genomics, Pleasanton, CA) was used to align sequencing reads to the zebrafish reference genome (dR10), which was constructed using the mkref command (78). Count data was imported to Seurat (version 4.0.4) for QC filtering, clustering, dimensionality reduction, visualization, and differential gene expression (79, 80). Each sample was filtered to cells containing at least 500 features with clusters requiring a minimum of 30 cells. Samples were merged prior to normalization and clustering (resolution 0.3). Differentially expressed genes between conditions for each cluster were identified using the “FindMarkers” function.

IPA

The functional pathways in each comparison were generated through the use of IPA (Qiagen Inc., https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis). Significant DEGs included in pathway analysis required a log2-fold change of ≥ 1 or ≤ −1, and P-value < 0.01.

Apoptosis assay

A separate cohort of 1-year-old adult male zebrafish under the same exposure scheme were euthanized in tricaine methanesulfonate, fixed in 10% Zn formalin, decalcified with Cal-ExII, bisected along the sagittal plane, dehydrated in a graded series of ethanol, cleared in xylene, and paraffin-embedded as described in Baker et al. (81). Immunohistochemistry was performed by the Wayne State University Biobank and Correlative Sciences Core. Formalin-fixed paraffin-embedded sections of bisected zebrafish were dewaxed and rehydrated in a xylene–ethanol–water series. Endogenous peroxides were removed by a methanol/1.2% hydrogen peroxide incubation at room temperature for 25 minutes. HIER antigen retrieval was done with a pH6 citrate buffer and the BIOCARE Decloaking Chamber (Concord, CA). A 40-minute blocking step with SuperBlock Blocking Buffer (Thermo Scientific, Waltham, MA) was performed prior to adding the primary antibody. A 1:100 dilution of Cleaved Caspase 3 (cC3) antibody (9664S; Cell Signaling, Danvers, MA) was used overnight at 4°C. Detection was obtained using GBI Labs (Bothell, WA) DAB Chromogen Kit and counterstained with Mayer's hematoxylin. Sections were then dehydrated through a series of ethanol to xylene washes and coverslipped with Permount (Fisher Scientific). The authors analyzed cC3 labeling to determine presence and/or extent of apoptosis (control fish: N = 3; TCDD-exposed: N = 3). We obtained up to three distinct images from replicate testes slides at 40x magnification, and manually quantified a quadrant of each image (control = 8 quadrants, TCDD = 9 quadrants). Significance of the % apoptotic cells per cell type between controls and TCDD images was measured via student's two-tailed t test.

Click here for additional data file.