Single-cell transcriptomic insights into endosulfan-induced liver injury: Key pathways and inflammatory responses.

To explore the molecular and cellular changes in the liver microenvironment after ENDO treatment, we first constructed a mouse model and performed scRNA-seq as described in the workflow chart ( Fig. 1 A). Histopathological analysis revealed significant inflammatory infiltration in the ENDO-treated group compared to the Ctrl group, indicating that ENDO induces liver damage ( Fig. 1 B). Using droplet-based scRNA-seq (10x Genomics platform) of liver tissues from Ctrl and ENDO-treated groups ( n = 3, respectively), we obtained 76,470 single cells ( n = 37,298 in Ctrl and n = 39,172 in the ENDO-treated group) with strict standards ( Supplemental Fig. 1A ) after quality control, in which the number of detected genes (nFeature) and UMIs (nCount) showed a strong positive correlation across each sample ( Supplemental Fig. 1B ). These cells were then integrated and subjected to principal component analysis (PCA) to dimensional reduction, resulting in 17 clusters ( Supplemental Fig. 1C ), which were identified as major cell types based on the canonical markers expression levels ( Fig. 1 C, Supplemental Fig. 1D ): hepatocyte (Hep, Ttr + ), cholangiocyte (Cho, Sox9 + ), hepatic stellate cell (HSC, Dcn + ), endothelial cell (Endothelial, Kdr + ) and immune cell (Immune, Ptprc + ). The DEGs of each cell type revealed distinct transcriptional features and enriched biological processes ( Fig. 1 D), for example, epithelial cell development and extracellular matrix organization enriched for Hep and HSC cells, respectively. The cellular composition of each cell type in the two groups ( Fig. 1 E) showed that Hep and HSC cells were increased, while Cho, endothelial, and immune cells were reduced after ENDO treatment, which showed slight individual heterogeneity in the individual samples ( Fig. 1 F, Supplemental Table 1 ). In addition, we investigated the transcriptional program alteration between Ctrl and ENDO-treated groups and found that Hep cells harbored the most DEGs ( n = 353 and n = 501 in up- and down-regulation, respectively) ( Fig. 1 G), suggesting the need to perform more in-depth exploration of the major cell types in the follow-up analysis. Fig. 1 Cellular composition of mice liver with ENDO treatment via scRNA-seq. (A) Schematic model of the experimental design of this study. (B) Histological examination staining of livers in Ctrl (left) and ENDO-treated groups (right), scale bar = 100 μm. (C) tSNE of the 76,470 profiled cells, colored by cell types. (D) Heatmap of cluster-specific DEGs in each cell type (left) and enriched biological processes (right). (E) tSNE plot of cell clustering in Ctrl (left) and ENDO-treated (right) groups. The relative proportion of each cell type in each group (labeled at right) is shown at the bottom by bar plots. (F) Bar plots of the relative proportion of each cell type in each sample. (G) Bar plots of DEG numbers in each cell type (labeled on the right). Abbreviations: ATP, adenosine triphosphate; Cho, cholangiocyte; Ctrl, control; DEG, differentially expressed gene; ENDO, endosulfan; Endothelial, endothelial cell; Hep, hepatocyte; HSC, hepatic stellate cell; Immune, immune cell; scRNA-seq, single-cell RNA sequencing; tSNE, t-distributed Stochastic Neighbor Embedding. Fig. 1 Cellular composition of mice liver with ENDO treatment via scRNA-seq. (A) Schematic model of the experimental design of this study. (B) Histological examination staining of livers in Ctrl (left) and ENDO-treated groups (right), scale bar = 100 μm. (C) tSNE of the 76,470 profiled cells, colored by cell types. (D) Heatmap of cluster-specific DEGs in each cell type (left) and enriched biological processes (right). (E) tSNE plot of cell clustering in Ctrl (left) and ENDO-treated (right) groups. The relative proportion of each cell type in each group (labeled at right) is shown at the bottom by bar plots. (F) Bar plots of the relative proportion of each cell type in each sample. (G) Bar plots of DEG numbers in each cell type (labeled on the right). Abbreviations: ATP, adenosine triphosphate; Cho, cholangiocyte; Ctrl, control; DEG, differentially expressed gene; ENDO, endosulfan; Endothelial, endothelial cell; Hep, hepatocyte; HSC, hepatic stellate cell; Immune, immune cell; scRNA-seq, single-cell RNA sequencing; tSNE, t-distributed Stochastic Neighbor Embedding. Due to the significant increase in the number of hepatocytes and the transcriptional alterations following ENDO treatment, we next reclustered 6710 Hep cells into 7 subtypes according to their changes in cellular composition and expression matrix ( Fig. 2 A, Supplemental Fig. 2A ). We found that Hep-c1, Hep-c2 and Hep-c6 subtypes, which were associated with fatty acid metabolic process, enriched more in ENDO-treated group, while Hep-c4 with blood coagulation, and Hep-c7 with humoral immune response decreased the abundance ( Fig. 2 B, Supplemental Fig. 2B ). Besides, Hep-c3 and Hep-c5, with RNA splicing regulation, showed increased and decreased, respectively. We then mapped liver zone-related markers (from central vein (CV) to portal node (PN)) to the distinct hepatocyte subtypes ( Supplemental Fig. 2C ) based on the previous studies. 41 , 48 , 49 It is found that CV markers were enriched in Hep-c1 cells, and PN markers were mainly located in Hep-c2, c6, and c7 cells; however, Hep-c3, c4, and c5 cells did not show any obvious location characteristics. Furthermore, we calculated the key pathways shaping liver zonation profiles ( Supplemental Fig. 2D ), including Wnt, RAS, and hypoxia. 47 , 48 , 49 , 50 This analysis clarified that Wnt signaling was activated and Ras signaling was inhibited in CV-related cells (Hep-c1), while hypoxia features were inhibited in PN-related cells ( e.g., Hep-c7). Thus, our results clustered hepatocytes into 7 subtypes, including CV-related Hep-c1 cells and PN-related Hep-c2, c6, and c7 cells, which exhibited zone-specific functions. Fig. 2 Hepatocyte-specific dynamic response to ENDO. (A) tSNE of 6710 hepatocytes, colored by subtypes, split by group type (right). (B) Dot plots of each hep subtype in Ctrl and ENDO-treated groups. (C) Strip chart of up-(top) and down-(bottom) regulated DEGs between ENDO-treated and Ctrl groups in each subtype. The enriched terms were listed. (D) Pathway scores comparison in the two groups of each subtype. (E) Heatmap showing expression levels of selected markers. (F) Dot plots indicating the top 3 regulons in each subtype. ∗ P < 0.05, ∗∗ P < 0.01, ∗∗∗ P < 0.001, ∗∗∗ P < 0.0001; ns, not significant. Abbreviations: Ctrl, control; ENDO, endosulfan; Hep, hepatocyte; tSNE, t-distributed Stochastic Neighbor Embedding. Fig. 2 Hepatocyte-specific dynamic response to ENDO. (A) tSNE of 6710 hepatocytes, colored by subtypes, split by group type (right). (B) Dot plots of each hep subtype in Ctrl and ENDO-treated groups. (C) Strip chart of up-(top) and down-(bottom) regulated DEGs between ENDO-treated and Ctrl groups in each subtype. The enriched terms were listed. (D) Pathway scores comparison in the two groups of each subtype. (E) Heatmap showing expression levels of selected markers. (F) Dot plots indicating the top 3 regulons in each subtype. ∗ P < 0.05, ∗∗ P < 0.01, ∗∗∗ P < 0.001, ∗∗∗ P < 0.0001; ns, not significant. Abbreviations: Ctrl, control; ENDO, endosulfan; Hep, hepatocyte; tSNE, t-distributed Stochastic Neighbor Embedding. Similarly, the transcriptional alterations in each subtype could jointly upregulate metabolic progress, while downregulating RNA splicing and epithelium migration ( Fig. 2 C). Thus, we assessed the metabolic pathway activity in each subtype of the Ctrl and ENDO-treated groups using the MSigDB database ( Fig. 2 D), focusing on pathways such as fatty acid and bile acid metabolism, which in agreement with the results that fatty acid metabolic related subtypes mainly increased the abundance after ENDO treatment ( Fig. 2 B, Supplemental Fig. 2B ). In addition, previous studies have demonstrated that ENDO exposure can lead to the excessive production of reactive oxygen species (ROS), resulting in oxidative stress and toxic effects on cells. 51 , 52 To investigate this, we calculated the ROS and peroxisome signature scores and observed the stronger activity after ENDO treatment ( Fig. 2 D), consistent with the toxic effects of ENDO that disrupt cellular redox balance and lead to oxidative damage. To further dissect these changes, we calculated the expression levels of several key genes involved in metabolism and ROS. We observed an increase in fatty acid metabolic genes, alongside a general decrease in several antioxidant genes, such as catalase ( Cat ), superoxide dismutase 1 ( Sod1 ), glutathione transferase ( Gsta1 , Gstm3 , Gsto2 , Gstp2 ) ( Fig. 2 E), although differed in certain subtypes. This weakening of the cell’s ability to scavenge ROS leads to the accumulation of ROS within the cell, triggering oxidative stress. These changes were predominantly enriched in Hep-c1 and Hep-c6 cells, suggesting that these subtypes may be particularly sensitive to ENDO-induced oxidative stress. To further identify the potentially driven TFs that regulate gene expression in Hep cells, we applied the SCENIC pipeline to each Hep subtype. Among the top 3 regulons identified ( Fig. 2 F), we distinguished the activation of interferon regulatory factor 1 (Irf1) in Hep-c1 and activating transcription factor 5 (Atf5) in Hep-c6 cells. Notably, Irf1 knockout could attenuate ROS production and Atf5 dysfunction might impair electron transport chain efficiency, leading to increased electron leakage and ROS generation. 53 , 54 These findings further support the hypothesis that ENDO-induced oxidative stress may be mediated, at least in part, by the Irf1 and Atf5 TFs, which regulate key genes involved in the cellular response to oxidative damage. HSCs showed increased abundance after ENDO treatment under liver injury conditions. Based on the transcription programs and locations of HSC cells, we divided them into three subtypes: one quiescent (HSC-c1: Rgs5 + Lrat + ) and two activated (HSC-c2: Fmod + Timp1 + , HSC-c3: Timp2 + Col1a1 + ) cells ( Fig. 3 A–C). Compared to the Ctrl group, ENDO treatment upregulated the expression levels of genes involved in extracellular matrix (ECM) (such as Lamb1 , Fbln5 , Smoc2 , Fbln1 ), ECM remodeling (such as Adamts5 , Adamts2 , Timp2 , Htra3 , Pcolce ), and proteoglycans (such as Srpx2 , Dcn , Fmod , Ogn ), illustrating the collagen deposition translation of HSCs ( Fig. 3 D). Subsequently, we depicted the trajectory of HSCs, which was rooted from HSC-c1 and ended by HSC-c2 ( Fig. 3 E). The epithelial-mesenchymal transition (EMT) signature score was enriched in HSC-c3 cells ( Fig. 3 E), and showed an increase in the abundance after ENDO treatment ( Fig. 3 C), thereby further confirming the activation of HSC cells. Along with the trajectory development, HSC cells were associated with cell junction assembly, extracellular matrix organization, and epithelial cell proliferation ( Fig. 3 F). Thereinto, genes in Fig. 3 D also showed a similar pattern of activation ( Fig. 3 G). Fig. 3 Identification of activated subtypes in HSC cells. (A) tSNE of 243 HSC cells, colored by subtype. (B) Heatmap of cluster-specific DEG expression in each subtype, labeled by representative markers. (C) tSNE plot of cell clustering in Ctrl (up) and ENDO-treated (bottom) groups. (D) Heatmap of selected genes expression level during ENDO treatment. (E) Pseudotime trajectory of HSC subtypes and EMT score distribution along the cell development. (F) Heatmap of DEGs associated with pseudotime change and enriched biological processes. (G) Expression of genes during pesudotime development. Abbreviations: Ctrl, control; DEGs, differentially expressed genes; EMT, epithelial-mesenchymal transition; ENDO, endosulfan; HSC, hepatic stellate cell; tSNE, t-distributed Stochastic Neighbor Embedding. Fig. 3 Identification of activated subtypes in HSC cells. (A) tSNE of 243 HSC cells, colored by subtype. (B) Heatmap of cluster-specific DEG expression in each subtype, labeled by representative markers. (C) tSNE plot of cell clustering in Ctrl (up) and ENDO-treated (bottom) groups. (D) Heatmap of selected genes expression level during ENDO treatment. (E) Pseudotime trajectory of HSC subtypes and EMT score distribution along the cell development. (F) Heatmap of DEGs associated with pseudotime change and enriched biological processes. (G) Expression of genes during pesudotime development. Abbreviations: Ctrl, control; DEGs, differentially expressed genes; EMT, epithelial-mesenchymal transition; ENDO, endosulfan; HSC, hepatic stellate cell; tSNE, t-distributed Stochastic Neighbor Embedding. Lymphocytes play a crucial role in immune regulation within the liver, especially under injury conditions. We further identified lymphocytes (T/NK: Cd3d + Ncr1 + , B: Cd79a + Jchain + ) in the immune population ( Fig. 4 A–B), and then reclustered into 19 subtypes ( Fig. 4 C). Among these, B-c7 was characterized as plasma cells with high expression of Jchain and Igha ( Fig. 4 D), as well as the highest cumulative distribution of plasma genes ( Fig. 4 E). Similarly, Cd4-c1 and Cd8-c1 were identified as naive cells ( Lef1 + Tcf7 + ), Cd4-c2 as regulatory T cells (Treg, Foxp3 + Il2ra + ), and Cd8-c2 as effector cells ( Cx3cr1 + Ccl5 + ) ( Fig. 4 D–E). 55 , 56 To investigate the immune effector induced by ENDO treatment, we quantified cytotoxic and cytokine signature scores ( Fig. 4 F), represented by several key genes ( Nkg7 , Gzmb , Gzmk , Fasl , Ifng , Ccl5 ) ( Supplemental Fig. 3A ). A significant increase in these scores was observed, especially in Cd4-c2, Cd4-c4, Cd8-c2, Cd8-c4, and gdT cells ( Fig. 4 F), we also noticed that the above subtypes were more abundance in the ENDO-treated group ( Fig. 4 C), revealing an enhanced immune response. In addition, ENDO treatment could reprogram immune cells, characterized by a well-defined developmental trajectory that started with Ctrl cells and ended with ENDO-treated cells ( Fig. 4 G, Supplemental Fig. 3B–C ). As for Cd8 + T cells, the naive subtype (Cd8-c1) developed into effector/memory subtypes (Cd8-c2 and Cd8-c4), along with a proliferation feature (Cd8-c5), characterized by lymphocyte chemotaxis ( Fig. 4 H). We identified DEGs in the three increased subtypes ( Fig. 4 I) and found the chemokines ( Ccl3 , Ccl4, Ccl5 ) were upregulated after ENDO treatment, consistent with the enrichment processes shown in Fig. 4 H. The relative expression levels of the key chemokines ( Ccl3 , Ccl4 , Ccl5 ) and cytotoxic genes ( Fasl , Klrc1 ) along with Cd8 + T cell’s trajectory were upregulated in ENDO-treated group, as presented in Fig. 4 J. Fig. 4 Lymphoid immune cells diversity and variation . (A) Lymphocyte distribution among immune cells. (B) Expression levels of lymphocyte markers. (C) tSNE of lymphocytes, colored by subtype, with the relative proportion of each subtype listed at the right bottom. (D) Heatmap of cluster-specific DEG expression. (E) Cumulative plot showing the distribution of naive B, plasma B, naive T, cytotoxic T, and regulatory state scores across lymphocyte subtypes. (F) Cytotoxic (top) and cytokine (bottom) gene scores in the two groups for each subtype. (G) Pseudotime trajectory of B (top) and CD8 + T cells (bottom), colored by group type. (H) Heatmap of gene dynamics along with the trajectory. Genes emerged at the end and enriched biological processes were listed on the right. (I) Fold changes of DEGs existed at Cd8-c2, Cd8-c4, and Cd8-c5 cells. (J) Expression of genes during pseudotime development. ∗ P < 0.05, ∗∗ P < 0.01, ∗∗∗ P < 0.001, ∗∗∗∗ P < 0.0001; ns, not significant. Abbreviations: CD, cluster of differentiation; Ctrl, control; DEGs, differentially expressed genes; ENDO, endosulfan; FCs, fold changes; HSC, hepatic stellate cell; tSNE, t-distributed Stochastic Neighbor Embedding. Fig. 4 Lymphoid immune cells diversity and variation . (A) Lymphocyte distribution among immune cells. (B) Expression levels of lymphocyte markers. (C) tSNE of lymphocytes, colored by subtype, with the relative proportion of each subtype listed at the right bottom. (D) Heatmap of cluster-specific DEG expression. (E) Cumulative plot showing the distribution of naive B, plasma B, naive T, cytotoxic T, and regulatory state scores across lymphocyte subtypes. (F) Cytotoxic (top) and cytokine (bottom) gene scores in the two groups for each subtype. (G) Pseudotime trajectory of B (top) and CD8 + T cells (bottom), colored by group type. (H) Heatmap of gene dynamics along with the trajectory. Genes emerged at the end and enriched biological processes were listed on the right. (I) Fold changes of DEGs existed at Cd8-c2, Cd8-c4, and Cd8-c5 cells. (J) Expression of genes during pseudotime development. ∗ P < 0.05, ∗∗ P < 0.01, ∗∗∗ P < 0.001, ∗∗∗∗ P < 0.0001; ns, not significant. Abbreviations: CD, cluster of differentiation; Ctrl, control; DEGs, differentially expressed genes; ENDO, endosulfan; FCs, fold changes; HSC, hepatic stellate cell; tSNE, t-distributed Stochastic Neighbor Embedding. In addition to lymphoid immune cells, we reclustered myeloid cells, including Kupffer cell (KC: C1qa + Clec4f + ), monocyte/macrophage (Mono/Macro: Ccr2 + Cx3cr1 + ), neutrophil (Neut: S100a8 + S100a9 + ), plasmacytoid dendritic cell (pDC: Siglech + Ccr9 + ), dendritic cell (DC: Naaa + Cst3 + ), and mast cell (Mast: Cpa3 + Mcpt8 + ) ( Fig. 5 A, Supplemental Fig. 4A ). Furthermore, 14 subtypes were identified based on cell composition and transcriptional programs ( Fig. 5 A–C). We characterized the main functions of each subtype and found that Macro-c2 ( Cx3cr1 + Trem2 + Cxcl16 + ), Macro-c3 ( Cd209a + Ifi30 + Malt1 + ), KC-c1 ( C1qa + Cd163 + Folr2 + , M2-like), and DC cells exhibited a stronger capacity for antigen processing and presentation; while Mono ( Cd14 + Lyz2 + ), KC-c1, and Macro-c4 ( Ear2 + Flna + Ace + ) cells exhibited significant phagocytosis capacity, as well as Mast, Neut-c2 ( Slpi + ), and Neut-c1 ( Fabp1 + ) showed elevated cytokine score ( Fig. 5 D). To explore the molecular basis of these functions, we calculated the relative expression levels of key genes from the signatures above, primarily enriched in KC-c1, KC-c2, Macro-c1, Macro-c2, and Neut-c1 cells ( Fig. 5 E). Therefore, we further explored the transcriptional changes of these key subtypes. DEGs between Ctrl and ENDO-treated groups were identified ( Supplemental Fig. 4B ), and pathway enrichment revealed dysregulation of antigen processing and presentation, interferon response, and cell chemotaxis ( Supplemental Fig. 4C ). There were 6 DEGs jointly differentially expressed in the 5 subtypes, in which Pdia3 and AY036118 were upregulated ( Fig. 5 F). Pdia3 has been reported to mediate anti-inflammatory effects, 57 , 58 thus we calculated the correlation coefficient between Pdia3 and inflammatory cytokines. We found a positive correlation with anti-inflammatory cytokines ( Il10 , Tgfb1 , Mrc1 , Cd163 ) and a negative correlation with proinflammatory cytokine ( Il1b ) ( Fig. 5 G). However, the biological mechanisms and significance of Pdia3 during ENDO treatment in mice liver require further investigation. Fig. 5 Clustering of distinct myeloid cells. (A) tSNE of myeloid cells including KC, Mono, Macro, neutrophil (Neut), pDC, DC, and mast cells. (B) Dot plots of each myeloid subtype in Ctrl and ENDO-treated groups. (C) Expression levels of representative markers identifying each subtype. (D) Module scores sorting in each subtype. (E) Different expression of cytokines in the selected subtypes. (F) Venn plot showing the overlapping DEGs of the key subtypes. (G) Heatmap illustrating the correlation between Pdia3 and pro- or anti-inflammatory cytokines. Correlation coefficients were labeled in the block. Abbreviations: Ctrl, control; DC, dendritic cell; DEGs, differentially expressed genes; ENDO, endosulfan; HSC, hepatic stellate cell; KC, Kupffer cell; Macro, macrophage; Mono, monocyte; Neut, neutrophil; pDC, plasmacytoid dendritic cell; tSNE, t-distributed Stochastic Neighbor Embedding. Fig. 5 Clustering of distinct myeloid cells. (A) tSNE of myeloid cells including KC, Mono, Macro, neutrophil (Neut), pDC, DC, and mast cells. (B) Dot plots of each myeloid subtype in Ctrl and ENDO-treated groups. (C) Expression levels of representative markers identifying each subtype. (D) Module scores sorting in each subtype. (E) Different expression of cytokines in the selected subtypes. (F) Venn plot showing the overlapping DEGs of the key subtypes. (G) Heatmap illustrating the correlation between Pdia3 and pro- or anti-inflammatory cytokines. Correlation coefficients were labeled in the block. Abbreviations: Ctrl, control; DC, dendritic cell; DEGs, differentially expressed genes; ENDO, endosulfan; HSC, hepatic stellate cell; KC, Kupffer cell; Macro, macrophage; Mono, monocyte; Neut, neutrophil; pDC, plasmacytoid dendritic cell; tSNE, t-distributed Stochastic Neighbor Embedding. To explore the intercellular interactions of ENDO-induced liver injury, we performed an analysis of cellular communication in each subtype between the Ctrl and ENDO-treated groups ( Fig. 6 A). The ENDO-treated group harbored more interactions ( n = 3356) than the Ctrl group ( n = 3109) ( Fig. 6 B). We conducted a differential expression analysis to identify upregulated pathways ( Fig. 6 C), including COLLAGEN, CCL, and MHC-I, which corresponded to the activation of HSC cells, enhanced chemotaxis, and antigen processing and presentation of immune cells. Accordingly, we confirmed that HSC cells were the major sender of COLLAGEN signaling, and CD8 + T cells both initiated CCL signaling and served as the primary receiver of MHC-I signaling ( Fig. 6 D–F). Thus, we evaluated the communication probabilities of the LR pairs from the three pathways, finding that H2-k1-Cd8b1, H2-d1-Cd8a (in MHC-I), Col1a2-(Itga1+Itgb1) (in COLLAGEN), Ccl5-Ccr5 and Ccl5-Ccr1 (in CCL) pairs were increased in the ENDO-treated group ( Fig. 6 G), further supporting our in-depth analysis. Fig. 6 ENDO changes inter cellular communication in the liver. (A) Heatmap plot showing the differences in interaction number between sender and receiver subtypes in the ENDO vs . Ctrl comparison. (B) The number of interactions in the cell-cell communication network of Ctrl and ENDO-treated groups. (C) The significantly enriched pathways in Ctrl and ENDO-treated groups. (D – F) The number of interactions of COLLAGEN (D) , CCL (E), and MHC-I (F) pathways in Ctrl and ENDO-treated groups. (G) Bobble plot showing the specific ligand-receptor pairs in MHC-I, COLLAGEN, and CCL pathways. Abbreviations: CCL, C–C motif chemokine ligand; Ctrl, control; ENDO, endosulfan; Hep, hepatocyte; MHC-I, major histocompatibility complex class I; tSNE, t-distributed Stochastic Neighbor Embedding. Fig. 6 ENDO changes inter cellular communication in the liver. (A) Heatmap plot showing the differences in interaction number between sender and receiver subtypes in the ENDO vs . Ctrl comparison. (B) The number of interactions in the cell-cell communication network of Ctrl and ENDO-treated groups. (C) The significantly enriched pathways in Ctrl and ENDO-treated groups. (D – F) The number of interactions of COLLAGEN (D) , CCL (E), and MHC-I (F) pathways in Ctrl and ENDO-treated groups. (G) Bobble plot showing the specific ligand-receptor pairs in MHC-I, COLLAGEN, and CCL pathways. Abbreviations: CCL, C–C motif chemokine ligand; Ctrl, control; ENDO, endosulfan; Hep, hepatocyte; MHC-I, major histocompatibility complex class I; tSNE, t-distributed Stochastic Neighbor Embedding.

Pan Huang: Investigation, Writing –original draft, review & editing. Yunmeng Bai: Methodology, Investigation, Data analysis, Methodology, Writing – review & editing. Chaohua Zhou: Conceptualization, Writing – review & editing, Resources. Xiaoru Zhong: Conceptualization, Writing – review & editing, Resources. Ashok Iyaswamy: Methodology, Writing – review & editing. Peng Chen: Methodology, Resources. Xu Wei: Methodology, Resources. Wei Zhang: Supervision, Conceptualization, Writing – review & editing, Finding. Chuanbin Yang: Supervision, Conceptualization, Writing – review & editing, Funding. Jigang Wang: Software, Supervision, Conceptualization, Writing – review & editing, Funding.

The protocol for animal experiments was ethically reviewed and approved by the Institutional Animal Care and Use Committee (IACUC) of Shenzhen People’s Hospital (No. AUP-211009-WJG-0001-01). All animal experiments complied with the ARRIVE (Animal Research: Reporting In Vivo Experiments) guidelines and were carried out in accordance with the National Research Council’s Guide for the Care and Use of Laboratory Animals. C57BL/6 mice (aged 4–6 weeks) were purchased from GemPharmatech (Guangdong, China) and randomly divided into vehicle control (Ctrl) and ENDO-treated groups ( n = 6 per group). Mice were housed under controlled environmental conditions (temperature: 18–23 °C, humidity: 40%–60%, 12-h light/dark cycle) with free access to standard laboratory diet and water. The ENDO-treated group received daily intraperitoneal (i.p.) injections of ENDO (TRC, Shanghai, China) with 10 mg/kg dissolved in corn oil for two weeks, 10 while the Ctrl group received equivalent volumes of corn oil. At the end of drug treatment, mice were euthanized by CO 2 asphyxiation (fill rate of 30%–70% of the chamber volume per minute). Liver tissues were immediately collected, with a portion used for scRNA-seq and histological analysis. Liver tissues of Ctrl and ENDO-treated mice were collected for scRNA-seq analysis according to previous studies. 31 , 37 In brief, liver tissues of Ctrl and ENDO-treated mice ( n = 3 in each group) were cut into small pieces and digested with a Liver Dissociation Kit for mice (Miltenyi Biotec, Bergisch Gladbach, Germany) according to the manufacturer’s instructions. Digested cells were then filtered through a 70 μm strainer and single-cell suspension was acquired by centrifugation at 300 g for 10 min. Then the Debris Removal Solution (Miltenyi Biotec, Bergisch Gladbach, Germany) was applied to eliminate cellular debris, which was followed by the Red Blood Cell Removal Solution (Miltenyi Biotec, Bergisch Gladbach, Germany) to clear erythrocytes, and cell suspension was used for subsequent analysis. According to the instructions of the Chromium Single Cell 3’ Reagent Kit (v3.1; 10x Genomics, Pleasanton, CA, USA), scRNA-seq libraries were constructed using the single-cell suspensions prepared as described in Section 2.3 . Raw sequencing data were processed using Cell Ranger (v6.0.1; 10x Genomics, Pleasanton, CA, USA) and subsequently input into the Seurat R package (v4.0.4; Satija Lab, New York, NY, USA). 38 Low-quality cells were filtered based on the following criteria: 39 , 40 (1) gene number 6000; (2) unique molecular identifier (UMI) number >20,000; (3) the percentage of mitochondrial genome >15%. Data from the six samples were normalized and scaled using the SCTransform function and then integrated into a combined matrix. The correlation between the number of detected genes (nFeature) and UMIs (nCount) in each sample was calculated by the FeatureScatter function. Cell clustering was performed via FindNeighbors and FindClusters functions, and visualized by t-distributed Stochastic Neighbor Embedding (t-SNE) and Uniform Manifold Approximation and Projection (UMAP) algorithms. Cell types were identified according to the expression levels of canonical markers. 30 , 41 Differentially expressed genes (DEGs) between Ctrl and ENDO-treated groups were identified using the Wilcoxon rank-sum test in Seurat (v4.0.4; Satija Lab, New York, NY, USA) with (1) |log2 fold change| > 0.25; (2) Bonferroni-adjusted P value < 0.05. 42 Next, the Bioconductor package clusterProfiler (v4.0) was used to perform gene ontology (GO) analysis with the default parameter. 43 Besides, gene set variation analysis (GSVA) was performed by the Bioconductor R package GSVA (v1.40.0) using the pathways from the Molecular Signatures Database (MSigDB). 44 Genes’ correlations were calculated by the cor.test function. Transcription factors (TFs) were identified by single-cell regulatory network inference and clustering (SCENIC) (v1.2.4) using the gene-motif ranking: mm10_refseq-r80_500bp_up_and_100bp_down_tss.mc9nr and mm10_refseq-r80_10kb_up_and_down_tss.mc9nr. 45 Firstly, the correlation among genes was calculated by the runCorrelation function, and the potential TF-target gene interactions were inferred using the runGenie3 function. Then coexpression modules were generated by the runSCENIC_1_coexNetwork2modules function, and filtered cis-regulatory motif from cisTargetDB by the runSCENIC_2_createRegulons function. Subsequently, the activity of each regulon was scored by the runSCENIC_3_scoreCells function. The R package monocle2 (v2.18.0) was used to explore cell-state transitions and development. 46 The expression matrix and cell meta-information were obtained, and a CellDataSet object was constructed using the newCellDataSet function. The data was then normalized by estimateSizeFactors and estimateDispersions functions. Genes expressed in fewer than 10 cells were removed and ordering genes that might determine the progress were selected by the differentialGeneTest function, then set into CellDataSet object by the setOrderingFilter function. Next, dimensionality reduction with the DDRTree algorithm was performed by the reduceDimension function , and finally, cells were ordered along the trajectory by the orderCells function and visualized by the plot_cell_trajectory function. We performed cellular communication analysis by R package CellChat (v1.1.3). 47 First, the normalized gene matrix and meta information of the Ctrl and ENDO-treated groups were input and merged to compare the numbers of pairs between the two groups. Differentially expressed pathways were filtered with the rankNet function. Finally, the key ligand-receptor (LR) pairs that are differentially expressed in Ctrl and ENDO-treated groups of the pathways above were displayed by the netVisual_bubble function.

In this study, we may present the first single-cell atlas of ENDO-associated toxicity in the liver, providing comprehensive heterogeneous gene expression data and pathways across different cell types in mouse livers. Apart from hepatocyte injury, our scRNA-seq analysis results showed that inflammatory responses, such as myeloid cell recruitment and activation, along with the proinflammatory microenvironment, play key roles in ENDO-induced liver injury. Thus, our findings not only serve as a valuable resource for understanding the toxicity of environmental toxins but also may provide novel therapeutic targets for treating liver diseases. Consistent with previous findings that oxidative stress is involved in ENDO-induced hepatocyte injury, 25 , 59 we showed that multiple antioxidant enzymes such as Cat , Sod1 , and glutathione transferase ( Gsta1 , Gstm3 , Gsto2 , Gstp2 ) were mostly decreased in hepatocytes of ENDO-treated mice. Interestingly, in the brain and several other cells, ENDO-induced toxicity was also associated with increased oxidative stress. 60 , 61 These results suggest that oxidative stress may be a common mechanism for ENDO-induced cell toxicity in different cells. Activated HSCs are critical for liver injury and liver fibrosis through multiple mechanisms, such as chemotaxis secretion, increased cytokine-mediated inflammation, and NKT cell proliferation. 62 Our scRNA-seq results revealed that HSCs in ENDO-treated mice exhibited increased signaling pathways associated with extracellular matrix remodeling and proteoglycan formation, suggesting a potential promoting fibrosis effect of ENDO. Intriguingly, ENDO has also been shown to increase profibrotic markers in kidney cells. 63 These findings indicate that active profibrotic signaling in HSCs may also be involved in ENDO-induced liver injury. Immune cells, particularly lymphocytes, and the inflammatory response are critical for liver diseases. 64 Our scRNA-seq results revealed that ENDO increased cytotoxic and cytokine gene signature scores in NK/T cells, indicating the activation of NK/T cells and proinflammatory response in these cells. Furthermore, we observed upregulation of chemokines such as Ccl3 , Ccl4 , and Ccl5 in CD8 + T cells. These findings indicate that activated lymphocytes may promote liver injury in response to ENDO treatment. Notably, ENDO has been reported to induce the production of proinflammatory cytokine in human lymphocytes, 65 highlighting the proinflammatory effects of ENDO. In addition to lymphocytes, the recruitment and/or activation of myeloid cells in the liver also play important roles in liver injury and various liver diseases. Under certain conditions, such as liver injury, Kupffer cells produce multiple proinflammatory cytokines, contributing to liver toxicity. 66 Our findings showed that ENDO-treated mice, especially in KCs and DCs, exhibited dysregulation of genes associated with antigen presentation, interferon response, and cell chemotaxis, suggesting that activation of KCs and DCs in the liver of ENDO-treated mice. This activation may then recruit other immune cells into the liver, create an excess proinflammatory environment, and finally lead to liver injury. Furthermore, ENDO has been shown to induce an inflammatory response in mouse macrophages, increasing the expression of COX-2 and tumor necrosis factor alpha. 67 , 68 Importantly, in addition to direct liver cell injury, our results showed that inflammatory response plays a key role in ENDO-induced liver injury, suggesting potential therapeutic opportunities for ENDO-related liver diseases. Future studies should investigate this possibility in animal models. Notably, ENDO treatment also resulted in changes in cell-cell communication, especially the interaction of activated HSC with other cells, particularly CD8 + T cells, highlighting the critical roles of liver microenvironmental for its toxicity effects. Finally, there are several limitations to this study. First, the mouse model may not fully recapitulate the complexity of human liver physiology and diseases. Second, the short exposure duration may not capture the long-term effects, such as chronic liver fibrosis or cirrhosis. It would explore prolonged ENDO exposure to better understand its role in chronic liver diseases in future studies. Third, as the study mainly focused on transcriptomic data, the integration of proteomic/metabolomic analyses and experimental verification could offer a more comprehensive understanding of ENDO-induced liver injury.

For the first time, this study provides a transcriptome map of ENDO-induced liver injury in mice. This study comprehensively reveals that ENDO induces oxidative stress in hepatocytes, which may then promote the infiltration of cytotoxic T cells and the activation of macrophages to induce inflammatory responses. Our study indicates that the immune response plays a critical role in ENDO-associated toxicity, and provides a novel insight into the molecular characteristics of ENDO-induced liver injury at single-cell resolution.

Endosulfan (ENDO) is an organochlorine insecticide that was widely used from the 1950s until its global phase-out to control pests on multiple crops, including vegetables, tea, fruits, coffee, and cotton. 1 , 2 Due to its bioaccumulation in the environment and potential toxic effects, ENDO was classified as a persistent organic pollutant (POP), and its use was globally phased out under the Stockholm Convention in 2011. 3 Despite the global restrictions, ENDO continued to be used in certain regions, such as India, because of its broad-spectrum efficacy against insects, high efficiency, and low cost. For instance, a recent study showed that the concentration of ENDO in freshwater from Yavtmal, Maharashtra, India, reached a concerning level of 0.78 μg/L. 4 Additionally, multiple pesticide residues, including ENDO, in vegetables from areas such as the North-Western Himalayan region of India exceeded the maximum permissible concentrations established by the European Commission. 5 Moreover, due to ENDO’s ability to travel long distances in the environment, it still poses potential global health risks. 6 The United States Environmental Protection Agency listed ENDO as a “highly acutely toxic” chemical, and the World Health Organization (WHO) has classified it as a “moderately hazardous” class II chemical based on its acute toxicity effects. 2 It is well-established that ENDO exerts multiple toxicity effects across various organs in both insects and humans. 7 ENDO’s acute neurotoxicity is associated with its binding to the type A gamma-aminobutyric acid (GABA(A)) receptor in animal models. 8 Furthermore, ENDO is recognized as an endocrine disruptor, affecting hormone production, causing reproductive toxicity, and interfering with steroid hormone synthesis. 9 , 10 In addition, ENDO exposure has been associated with delayed sexual maturation in boys, probably due to reduced testosterone levels. 11 ENDO has also been shown to induce cardiotoxicity by triggering mitochondrial-mediated apoptosis, 12 and increasing the production of free radicals. 13 Apart from its effects on the brain, testes, and heart, ENDO also causes kidney toxicity, likely through apoptosis induction and oxidative stress. 14 , 15 Given the broad spectrum of organ toxicity associated with ENDO, it is important to consider its potential impact on other vital organs. Liver diseases, such as metabolic dysfunction-associated fatty liver disease (MAFLD) and hepatitis, are a major global health burden. 16 , 17 It is estimated that more than 20% of adults worldwide are affected by MAFLD. 16 , 17 Recently, environmental toxins, particularly pesticides, have raised public awareness due to their association with liver health disorders. 18 , 19 Exposure to environmental toxins may induce liver injury and is a key factor in the pathogenesis of many liver diseases. 20 Notably, as a metabolic organ, the liver is highly susceptible to ENDO-induced toxicity, which has been demonstrated in multiple cells and animal models, as well as in humans. For instance, ENDO exposure was reported to cause acute liver injury in a woman who ingested ENDO-contaminated food, characterized by elevated serum alanine aminotransferase (ALT) and aspartate aminotransferase (AST) levels. 21 Additionally, an autopsy of a man who ingested ENDO showed the characteristics of liver congestion. 22 In rats, ENDO exposure was shown to induce liver injury, increase key enzymes ( e.g., serum ALT and AST), and cause histopathological changes, inducing sinusoidal dilation and pyknotic nuclei. 23 , 24 Moreover, inhibition of nicotinamide adenine dinucleotide dehydrogenase and increased oxidative stress have been linked to ENDO-induced liver injury in rats. 25 However, the underlying mechanisms of ENDO-induced liver injury are largely unclear. Importantly, beyond hepatocytes, immune cells and inflammatory responses play crucial roles in the pathogenesis of liver diseases. 26 , 27 Yet, how ENDO affects the liver microenvironment is unclear. Recently, the use of advanced single-cell RNA sequencing (scRNA-seq) has expanded our understanding of the pathogenesis of liver diseases and the underlying mechanisms of toxin-induced liver injury. 28 , 29 , 30 , 31 , 32 Compared with traditional bulk sequencing, such as RNA sequencing, scRNA-seq enables detailed analysis of gene expression changes in different cell populations within a given organ, advancing our understanding of cell-type-specific gene expression and tissue microenvironmental changes. 28 , 29 , 33 , 34 , 35 , 36 Given the highly heterogeneous nature of the liver, especially its diverse immune cell subsets with distinct functions, the responses of different cell populations within the liver to environmental toxins such as ENDO remain unclear. Here, we collected the liver tissues from mice exposed to ENDO and performed scRNA-seq to explore gene expression changes in hepatocytes and intrahepatic immune cell types at single-cell resolution. Our results are expected to present ENDO-associated characteristics, especially the immune response of different cell subtypes in the liver, and offer a map of the mouse hepatic immune microenvironment in response to environmental toxins such as ENDO.

The authors declare that they have no conflict of interest.

Data are contained within the article or supplementary material. The data in the current study are available from the corresponding author upon request. The sequencing data reported in this study 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/release/OMIX004105 ).