GLY induces liver damage in mice

To investigate the potential hepatotoxicity of GLY at a single-cell level, C57BL/6 male mice were administered GLY for consecutive 7 days (Fig. 1A). Male mice were chosen since male animals suffered more acutely than females from liver damage [34]. Histopathological observation revealed fibroblast proliferation in the Glisson’s capsule (black arrow), exudation (green arrow), slight inflammatory cell infiltration (yellow arrow), and focal necrosis (blue arrow) in GLY-treated mice, suggesting that GLY can induce liver damage (Fig. 1B). Using the criteria in Additional file 1: Fig. S1A, the degree of necrosis of the liver tissue of each mouse from the CON- and GLY-treated group was examined. We found single-cell necrosis (numerical assessment number 1) in two of the three CON mice, while all three mice in GLY group showed more severe necrosis (numerical assessment numbers 3, 2, and 4) (Additional file 1: Fig. S1A). The reported pathological features of the liver after GLY exposure were cytoplasmic vacuolation, degeneration of hepatocytes, steatosis, multifocal necrosis, leukocyte infiltration, massive deposition of reticular fibers, and cytoplasmic glycogen deposition in hepatocytes [16, 34,35,36]. The histopathological effects of GLY on liver tissue in our study are in accordance with these studies using laboratory rats.

Fig. 1

The global cell types and changes in mice liver tissues delineated by scRNA-seq. A An overview of the study design. B The H&E staining of liver tissue (scale bar, 100 μm). Note the fibroblast proliferation in the Glisson’s capsule (black arrow), exudation (green arrow), slight inflammatory cell infiltration (yellow arrow), and focal necrosis (blue arrow) in GLY-treated mice. C UMAP plot displaying 13 major cell clusters based on 60,348 single-cell transcriptomes. D Violin plots showing the expression of two well-established marker genes for each cell type. E Distribution comparison of clusters from control (CON) and GLY groups. F Percent contribution of control (blue) and GLY (orange) mouse liver cells for each cell type. Note that the relative percentages of HSCs and LCMs were significantly increased by GLY. G The distributions of upregulated and downregulated genes for each cell type after GLY treatment. Note that HSCs and LCMs contained the largest number of DEGs. H The number of upregulated and downregulated genes for all cells after GLY treatment. Note that HSCs and LCMs contained the largest number of DEGs

Construction of mouse liver cell atlas upon GLY treatment

We performed droplet-based scRNA-seq (10x Genomics platform) of liver tissues from three control (CON) and three GLY-treated mice. After rigorous quality control, we obtained a total of 20,786 mouse genes from 60,348 qualified cells for subsequent analysis (Additional file 1: Fig. S1B). Of these, 27,733 cells (46%) were originated from GLY-treated samples and 32,615 cells (54%) were from control samples. The unsupervised clustering analysis partitioned all cells into 13 major cell lineages as displayed by uniform manifold approximation and projection (UMAP) plot (Fig. 1C and Additional file 1: Fig. S2A, B). The cell lineages were annotated by the expression of well-established markers, including endothelial cell (EC, expressing Oit3 and Clec4g), cholangiocyte (Chol, expressing Epcam and Sox9), hepatocyte (Hep, expressing Tat and Bhmt), hepatic stellate cell (HSC, expressing Dcn and Col3a1), Kupffer cell (KC, expressing Clec4f and Marco), liver capsular macrophage (LCM, expressing S100a4 and Cx3cr1), NK cell (NK, expressing Gzma and Nkg7), B lymphocyte (B lymph, expressing Cd79a and Ms4a1), T lymphocyte (T lymph, expressing Cd3d and Trbc2), plasmacytoid dendritic cell (pDC, expressing Siglech and Ccr9), conventional dendritic cell (cDC, expressing Xcr1 and Flt3), basophil cell (Baso, expressing Ms4a2 and Cap3), and neutrophil (Neutro, expressing Retnlg and S100a8) (Fig. 1D, Additional file 1: Fig. S2C and Additional file 2: Table S1).

Next, we quantified and compared the proportions of each cell type in control and GLY groups to reveal the effect of GLY treatment on cell composition. Overall, the healthy and GLY-treated mouse livers both contain the 13 cell clusters (Additional file 1: Fig. S2D). The relative percentages of HSCs and LCMs were significantly increased by GLY, whereas the proportions of ECs, hepatocytes, and lymphocytes were decreased compared with those in healthy liver tissues (Fig. 1E, F, Additional file 1: Fig. S2E and Additional file 2: Table S2). Transcriptomic analysis identified differentially expressed genes (DEGs) associated with GLY treatment in various cell types (Fig. 1G and Additional file 2: Table S3). Remarkably, HSCs and LCMs contained the largest number of DEGs. The number of DEGs from scRNA-seq analysis in HSCs and LCMs (Fig. 1G) is significantly higher than that of traditional analysis (Fig. 1H). These observations strongly suggest that GLY exerts hepatotoxicity effects by affecting specific cell (sub)types and corresponding genes.

GLY exposure induces HSC activation and proliferation

HSCs reside in the subendothelial space of Disse, interposed between hepatocytes and LSECs. Upon liver injury, nonproliferative quiescent HSCs transdifferentiate into proliferative, inflammatory, and contractile myofibroblasts with enhanced extracellular matrix (ECM) production (also known as “cell activation”) [37]. GLY treatment dramatically increased the relative proportion and absolute numbers of HSCs (Figs. 1E, 2E and Additional file 1: Fig. S3A–C), along with upregulation of key genes involved in extracellular matrix organization, collagen metabolic process, collagen fibril organization, and collagen biosynthetic process (Fig. 2A, B and Additional file 2: Table S5, 6).

Fig. 2

Transcriptomic changes of HSC activation after GLY stimulation. A Volcano plot of genes differentially expressed between control and GLY-treated HSCs. B GO analysis for upregulated DEGs of total HSCs after GLY treatment. C UMAP plot showing the cell subtypes of HSCs. D The expression of specific marker genes identified HSC subtypes. E Merged UMAP plot for control and GLY groups. F Violin plots of the expression for three genes associated with ECM. G Pseudotime trajectory analysis of HSC subtypes. H Proportion changes of each cell type for different states during the pseudotime analysis. I Typical images and bar graph showing the increased area percentage of Desmin after GLY treatment. Scale bar, 100 μm; Error bar, mean ± standard deviation (SD); P value is from Student’s t test (unpaired and two-tailed); ****P ≤ 0.0001; the quantification was based on six random areas from each group

Clustering analysis of HSCs in both control and GLY mice identified five distinct subpopulations (Fig. 2C)—one quiescent subtype qHSC and four activated subtypes composed of aHSC1-3 and aHSC_prof—using the markers shown in Fig. 2D. The qHSC cluster accounted for the largest proportion (~ 44%) in control mice but dropped to less than 10% after GLY treatment, whereas the relative proportions of activated subtypes, particularly aHSC_prof and aHSC3, were increased by GLY treatment (Fig. 2E, Additional file 1: Fig. S3C and Additional file 2: Table S4). The aHSC_prof cluster was defined by the high expression of proliferation markers such as Mki67 (Fig. 2D and Additional file 1: Fig. S3F), and aHSC3 was inflammation-related with strong expression of chemokines, including Ccl6, Ccl24, Cxcl9, and C3ar1 (Fig. 2D). DEG analysis identified that genes related to ECM (Col1a1, Col1a2, and Col3a1) were among the top ten upregulated ones after GLY exposure (Fig. 2F, and Additional file 2: Table S5). This finding suggested that GLY induces excessive deposition of fibrillary collagen, which may further activate HSCs and lead to liver fibrosis and other liver diseases [38].

To gain insight into the transdifferentiation fate of HSCs, pseudotime analysis was conducted to acquire the transcriptomic changes of HSCs during GLY exposure. We observed that the branch 1, primarily composed of quiescent HSCs from control group, developed into aHSC2-enriched branch 2 (cell fate 1) or aHSC3 and aHSC_prof-enriched branch 3 (cell fate 2), accompanied by other cell types (Fig. 2G, H and Additional file 1: Fig. S3D, E). Of note, nearly all the cells in branches 2 and 3 were activated HSCs from GLY-treated mice (Fig. S3D, E). The activation of HSCs is now well established as a key driver of fibrosis upon liver injury [39], and our data implied that GLY could induce HSC proliferation and activation. We next examined the expression of Desmin, a marker of HSC [40], in CON- and GLY-treated samples, and observed a significantly increased number of cells expressing Desmin (presented in area percentage, from 3.8% to 6.6%) upon GLY exposure which is consistent with our finding that GLY induced HSC proliferation (Figs. 1E, 2E, I).

Taken together, we observed activation of quiescent HSCs, including increase of proliferative and inflammatory HSCs, and enhancement of collagen genes, implying that GLY could induce liver damage.

Dysregulation of hepatocyte signaling induced by GLY

As GLY treatment reduced the proportion of hepatocytes (Fig. 1F and Additional file 1: Fig. S4A, B), we examined the subpopulations of this cell type. Hepatocytes were divided into three groups, namely Hep1, Hep2, and Hep3 (Fig. 3A, Additional file 1: Fig. S4C and Additional file 2: Table S7), using the specific genes listed in Fig. 3B. Functional enrichment revealed that different subtypes were highly involved in various processes (Fig. 3C and Additional file 2: Table S8). Hep1 was mainly involved in lipid metabolism such as lipid localization, steroid metabolic process, cholesterol metabolic process, and regulation of plasma lipoprotein particle levels. In Hep2, the enriched functions were focused on some catabolic processes such as organic acid catabolic process and carboxylic acid catabolic process. Hep3 participated in processes of vasculature development, tissue migration, and angiogenesis.

Fig. 3

Specific changes of hepatocyte subsets between control and GLY groups. A UMAP plot for the three subtypes of hepatocytes termed Hep1, Hep2, and Hep3. B The expression of specific marker genes in hepatocyte subtypes. C Heatmap of DEGs for the three cell subtypes and their corresponding functional enrichment analysis. D Pie plots showing proportional changes of each subtype within the control and GLY groups. E GO analysis for upregulated DEGs of Hep1 and Hep2. F Boxplots of ROS score for each cell subtype. The middle lines indicate median values; boxes range from the 25th to 75th percentile. P value is from wilcox.test (unpaired and two-tailed); ***P < 0.001, NS, not significant. G Heatmaps of upregulated genes for Hep1 and Hep2 after GLY treatment

After GLY exposure, Hep1 displayed an increased proportion while a decreased trend was noted for Hep2 (Fig. 3D and Additional file 1: Fig. S4D). Hep3 accounted for only 2% of the total population, and no substantial difference in its composition was observed after GLY treatment (Fig. 3D and Additional file 1: Fig. S4D). Next, we examined the DEGs in the Hep1 and Hep2 (Fig. 3G and Additional file 2: Table S9). ROS metabolic process and ROS score were elevated in Hep1 of GLY-treated mice, which was the major contribution to the increased ROS score observed for all the hepatocytes (Fig. 3E, F and Additional file 1: Fig. S4E), suggesting that GLY induced oxidative stress in this subtype. A similar effect of GLY in causing oxidative stress has been reported in the SH-SY5Y cell line [41] and rat liver [17]. Mesnage et al. observed activated oxidative stress (reflected by Srxn1 and Blvrb) in mammalian stem cell-based ToxTracker system by glyphosate-based herbicides but not glyphosate itself [23], In our study, the expression level of Srxn1 and Blvrb (Fig. 3G and Additional file 2: Table S9) was improved after GLY exposure. The contradictory observation could be caused by different testing models (in vitro cell line versus in vivo animal model) or administered dose variance (≤ 5 mM versus 200 mg/kg). The dose-dependent toxicity of glyphosate has been reported by some researchers [16, 17, 42, 43]. One of the main events, induced by ROS and oxidative stress, is lipid peroxidation or even oxidative damage of lipids [44, 45]. Thus, we observed that the genes involved in fatty acid metabolic process were markedly upregulated in Hep1 (Fig. 3E and Additional file 2: Table S10). Besides, the expression of genes related to inflammatory response also increased after GLY treatment (Fig. 3E and Additional file 2: Table S10), displaying consistency with former statements that GLY could induce inflammatory response in multiple organs, including liver, lung, and small intestine [17, 46, 47]. In addition, upregulation of genes involved in ribonucleoprotein complex assembly and ribonucleoprotein complex subunit organization in both Hep1 and Hep2 suggests potential DNA damage of GLY (Fig. 3E), which has been reported previously [23]. Overall, GLY may cause oxidative stress, abnormal lipid metabolism, inflammatory response, and DNA damage, suggesting GLY-induced dysregulation of hepatocytes.

Gene expression heterogeneity of KCs in GLY-treated mouse livers

KCs, the resident macrophages of the liver, reside within the hepatic sinusoid together with other innate immune cells, such as natural killer cells, natural killer T cells and dendritic cells [48]. On the basis of the specific marker genes, KCs were further clustered into five subtypes: KC1 (expressing Csf2rb, Cd14, and Fpr1), KC2 (Fcer1g, Lyz2, and Psme1), KC3 (Sdc3, Grn, and Jund), KC4 (Tpt1, S100a4, and Spp1) and KC5 (Itga1, Tek, and Pecam1) (Fig. 4A, B and Additional file 1: Fig. S5A–C), and GLY treatment resulted in various proportional changes in these subtypes (Fig. 4C, Additional file 1: Fig. S5D, E and Additional file 2: Table S11). Next, we analyzed the GLY-induced DEGs in each subtype (Fig. 4D). KCs, especially KC1, KC3, and KC5, displayed enhancements in ROS level and oxidative stress response after GLY exposure (Fig. 4D, E and Additional file 1: Fig. S5G). We also noticed that genes associated with phagocytosis, such as Macro and Cyba, were upregulated in KCs (Fig. 4D, F and Additional file 1: Fig. S5F), suggesting that KCs might be activated by GLY to maintain the immunological homeostasis. It has been reported that the expression of CD68 (a macrophage marker) was enhanced after macrophage activation [49], and the increased CD68 expression after GLY treatment was also observed in our work (Fig. 4H). Since CD68 is a general macrophage marker, we checked its expression in the 13 major cell clusters identified in our study and noted that it was highly expressed in KCs, LCMs, pDCs, and cDCs (Additional file 1: Fig. S5H). The improved average intensity of CD68 should be more correlated with KCs since we excluded LCMs by quantifying non-liver capsular areas, and no significant difference of CD68 expression level in pDCs and cDCs was observed after GLY treatment and their numbers were much less than KCs. Besides, GLY induced some KCs polarized into a more circular morphology, an anti-inflammatory M2 state that is generally believed to promote tumorigenesis and tumor progression (Fig. 4H) [50, 51]. Collectively, these results further indicated that GLY might activate macrophages, including KCs.

Fig. 4

Transcriptomic changes of KCs after GLY treatment. A UMAP visualization of the five subtypes for KCs. B Heatmap of the expression levels of representative marker genes in KC subtypes. C Histogram showing proportional changes of each KC subtype for two groups. D GO analysis for upregulated DEGs of different cell subtypes. E Boxplots of ROS score for each cell subtype. The middle lines indicate median values, boxes range from the 25th to 75th percentile. P value is from wilcox.test (unpaired and two-tailed); ***P < 0.001, *P < 0.05. F Split violin plots of the expression of four genes associated with phagocytosis. G Violin plots comparing the expression level of Spp1 in each subtype between control and GLY groups (left). GO analysis for upregulated DEGs of KC4 involving Spp1 gene (right). H Typical images and bar graph showing the increased expression level of CD68 after GLY treatment. Boxed regions are amplified on the right; scale bar, 100 μm; error bar, mean ± SD; P value is from Student’s t-test (unpaired and two-tailed); ****P ≤ 0.0001; the quantification was based on six random areas from each group. I Heatmap of upregulated and downregulated genes for KC4 in control and GLY samples. J Differentially expressed cytokines for each KC subtype after GLY treatment. K Heatmap for common upregulated and downregulated genes of all KC clusters

By scrutinizing the DEGs in KCs, we noted that the GLY treatment remarkably increased the expression level of Spp1 in subtype KC4 compared with the control group (Fig. 4G, I and Additional file 2: Table S12). Spp1, which encodes the secreted fibrogenic factor osteopontin (OPN), has been reported to be upregulated in liver fibrosis [52]. OPN has been found to activate HSCs [53,54,55] and be involved in some physiological cellular functions, including migration, adhesion, and secretion (Fig. 4G) [56, 57]. The expression of cytokines that were significantly affected by GLY is shown in Fig. 4J (Additional file 2: Table S13). Moreover, other liver-fibrosis-related genes, such as Trem2 [58] and Tgfbi [59], also exhibited elevated expression levels (Fig. 4K and Additional file 2: Table S14). Thus, the results implied that GLY may potentially induce liver fibrosis.

Activation of LCMs with GLY treatment

Three LCM subsets were revealed by unique transcriptomic signatures and visualized in UMAP plot, termed LCM1, LCM2, and LCM3 (Fig. 5A and Additional file 1: Fig. S6A, B). We annotated the subclusters with multiple DEGs and analyzed functional enrichment of these genes (Fig. 5B and Additional file 1: Fig. S6D). For instance, highly expressed genes in LCM1, including Ly6c2, Clec4e, and Vcan, were involved in positive regulation of cytokine production, regulation of inflammatory response, and phagocytosis. LCM2 strongly expressed Adgre1, Cd63, Mrc1, and Trem2, genes involved in cell chemotaxis, macrophage migration, and activation. Antigen processing and presentation genes (H2-Eb1, H2-Ab1, and H2-Aa) were prominently expressed by LCM3. Overall, the number of LCMs was significantly increased after GLY treatment (Fig. 5C and Additional file 2: Table S15). After GLY exposure, LCM2 was the most increased cluster among the three, as shown by either absolute number or relative percentage of cells (Fig. 5D and Additional file 1: Fig. S6C). The increased number of LCMs after GLY treatment was confirmed by CD68 staining (Fig. 5J). Since the LCMs are located in the hepatic capsule, we co-stained the LCM marker Cx3cr1 with CD68 to differentiate LCMs from other macrophages and found that the number of LCMs (displayed by area percentage) was significantly increased from 3.8% to 32.4% upon GLY exposure (Fig. 5J). GLY treatment induced upregulation of unique genes tightly related to oxidative phosphorylation, phagosome, endocytosis, lysosome, and ROS in LCM2 (Fig. 5F and Additional file 1: Fig. S6E). There are also several genes and pathways commonly dysregulated after GLY treatment across different subclusters (Fig. 5E and Additional file 2: Table S16). Notably, LCM1 and LCM2 showed increased expression of genes (Trem2, Fcgr1, and Fcgr3) related to inflammatory response in GLY group (Fig. 5G and Additional file 1: Fig. S6G). The expression levels of ROS metabolic process genes, including Ncf1, Ncf4, Gpx1, and Gpx4, were enhanced by GLY in all subclusters (Fig. 5H and Additional file 1: Fig. S6H). More importantly, GLY treatment induced the high expression of Spp1, which was closely associated with liver fibrosis and nonalcoholic steatohepatitis (NASH) (Fig. 5I and Additional file 1: Fig. S6F). Tgfbi, the gene encoding a potent profibrotic cytokine, was also overexpressed after GLY treatment among LCM subclusters (Fig. 5I). Meanwhile, the expression of several chemotactic cytokines was significantly increased by GLY in all LCM subsets (Additional file 1: Fig. S6I).

Fig. 5

Inflammatory response and pro-fibrosis of LCMs after GLY treatment. A UMAP visualization of the three subtypes for LCMs. B The expression of specific marker genes identified the subtypes of LCM. C The distribution of each cell type in control and GLY groups. D Pie plots showing proportional changes of each cell type for two groups. E GO analysis for upregulated DEGs of different cell subtypes. F Venn diagram displaying the numbers of upregulated DEGs in the three subtypes (left). GO analysis for unique upregulated DEGs of LCM2 after GLY treatment (right). G Violin plots showing the expression levels of inflammatory response-related genes in each subtype. H Violin plots showing the expression levels of ROS metabolic process-related genes in each subtype. I Violin plots showing the expression levels of Spp1 and Tgfbi genes in each subtype. P value is from wilcox.test (unpaired and two-tailed); **P < 0.001; NS, not significant. J Typical images and bar graph showing the increased area percentage of CD68 in the liver capsule after GLY treatment. The liver capsule outlined by yellow-dotted lines was denoted by Cx3cr1 staining; scale bar, 100 μm; error bar, mean ± SD; P value is from Student’s t-test (unpaired and two-tailed); ****P ≤ 0.0001; the quantification was based on eight random areas from each group

Intercellular communication in crucial subtypes

The above findings have demonstrated that certain cell types in liver tissue underwent significant proportional and transcriptional alterations following GLY treatment, particularly LCMs, KCs, and HSCs. These aberrations may induce massive and complex changes of cell–cell communication networks in response to GLY stimulation. To disclose the underlying signaling patterns between cells, the statistical and biological communication networks were constructed for significant ligand–receptor pairs using CellChat program (Additional file 1: Fig. S7A). By comparing the detailed signaling features of the two groups, we identified that GLY treatment caused increased counts of interaction, as well as enhancement of interaction strength (Additional file 1: Figs. S7A, 8A). It is prominent that the interaction strength between LCMs and KCs increased significantly in the GLY group (Additional file 1: Figs. S7B, 8B). Our transcriptional analysis of LCMs has found that GLY increased the expression of Spp1 (Additional file 1: Fig. S6I), which could induce the activation and infiltration of macrophages and generate pro-inflammatory environment through binding to integrin [60, 61]. Further ligand–receptor pair analysis on this pathway demonstrated that LCMs are the source of SPP1 ligands that act on LCMs, DCs, and KCs after GLY exposure (Additional file 1