Single-cell metabolic profiling reveals subgroups of primary human hepatocytes with heterogeneous responses to drug challenge

Abstract

Background

Xenobiotics are primarily metabolized by hepatocytes in the liver, and primary human hepatocytes are the gold standard model for the assessment of drug efficacy, safety, and toxicity in the early phases of drug development. Recent advances in single-cell genomics demonstrate liver zonation and ploidy as main drivers of cellular heterogeneity. However, little is known about the impact of hepatocyte specialization on liver function upon metabolic challenge, including hepatic metabolism, detoxification, and protein synthesis.

Results

Here, we investigate the metabolic capacity of individual human hepatocytes in vitro. We assess how chronic accumulation of lipids enhances cellular heterogeneity and impairs the metabolisms of drugs. Using a phenotyping five-probe cocktail, we identify four functional subgroups of hepatocytes responding differently to drug challenge and fatty acid accumulation. These four subgroups display differential gene expression profiles upon cocktail treatment and xenobiotic metabolism-related specialization. Notably, intracellular fat accumulation leads to increased transcriptional variability and diminishes the drug-related metabolic capacity of hepatocytes.

Conclusions

Our results demonstrate that, upon a metabolic challenge such as exposure to drugs or intracellular fat accumulation, hepatocyte subgroups display different and heterogeneous transcriptional responses.

Background

Single-cell sequencing technologies have revealed a wealth of information regarding cellular heterogeneity in health and disease in multiple tissues [1,2,3,4,5,6,7,8]. In particular, single-cell RNA sequencing enables the identification of groups of cells in a population with similar transcriptomic profiles, which is generally associated with similar functionality [9, 10]. In the liver, single-cell transcriptomic analyses have shown novel cell types and states involved in the development and progression of liver disease [1, 7, 11], as well as in the transcriptomic responses to xenobiotics [12]. However, drug toxicity and safety are generally assessed in primary human hepatocytes (PHH) as the gold standard model to study the drug metabolism in humans, which are widely considered a seemingly homogeneous population of cells [13, 14]. Whether all hepatocytes share the same functional molecular phenotype in vitro remains unknown.

In the liver, drug and xenobiotic metabolism is mainly performed by hepatocytes, the major and predominant cell type of the parenchyma [15,16,17]. The hepatic metabolism of drugs occurs in three phases. In phase I, oxidation, hydrolysis, reduction, and cyclization reactions are catabolized mainly by the cytochrome P450 (CYP450) superfamily of monooxygenase enzymes. The main members are the isoforms CYP1A2, 2C9, 2C19, 2D6, and 3A4, which account for the metabolism of around 70–80% of the clinically available drugs [18,19,20]. The expression of these cytochromes is induced by the presence of their substrate compounds, which are indirectly used as a measure of liver metabolic capacity [10, 21]. These substrates act as hepatocyte probes when administered to phenotype and monitor the cytochrome enzymatic activity, formally known as the “phenotyping cocktail approach” [22,23,24,25,26,27,28]. Phase II comprises conjugation reactions catabolized mainly by transferase enzymes and are purposed to hydrophilize compounds for their elimination [29]. During phase III of xenobiotic metabolism, conjugated compounds are excreted out of the cell through active transmembrane transporters [30]. The cytochrome superfamily consists of nearly 60 members in humans (Human Genome Project 2013), which may be expressed differently in individual cells when challenged by xenobiotics, leading to cellular heterogeneity that remains concealed in a bulk analysis.

An additional potential source of cellular heterogeneity is the intracellular accumulation of triglycerides in hepatocytes, known as hepatic steatosis [31, 32]. This is a hallmark of non-alcoholic fatty liver disease (NAFLD) and is generally associated with metabolic dysfunction, inflammation, and risk of fibrosis [33, 34]. In culture, intracytoplasmic fat accumulation can be modeled by incubating primary human hepatocytes with free fatty acids (FFA) to mimic benign chronic steatosis [35, 36]. However, lipid accumulation is heterogenous in regards to the number and size of lipid droplets [37] and it still remains unclear whether all cells respond to lipid accumulation in a coordinated manner.

Moreover, the co-administration of five or more drugs (polypharmacy) is associated with a higher incidence of adverse drug reactions (ADR) and drug-induced liver injury (DILI) [38,39,40]. Importantly, a higher incidence of DILI has been reported in patients suffering from NAFLD [41,42,43]. Therefore, the precise molecular pathways commonly dysregulated between chronic accumulation of lipid and drug metabolism remain unexplored at cellular resolution.

Here, we report concealed cellular heterogeneity in primary human hepatocytes which are classically considered a seemingly homogenous population of cells and a “gold standard” for in vitro hepatic assays. The comparison between pseudobulk and single-cell transcriptomics reveals four distinct hepatocyte subgroups that can be found also in fourteen human donors from publicly available data sets. We present evidence that chronic accumulation of lipids and xenobiotics leads to transcriptional variability and the impairment of multiple metabolic pathways in a subgroup-specific manner. Our results suggest that liver steatosis combined with multiple drug intake differentially affects the transcriptional profiles of individual cells decreasing the expression of key cytochrome P450 genes.

Results

Single-cell RNA-seq reveals four major subgroups of hepatocytes showing cellular heterogeneity and functional specialization in primary human hepatocytes

Liver function is compartmentalized along the liver lobule, and the maintenance of liver function is regulated by liver-enriched transcription factors [44, 45]. Primary human hepatocytes are considered the gold standard model to predict drug responses in humans, but they are characterized by phenotypic instability in culture and rapid loss of the hepatic phenotype [13, 46]. Therefore, we aimed to investigate the level of cellular heterogeneity that remains in a seemingly homogeneous population of primary human hepatocytes (PHHs), and how this heterogeneity might affect pharmaco-toxicological studies [47].

Cryopreserved and metabolically competent PHHs from four donors were used as in vitro model to characterize the metabolic profile of individual hepatocytes in response to a drug challenge and chronic accumulation of fat (Fig. 1A, “Methods”). The hepatocytes were plated for 6 h and incubated for 66 h with either vehicle (DMSO) or a five-drug cocktail. Chronic accumulation of fat was achieved by incubation with free fatty acids (FFA) corresponding to a 200 µM mixture of a 2:1 ratio of oleic acid to palmitic acid for 72 h. In brief, four different conditions were studied: (i) vehicle (DMSO 0.5% v/v); (ii) Cocktail (66 h five-drug cocktail incubation); (iii) FFA (2:1 oleic:palmitic acid), and (iv) FFA + Cocktail (combination of FFA incubation and five-drug cocktail) (Fig. 1A). Donors were healthy males between 18 and 57 years, with a body mass index (BMI) classified as normal, non-diabetic, and representing the most common age range commercially available (Additional file 2: Table S1). After isolation of viable cells by a quick and non-harsh dead cell removal step, single-cell RNA-seq was immediately performed using 10x  Genomics (Methods).

Fig. 1

Transcriptomic profiling reveals four subgroups of PHHs independently of donor and treatment condition. A Overview of experimental design. Purified cryopreserved PHHs from four human donors were plated, incubated with or without free fatty acids to model hepatic fat accumulation and with or without a phenotypic 5-drug cocktail (Sanofi-Aventis). B UMAP colored by subgroup, treatment, and donor showing that the annotated subgroups are found throughout all donors and conditions. C UMAP colored by expression levels of—from left to right: (i) key hepatocyte marker genes, (ii) bile acid and sterol metabolism; (iii) carbohydrate and phase II metabolism, (iv) lipid and phase III metabolism marker genes in four subgroups of hepatocytes. D Boxplot showing the percentage of cells in which a given gene is expressed, colored by the identified subgroups. E Bar plot showing the percentage of cells per subgroup assigned to phases G1, S, and G2M by performing cell cycle analysis using cyclone. F Dot plot highlighting marker gene expression in four subgroups of hepatocytes identified in vitro (top) and in vivo [1]. Subgroup marker gene expression was grouped by aggregated Louvain clusters (dot size: fraction of cells in the group; color scale: mean expression in the group)

A total of 38,232 high-quality hepatocytes with an average of 2550 genes per cell were analyzed. Lenient doublet filtering was performed to account for polyploid hepatocytes (Methods). The variation in the number of cells profiled per donor was attributed in part to sample viability, technical differences in cell capture rates in each single-cell RNA-sequencing (scRNA-seq) run, and sequencing depth (Additional file 1: Fig. S1, Methods).

Across all four donors, Louvain clustering led to the detection of heterogeneity separating four subgroups of hepatocytes, independently of the treatment (Fig. 1B, Additional file 1: Fig. S1B and C, Methods). Harmony [48] was used to correct for batch effects, and Louvain clustering was performed on the combined data set (Additional file 1: Fig. S2A). Overall, hepatocyte purity was confirmed by the expression of the hepatocyte-specific genes ALB, SERPINA1, and TTR, among others (Fig. 1C, and Additional file 3: Table S2). Differential expression analysis was used to annotate the four subgroups of hepatocytes present in all four donors and treatment conditions (Fig. 1B, Additional file 3: Table S2, Methods). Representative marker genes illustrate functional specialization. For instance, ATF6, CYP8B1, and HMGCS2 for sterol and bile acid synthesis were up-regulated in subgroup I. Subgroup II was represented by LDHA, GSTZ1, and GSTO1, involved in the carbohydrate and phase II metabolism. Subgroup III was characterized by ABCC2, ABCC3, and APOB over-expression, among other genes (Additional file 1: Fig. S2B) involved in lipid and phase III metabolism (Fig. 1C). Lastly, subgroup IV was characterized by the loss of gene expression, which was 11-fold reduced compared to subgroup I (Fig. 1D, “Methods”). The instability of PHHs has been previously characterized for the loss of expression of liver-specific transcription factors and downstream target genes [49, 50].

Across all conditions, 38.2% of the cells were specialized in the metabolism of bile acids and sterols (subgroup I) [51]; 38.7% of the cells were enriched in genes involved in carbohydrate metabolism and phase II enzymes (subgroup II); 5.4% of the cells were responsible for lipid metabolism and expression of phase III enzymes (subgroup III) (Additional file 1: Fig. S2C); and 17.7% of the cells were losing expression (subgroup IV) (Additional file 1: Fig. S2C). To further investigate the metabolic specialization of the PHH subgroups, we focused our analysis on metabolic marker genes, taken from metabolicatlas.org (Methods). We found that clustering analysis using metabolic markers alone led to the identification of the same subgroups of hepatocytes as when using all genes (Additional file 1: Fig. S2E and 2F). These results suggest that differences between subgroups of hepatocytes might be driven by their metabolic function.

Interestingly, in subgroup III, the majority of the cells (73.2%) were assigned to S-phase according to the cell cycle phase classification using the tool cyclone (Fig. 1E, Additional file 1: Fig. S2A and S3) [52]. Moreover, subgroup II showed a higher expression level of stress marker genes such as RSP19, PRDX1, BAX, GSTA1, LGALS1, MTH1, and MTHM (Additional file 1: Fig. S2D, Additional file 3: Table S2). Likewise, subgroup II showed a low percentage of cells in which a gene is expressed, suggesting that these cells might be prone to lose their characteristic hepatocyte-like expression profile in culture.

To further characterize underlying upstream molecular events, we used the ChIP-X Enrichment Analysis tool, ChEA3 [53] to reconstruct the network of putative transcription factors (TFs) regulating gene expression. We focused our analysis on the untreated cells (i.e., DMSO-treated cells) across all donors, to prevent potential effects due to treatment conditions. Considering the top 500 differentially expressed genes per subgroup, key hepatic TFs, such as HNF4A [54, 55] and MLXIPL [56] were found among the top 25 predicted TFs. Only in subgroup IV a decrease in the expression levels of key transcription factors was detected (Additional file 1: Fig. S4A, Additional file 4:Table S3), while subgroups I, II, and III were defined by a unique combination of master regulators (Additional file 1: Fig. S4A, Additional file 4:Table S3). To validate the specific TF activity per subgroup, we performed single-cell assay for transposase-accessible chromatin sequencing (scATAC-seq) on PHH incubated only with DMSO (Additional file 1: Fig. S4B). After scATAC-seq in one donor, we found that the promoter regions of subgroup-defining marker genes were open and accessible to the TFs predicted by ChEA3 in DMSO-treated cell (Additional file 1: Fig. S4, Methods).

Primary human hepatocytes retain functional specialization in vitro in the absence of liver zonation

Among the four subgroups of hepatocytes identified in vitro, only subgroups I, II, and III were considered functional and metabolically active. We then further investigated whether these three hepatocyte subgroups could also be identified in vivo, using two distinct publicly available datasets [1, 5].

First, we investigated nine human livers described in Aizarani et al. in 2019 [1]. Hepatocytes were extracted from the data set based on the expression of mature hepatocyte maker genes ALB, TTR, and HNF4A. After performing Louvain clustering on the in vivo hepatocytes (resolution 0.2), we used the transcriptional profiles of our subgroups (I, II, and III) to identify hepatocyte specialization per cluster. We observed that the transcriptional profiles defining our subgroups are similarly expressed in hepatocytes in vivo with a correlation up to 0.94 for the top ten differentially expressed genes (DEGs) per subgroup (Fig. 1F and Additional file 1: Fig. S5).

In vivo, expression profiles of hepatocytes are affected by their position in the liver lobule [7, 57, 58]. Therefore, we investigated whether the three hepatocyte subgroups identified as subgroups I, II, and III in vivo were related to liver zonation patterns (Additional file 1: Fig. S5A). Based on their expression profile, cells were scored for zonation marker genes (Additional file 1: Fig. S5B, Methods) and assigned to three zones: pericentral, mid-zone, and periportal (Additional file 1: Fig. S5A and 5B). In subgroup I, 64.6% of cells were assigned to the periportal area, 6.6% to the midzone, and 46.5% to the pericentral area. The majority of cells in subgroup II, 52.3%, were assigned to the midzonal area, with 12.1% of cells in the pericentral area and 27.9% in the periportal area. Subgroup III showed an enrichment in mid- and pericentral expression profiles, 41.1% and 41.4%, respectively, and 7.5% in the periportal area (Additional file 1: Fig. S5B and 5C, Methods).

Additionally, hepatocyte subgroups and zonation were found to be intertwined in vivo (Additional file 1: Fig. S5D, left). For instance, CYP27A1 in subgroup I was pericentrally enriched [59,60,61], while HSD11B1, involved in bile synthesis [62], was distributed periportally (Additional file 2: Table S1; Additional file 1: Fig. S5, left). Remarkably, in vitro, PHH retained specialization into the three functional subgroups (I, II, and III) in the absence of zonation (Additional file 2: Table S1; Additional file 1: Fig. S5, right). In subgroup I, the CYP27A1 and HSD11B1 had similar expression levels in all the cells. Similarly, zonation marker genes in subgroups II and III were evenly expressed in all cells studied.

These findings were extended to additional publicly available data sets for five human donors [5] (Additional file 1: Fig. S5E and 5F, Methods). In these five donors, MacParland et al. annotated six hepatocyte clusters, named Hep 1 to Hep 6 [5]. Based on marker gene overlaps, we found similarities in the gene expression profiles between those clusters and our hepatocyte subgroup in vitro.

In order to increase the power of the in vivo and in vitro comparisons, we computationally integrated the two in vivo datasets, resulting in a human cohort of 14 donors (Additional file 1: Fig. S5F). The identified metabolic subgroups showed high gene expression correlation between in vivo and in vitro (Additional file 1: Fig. S5G). In summary, our results using publicly available data sets from human donors indicate that PHH in vitro retained functional specialization in the absence of liver zonation.

Phenotyping cocktail used to assess the induction of cytochrome P450 shows differential metabolic profiles among hepatocyte subgroups

In order to study the impact of cellular heterogeneity on liver function, we further characterized the hepatocyte subgroups by assessing their response to a phenotyping cocktail. The Sanofi-Aventis five-drug cocktail was used to simultaneously monitor the expression levels of the main five cytochrome P450 genes as a readout of the metabolic capacity of primary human hepatocytes [25, 63]. Therefore, this cocktail, consisting of a mixture of individual selective substrates of CYP2D6 (metoprolol), CYP2C19 (omeprazole), CYP2C9 (S-warfarin), CYP3A (midazolam), and CYP1A2 (caffeine), was used to monitor CYP induction [64] and incubated for 66 h in primary human hepatocytes (Methods).

Incubation with the five-drug cocktail did not lead to substantial differences in the number of captured cells in each condition per subgroup (Additional file 1: Fig. S6A). Upon incubation, the induction of the mRNA of the five cytochrome P450 genes (CYPs) involved in the metabolism of the five-drug cocktail was monitored showing their up-regulation in pseudobulk (Fig. 2A left, Additional file 1: Fig. S6B). However, differences in the basal transcriptional levels were observed between subgroups showing transcriptional heterogeneity in a seemingly homogenous population of PHH (Fig. 2A, right). While the pseudobulk shows the expected mRNA basal levels of CYP2C9 compared to CYP1A2 and CYP2D6 [65], scRNA-seq identified subgroups of hepatocytes with different basal mRNA levels of CYPs. Additionally, a subgroup-specific upregulation of CYPs upon cocktail incubation was also detected. For instance, in pseudobulk, CYP2C9 was up-regulated 1.7-fold in Cocktail, while in subgroup III only a 1.1-fold change was detected. Likewise, CYP3A4 was not significantly up-regulated in subgroup III, but a 4-fold increase was detected in subgroup I (Fig. 2A, right).

Fig. 2

Subgroups of PHHs show different metabolic profile in response to 5-drug cocktail. A Violin plots depicting expression levels between Cocktail and DMSO of the five CYPs involved in the metabolism of the five-drug cocktail in pseudobulk (left), and in each hepatocyte subgroup (right) (* = p-value < 0.05 and |log₂-fold change|> 1, t-test). B Volcano plot depicting the differential expression between Cocktail (green) and DMSO (purple) in pseudobulk. Text highlights the genes identified as DEGs in all subgroups. C Dot plot showing log₂-fold change (color scale) and p-value (dot size) between Cocktail and DMSO for genes significantly up-regulated in all subgroups (left); representative subgroup-specific up-regulated genes (middle); and genes significantly down-regulated in all subgroups (right). D Venn diagram showing overlaps of significantly up-regulated genes upon cocktail treatment between the subgroups. E Scatter plot depicting enrichment of the genes specifically up-regulated in the metabolically active hepatocyte subgroups I, II, and III in pathways known to be involved in the metabolism of given chemical compounds (Drug.CTD database). The size of the dot corresponds to the number of overlapping genes in a given pathway

Thus, we wondered whether the subgroup of hepatocytes shared global transcriptional signatures in response to xenobiotics. We first analyzed global changes upon five-drug cocktail incubation, detecting 161 significantly up-regulated genes compared to DMSO (Fig. 2B and Additional file 1: Fig. S6C). All genes with a corrected p-value below 0.05 and a log₂ fold change greater than 1 were designated as significantly up-regulated in cocktail-treated cells. Genes with log₂ fold change below − 1 were considered significantly up-regulated in DMSO-treated cells (Fig. 2B and C, Additional file 1: Fig. S6D, Additional file 6: Table S5, Methods). Additionally, removing hepatocyte subgroup IV (characterized by the loss of expression) leads to the detection of 205 significantly up-regulated genes (Additional file 1: Fig. S6C, Additional file 5: Table S4).

Investigating the shared and subgroup-specific signatures in response to five-drug cocktail, only eight genes were significantly up-regulated upon cocktail treatment in all four subgroups: CYP1B1, POR, CYP1A1, CYP1A2, RGS9, GDF15, CYP2A7, and CYP2B6 (Fig. 2B and C). These genes were also detected as significantly up-regulated in the pseudobulk analysis and correspond to genes involved in the phase I metabolism of xenobiotics [20, 21].

We focused our attention on subgroup-specific DEGs that might be concealed by studies performed in bulk (Fig. 2C middle). We found subgroup-specific upregulated genes that were not detected as significantly up-regulated in the pseudobulk analysis. For example, upon five-drug cocktail incubation, hepatocytes in subgroup I (bile acid and sterol metabolism) specifically up-regulated ATF3 and SRD5A2; subgroup II (carbohydrate and phase II metabolism) specifically up-regulated CYP2U1 and SLC4A7; and subgroup III (lipids and phase III metabolism) specifically up-regulated PLIN2 and OSGIN (Fig. 2C middle). A similar number of genes were specifically up-regulated in every metabolically active hepatocyte subgroup: 122 genes in subgroup I, 102 genes in subgroup II, and 126 genes in subgroup III. Meanwhile, only 64 genes were specifically up-regulated in subgroup IV (Fig. 2D). These results revealed subgroup-specific transcriptional signatures in response to xenobiotics that might be undetected in bulk studies.

Additionally, gene ontology (GO) analysis of subgroup-specific DEGs was performed using the Comparative Toxicogenomics Database (CTD) to explore toxicological interactions. This database is particularly suited for drug-disease or drug-phenotype interactions [66]. GO analysis using CTD showed that each hepatocyte subgroup specialized in the metabolism of certain xenobiotics, based on their differential transcriptomic profile (Fig. 2E, Additional file 1: Fig. S6E). For instance, the metabolic pathway of abrine was up-regulated across all hepatocyte subgroups I, II, and III. However, some compounds were only enriched in two subgroups, like cisplatin, in subgroups I and II. Finally, the pathways for the metabolism of other compounds were only enriched in one subgroup, like of ciglitazone and aflatoxin B1, enriched in subgroups II and I, respectively.

In summary, these results indicate that upon five-drug cocktail treatment, hepatocyte subgroups showed differential transcriptional responses, characterized by shared metabolic pathways and subgroup-specific transcriptional profiles associated with different potentials for metabolizing endobiotic (endogenous) and xenobiotic (exogenous) chemical compounds.

Intracellular lipid accumulation leads to differential transcriptional variability among hepatocyte subgroups

Recently, it has been shown that changes in cytochrome P450 activity correlate with altered lipid metabolism [67,68,69]. For example, hepatic steatosis affects the transcriptomic profile of parenchymal and non-parenchymal cells as well as the cellular composition in the liver [32, 70, 71]. Particularly in hepatocytes, the lipid metabolism is disrupted upon fat accumulation through an alteration of key enzymes in the lipid synthesis, storage, and clearance pathways [72]. Furthermore, an increase in the production of chemokines has been observed, associated to the inflammation occurring in NAFLD [73].

To shed light on the effect of lipid accumulation on the metabolic capacity of functional subgroups of PHH, hepatic steatosis in vitro was modeled by incubating the cells with a 200-µM mixture of a 2:1 ratio of oleic acid to palmitic acid [35] (Methods). This mixture has previously been shown to mimic benign chronic steatosis with minor lipotoxic and apoptotic effects [35, 36]. First, we investigated whether — and to what extent — fat accumulation triggers a coordinated transcriptional response in PHHs. Thus, the coefficient of variation for DMSO- and FFA-treated cells was calculated per subgroup. Overall, cells losing expression (subgroup IV) showed the highest transcriptional variability (Fig. 3A, Additional file 1: Fig. S7A). Moreover, lipid accumulation significantly increased the variability in functional subgroups I and II, but decreased it in subgroup III (Fig. 3A). This indicated that subgroup III showed a more coordinated response towards accumulation of lipids.

Fig. 3

Intracellular lipid accumulation increases loss of expression and transcriptional variability. A Box plots depicting the coefficient of variation per gene in cells treated with DMSO or with FFA per subgroup (* = p-value < 0.05, Mann–Whitney U). B Bar plot showing the percentage of cells treated with DMSO or with FFA in every subgroup. C Stacked violin plots of the genes used in Fig. 1 to identify the subgroups, split into DMSO- and FFA-treated cells (top) and top 5 up-regulated genes upon fat accumulation per subgroup (bottom) in DMSO and FFA treatments (* = p-value < 0.05 and |log₂-fold change|> 0.75, t-test). D Heatmap depicting the logarithmic mean expression of genes related to lipid metabolism, lipid storage, NAFLD-related genes, inflammation, and stress response in DMSO- and FFA-treated cells per subgroup (↑ indicates up-regulation towards DMSO; ↓ indicates down-regulation towards DMSO, t-test). E Scatter plot of the gene ontology (GO) analyses using gprofiler of the genes up-regulated upon fat accumulation in each of the subgroups. The top 7 most significantly enriched terms are depicted

Secondly, we observed that the proportion of cells treated with FFA differed between the functional hepatocyte subgroups. More than 74% of the cells in subgroup IV were FFA-treated cells, while in subgroup III 24% of the cells were FFA-treated cells. A similar percentage of FFA-treated cells were found in subgroups I and II (Fig. 3B, Additional file 1: Fig. S7B).

To explore global changes triggered by lipid accumulation, differential expression analysis was performed between FFA- and DMSO-treated cells. Compared to cocktail-treated cells, fewer up-regulated genes (3.5 times) were detected in FFA-treated cells over DMSO expression level (Additional file 1: Fig. S7C). Thus, genes with a corrected p-value below 0.05 and a log₂ fold change greater than 0.75 were designated as significantly up-regulated in FFA-treated cells. Genes with log₂ fold change below − 0.75 were considered significantly up-regulated in DMSO-treated cells (Additional file 1: Fig. S7C, Methods).

Thirdly, to track down the impact of intracellular fat accumulation on cellular identity, transcriptional changes in the previously selected marker genes — defining the hepatocyte subgroups — were investigated (Figs. 1F and 3C top). No significant changes in mean expression were observed for most of these marker genes upon accumulation of fat. Only HMGCS2, a key enzyme responsible for the synthesis of ketone bodies [74, 75] was significantly up-regulated in subgroup I upon fat accumulation (Fig. 3C, top).

In the metabolically active subgroups (I–III), the top five significantly up-regulated genes under fat accumulation were known markers of lipid droplet formation and lipid metabolism (Fig. 3C, bottom). For instance, the lipid droplet-associated perilipin protein PLIN2 was significantly up-regulated in all subgroups, which has previously been shown to be relevant in diet-induced NAFLD [76, 77]. The inflammation marker TNFAIP3, also discovered to ameliorate NAFLD and be protective against its progression [78, 79], was up-regulated among the five top DEG genes in subgroup I (Fig. 3C).

Subsequent analysis of fat-metabolism-related pathways and genes involved in stress response [80] showed subgroup-specific signatures upon intracellular fat accumulation (Fig. 3D, and Additional file 3: Table S2 and Additional file 5: Table S4). For example, CYP4A11, involved in NAFLD progression by inducing ROS-related lipid peroxidation and inflammation [81, 82], was significantly up-regulated in subgroups II and III; and CIDEC, a promotor of triglyceride accumulation [83, 84], was significantly up-regulated in subgroups I, II, and III. Cells losing expression (IV), up-regulated LGALS1 [85, 86], and HSPB1 [87, 88] (Fig. 3D and Additional file 1: Fig. S6D), which are known markers in the gene ontology pathway GO:0006950: “response to stress.”

To identify the main biological processes associated to each metabolically active subgroup upon fat accumulation, gene ontology (GO) analyses were performed using all the significantly up-regulated genes per subgroup. Subgroup I showed gene overlaps in pathways related to cellular response to lipids, together with lipopolysaccharide and chemokine metabolism [89] (Fig. 3E). For instance, chemokines CXCL8, CXCL1, CXCL10, and CXCL11 were up-regulated in FFA condition (Additional file 8: Table S7), suggesting that lipid accumulation could lead to increased inflammation through chemokine signaling [90, 91]. Moreover, subgroup II exhibited a high overlap of genes involved in the regulation of triglyceride metabolic processes, as well as in the acylglycerol catabolic process, denoting that these cells were rather involved in the clearance of neutral lipids [92, 93]. Finally, subgroup III cells showed enrichment in lipid, monocarboxylic acid, and fatty acid-related metabolic processes and lower transcriptional variability, most likely due to their coordinated response to fat accumulation [94,95,96].

Taken together, in metabolically active hepatocyte subgroups I and II, intracellular lipid accumulation increased transcriptional variability, thus affecting the fine-tuned regulation of lipid metabolism. Importantly, in subgroup III, specialized in the metabolism of lipids, transcriptional variability was reduced upon fat accumulation, suggesting a robust and tight coordinated response to chronic accumulation of lipids.

Intracellular lipid accumulation impairs drug metabolism phases I, II, and III, with concomitant up-regulation of stress-related pathways

Large inter-individual variability has been observed in the metabolism of CYP substrates in vivo. This variability might be caused by genetic, epigenetic, and environmental factors. For instance, the simultaneous administration of five or more drugs, known as polypharmacy, is highly common in the clinical practice [38]. The co-administration of drugs increases the risk for developing hepatotoxicity and adverse drug reactions (ADRs), such as drug-induced liver injury (DILI) [97, 98]. Additionally, a higher incidence of DILI has been reported in patients suffering from NAFLD [41,42,43]. Therefore, we assessed the impact of fat accumulation in the phase I metabolism of drugs, using the previously characterized phenotyping cocktail (Sanofi-Aventis) (Fig. 2). The expression of the individual CYPs targeted by the drug cocktail was analyzed, comparing Cocktail- and FFA + Cocktail-treated cells (Fig. 4A). In all subgroups, fat accumulation decreased the expression levels of the five targeted cytochromes (Fig. 4A). For instance, in subgroups I, II, and IV, CYP3A4 was significantly up-regulated upon cocktail treatment, but its induction was attenuated upon chronic lipid exposure. In subgroup III, the largest magnitude change of CYP expression was observed. Notably, CYP2D6, CYP2C19, CYP2C9, and CYP3A4 were significantly down-regulated in comparison to baseline DMSO levels upon FFA + Cocktail treatment.

Fig. 4

Intracellular lipid accumulation impairs drug metabolism phases I, II, and III. A Scatter plot depicting the log₂-fold change of the 5 induced cytochromes between cocktail and DMSO (green), and between FFA + Cocktail and DMSO (red) in each subgroup (* = p-value < 0.05 and |log₂-fold change|> 1, t-test). Dot size corresponds to the number of cells in which the gene is expressed. B Venn diagram showing the overlap of DEGs between Cocktail- vs. DMSO-treated cells (green) and FFA + Cocktail vs. DMSO-treated cells (red). C Scatter plot of the gene ontology (GO) analyses using gprofiler of the genes up-regulated specifically upon cocktail treatment (green); up-regulated specifically upon FFA + Cocktail treatment (magenta), and in both Cocktail and FFA + Cocktail (beige). The top 5 most significantly enriched terms are depicted. D Bar plot showing for each subgroup the percentages of genes specific to (a) Cocktail vs. DMSO (green); (b) genes up-regulated in both, Cocktail vs. DMSO and FFA + Cocktail vs. DMSO (beige); (c) specific to FFA + Cocktail vs. DMSO (magenta). E GSEA plot for FFA + Cocktail vs. Cocktail on the pathway of “Metabolism of xenobiotics by CYP450,” enriched in the Cocktail vs. DMSO-specific genes. F Heatmap depicting log₂-fold change to DMSO level of genes involved in drug-metabolism related pathways that were enriched in Cocktail vs. FFA + Cocktail treatment

Additionally, to identify global transcriptional changes between Cocktail and FFA + Cocktail, differential expression analysis was performed in pseudobulk (Methods). Using DMSO as a baseline expression, 264 genes were up-regulated uniquely under Cocktail treatment; 234 genes were commonly up-regulated in both Cocktail treatment and FFA + Cocktail; and 602 genes were up-regulated solely in FFA + Cocktail, suggesting that more biological processes are affected (Fig. 4B, Additional file 1: Fig. S6D and S8A, Additional file 6: Table S5 and Additional file 9: Table S8). To further investigate the affected pathways, gene ontology (GO) analyses were performed on these genes (Fig. 4C). Specifically, upon Cocktail treatment, we observed an evident enrichment in pathways responsible for the metabolism of xenobiotic compounds (Fig. 4C, top, green). Genes commonly up-regulated upon Cocktail and FFA + Cocktail were less specific for drug metabolism and showed an enrichment in pathways for general stimulus responses (Fig. 4C, middle, beige). Finally, the genes specifically up-regulated in FFA + Cocktail showed an enrichment in stress-related pathways (Fig. 4C, bottom, magenta). The percentage of genes in each category was similar in all four subgroups of hepatocytes, suggesting that drug metabolism was similarly affected by the accumulation of fat (Fig. 4D). In addition, a down-regulation of key hepatic marker genes was observed in FFA + Cocktail (Additional file 1: Fig. S8B). For instance, CYP2A7 was up-regulated in subgroups I, III, and IV upon Cocktail treatment, but the induction was impaired in the presence of fat (Additional file 1: Fig. S8B).

With the aim of dissecting the impact of chronic accumulation of fat on drug metabolism, we compared FFA + Cocktail and Cocktail treated cells in a gene set enrichment analysis (GSEA) (Methods) [99]. Among the top significant enriched pathways, the metabolism of xenobiotics by the CYP450 family was down-regulated (Fig. 4E), while the insulin resistance pathway was up-regulated in the presence of fat (Additional file 1: Fig. S8C), indicating a dysregulation of multiple metabolic pathways. For example, while chronic exposure to lipids did not trigger a consistent up-regulation of the insulin resistance pathway, the combination of FFA + Cocktail led to the up-regulation of NFKBIA, TRIB3, CREB3L3, IRS2, SLC2A1, SOCS3, PIK3CD, CREB5, RELA, and CPT1B in all subgroups (Additional file 1: Fig. S8D). Moreover, upon incubation with Cocktail, drug metabolism-related genes were up-regulated consistently in all functional subgroups of hepatocytes (I, II, II) (Fig. 4F). However, in the FFA + Cocktail treatment, subgroup III showed a drastic down-regulation in multiple phase I CYPs, e.g., CYP2C19, CYP2C8, CYP2C9, CYP2D6, CYP3A4, and CYP3A5; in phase II enzymes, GSTA1, GSTA2, and SULT2A1; and in phase III enzymes, SLCO1B1 and ABCG5, exhibiting an overall impairment in all three phases of drug metabolism (Additional file 1: Fig. S9).

In summary, we observed that intracellular lipid accumulation led to an impairment of drug metabolism characteristics for each subgroup, in which multiple metabolic pathways are simultaneously affected. Therefore, primary human hepatocytes lose their drug-related metabolic specificity in the presence of chronic accumulation of lipids.

Discussion

The application of single-cell genomics technologies allows the dissection of subtle differences within a seemingly homogeneous population of cells, which has been demonstrated in a plethora of organs, tissues, and cell types, including the liver [100,101,102]. The assessment of safety, toxicity, and efficacy of drugs is generally performed in bulk analyses, representing average features of the most abundant transcripts or the most abundant cell type [35, 50]. However, this approach hinders the assessment of cellular heterogeneity and the identification of cellular phenotypes that might be rare or display opposite transcriptional responses upon exposure to xenobiotics.

To circumvent this limitation, we have used single-cell transcriptomics to assess the metabolic profiles of individual PHHs, a gold standard human in vitro liver model, to study drug-related metabolic capacity in healthy condition and in response to environmental factors. Across all donors and treatment conditions, four major subgroups of hepatocytes were identified. Consistently with previous literature on PHHs [103,104,105], we identified a subgroup of PHHs losing the characteristic mature hepatocyte signature after 72 h in culture (subgroup IV, Fig. 1, Additional file 1: Fig. S2). For instance, in this subgroup, we observed a down-regulation of upstream liver-enriched transcription factors such as MLXIPL (ChREBP), RXRA, NH1H4 (FXR), PPARA, HNF4A, and CEBPA [106, 107] (Additional file 1: Fig. S7B). Further analysis of metabolic profiles showed that subgroups I, II, and III were specialized in sterol and bile acid, carbohydrate and phase II, and lipids and phase III metabolism, respectively. Thus, we defined subgroups I, II, and III as metabolically active based on their respective transcriptional profiles (Fig. 1). Interestingly, in subgroup II, a high expression of stress markers was observed, suggesting that this group is prone to eventually lose mature hepatocyte-like expression (Additional file 1: Fig. S2) [108, 109]. This observation is supported by the high correlation between subgroup II and subgroup IV in vitro (Additional file 1: Fig. S4F). Moreover, the percentage of cells in each subgroup changed upon treatment, suggesting that fat accumulation and xenobiotics might change the number of cells per hepatocyte subgroup and the functional specialization of the liver tissue.

Another source of functional specialization is the well-known hepatic zonation in the liver. The zonation patterns found in vivo due to the oxygen and nutrient gradient along the lobule axis are not conserved in 2D in vitro models [91, 110,111,112,113,114]. To deeply investigate if hepatocyte subgroups are related to reminiscent liver zonation, we have compared our findings in vitro with two publicly available in vivo data sets comprising a cohort of 14 human donors [1, 5]. We have revealed the presence of the three metabolically active subgroups in vivo, mostly independent of zonation. Therefore, subgroup specialization and hepatocyte zonation might play both a key role in liver function. Since liver zonation has a profound impact in transcriptional regulation in individual cells [10], it is technically challenging to dissect the role of hepatocyte subgroups in vivo.

In vitro, we were able to assess the drug metabolic capacity of the hepatocyte subgroups by means of a well-known phenotyping cocktail (Fig. 2). The drug metabolism in vitro can be defined by three major phases. In phase I, oxidation, hydrolysis, reduction, and cyclization reactions are catabolized mainly by the cytochrome P450 (CYP450) superfamily of monooxygenase enzymes. The main members are the isoforms CYP1A2, 2C9, 2C19, 2D6, and 3A4, which account for the metabolism of around 70–80% of the clinically available drugs [19, 20]. This strategy, formally known as the “cocktail approach,” has been used to monitor the cytochrome enzymatic activity and changes in the induction of their corresponding mRNA levels in bulk [22,23,24,25,26,27,28]. We have selected the Sanofi-Aventis cocktail to dissect changes in the transcriptome of individual human hepatocytes. The Sanofi-Aventis cocktail can be used in different species, including mouse [63, 115], primates [116, 117], dogs and minipigs [116, 118], and humans, and in both in vivo [25, 63, 117] and in vitro studies [116].

With the Sanofi-Aventis cocktail, we have revealed transcriptomic responses in individual cells that were otherwise concealed in bulk studies. For instance, a high percentage of cells in subgroup IV might lead to an underestimated effect of drug induction, which is highly relevant in the safety evaluation of xenobiotics [119]. Additionally, exploration of toxicological interactions (i.e., CTD) showed that each hepatocyte subgroup is specialized in the metabolism of certain xenobiotics, which suggests that certain subgroups of hepatocytes could be more susceptible to develop adverse drug reactions (ADRs) and toxic metabolites when challenged by a specific compound (Fig. 2E). Furthermore, we anticipate that chronic liver diseases might also affect the percentage of cells in each hepatocyte subgroup and therefore have major implications in the assessment of drug efficacy, safety, and toxicity in the early phases of drug development.

To extend the impact of cellular heterogeneity in the human liver, we have mimicked hepatic steatosis and early stages of non-alcoholic fatty liver disease (NAFLD) by loading the cells with intracellular lipids (Fig. 3) [35, 120]. Chronic accumulation of fat led to an increase in transcriptional variability in subgroups I and II, indicating transcriptional noise and random fluctuation in the mRNA level in individual cells. Importantly, in subgroup III, specialized in the metabolism of lipids, transcriptional variability was reduced upon fat accumulation, suggesting a robust and tight coordinated response to chronic accumulation of lipids. In fact, it is well-known that after two-third partial hepatectomy, there is a prominent hepatic fat accumulation in the first round of DNA synthesis before mitotic activity [121, 122]. During partial hepatectomy-induced regeneration, four continuous waves of hepatic replication have been described as associated to a tight control on the number of hepatocytes entering into cell cycle. Interestingly, three waves of hepatic fat accumulation were coupled to hepatocyte replication [123], suggesting that lipid accumulation and cell cycle in hepatocytes are co-regulated processes. However, the precise molecular mechanism underlaying changes in cell cycle and fat accumulation still remains unknown [124].

It is also possible, that the lipid accumulation leads to increased inflammation in a specific cell population. Upon fat accumulation, we found subgroup-specific signatures in which lipopolysaccharide and chemokine metabolism were upregulated in subgroup I (Fig. 3). Subgroup I is the largest in our studied primary hepatocytes. These results reinforce the notion that changes in the proportion of hepatocyte subgroups might determine the functional outputs in response to environmental or dietary factors.

For instance, hepatic fat accumulation also occurs during heathy aging [125], and an age-related increase in transcriptional variability has been reported in several tissues and cell types [9, 94, 96, 126]. Additionally, aging affects the hepatocyte function [127] and cytochrome P450 enzymes [40, 128,129,130]. Therefore, further investigations on the age-associated cellular heterogeneity and drug metabolism could anticipate the unexpected adverse drug reactions or drug-induced liver toxicity in the elderly.

Additionally, in the elderly, the co-administration of drugs known as polypharmacy is highly common for the treatment of age-related comorbidities, increasing the risk for the development of ADRs and, more specifically, DILI [38, 39]. In all subgroups of hepatocytes, intracellular lipid accumulation diminished their capacity to metabolize the five-drug cocktail (Fig. 4). Especially, subgroup III, responsible for lipid and phase III metabolism, showed the most prominent decrease in the five targeted cytochromes simultaneously (Fig. 4).

Conclusions

In summary, cellular heterogeneity is affected by intrinsic and extrinsic factors which might lead to aberrant metabolism, producing toxic metabolites and associated stress-related responses and inflammation. Intracellular lipid accumulation increases transcriptional variability and loss of expression in vitro, diminishing the capacity for drug metabolism in individual cells. Extrapolating these results to in vivo human biology, our findings would suggest that fat accumulation could further accelerate age-related processes by means of increased transcriptional noise and susceptibility to adverse drug reactions. Assessing the impact of cellular heterogeneity on tissue function will shed light on novel molecular mechanisms underlaying chronic or age-related pathologies.

Methods

An outline of the experimental design is shown in Fig. 1. Briefly, commercially plateable and interaction-qualified cryopreserved human hepatocytes (Lonza, Walkersville, MD, USA) were purchased from Lonza from four different donors: HUM180812 (male, 57 years old, Hispanic) and HUM4152 (male, 18 years old, Caucasian), HUM181641 (male, 56 years old, Caucasian) and HUM4190 (male, 26 years old, Caucasian) (Additional file 2: Table S1). All donors had a body mass index in the normal range and were not diabetic (Additional file 2: Table S1). Cell lines were not authenticated and were not checked for mycoplasma contamination. The first two donors were processed simultaneously in a first batch, and the second two in a second batch, aiming for a maximum of eight samples processed at a time.

Cell culture

Each cryovial of PHH was thawed and plated according to Lonza’s “Suspension and plateable cryopreserved hepatocytes: technical information and instructions.” the protocol was followed stepwise minutely, using the recommended thawing and plating media (Lonza, MCHT50, and MP250, respectively). The cells were dispensed and mixed using only wide orifice tips (Rainin, Ref. 17014297). For efficient cell seeding densities and attachment, cells were counted using the Trypan Blue Exclusion Method and seeded into Collagen-I coated plates at a density of approximately 100,000 cells/cm² following Lonza’s instructions (Lonza, “Suspension and Plateable Cryopreserved Hepatocytes Technical Information and Instructions”). Six hours post-seeding, cells were washed with 1 mL of pre-warmed Maintenance Medium (Lonza, MCHT50) before the addition of treatment media. The treatment medium was renewed every 24 h for a total incubation period of 72 h post-seeding. Free fatty acids (FFA) consisting of a 200 µM mixture of a 2:1 ratio of the unsaturated oleic acid to the saturated palmitic acid were added to the maintenance media, to mimic the levels in human steatosis [35]. In order to facilitate FFA uptake, pre-bounding of free fatty acids to 1% bovine serum albumin in a 1:5 molar ratio (Sigma-Aldrich) was performed by heating the mixture at 40 °C for 2 h [36].

Drug cocktail preparation and storage

The individual components of the 5-drug cocktail [24] were dissolved in sterile DMSO, filtered through a 0.2-µM syringe filter (Merck, SLGVV255F), and stored at − 80 °C for a maximum of six months (“compound stock concentration”, Additional file 11: Table S10). The individual drugs were mixed at 200 × concentration (“working concentration”, Additional file 11: Table S10) and added to the cells to a final concentration of 80 µM Caffeine (Sigma-Aldrich, Ref. 56396-100MG), 5 µM Midazolam (LGC Chemicals, Ref. LGCFOR1106.00), 17 µM Omeprazole (TRC Chemicals, Ref. 0635000), 20 µM S-Warfarin (Sigma-Aldrich, Ref. UC214-5MG) and 23 µM Metoprolol (TRC Chemicals, Ref. M338815). The final DMSO concentration used on the cells was 0.5% (v/v), in all conditions.

Single-cell RNA-seq sample preparation and sequencing

After a total of 72-h incubation in a treatment culture medium, cells were detached with pre-warmed 0.25% Trypsin (Gibco, 25200056) for 7 min, followed by trypsin inactivation. The dissociated cells were then collected to pellet by centrifugation at 100 × g for 5 min at room temperature (RT). Cells were washed twice with 1 mL of pre-warmed 1 × PBS pH 7.4 (Gibco, 10010023), followed by cells pelleting at 100 × g for 5 min at RT. A live-cell selection was performed using a Dead Cell Removal Kit (Miltenyi Biotec, 130–090-101) as follows: cells were pelleted at 100 × g for 5 min at RT, resuspended in 100 µL of dead-cell removal microbeads and incubated for 15 min at RT. Dead cells were positively selected by passing the cell suspension through a MS column and performing a wash with a total of 2 mL of binding buffer. Living cells were eluted from the column and collected in 2-mL Eppendorf tubes. After pelleting by centrifugation at 100 × g for 5 min at RT, cells were resuspended in 1xPBS pH 7.4 supplemented with 0.04% BSA (Miltenyi Biotech, 130–091-376), stained with trypan blue to assess viability, and counted in a hemocytometer. A single-cell suspension was obtained by dissociating cells with wide orifice pipette tips and preparing the target cell stock concentration for loading the 10 × chip.

Single-cell RNA-seq libraries were prepared from each sample following the 10x Genomics Single Cell 3′ Reagent Kit User Guide (v3 or v3.1, manual CG000183 and CG000204, respectively) in the Chromium Controller (10x Genomics). The quality control of cDNA and obtained final libraries was performed using a Bioanalyzer High Sensitivity DNA Analysis assay (Agilent). Library quantification was performed using the Collibri™ Library Quantification Kit (Thermo Fischer Scientific, A38524500) in a QuantStudio™ 6 Flex Real-Time PCR System (Thermo Fisher Scientific). Both batches were sequenced in a NovaSeq6000 sequencer (Illumina) at the HMGU Core Facility for NGS Sequencing. The first batch (HUM180812 and HUM4152) was sequenced in a S2 flowcell at a depth of 250,000 reads per cell. The second batch (HUM181641 and HUM4190) was sequenced in a SP flowcell at a sequencing depth of 20,000 reads per cell. The sequencing length was set as indicated by 10x Genomics: 28/8/0/91.

Single-cell ATAC-seq sample preparation and nuclei isolation

After a total of 72-h incubation in treatment culture medium (DMSO), cells were trypsinated with pre-warmed 0.25% Trypsin (Gibco, 25200056) for 7 min, followed by inactivation. Approximately 500,000 cells were pelleted at 100 × g for 5 min at RT. Cells were washed twice with 1 mL of cold 1 × PBS pH 7.4 (Gibco, 10010023), centrifuged down at 500 rcf for 5 min at 4 °C and resuspended in 1 mL of homogenization buffer (5 mM MgCl₂, 25 mM KCl, 10 mM Tris buffer pH 8.0, 250 mM sucrose, 1% DTT, protease-inhibitor (Life Technologies, 1187358001), 0.2% NP-40 and 0.3% Triton-X). Cells were transferred to a 2-mL homogenizer and 5 strokes were performed with the lose pestle (A) followed by 10 min incubation on ice and 25 strokes with the tight pestle (B). Nuclei were then transferred to a 2-mL Eppendorf tube. The douncer was rinsed using additional 400 µL of HB and added to the tube, which was then centrifuged at 500 rcf for 5 min at 4 °C. After discarding the supernatant, the pellet was resuspended in 1 mL of swelling buffer (10 mM Tris–HCl pH7.5, 2 mM MgCl₂, 3 mM CaCl₂) and 1 mL of prechilled swelling buffer with 10% glycerol was added drop-wise and then pipette-mixed 5 times [131]. The mixture was incubated for 10 min on ice with occasional flicking. Then, the nuclei were pelleted and resuspended in 1 mL of 10 × Genomics Wash Buffer (10 mM Tris–HCl pH 7.4, 10 mM NaCl, MgCl₂, 1% BSA, 0.1% Tween-20 and nuclease-free water), counted nuclei per µL and incubated for 15 min on ice. Proceeding to nuclei lysis, nuclei were centrifuged at 500 rcf at 4 °C for 5 min. Then, the supernatant discarded and nuclei resuspended in 200 µL of 10x Genomics Lysis Buffer with 0.2% NP-40 instead of the original 0.1% concentration (10 mM Tris–HCl pH 7.4, 10 mM NaCl, MgCl₂, 1% BSA, 0.1%Tween-20, 0.01% Digitonin, as a minor 0.2% NP-40 and nuclease-free water). Nuclei were lyzed for 10 min on ice and then 1 mL of wash buffer was added followed by rocking the tube gently. Nuclei were pelleted at 500 rcf for 5 min at 4 °C, then 500 µL of PBS added without resuspending and nuclei spun down at 500 rcf for 5 min at 4 °C.

Tagmentation, library preparation, and sequencing

Tagmentation reaction was prepared following the 10 × Genomics demonstrated protocol [131]. In brief, isolated nuclei were resuspended in 500 µL of 1X Nuclei Storage Buffer, pelleted at 500 rcf for 5 min at 4 °C, and discarding the supernatant. The transposition master mix was prepared by mixing 7 µL of ATAC Buffer and 3 µL of ATAC enzyme [131]. Then, using the Nuclei Concentration Guidelines calculator, 5 µL of the mixture of nuclei (2.9 µL) and diluted Nuclei Storage Buffer (2.1 µL) were incubated for 60 min at 37 °C in a thermocycler.

Single-cell ATAC-seq library was prepared following the 10x Genomics Chromium Next GEM Reagents Kit and Reagent Kit User Guide (v1.1, manual CG000209, RevD) [131] in the Chromium Controller (10x Genomics). The cDNA and library’s quality were assessed in a Bioanalyzer 2100 using the High Sensitivity DNA Analysis assay (Agilent). Final libraries’ quantification was performed using the Collibri™ Library Quantification Kit (Thermo Fischer Scientific, A38524500) in a QuantStudio™ 6 Flex Real-Time PCR System (Thermo Fisher Scientific). Sequencing was performed in a SP 100 flowcell in a NovaSeq 6000 sequencer (Illumina) at the HMGU Core Facility for NGS Sequencing with the sequencing length recommended by 10x Genomics (50, 8, 16, 50). Approximately, 50,000 reads per nucleus were yielded.

Read alignment, counting, and filtering of the combined batches

Reads were aligned to GRCh38 and counted using 10x Genomics Cellranger 4.0.0 with standard parameters for each batch individually (Additional file 3: Table S2). The resulting count matrices were combined into a count matrix of 63,527 cells times 19,971 genes. Cells with at least 1000 counts and 500 genes were kept. Genes were kept if they were present in at least 5 cells and had fewer than 5 million reads. Scrublet [132] was used to identify doublets. Due to hepatocytes being subject to polyploidyzation, a lenient cutoff of 0.15 was used to avoid unwanted removal of tetraploid hepatocytes. This led to 1.7% of cells being annotated as doublets and subsequently being removed. Lastly, cells with more than 1% mitochondrial reads were removed, resulting in a filtered matrix of 49,378 cells times 16,256 genes.

Normalization and initial clustering

Scran was used for library size normalization with parameter min.mean = 0.05. After normalization, cells with more than 20,000 normalized counts were removed. Scanpy functions were used for calculating principal components (PCs) and clustering. Shortly, the top 50 PCs were used to construct a neighboring graph before calculating a UMAP embedding and Louvain clusters. For two samples, encapsulation of the cells in droplets partially failed during the 10X library preparation (wetting failure). Therefore, an initial Louvain resolution of 0.5 was used to verify that cells stemming from these samples clustered separately. Indeed, some cells from one of the failed samples clustered apart from the rest of the cells, and another two clusters showed fractal structures that were not present in the samples where encapsulation worked correctly. Hence, based on the knowledge of the wetting failure, these three clusters were removed, resulting in a final matrix of 38,232 cells times 16,256 genes. To remove batch effects in downstream analysis, the two batches were integrated using Harmony [48].

Subgroup annotation based on individual donor analysis and clustering

In each of the four donors, the resolution of the Louvain clustering was selected so that every cluster contained every treatment condition (Additional file 1: Fig. S2A). Then, the top 1000 genes per Louvain cluster were calculated using sc.tl.rank_genes_groups with n_genes = 1000. To find out which clusters were shared between donors, the percentage of overlaps between the top 1000 genes per cluster was calculated and hierarchical clustering was performed. By this approach, three distinct groups of similar Louvain clusters were detected between donors (Additional file 1: Fig. S1C). Two of the donor-specific clusters did not group together with clusters from other donors. Therefore, they were labeled as individual clusters to inspect where they would fall on the integrated data set (Additional file 1: Fig. S1C). Due to differences in sequencing depth, one group was identified in the first batch which was assigned to the cells losing expression. After integration and putting the group labels from the individual donor analysis on top of the combined UMAP, Louvain clustering was performed, showing that cells from the second batch that we had assigned to subgroup II were clustering with the cells losing their expression from the first batch. Therefore, those cells we re-labeled as losing their expression, subgroup IV (Additional file 1: Fig. S2A). One of the two donor-specific clusters that could initially not be assigned to a shared cluster was split into two clusters on the combined data set, which were assigned to subgroup I and subgroup II, respectively. The other individual cluster could be associated to the cells losing their expression. Marker genes were taken from reported literature to assign metabolic preferences to the functional subgroups I, II, and III (Additional file 3: Table S2).

Comparison to publicly available in vivo data

Data were obtained from GEO (Accession number GSE124395) [1]. Genes with zero expression across all cells, and cells that did not express any genes, were removed. During further filtering, cells with 100 to 6000 genes and 800 to 30,000 reads were kept and genes were kept if they were covered in at least 10 cells, resulting in a count matrix of 11,059 cells times 19,416 genes. To keep the processing comparable to our data, normalization was performed using scran with parameter min.mean = 0.05. After normalization, cells with more than 20,000 normalized counts were further removed, resulting in a count matrix of 11,043 cells times 19,416 genes. Clustering was performed using scanpy functions as described above. Louvain clustering was performed at a resolution of 0.08 to computationally separate the cell types in this data set from each other. Expression of mature hepatocyte markers, such as ALB, HNF4A, and TTR was used to identify the hepatocyte cluster. Louvain clustering with a resolution of 0.2 was then performed on the hepatocytes to identify hepatocyte subgroups in vivo. The expression of marker genes used for subgroup identification in our data was investigated on the Louvain clusters. Based on this marker gene expression, Louvain clusters were assigned to subgroups I, II, and III (Fig. 1D, Additional file 1: Fig. S5A).

As zonation impacts gene expression along the pericentral-periportal axis in vivo, we investigated its connection to our subgroups. Zonation markers were taken from Aizarani et al. (2019) [1]. Based on the 35 zones reported in their study, the genes were grouped into three zones, pericentral, mid, and periportal through binning. Then sc.tl.score_genes was used to calculate scores for each of the three zones (Additional file 1: Fig. S5C). Cells were assigned to pericentral [133] if they had a central vein (CV) score > 0.45 and a periportal (PV) score < 1. If cells had a PV score >  = 0.65 and a CV score <  = 0.45, they were assigned to periportal. The rest of the cells were assigned to mid-zone (Additional file 1: Fig. S5B and C). For each of the three subgroups, the percentage of cells in CV, mid, and PV was calculated subsequently (Additional file 1: Fig. S5B). For visualization purposes, the subgroups were separated and zonation marker genes for each subgroup were depicted on a UMAP for the in vivo data set and for our in vitro data (Additional file 1: Fig. S5D).

To further investigate the presence of our four hepatocyte subgroups in vivo, an additional human data set was downloaded from GEO (Accession number GSE115469) [5] filtered and normalized as before. In order to reduce potential noise, genes were removed if they were not present in at least three cells and counts were log-transformed for better comparability between datasets. Hepatocytes were isolated based on the authors’ annotation available at GSE115469. Scanpy functions were used as described above to perform clustering with a Louvain resolution of 0.5 to separate subgroups of hepatocytes, leading to six Louvain clusters highly overlapping with the clusters reported in the initial study. Marker gene expression of the three identified subgroups in our data was again used to assign the six Louvain clusters to the three subgroups. To increase power, the two in vivo data sets were integrated using scGen [134] (Additional file 1: Fig. S4E). The top 10 DEGs per our in vitro subgroups were calculated and their scaled mean expression was correlated to the combined in vivo scaled mean expression. This showed a high correlation of marker genes between identified subgroups in vitro and in vivo (Additional file 1: Fig. S4F).

Analysis of upstream regulators using ChEA3 and single-cell ATAC-seq

After defining metabolically functional subgroups in our in vitro data, we first calculated the top 500 DEGs per subgroup under DMSO in comparison to each other. To assess what drives the basal functional specialization of the subgroups, these top 500 genes per subgroup were taken as input for the online tool ChEA3 [53] which predict transcription factors regulating the gene expression per subgroup. For the purpose of visualization, five transcription factors among the top 25 per subgroup were depicted as a stacked violin plot (Additional file 1: Fig. S2E).

Raw reads were aligned to GRCh38 using cellranger-atac 2.1.0. Transcription factor binding site (TFBS) annotation was downloaded from ReMap2022 (https://remap2022.univ-amu.fr/) and filtered to only retain the binding sites of the top 25 transcription factors per subgroup that were predicted by ChEA3 (Additional file 4: Table S3). Using the fragment file obtained as output from cellranger-arc, we constructed the count matrix for these binding sites using epiScanpy’s function epi.ct.peaks_mtx(). This resulted in a matrix of 360,878 barcodes × 4,029,591 TFBS. Cells were removed if they had fewer than 2000 TFBS covered, a nucleosome signal score of > 5, and a TSS-enrichment score of < 2. TFBS were removed if they were present in fewer than 30 cells. Furthermore, highly variable TFBS were selected by running epi.pp.highly_variable() with parameter min_score = 0.515. To only consider binding sites at close proximity of the subgroup-specific genes, TFBS were annotated by the genes in their proximity using epi.tl.find_genes() with parameters upstream = 1000 and downstream = 100. We used this annotation to select only the TFBS falling within 1000 bp upstream and 100 bp downstream of the TSS of any of the top 500 DEGs per subgroup. This resulted in a final matrix of 3582 cells × 16,488 TFBS.

To confirm the presence of characteristic chromatin landscapes for our identified subgroups, we investigated the co-accessibility of the TFBS in the promoter regions of 5 genes per subgroup (subgroup I: IQGAP2, LARP1, CTSB, DIAPH1, SPTBN1; subgroup II: RAC1, UGP2, RPL5, RPS12, ZNF706; subgroup III: OGT, ALB, CCNL1, NRBP2, PPP1R12B). Briefly, the accessibility of each binding site across all cells was correlated to all others. Then, per TFBS, we calculated its average correlation to all the binding sites in the promoter of every selected gene. Clustering the TFBS based on these average correlations showed that TFBS at the promoters of subgroup-specific genes were clustering together. This indicates their co-accessibility in the cells (Additional file 1: Fig. S5A).

Additionally, Louvain clustering was applied using the above count matrix, to group cells showing similar openness profiles at the TFBS regions in front of subgroup-relevant genes. Cells that could not be assigned to one of the subgroups were disregarded. Finally, for the 1176 annotated cells, we visualized chromatin openness per subgroup at the binding sites of the ChEA3-predicted transcription factors in proximity to subgroup-defining genes (Additional file 1: Fig. S5B-D) [http://andrewjohnhill.com/blog/2019/04/12/streamlining-scatac-seq-visualization-and-analysis/)].

Differential expression analysis between conditions in the subgroups

After defining metabolically functional subgroups in our in vitro data, we first calculated the top 500 DEGs per subgroup under DMSO in comparison to each other. To assess what drives the basal functional specialization of the subgroups, these top 500 genes per subgroup were taken as input for the online tool ChEA3 [53] which predict transcription factors regulating the gene expression per subgroup. For the purpose of visualization, five transcription factors among the top 25 per subgroup were depicted as a stacked violin plot (Additional file 1: Fig. S2E). To explore the implications of the functional specialization on the metabolic capacity of hepatocytes, we then performed differential expression analysis between Cocktail- and DMSO-treated cells in each of the functional subgroups. Genes were defined significantly up-regulated if they had a log₂-fold change of greater than 1 and a Bonferroni-adjusted p-value below 0.05. We used the venn package to visualize the overlap of significantly up-regulated genes between the subgroups. Furthermore, ShinyGO [135] was used to investigate the enrichment of the subgroup-specific up-regulated genes in pathways known to be relevant for drug metabolism (Additional file 5: Table S4).

When comparing FFA- with DMSO-treated cells, we observed that there were on average 3.5 times fewer genes with a positive log₂-fold change than when we compared Cocktail and DMSO. Therefore, to capture the subtler effects of fat accumulation on the cells, a gene was identified as significantly up-regulated if it had a Bonferroni-adjusted p-value below 0.05 and a log₂-fold change greater than 0.75. We used the same marker genes as described above (Additional file 3: Table S2) to count overlaps between groups of marker genes and significantly up-regulated genes in the subgroups. To investigate in which pathways related to biological processed the genes up-regulated upon FFA-treatment were enriched in every subgroup, gprofiler was used (https://pypi.org/project/gprofiler-official/).

Assessing transcriptional variability through the coefficient of variation

It is generally assumed that lowly expressed genes have an inflated transcriptional variability. Therefore, genes with a mean normalized log-transformed expression smaller than 0.25 were removed, and the coefficient of variation per condition in each of the subgroups was calculated on 3434 genes. To obtain the coefficient of variation on normalized, log-transformed counts, we used the formula described in Canchola et al. [136].

$$\mathrm{CV}=\sqrt{{e}^{{\sigma }^{2}}}-1,$$

where σ² is the variation of gene j in the group of interest.

In every subgroup, a Mann–Whitney U test was performed to check if the coefficient of variation differed significantly between treatment condition.