Characterization of 2,2',4,4'-tetrabromodiphenyl ether (BDE47)-induced testicular toxicity via single-cell RNA-sequencing.

Background: The growing male reproductive diseases have been linked to higher exposure to certain environmental compounds such as 2,2',4,4'-tetrabromodiphenyl ether (BDE47) that are widely distributed in the food chain. However, the specific underlying molecular mechanisms for BDE47-induced male reproductive toxicity are not completely understood.
Methods: Here, for the first time, advanced single-cell RNA sequencing (ScRNA-seq) was employed to dissect BDE47-induced prepubertal testicular toxicity in mice from a pool of 76 859 cells.
Results: Our ScRNA-seq results revealed shared and heterogeneous information of differentially expressed genes, signaling pathways, transcription factors, and ligands-receptors in major testicular cell types in mice upon BDE47 treatment. Apart from disruption of hormone homeostasis, BDE47 was discovered to downregulate multiple previously unappreciated pathways such as double-strand break repair and cytokinesis pathways, indicative of their potential roles involved in BDE47-induced testicular injury. Interestingly, transcription factors analysis of ScRNA-seq results revealed that Kdm5b (lysine-specific demethylase 5B), a key transcription factor required for spermatogenesis, was downregulated in all germ cells as well as in Sertoli and telocyte cells in BDE47-treated testes of mice, suggesting its contribution to BDE47-induced impairment of spermatogenesis. Conclusions: Overall, for the first time, we established the molecular cell atlas of mice testes to define BDE47-induced prepubertal testicular toxicity using the ScRNA-seq approach, providing novel insight into our understanding of the underlying mechanisms and pathways involved in BDE47-associated testicular injury at a single-cell resolution. Our results can serve as an important resource to further dissect the potential roles of BDE47, and other relevant endocrine-disrupting chemicals, in inducing male reproductive toxicity.

Introduction

Polybrominated diphenyl ethers (PBDEs) are widely used as brominated flame retardants in multiple highly flammable consumer products, including electrical goods, furniture, electronic goods, and textiles. As environmental pollutants, PBDEs are widely present in the environment such as in the atmosphere, soil, and water, and they can be detected in multiple animals and human fluids.[3] Increasing public concern regarding the toxicity of PBDEs has been aroused. Particularly, PBDEs show developmental toxicity, neurotoxicity, liver toxicity, and endocrine disruption-mediated reproductive toxicity.[4] Therefore, PBDEs were restricted under the Stockholm Convention in 2009.[5] However, PBDEs are still being released from consumer products as a result of e-waste recycling, or via the debromination of higher PBDE congeners, and recent studies have shown that humans are still being exposed to PBDEs.[1,3], As such, exposure to PDBEs is likely to persist,[1] and the potential health risks are still a major concern.

2,2',4,4'-Tetrabromodiphenyl ether (BDE47) is a predominant PBDE congener that is distributed in the environment, in animals, and in human fluids. BDE47 exhibits relatively high toxicity among other major PBDE congeners.[8,9] In animals, BDE47 shows developmental, reproductive, and neurological toxicities.[4,10] Exposure to BDE47 is negatively correlated with human reproduction system disorders such as infertility, poor sperm quality, and testicular cancer, possibly via disruption of thyroid homeostasis and other mechanisms.[4,11,12] The effects of BDE47 on disruption of spermatogenesis, impairment of sperm motility, and reduction of testes size in animal models have been well-documented.[6], Prepuberty is a key period for the development of the male reproductive system. Adversely affected at this stage by a variety of endocrine-disrupting environmental chemicals (e.g. BDE47, and BPA (bisphenol A)) may affect male fertility even in later adult life.[18] Notably, accumulating evidence has suggested that infants and young children show a much higher amount of PBDEs exposure than adults. A variety of mechanisms such as induction of apoptosis, increased reactive oxygen species, and disruption of homeostasis hormones have been proposed for BDE47-induced testicular injuries.[6],,[23] However, previous studies mainly focus on certain cell types within the testis, and the underlying mechanism of BDE47-induced testicular injury is still unclear. Thus, comprehensively elucidating the underlying mechanisms of BDE47-induced testicular injury is an interesting topic.

Single-cell RNA-sequencing (ScRNA-seq) has provided an unprecedented solution to characterize complex biological processes by a comprehensive and unbiased analysis of gene expression at the cellular level. Though previous studies successfully identified major different regulated genes and pathways associated with BDE47-induced testicular injury,[24] these studies mainly used aggregated RNA from multiple cell populations[25] that cannot resolve BDE47-related gene expression across different cell types in testicular tissue. Thus, to characterize how different cell types are affected by BDE47 in the testis and to clarify if different cells in the testis follow a similar blueprint or whether certain cell types display unique transcriptional changes upon BDE47 treatment are urgently needed. Here, male C57BL/6J mice were exposed to BDE47 during the prepubertal period, and ScRNA-seq was performed in testicular tissues. We aimed to comprehensively understand the underlying mechanism of reproductive toxicity of prepubertal BDE47 exposure in a mouse model. Our results provide novel insight into our understanding of the underlying mechanisms of BDE47-associated male reproductive toxicity, which may also serve as an important resource for future dissection of the role of BDE47 and other relevant endocrine-disrupting chemicals in inducing testicular toxicity.

Methods

Animals and drug treatments

C57BL/6J mice were obtained from Charles River, China. The mice were housed in a pathogen-free environment with 12 h light and dark cycles. All animal experiments were approved and followed the guidelines and regulations of our institution. For drug treatment, mice were given BDE47 at a concentration of 10 mg/kg/day by intraperitoneal injection once daily for three weeks. Dosage was selected based on environmental relevance, human concentrations, and previous toxicity studies.[4,8,26,27]

Hematoxylin and eosin staining

Following fixation, the testicular tissues from untreated and BDE47-treated mice were embedded and then serially sliced into 4 μm sections followed by hematoxylin and eosin (H&E) staining. Images were observed and captured by a microscope.

Immunostaining and western blotting

For immunostaining, after fixation, the testicular tissues were embedded and then serially sliced into small sections and subjected to immunostaining with the indicated primary antibodies, followed by incubation with 3,3-N-diaminobenzidine tertrahydrochloride or fluorescent secondary antibodies as described previously.[28,29] Western blotting was performed as previously described.[30,31] Briefly, tissue sections were lysed in radioimmunoprecipitation assay (RIPA) buffer and extracted proteins were resolved by SDS-PAGE separation and then incubated with the indicated primary and secondary antibodies.

Reverse transcriptase-quantitative PCR (RT-qPCR)

RT-qPCR was used as described previously.[32] 10 mg of testicular tissue were used for RNA isolation using the RNAprep Pure Tissue Kit (CAT# DP431) according to the product manual. A total of 1 µg of RNA was used for cDNA synthesis using the PrimeScript™ RT reagent Kit with gDNA Eraser (Cat. # RR047A). cDNA was analyzed by RT-qPCR using primers specific for the mouse testicular genes Kdm5b, Hsd3b6, Star, Asah1, Dhcr24, Cyb5r3, and reference gene Gapdh (glyceraldehyde-3-phosphate dehydrogenase) (primer sequences are shown in Supplementary Table 6). PCR was carried out using the TB Green® Premix Ex Taq™ II (Tli RNaseH Plus) (Cat. # RR820A).

For each gene, relative expression levels were calculated by the comparative CT method (StepOne software v2.3, from Applied Biosystems) obtained by qPCR assays of cDNA samples. Finally, Kdm5b, Hsd3b6, Star, Asah1, Dhcr24, and Cyb5r3 expression levels were normalized to that of Gapdh.

Single-cell suspension preparation

Single-cell suspension was prepared according to reported protocols.[33] Briefly, testes from untreated and BDE47-treated mice were excised. Seminiferous tubules were digested with collagenase IA, DNase I, and trypsin. Single cells were then filtered through a 100 µm strainer, washed with phosphate-buffered saline (PBS), and resuspended in MACS buffer, and single-cell suspensions were used directly for Droplet-Sequencing (Drop-Seq).

Single-cell libraries establishment and sequencing

Drop-seq libraries were established according to a previously described protocol.[34] Briefly, single cell suspension and beads were mixed, and individual droplets were harvested for reverse transcription, and then cDNA amplification. A 3ʹ gene expression library was then established with a 10x Genomics Chromium Single Cell system with v3 chemistry. The NovaSeq 6000 (Illumina) was used for sequencing by Novogene (China).

Read alignment and gene expression quantification

Reads to the mm10 reference genome were aligned and obtained from Ensembl. The GENCODEvM20 was employed for gene annotations. 10x Genomics Cell Ranger 6.1.1 was used for gene alignment, demultiplexing unique molecular identifiers (UMIs), and identification of cell barcode to generate a cells x genes count matrix.

Quality control

A single-cell expression matrix was calculated using the Seurat package (version 3.6.1) for filtering, normalization of data, reduction of dimensionality, clustering, and analysis of differentially expressed genes (DEGs). We retained the high-quality single-cell data that met the following three criteria: (1) cells with fewer than 300 genes detected or a mitochondrial gene ratio of greater than 20% were excluded; (2) the number of detected genes was greater than 500; and (3) the number of detected UMI was greater than 3000. The DoubletFinder package (version 2.0.2) was used to detect doublets.[35] The mean-variance-normalized bimodality coefficient of individual samples was calculated to detect the neighborhood size (pK), and the number of artificial doublets was set to a value of 0.25. Doublets were then determined. Variable genes were selected based on the default parameters. Variable genes were projected onto a low-dimensional subspace by using canonical correlation analysis across different samples to correct batch effects.

Identification of cell types and cell-type-specific markers

We normalized the cells x genes matrix by the total number of UMIs in each cell and scaled by 106 to yield counts per million. We transformed a natural log of this matrix after the addition of a pseudocount of 1 to avoid undefined values. To enable unsupervised clustering and cell-type identification, we performed dimensionality reduction using principal component analysis (PCA) on the combined set of samples after the selection of high variable genes. Once embedded in this PCA space, we constructed a nearest neighbor graph identifying the k = 45 nearest neighbors for each cell. Cells were visualized using the projection for dimension reduction (UMAP) algorithm and were clustered at an appropriate resolution. Cell types were assigned to each cluster using the abundance of known marker genes. Through the above pipeline, we processed the ScRNA-seq data from six testes of mice.

Differential gene expression analysis

Differential gene expression analysis for each cell type between untreated and BDE47-treated mice was performed with the Wilcoxon rank-sum test. Before conducting the differential gene expression analysis, we excluded the cell types that were missing or had fewer than 3 cells in the comparison groups. Differential expression analysis between the groups was performed to generate |Log2Fold Change (FC)| > 0.1, P value < 0.05.

Gene ontology enrichment analysis

Gene ontology (GO) analysis of DEGs was performed by ClusterProfiler[36] and the results were visualized using the ggplot2 R package (https://ggplot2.tidyverse.org/) (version 3.2.1). Representative terms selected from the top-ranked GO terms or pathways (P < 0.05) were displayed.

Transcriptional regulatory network analysis

To analyze the transcriptional regulatory network, SCENIC workflow (version 1.1.2.2) was used according to default parameters with the mm10 database from RcisTarget (version 1.6.0).[37] The SCENIC workflow consists of three steps. In the first step, sets of genes that are co-expressed with transcription factors (TFs) are identified using GENIE3. Being based only on co-expression, these modules may include many false positives and indirect targets. To identify putative direct-binding targets, each co-expression module is analyzed using cis-regulatory motif analyses using RcisTarget. Only modules with significant motif enrichment of the correct upstream regulator are retained, and pruned to remove indirect target genes without motif support. For the BDE47-related transcriptional regulatory network, only BDE47-related DEGs were used as input, and all selected cell types were calculated together. The transcriptional regulatory network analysis was performed by ggraph R package.

Cell-cell communication analysis

Cell-cell communication analysis was conducted from the data using the CellPhoneDB[38] software (version 1.1.0) (www.cellphonedb.org) and CellChat (version 1.1.0) as described previously.[39] Only receptors and ligands expressed in >10% of cells from either untreated or BDE47-treated samples were further detected, while cell-cell communication was considered nonexistent if the ligand or the receptor was unmeasurable. Averaged expression of each ligand-receptor pair was analyzed between various cell types and only those with a P value < 0.05 were employed for further prediction of cell-cell communication between any two cell types. For each receptor-ligand pair in each pairwise comparison between two cell types, this generates a null distribution. By calculating the proportion of the means that are as high or higher than the actual mean, we obtained a P value for the likelihood of cell-type specificity of a given receptor-ligand complex. Only those with a P value < 0.01 were used for the prediction of cell-cell communication between any two cell types.

Data availability and statistical analyses

All data were statistically analyzed by a two-tailed t-test, assuming equal variance with R package. A P value < 0.05 was considered statistically significant. The data reported in this paper have been deposited in the OMIX, China National Center for Bioinformation /Beijing Institute of Genomics, Chinese Academy of Sciences (https://ngdc.cncb.ac.cn/omix: accession No. OMIX792).

Results

BDE47 exposure induces testicular injury in mice

To determine the testicular toxicity of BDE47, a prevalent PBDE congener, during the prepubertal exposure, C57BL/6J mice were administrated with BDE47 at a concentration of 10 mg/kg/day via intraperitoneal injection starting at postnatal day 22 for 3 weeks. The experimental design is shown in the schematic model (Fig. 1A). H&E staining results showed that prepubertal BDE47 exposure induced testicular toxicity in mice as reflected by mice testis degeneration, which was accompanied by reduced seminiferous epithelium thinning, wrinkling on the boundaries of seminiferous tubules, and intensified atrophy (Fig.1B). The role of BDE47-induced testicular injury was further confirmed by reducing the testes index (Fig.1C). Our results are consistent with previous findings that BDE47 impairs the spermatogenetic process in multiple animal models, and prepubertal exposure to BDE47 results in male reproductive toxicity.[15,17,23,25]

Figure 1.

Single-cell transcriptome analysis for identification of cell types in the testis upon BDE47 treatment. (A) Schematic model showing the experimental design of our study. C57BL/6J (day 22) mice were given BDE47 for 3 weeks, and testes were harvested for single-cell RNA sequencing and bioinformatics analysis. (B) Representative H&E staining images from five mice in each group showing BDE47-induced testicular injury. Wrinkling on the boundaries of seminiferous tubules is indicated with a black arrow; intensified atrophy and loss of seminiferous tubules are shown with a blue arrow. (C) BDE47 reduces testis/body weight, *P< 0.05. (D) UMAP plot of cells from untreated and BDE47-treated testicular tissues showing 10 major colored clusters of cell types. (E) Violin plot showing average gene numbers in different cell types. (F) Plot of canonical cell-type markers in different cell types.

Single-cell transcriptome profiling identifies cell types in the testis upon BDE47 treatment

To comprehensively understand BDE47-induced testicular injury, we performed ScRNA-seq to characterize the underlying mechanism of how prepubertal BDE47 exposure may contribute to testicular toxicity in mice. Single-cell suspensions from mice testes were obtained and loaded into a droplet-forming microfluidic device. To define different cell types, quality control of the sequencing data was first processed with the DoubletFinder[35] and Seurat R packages. After normalization, corrected batch effect, and clustering, each cell type was annotated in the testes from 76 859 qualified cells based on the expression of canonical cell-type-specific markers as described previously.[40] We obtained all major germ cell and somatic cell populations covering the full testis developmental spectrum, including spermatogonia (SPG), meiotic spermatocytes, pachytene, acrosomal, post-meiotic haploid round spermatids, elongating spermatids, innate lymph, telocytes, Leydig, and Sertoli cells (Fig. 1D, Fig. S1B, see online supplementary material). Average gene counts and gene numbers detected in individual cell types were largely consistent with previous reports (Fig.1E and Fig. S1A).[40] Detailed marker genes for testis are given in Fig. 1F. Overall, we generated a single-cell atlas in mouse testis upon prepubertal BDE47 treatment, providing a cellular roadmap and resource for further dissection of BDE47-associated testicular toxicity in mice.

Characterization of BDE47-associated heterogeneous changes in gene expression

Our ScRNA-seq results were generated from both untreated and BDE47-treated mice showing each cell type with comparable UMI numbers and genes detected in testis. Our results showed that the cellular composition of multiple cell types was largely consistent between untreated and BDE47-treated testes, and cell identity was also largely preserved in BDE47-treated testis as reflected by unbiased clustering where all clusters represent cells of all mice (Fig. S1B and C). However, BDE47 treatment induced a decrease in SPG and increased innate lymph cells (Fig. S1C), indicating that the changes in these two cell types may play a major role in BDE47-induced testicular injury. Moreover, we examined the breadth of transcriptional changes that happen in the mouse testis treated with BDE47 by performing DEG analysis between all 10 major cell types in testis. Our results revealed that 3652 genes were significantly upregulated or downregulated upon BDE47 treatment in at least one cell type of testis (P value < 0.05, Fig. 2A, Table S1, see online supplementary material). Notably, of those, 2443 genes showed the same directionality regardless of the cell type (1060 upregulated and 1383 downregulated), while the direction of change in the expression of 1209 genes was different across the major identified cell populations (Fig.2A, Table S1). Among all the cell types, pachytene showed the greatest numbers of differentially expressed genes in response to BDE47 (Fig.2A).

Figure 2.

BDE47 induces the heterogeneous changes in gene expression in different cell types in the testis. (A) Strip chart showing the changes in all detected genes (dots) across ten cell types. Genes with colored dots were significantly upregulated or downregulated upon BDE47 treatment. (B) Top five cell-type-specific upregulated and downregulated genes in different cell types upon BDE47 treatment. (C) Heatmap showing bidirectional DEGs in different cell types upon BDE47 treatment. (D) Heatmap showing significantly upregulated genes in multiple cell types upon BDE47 treatment. Genes in red were significantly upregulated and genes in blue were significantly downregulated upon BDE47 treatment. Black genes were not significantly changed upon treatment ( False Discovery Rate (FDR)  <  0.05 and  Fold Change (FC)  > 10%).

Additionally, we found that multiple genes were dramatically downregulated in a cell-type-specific manner and the top five upregulated and downregulated genes in different cell types are listed in Fig.2B. From this list, it can be found that a variety of genes that play a major role in regulating spermatogenesis were shown to be downregulated. For instance, Atr (ataxia telangiectasia and Rad3 related) is significantly downregulated in pachytene cells upon BDE47 treatment. ATR is a serine/threonine-protein kinase that is critical for sensing DNA damage to preserve genome stability. It is also a key regulator of male mouse meiosis,[41] whose deficiency induces mid-pachytene germ cell elimination.[41]Akap9 (A-kinase anchoring protein 9), a key regulator of Sertoli cell maturation and spermatogenesis,[42] was found to be dramatically decreased in acrosomal cells in response to BDE47. The downregulation of Hsd3b6 (3 beta-hydroxysteroid dehydrogenase), a gene that is key for the production of progesterone,[43] was found in Leydig cells, which is consistent with previous findings that BDE47 inhibits progesterone synthesis.[44] These results indicate that multiple downregulated genes in different cell types may be involved in BDE47-induced prepubertal testicular injury. Moreover, the genes that are most significantly upregulated in a cell-type-specific manner are shown in Fig.2B. For instance, Dnajb8, which encodes a protein that is critical for regulating chaperone activity,[45] was only significantly increased in pachytene. As another example, Cdc34, which encodes the protein ubiquitin-conjugating enzyme that plays an important role in regulating the cell cycle of G1 regulators and DNA replication, was found to be upregulated only in meiotic division cells. These results suggest that BDE47-induced prepubertal injury may be via upregulating the expression of multiple genes in a cell-type-specific manner.

Our ScRNA-seq results also identified bidirectional regulated genes in different cell types upon BDE47 treatment (Fig.2C). For instance, Smc6, which encodes a protein involved in the structural maintenance of chromosomes protein complex that is crucial in maintaining genome stability and spermatogenesis,[46] was revealed to be significantly decreased in pachytene, meiotic division, Sertoli, and telocyte cells, but dramatically increased in Leydig cells. As another example, Fndc3a, which encodes fibronectin type-III domain-containing protein 3A that plays an important role in mediating the adhesion of spermatids and Sertoli cells during spermatogenesis,[47] was significantly decreased in pachytene, meiotic division, acrosomal, and Sertoli cells, but was significantly increased in Leydig cells. These results indicate that certain changes may be masked by using traditional bulk RNA sequencing methods, and ScRNA-seq results provide deep information for understanding BDE47-induced changes in gene expression and facilitating the understanding of the underlying mechanism of BDE47-induced injury in testis.

In addition, our ScRNA-seq results showed that multiple genes were significantly increased (Fig.2D) or decreased (Fig. S2, see online supplementary material) in several different cell types upon BDE47 treatment. For instance, Hspb1 was significantly increased in SPG, acrosomal, round spermatids (STids), Sertoli, and telocyte cells. Hspb1, which encodes heat shock factor binding protein 1, is a negative regulator of the heat shock response by inhibiting HSF1’s DNA-binding activity.[48] Since HSF1’s transcriptional activity is crucial for spermatogenesis and male fertility,[49] it has been postulated that the upregulated expression of Hspb1 in several cell types in testis upon BDE47 treatment may be involved in BDE47-mediated testicular injury. Notably, multiple genes dramatically downregulated in different cell types upon BDE47 treatment were also identified (Fig. S2). For instance, several genes encoding deubiquitinating enzymes including Usp1, Usp8, Usp16, and Usp47 were significantly downregulated in several cell types. These deubiquitinating enzymes are critical for spermatogenesis. For instance, ubiquitin-specific peptidase 1 (USP1) is a protein involved in regulating DNA repair and its deficiency caused male infertility [27]. Several cilia- and flagella-associated proteins encoding genes including Cfap43, Cfap69, and Cfap97 were also significantly decreased in several cell types in the testis. These genes are crucial for spermatogenesis, and it has been shown that mutations in Cfap43 induces flagellum defects and male infertility.[50] Similarly, CFAP69 is also crucial for flagellum assembly/stability, and spermatogenesis.[51] These results suggest that BDE47-induced prepubertal injury may be via regulating the expression of multiple genes involved in spermatogenesis in different cell types.

Characterization of BDE47-induced testicular injury in somatic cells

Spermatogenesis is a conserved complex program that is driven by interactions of multiple germ cells and somatic cells. This process involves mitotic, meiotic, and differentiation events in the seminiferous tubules within the testis. To better understand BDE47-induced testicular toxicity, we then investigated BDE47-induced changes in cellular pathways and processes in somatic cells. By analyzing ScRNA-seq data, we dissected the shared and cell type-specific DEGs in different somatic cells in the testis and revealed that most DEGs were cell type-specific (Fig. 3A). Next, we investigated the changes in cellular pathways and processes upon BDE47 treatment by performing gene functional enrichment analysis. Our results showed that the most common up-regulated pathways in somatic cells in response to BDE47 treatment included ribosome biogenesis, purine nucleotide metabolic process, lactate metabolic process, and energy deprivation by oxidation of organic compounds (Fig. 3B, Table S2, see online supplementary material). GO analysis showed that BDE47-associated down-regulated pathways in somatic cells were steroid hormone-mediated signaling pathways, regulation of cellular amide metabolic process, histone modification, and double-strand break repair (Fig. 3C, Table S3, see online supplementary material). These results indicate that BDE47 may induce male reproductive toxicity via regulating these pathways.

Figure 3.

Characterization of BDE47-induced changes in cellular pathways and TFs in somatic cells. (A) Venn diagram showing the distribution of DEGs across all the somatic cell types upon BDE47 treatment. Color intensity indicates the changed levels (P value < 0.05 and FC > 10%). (B and C) Representative top GO terms of up-regulated (B) and down-regulated (C) genes in different somatic cell types in the testis upon BDE47 treatment. Circle size indicates the gene count. (D) Violin plots showing gene set scores of steroid metabolic process. **P < 0.01; ***P<0.001; ****P < 0.0001. (E) BDE47 decreases testosterone levels in testicular tissues of BDE47-treated mice as reflected by ELISA assay. **P<0.01. (F) BDE47 decreases mRNA levels of key genes (Asah1, Cyb5r3, Dhrc24, Hsd3b6, and Star) that are involved in testosterone synthesis compared with the Ctrl groupas as reflected by qPCR analysis. **P<0.01. (G and H) Violin plots showing gene set scores of cytokinesis steroid metabolic process (G) and ER stress (H) in different somatic cells upon BDE47 treatment. **P < 0.01; ***P<0.001; ****P < 0.0001. (I) Immunofluorescence showing increased caspase 3 upon BDE47 treatment. **P < 0.01; ***P<0.001; ****P < 0.0001. (J) Dot plot showing significantly changed TFs in different somatic cells. (K) Network plot showing the differentially expressed TFs and corresponding differentially expressed target genes in different somatic cells. (L) Bar plot showing the siginifacant enrichment of GO terms of TF CREB1 target genes, Y axis represents the Log 10(P value). Data were presented as mean ± SEM of three control and three BDE47-treated   testes.

Interestingly, gene set scores results revealed that BDE47 decreased the expression of genes associated with steroid metabolic processes in Leydig and Sertoli cells (Fig.3D), which is consistent with previous findings that BDE47 inhibits progesterone synthesis[44]. Interestingly, BDE47 reduced testicular testosterone levels as reflected by our enzyme-linked immunosorbent assay (ELISA) results (Fig. 3E). Consistent with ScRNA-seq analysis results, the mRNA levels of several key genes such as Hsd3b6, Star, Asah1, Dhcr24, and Cyb5r3 were significantly reduced in BDE47-treated mice (Fig.3F), which may explain how BDE47 inhibits testosterone synthesis. Gene set scores also showed that BDE47 induced the downregulation of genes associated with cytokinesis in somatic cells within testis (Fig. 3G). Cytokinesis is a conserved process for separating the cytoplasm and the genome into two daughter cells during cell division, and deregulation of this process is linked to multiple diseases and affects spermatogenesis.[52] Moreover, BDE47 specifically increased gene set scores associated with ER stress in Leydig cells (Fig.3H). Because chronic ER stress can induce apoptosis, we then found that BDE-47 indeed increased caspase 3 levels (Fig.3I), as reflected by immunostaining analysis and western blotting analysis (Fig. S4A, see online supplementary material). TUNEL (terminal deoxynucleotidyl transferase dUTP nick end labeling) staining further confirmed that BDE-47 induced apoptosis (Fig. S4B). These results highlight the critical roles of steroid metabolic, cytokinesis, and ER stress pathways in BDE47-induced testicular injury.

To better understand BDE47-induced testicular toxicity, we characterized BDE47-associated changes in key TFs using gene module co-expression analysis and known DNA-binding sites analysis from the RcisTarget database, and identified multiple shared or cell-type-specific TFs that were changed (Fig. S3, see online supplementary material). In somatic cells, we discovered that Zfp10, Kdm5b, and Creb1 were major changed TFs after BDE47 treatment (Fig.3J). We further analyzed the differentially regulated genes of these TFs in somatic cells of testis after BDE47 treatment (Fig.3K). CREB TFs are critical for regulating spermatogenesis[53] and GO analysis showed that Creb1 target genes were involved in multiple key processes of spermatogenesis such as response to growth hormone, regulation of stem cell differentiation, cell junction organization, and cell cycle checkpoint (Fig. 3L). These results indicate that BDE47 impairs spermatogenesis possibly via regulating several critical TFs in somatic cells of the testis.

Characterization of BDE47-induced testicular injury in germ cells

To further dissect BDE47-associated changes in cellular pathways in different germ cells, we first analyzed DEGs in germ cells and found multiple shared and cell type-specific DEGs upon BDE47 treatment (Fig. 4A). GO analysis showed that BDE47-associated top upregulated pathways include macrophage chemotaxis, ketone biosynthetic processes, endoplasmic reticulum unfolded protein response, and ATP metabolic processes (Fig.4B, Table S4, see online supplementary material). In contrast, GO analysis showed that the top downregulated pathways assiciated with BDE47 included multiple cellular processes regulating spermatogenesis such as histone H2A ubiquitination, microtubule cytoskeleton organization involved in mitosis, chromatin remodeling, centriole replication, and cellular response to steroid stimulus (Fig. 4C, Table S5, see online supplementary material). These results highlight the critical functions of the DEG-associated pathways mentioned above in response to BDE47-associated testicular injury in germ cells. Gene set scores showed that double-strand break repair was significantly decreased in response to BDE47 treatment (Fig.4D). This result was further supported by immunostaining showing an increase in γ-H2AX (Fig.4E), a marker for DNA damage.[54] Our results also showed decreased gene set scores associated with chromatin remodeling (Fig.4F). Collectively, these results suggest critical roles of DNA damage response and chromatin remodeling activity in BDE47-mediated testicular injury in germ cells.

Figure 4.

Characterization of BDE47-induced changes in cellular pathways and TFs in germ cells. (A) Venn diagram showing the distribution of DEGs across all germ cells upon BDE47 treatment. Color intensity indicates the changed levels (P value < 0.05 and FC > 10%). (B and C) Representative top GO terms of up-regulated (B) and down-regulated (C) genes in different germ cells upon BDE47 treatment. Circle size indicates the gene count. (D) Violin plots showing gene set scores of double-strand breaks in different germ cells upon BDE47 treatment. Data were presented as mean ± SEM of three control and three BDE47-treated testis; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001. (E) Immunochemistry showing BDE47-increased γ-H2AX, a marker for DNA damage. (F) Violin plots showing gene set scores of chromatin remodeling in different germ cells upon BDE47 treatment. Data were presented as mean ± SEM of three control and three BDE47-treated testis; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001. (G) Dot plot showing significantly changed TFs in different germ cells upon BDE47 treatment. (H) Network plot showing the differentially expressed TFs and corresponding differentially expressed target genes in different somatic cells. (I) Bar plot showing the siginifacant enrichment of GO terms of Kdm5b target genes, Y axis represents the Log10(P value). (J) Netwok plot showing Kdm5b-regulated DEGs in response to BDE47 in testes. Color intensity indicates the changed levels (P value < 0.05 and FC > 10%). (K) BDE47 decreases Kdm5b mRNA levels as reflected by qPCR analysis.**P<0.01. (L) BDE47 decreases KDM5B levels as reflected by immunostaining analysis  .

Moreover, we determined the changes in TFs and their target genes that were changed in germ cells upon BDE47 treatment (Fig. 4G and H, and Fig. S3). We identified multiple common or cell type-specific TFs that are changed upon BDE47 treatment (Fig. S3). Several of these TFs play critical roles in regulating multiple stages in spermatogenesis. For example, Rfx1 (regulatory factor X 1) was significantly decreased in SPG though it has also been reported to be expressed in other cells such as Sertoli cells.[55]Rfx2, a key regulator of spermiogenesis,[56] has been identified to be specifically downregulated in pachytene though it is expressed in multiple germ cells.[56] Since the regulatory factor for X-box (RFX) family of TFs is involved in gene regulation during spermatogenesis, our results suggest that cell-type-specific downregulation of Rfx1 and Rfx2 may be involved in BDE47-induced testicular injury that affects spermatogenesis. As another example, significant downregulation of Sox30 was identified in round STids. Interestingly, deficiency of Sox30 induces the arrest of meiotic germ cells at the post-meiotic round spermatid period,[57] suggesting that SOX30 may play a role in BDE47-induced impairment of spermatogenesis. Furthermore, our results revealed that the downregulation of Smarca4 was detected in pachytene, meiotic division, and innate lymph cells. SMARCA4/BRG1, which encodes the catalytic subunit of the switch/sucrose‐non‐fermentable (SWI/SNF) chromatin remodeling complex, is critical for regulating spermatogenic transcription and meiotic progression.[58]

Interestingly, we found that Kdm5b (lysine-specific demethylase 5B), which plays an important role in post-meiotic spermatid cells,[59] was significantly decreased in all germ cells. GO analysis showed that Kdm5b affected multiple processes involved in spermatogenesis, such as steroid hormone-mediated signaling pathway, spermatid differentiation, RNA splicing, mRNA processing, DNA repair, and chromatin remodeling (Fig.4I). We also mapped multiple Kdm5b-regulated DEGs in reposne to BDE47 in testes (Fig.4J). Interestingly, downregulation of Kdm5b was also found in Sertoli and telocyte cells, further highlighting the critical roles of Kdm5b in BDE47-associated testicular injury. The reduced KDM5B at both mRNA (Fig. 4K) and protein (Fig.4L) levels in response to BDE47 were further verified by qPCR and immunostaining analysis, respectively. Overall, our results reveal the critical roles of multiple TFs involved in BDE47-induced impairment of spermatogenesis. Future studies to dissect the roles of these TFs in BDE47-induced testicular injury are an interesting topic. Taken together, our ScRNA-seq results reveal shared and cell type-specific genes, pathways, and TFs that are associated with spermatogenesis in response to BDE47 in germ cells within testis, shedding light on how BDE47 affects spermatogenesis.

Characterization of BDE47-related changes in intercellular communication

Characterization of ligand-receptor pairs can provide experimental information for understanding central cellular components that affect tissue fate. The interaction of somatic cells and germ cells is critical for spermatogenesis.[60] To understand the potential ligand-receptor interactions between somatic cells and germ cells upon BDE47 treatment, we further examined the changes in receptors and ligands in different cell types. Our results showed that BDE47 reduced both the overall interaction numbers and interaction strength across major cell types in testis (Fig. 5A-C). To systematically dissect cellular interactions between different cells within testis upon BDE47 treatment, we further constructed a changed comprehensive intercellular network of potential ligand–receptor interactions upon BDE47 treatment using a previously established database of ligand-receptor pairs.[61] Our result showed that multiple receptor-ligand pairs were significantly decreased upon BDE47 treatment (Fig. 5D). For instance, Nectin 3-Nectin 2 interaction was downregulated in multiple cells. Nectins are Ca2+-independent cell adhesion molecules at adherens junctions that play a critical role in spermatid development.[62] These findings suggest that impairment of the interaction of Nectin 3-Nectin 2 may contribute to BDE47-induced testicular injury. We also showed an increase in several receptors-ligands upon BDE47 treatment (Fig. 5E). For instance, the interaction of ectonucleoside triphosphate diphosphohydrolase 1 (ENTPD1)-transmembrane and immunoglobulin domain containing 3 (TMIGD3) was increased in several cell types upon BDE47 treatment. ENTPD1 metabolizes extracellular ATP and ADP, and is critical for regulating the inflammatory response.[63] The increased ENTPD1-TMIGD3 may indicate the increased inflammatory response upon BDE47 treatment in multiple cell types (such as pachytene and innate lymph cells).

Figure 5.

Characterizations of intercellular communication upon BDE47 treatment. (A-C) Heatmap plot (A) and quantification data (B and C) showing interaction numbers and strengths of receptors and ligands between different cell types upon BDE47 treatment. (D) Bubble plot showing significantly downregulated ligand-receptor pairs upon BDE47 treatment. (E) Bubble plot showing significantly increased ligand-receptor pairs upon BDE47 treatment.

Conclusion and discussion

Here, for the first time, we establish a comprehensive ScRNA-seq atlas of BDE47-associated changes in genes, pathways, TFs, and receptor-ligand pairs, which helps to understand shared and cell type-specific changes in genes and relevant biological processes that are associated with BDE47-induced testicular injury at unprecedented resolution. Our ScRNA-seq atlas not only provides novel insight into our understanding of the underlying mechanisms and pathways for BDE47-induced testicular injury but also serves as an important resource for future characterization of endocrine-disrupting chemicals such as BDE47-induced testicular toxicity.

Our ScRNA-seq comprehensively detected the dynamic gene expression profiles in different cell types within testis upon BDE47 treatment. Since spermatogenesis is a complex process and transcriptomic changes in different cell types and even in individual cells within the same cell type vary greatly, it is much difficult to distinguish such changes by the traditional RNA-seq method. Thus, our results provide additional information and act as a rich resource for understating and dissecting the molecular mechanisms of BDE47-induced testicular toxicity in mice. By using these genes, GO analysis was further used to dissect relevant pathways that may be associated with BDE47-induced testicular injury in somatic and germ cells. Consistent with previous findings that several pathways such as hormone homeostasis, inflammation response, and ER stress were enriched in response to BDE47,[23,25,64] our ScRNA-seq results further showed multiple shared or cell-type-specific enrichment of pathways such as double-strand break repair, cytokinesis, and histone modification in different cell types, suggesting the critical roles of these pathways in BDE47-induced toxicity. Interestingly, our ScRNA-seq results revealed an increase of innate lymph cells upon BDE47 treatment, which may reflect the increased inflammatory response in testicular tissue possibly through disruption of the blood–testis barrier upon BDE47 treatment in mice testis. However, the exact underlying mechanism needs to be further clarified.

Leydig cells are crucial for the secretion of testosterone in response to luteinizing hormone and they play critical roles in spermatogenesis. Previous studies regarding BDE47-associated testicular toxicity have mainly focused on Leydig cells. For instance, BDE47 has been reported to act as an endocrine disruptor chemical by interfering with thyroid homeostasis and impairing testicular steroidogenesis.[65] Indeed, GO analysis showed that BDE47 inhibited gene sets associated with the steroid metabolic process (Fig.3D). Our ScRNA-seq results showed that Hsd3b6 is one of the most dramatically downregulated genes in Leydig cells upon BDE47 treatment. Hsd3b6 is a gene that encodes hydroxy-δ-5-steroid dehydrogenase, or 3β- and steroid delta-isomerase 6. The 3-β-hydroxysteroid dehydrogenase (HSD) enzymatic system is critical for the biosynthesis of all classes of hormonal steroids and is involved in the local production of progesterone. The downregulation of Hsd3b6 in Leydig cells indicates that BDE47 may inhibit the production of thyroid hormone and progesterone, which is consistent with previous findings that BDE47 inhibited progesterone synthesis[44] and decreased thyroid hormone levels in female mice.[66] Thyroid hormone controls androgen biosynthesis and signaling through direct and indirect regulation of steroidogenic enzyme expression and activity such as via induction of Star expression. STAR encodes a steroidogenic acute regulatory protein, a rate-limiting enzyme that is critical for the modulation of steroid hormone synthesis by promoting the conversion of cholesterol into pregnenolone.[67] Pregnenolone is then converted by a series of enzymes into various steroid hormones in specific tissues.[67] The downregulated Star suggests that BDE47 may inhibit steroidogenesis. This result is consistent with a previous report showing that BDE47 reduced Star levels in Leydig cells.[68] Furthermore, the downregulation of Asah1 further suggests that BDE47 inhibits steroidogenesis. Acid ceramidase (ASAH1) is a key regulator of steroidogenic capacity by modulating the expression of multiple steroidogenic genes.[69] Interestingly, BDE47 also significantly decreased the expression of key genes in cholesterol biosynthesis such as DHCR24 (3β-hydroxysterol Δ[24]-reductase) gene, and CYBR3 (cytochrome B5 Reductase 3) gene. These results further revealed that BDE47 inhibits testosterone synthesis since cholesterol is the precursor to all steroid hormones. Overall, our Sc-RNA results reveal that BDE47 downregulates multiple genes responsible for steroid hormone synthesis such as Hsd3b6, Star, Asah1, Dhcr24, and Cyb5r3 in Leydig cells,[70] indicating that BDE47 may inhibit spermatogenesis via inhibition of testosterone production.

Sertoli cells are a kind of somatic cells within the testis that are crucial for testis formation and spermatogenesis in response to follicle stimulating hormone (FSH) and testosterone.[71] Here, we show that several genes involved in the testosterone receptor signaling pathway, such as the nuclear receptor corepressor 1 (NCoR 1) 1 and Kdm3a were significantly decreased in Sertoli cells. NCoR 1 is a corepressor of the unliganded thyroid receptor and a physiological regulator of the androgen receptor.[72] The histone lysine demethylase 3A (KDM3A) is a key factor for modulating the transcriptional program of the androgen receptor and spermatogenesis.[73] Our ScRNA-seq results revealed significant downregulation of Ncor1 and Kdm3a, suggesting that BDE47-induced testicular injury may be via regulating androgen receptor signaling in Sertoli cells. Moreover, several other key pathways such as cilium organization, cilium assembly, and regulation of chromosome organization were also significantly enriched in Sertoli cells upon BDE47 treatment.

Our ScRNA-seq results also showed multiple pathways that were unappreciated previously, such as downregulation of double-strand break repair, cytokinesis, histone modification, chromatin remodeling, and others in different subtypes of the testis, which may be involved in BDE47-induced testicular injury. Future studies are required to dissect the functions of these pathways in BDE47-induced reproductive toxicities in mice.

We also uncovered multiple key TFs (e.g.Kdm5b) that may be associated with BDE47-induced toxicity by analyzing ScRNA-seq data. These findings provided key evidence for our further understanding of how TFs function in different cell types that affects spermatogenesis in response to BDE47. Though the majority of TFs identified are associated with spermatogenesis, how these TFs (e.g. Tcfl5) contribute to BDE47-induced testicular injury is largely unclear. Future studies aiming to dissect the specific roles of the above TFs identified in BDE47-induced reproductive toxicity are an important topic.

Finally, we showed the changes in receptor/ligand upon BDE47 treatment, providing another angle for understanding BDE47-induced prepubertal testicular injury in mice. Certain factors secreted from either somatic cells or germ cells may affect multiple biological processes to impair spermatogenesis in mice. Interestingly, several changes in the interaction of multiple receptor/ligands upon BDE47 treatment that we discovered have not been reported previously. Therefore, our results provide a novel resource to better understand how intercellular communication impairs spermatogenesis in mice.

Notably, 10 mg/kg/day of BDE47 treatment of mice was used in our study. The dosage is relevant and consistent with previous studies of laboratory exposures in mammalian animal models[4,8,26,27]. However, future studies using lower concentrations and longer exposure times to BDE47 to systematically characterize BDE47-induced testicular injuries are warranted.

Overall, for the first time, to our knowledge, we have constructed a ScRNA-seq atlas of BDE47-associated testicular injury in mice. These results not only provide a novel insight into the underlying mechanism but also serve as a novel resource for future characterization of BDE47-induced male reproductive toxicity.

Click here for additional data file.