Single-cell RNA sequencing of human liver reveals hepatic stellate cell heterogeneity.

BACKGROUND & AIMS: The multiple vital functions of the human liver are performed by highly specialised parenchymal and non-parenchymal cells organised in complex collaborative sinusoidal units. Although crucial for homeostasis, the cellular make-up of the human liver remains to be fully elucidated. Here, single-cell RNA-sequencing was used to unravel the heterogeneity of human liver cells, in particular of hepatocytes (HEPs) and hepatic stellate cells (HSCs).
METHOD: The transcriptome of ~25,000 freshly isolated human liver cells was profiled using droplet-based RNA-sequencing. Recently published data sets and RNA in situ hybridisation were integrated to validate and locate newly identified cell populations.
RESULTS: In total, 22 cell populations were annotated that reflected the heterogeneity of human parenchymal and non-parenchymal liver cells. More than 20,000 HEPs were ordered along the portocentral axis to confirm known, and reveal previously undescribed, zonated liver functions. The existence of 2 subpopulations of human HSCs with unique gene expression signatures and distinct intralobular localisation was revealed (i.e. portal and central vein-concentrated GPC3 + HSCs and perisinusoidally located DBH + HSCs). In particular, these data suggest that, although both subpopulations collaborate in the production and organisation of extracellular matrix, GPC3 + HSCs specifically express genes involved in the metabolism of glycosaminoglycans, whereas DBH + HSCs display a gene signature that is reminiscent of antigen-presenting cells.
CONCLUSIONS: This study highlights metabolic zonation as a key determinant of HEP transcriptomic heterogeneity and, for the first time, outlines the existence of heterogeneous HSC subpopulations in the human liver. These findings call for further research on the functional implications of liver cell heterogeneity in health and disease. LAY
SUMMARY: This study resolves the cellular landscape of the human liver in an unbiased manner and at high resolution to provide new insights into human liver cell biology. The results highlight the physiological heterogeneity of human hepatic stellate cells.

Introduction

Establishing a more comprehensive blueprint of the human liver cellular landscape is crucial to improve understanding of the physiology of this organ and to inform the interpretation of disease-associated changes. Important limiting factors have always been the lack of adequate methods and reliable cell markers, combined with the difficulties involved in sourcing human liver tissue for research. As a result, many low-abundance cell types of the human liver remain poorly characterised. One such cell type is the hepatic stellate cell (HSC). Although extensively studied in the context of liver fibrosis, little is known about the heterogeneity of human HSCs under physiological conditions. This is the result, at least in part, of the gold-standard methods of isolation relying on uniform physicochemical properties, such as low cell density and high autofluorescence. These properties, conferred by the intracytoplasmic storage of retinyl-ester-loaded lipid droplets, are considered to be functional hallmarks of quiescent HSCs.[1], [2], [3] However, in healthy mice, HSCs have been shown to be at least heterogeneous in their capacity for retinoid and lipid storage.

Recent advances in single-cell analyses have enabled whole-organ heterogeneous cell populations to be studied without the need to separate or expand the cells in culture, opening avenues to build comprehensive inventories of the different functional subpopulations that comprise organs.[6], [7], [8] In mice, single-cell RNA-sequencing (scRNA-seq) showed that HSCs are spatially and functionally zonated during centrilobular injury-induced fibrogenesis. In human liver, the technique was applied to reveal distinct hepatic macrophage (MP) subpopulations and bipotent epithelial progenitors,, and to study the fibrotic niche., Although some of these studies reported the first genome-wide, single-cell resolved human HSC transcriptomes, they did not enable a comprehensive characterisation of this important cell population.

To deepen the characterisation of the heterogeneity of human liver cells, in particular of HSCs, in the current study, a streamlined tissue dissociation protocol was applied to freshly explanted human livers to recover sufficient numbers of high-integrity, low-abundance liver cell types, void of any preselection bias. The duration of post-dissociation separation steps was reduced to limit transcriptomic alterations in response to the loss of cell–cell and cell–scaffold interactions; subsequently, the transcriptome of ~25,000 individual liver cells was profiled using droplet-based RNA-seq.

In total, 22 distinct cell populations were annotated that pertain to the liver epithelial, endothelial, mural, and immune cell compartments. Despite the implementation of centrifugation steps to balance the representation of parenchymal and non-parenchymal cells, hepatocytes (HEPs) represented ~80% of sequenced cells in the data set. Using previously described landmark genes of zonated HEPs, the cells were classified along the portocentral axis to comprehensively map zonated liver functions. This analysis, for the first time, identified 2 subpopulations of human HSCs, characterised by different spatial distribution in the native tissue and gene expression profiles suggesting functional disparity. Collectively, this study sheds new light on the cellular compartments that underlie the physiology of the human liver and represents a reliable reference for studies on human liver cell biology.

Material and methods

Sourcing of fresh human liver tissue

The right liver lobe of donor 156 (obtained after the liver was processed for reduced-size liver transplantation of the left lobe) and the whole liver of donor 158 (Crigler–Najjar syndrome; not suitable for transplantation) as well as paraffin-embedded liver fragments (Table S1) were provided by the Hepatic Biobank of the Cliniques Universitaires Saint-Luc (CUSL). All protocols and experiments were approved by the ethics committees of the CUSL and Université catholique de Louvain (Agreements 2015/03NOV/585 and 2019/04MAR/100). In accordance with Belgian law, the tissue donors and/or their families received the necessary information and provided active or passive consent for the use of human residual material for research purposes.

Isolation of human primary liver cells

The freshly explanted human livers were kept on ice in IGL-1 solution (Institut Georges Lopez, Lissieu, France). Total liver cell suspensions were obtained by a 2-step perfusion method optimised for cell therapy applications, as previously described,, between 6 h and 12 h after surgery. In brief, the liver tissue was sequentially perfused with 37°C prewarmed Earle’s balanced salt solution containing 440 μM EGTA and 1400 UI/L collagenase (Sigma-Aldrich, St. Louis, MO, USA). Given that HEPs represent most liver cells, and to guarantee cells from the non-parenchymal fraction would be represented in the data set, part of the resulting whole-liver cell suspension was further separated into parenchymal and non-parenchymal fractions by 2 low-speed centrifugation steps (160 g, 10 min each at 4°C). Samples from the whole fraction, the parenchymal fraction and the non-parenchymal fraction were transferred to DNA LoBind vials (Eppendorf, Hamburg, Germany). A total of 4 samples per donor were prepared for loading in the Chromium single-cell platform (10X Genomics, Pleasanton, CA, USA), according to the manufacturer's instructions. The cells were washed with 0.04% bovine serum albumin (BSA)-supplemented PBS and subsequently filtered on Flowmi cell strainers (Sigma-Aldrich). The cell concentration was adjusted to 106 total cells/ml and cell viability was evaluated by a Trypan blue exclusion assay (Table S2).

Library preparation and sequencing

Samples were further processed for droplet-based RNA-seq using the Chromium single cell 3' reagent kits version 2, chip and controller (10X Genomics), following the manufacturer's instructions. Post-cDNA amplification reaction and post-library construction quality controls were performed with an Agilent Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The Kapa DNA quantification kit (Sigma-Aldrich) was used for post-library construction quantification. Libraries were sequenced at an average depth of 65,000 reads/cell using a NextSeq 500 system (Illumina Inc., San Diego, CA, USA) with the following parameters: read 1, 26 cycles; i7 index, 8 cycles; read 2, 58 cycles.

Quality control and treatment of new data sets

Raw sequencing data files (.bcl files) were converted to FASTQ format using the mkfastq pipeline from CellRanger version 2.1.1 (10X Genomics). Reads were aligned to the human reference genome GRCh38 release 97 from Ensembl (www.ensembl.org) using the CellRanger count pipeline. The default filtering parameters of CellRanger were applied to obtain a gene expression matrix (unique molecular identifier [UMI] counts per gene per cell) for each library (CellRanger reports are provided in the supplementary material). Two additional exclusion filters were applied: (i) cells with a library size ≤2,000 UMIs; and (ii) cells with a mitochondrial transcript detection ratio ≥25%. Genes detected in fewer than 3 cells were labelled as undetected. Data sets of samples were combined per donor, and each merged data set was first analysed separately. After annotation of the cell clusters, the data sets of both donors were pooled to form a combined data set using the Seurat method for integration. Cell clusters presenting aberrant co-expression of distinct canonical cell-type markers were labelled as potential doublets and removed. All raw sequencing data were deposited in the Gene Expression Omnibus (GEO) under accession GSE158723.

Statistical analysis

Analyses of scRNA-seq data were performed using Seurat R package version 3.1.0 (http://satijalab.org/seurat/)., Differential expression analyses were performed using the Seurat FindAllMarkers function (with Wilcoxon rank sum test, min.pct = 0.25, logfc.threshold = 0.25 and return.thresh = 0.001). Cell cycle scoring was performed using the Seurat CellCycleScoring function and the list of cell cycle genes from Kowalczyk et al. The expression of different combinations of genes was used to define scores and signatures using the Seurat PercentageFeatureSet function (Table S3).

Gene ontology analysis

The lists of significantly upregulated genes (as defined by the Seurat FindAllMarkers function, Tables S4–S7) and total detected genes within a given liver cell population were introduced in GOrilla software (http://cbl-gorilla.cs.technion.ac.il/), according to the developer's instructions., All gene ontology (GO) terms (processes, functions, and components) and a p value threshold of 0.001 were considered for the analysis. Shown significant enrichments (FDR q-value <0.1, full list in Tables S8 and S9 as an output of GOrilla) were selected according to their biological relevance. DEGs associated with a given GO term were further used to define a score using the Seurat FindAllMarkers function (Table S3).

Import of previously published data sets

The raw scRNA-seq data from 5 liver donors recently published and deposited in the GEO by Ramachandran et al. (GSE136103) were downloaded. The following exclusion filters were applied: (i) cells with a library size <300 genes; and (ii) cells with a mitochondrial transcript detection ratio >30%.

Histology

Immunofluorescence was performed on 5 μm-thick sections of formaldehyde-fixed paraffin embedded (FFPE) liver fragments. Endogenous peroxidases were inhibited for 20 min with 3% hydrogen peroxide in methanol. Antigen was retrieved in 10 mM citrate buffer pH 5.7 and aspecific binding sites were blocked with Tris buffered saline (TBS) containing 5% BSA and 0.1% Tween20. An anti-human leukocyte antigen (HLA)-DR primary antibody (clone TAL.1B5, Agilent Technologies) was incubated for 1 h at room temperature at a dilution of 1/100 in TBS containing 1% BSA and 0.1% Tween20 and detected by anti-mouse horseradish peroxidase (HRP)-conjugated polymer secondary antibodies (Agilent Technologies) for 1 h at room temperature. HRP was then visualised by tyramide signal amplification (TSA) using AlexaFluor488-conjugated tyramides (Thermo Fisher Scientific, Waltham, MA, USA). After a new citrate buffer incubation step, the same protocol was applied with anti-vimentin (VIM) primary antibody (clone V9, Agilent Technologies, dilution 1/70) or anti-CD45 primary antibody (clone 2B11+PD7/26, Agilent Technologies, dilution 1/200) and revealed with AlexaFluor647-conjugated tyramide (Thermo Fisher Scientific). After a washing step in PBS, nuclei were stained with Hoechst 33342 (Thermo Fisher Scientific) diluted in TBS containing 10% BSA and 0.1% Tween20, washed in TBS containing 0.1% Tween20, and mounted with Dako fluorescence mounting medium (Agilent Technologies). Slides were digitalised using a Pannoramic 250 FlashIII scanner (3DHistech) at 20× magnification. Other immunohistochemistry images were obtained from the Human Protein Atlas database (www.proteinatlas.org). The URLs of used images are listed in Table S10. Haematoxylin and eosin, Sirius Red (collagen), Miller (elastin), and Alcian blue (glycosaminoglycan; GAG) stainings were performed on 4% FFPE liver tissue sections of 5-μm thickness.

RNAscope 2.5 HD duplex and multiplex fluorescent v2 assays (ACDBio, CA, USA) were manually performed according to the manufacturer's instructions on 5 μm-thick sections of FFPE liver fragments. Target retrieval was performed for 15 min, protease was applied for 30 min, and 1× probes were incubated for 2 h at 40°C. The probes used were: Hs-HHIP 464811-C1, Hs-DBH 545791-C1, Hs-GPC3 418091-C2, and Hs-NGFR 406331-C3 (ACDBio). In the fluorescent assay, C1, C2, and C3 probes were revealed with Opal 520, Opal 570, and Opal 690 fluorophores, respectively (Akoya Biosciences, Menlo Park, CA, USA), diluted 1:750. Images were acquired with a SCN400 slide scanner (Leica, Wetzlar, Germany) or with a Pannoramic 250 FlashIII scanner at 40× magnification.

Results

Dissociation and annotation of human liver cell populations

Two freshly explanted human livers (156 and 158) were digested by a 2-step collagenase perfusion procedure to collect primary liver cell suspensions for droplet-based scRNA-seq analysis (Figure 1A). The livers were unsuitable for transplantation but showed no significant steatosis, pathological immune infiltrates or fibrosis (Figure S1A). A total of 25,325 individual liver cells were sequenced at an average depth of 65,000 reads per cell.

Fig. 1

Clustering and annotation of healthy human liver cell transcriptomes.

(A) Whole-cell suspensions extracted from 2 freshly explanted human livers were centrifuged and prepared for droplet-based scRNA-seq using the 10X Chromium technology. (B) UMAP visualisation of 25,325 human liver cells clustered into 22 cell subpopulations. (C) Cell-type annotation of each subpopulation based on the differential expression of liver cell type-specific genes. (D) Heatmap displaying the expression level of established liver cell type-specific markers in each cell subpopulation. (E) Top 5 genes most highly expressed in each cell type. scRNA-seq, single cell RNA sequencing; UMAP, uniform manifold approximation and projection.

The individual liver cell transcriptomes were projected on a uniform manifold approximation and projection (UMAP) plot and organised into 22 distinct cell clusters (Figure 1B). A comparative analysis of the expression of 30 specific markers found each cell cluster to pertain to 1 of the 6 main cell types of the liver (i.e. HEPs, cholangiocytes [CCs], endothelial cells [ECs], HSCs, MPs, and lymphoid cells [LPs]) (Figure 1C,D). HEPs represented ~80% of the total analysed cells, whereas CCs, ECs, HSCs, MPs, and LPs represented 3.8%, 2.5%, 1.0%, 9.8%, and 2.6% of the total analysed cells, respectively (Table S11). Importantly, neither the size of the library nor the cell cycle phase (Figure S1B,C) drove the clustering of the 22 subpopulations, with the notable exception of proliferating MPs and LPs. Despite the age and gender differences between both donors (Table S1), there was good overlap in the expression profiles and cells originating from both donors were identified in each of the 6 cell types (Figure S1D and Table S11). The identity of each cell type was further confirmed by an analysis of the top DEGs (Figure 1E), most of which encoded proteins involved in well-documented functions of the different liver cell populations, such as secreted acute phase proteins (ORM1, ORM2, HP, and SAAI) in HEPs, keratins (KRT7 and KRT19) in CCs, FCN2 and FCN3 lectins in ECs, DCN proteoglycan in HSCs, HLA class II molecules (HLA-DPA1 and HLA-DRA) in MPs, and immunoglobulin chains (IGKC, IGHA2, and IGHM) in LPs.

Comparative analysis of the CC subpopulations identified CC2 as potential progenitor cells,, characterised by the co-expression of both canonical CC and HEP markers (Figure S2A–C). However, in the absence of functional studies, one cannot exclude that this subpopulation represents artefactual libraries caused by the capture of more than 1 cell per droplet or by abundant ambient HEP mRNA. Cells of the EC, MP, and LP compartment comprised different subpopulations recently described by scRNA-seq, including macrovascular ECs, as well as periportal (PP) and pericentral (PC) liver sinusoidal ECs (Figure S2D–F), MARCOhighLILRB5high immunoregulatory MPs, killer lectin-like receptor (KLR)-expressing natural killer cells, CD19high B LPs, and Ig-producing plasma cells (Figure S2G–J).,,, Regulatory MPs and, to a lesser extent, intermediate MPs showed a good match with the gene signature of the liver-resident Kupffer cells (i.e. ARL4C, CD163, CD5L, MERTK, NR1H3, SIGLEC1, TIMD4, and VCAM1) (Figure S2K).

Liver zonation drives the clustering of human hepatocytes

To allow for a more comprehensive understanding of functional HEP zonation across the human liver lobule, the differences in gene expression were mapped between PP, midzonal (MZ), and PC HEPs. To do so, HEPs were ordered according to the expression level of genes associated with well-documented zonated functions and previously described as landmark genes for PC and PP HEPs in the murine liver (i.e. GLUL, CYP2E1, ASS1, ASL, and ALB).,, The distribution of the corresponding proteins in human liver tissue confirmed the validity of these markers for human PC and PP HEPs (Figure 2A). For robustness of the analysis, only HEPs expressing at least 4 of the aforementioned markers (12%) were selected for an intermediate classification. This first pseudo-ordering demonstrated that the correlation between gene and protein expression for these markers was precise (Figure 2B) and that metabolic zonation clearly determined the polarisation of the HEP cluster (Figure 2C). Based on this first classification, a broader PC- and PP-specific gene signature was determined (Figure 2D and Table S4) and applied for the ordering along the portocentral axis of all 20,352 HEPs in the data set (Figure 2E,F). In addition to confirming well-documented zonated functions, such as xenobiotic metabolism (PC),, fatty acid biosynthesis (PC), and the secretion of plasma proteins (PP), GO analysis of DEGs between PC and PP HEPs (Tables S5 and S8) also revealed gene networks associated with less-characterised liver functions. These included retinoid metabolism, Ig transcytosis, components of HLA class I complex, and AP-1 transcription factor subunits by PC HEPs. By contrast, PP HEPs were enriched in the expression of genes associated with amino acid catabolism-related ammonia lyase activity and iron homeostasis (Figure 2G). The relevance of several DEGs from the PC and PP signatures was further supported by the spatial distribution of the corresponding proteins in human liver tissue (Figure 2H).

Fig. 2

Liver zonation drives the clustering of human hepatocytes.

(A) Immunostaining for GLUL, CYP2E1, ASS1 and ASL (images from the Human Protein Atlas). (B) HEPs ordered along the portocentral axis according to the expression of GLUL, CYP2E1, ASS1, ASL, and ALB, per donor. (C) Annotation of HEPs as PC, MZ, or PP HEPs based on the expression of ≥4 of the zonation markers GLUL, CYP2E1, ASS1, ASL, and ALB. (D) Distribution of the expression of genes significantly upregulated in (i) PC or (ii) PP HEPs. (E) HEPs ordered along the portocentral axis according to the expression of genes significantly upregulated in PC or PP HEPs. (F) Annotation of all HEPs as PC, MZ, or PP based on the expression of genes significantly upregulated in PC or PP HEPs. (G) Distribution of the expression of genes associated with zonated functions in HEPs. (H) Protein-level validation of a selection of differentially expressed genes from the PP and PC HEP gene signatures (images from the Human Protein Atlas). HEP, hepatocytes; MZ, midzonal; PC, pericentral; PP: periportal.

Human HSCs are transcriptionally heterogeneous and spatially zonated

The above analysis enabled the recovery of 246 HSC transcriptomes that clustered into 2 clearly distinct subpopulations (i.e. HSC1 and HSC2) (Figure 3A). The HSC identity of both subpopulations was confirmed by the specific expression of genes encoding well-established HSC markers, including RGS5, PTN, NGFR, LRAT, FBLN5, DPT, DCN, CYGB, and COLEC11 (Figure 3B and Table S6).,,,,[37], [38], [39], [40], [41], [42], [43] These data confirmed recently proposed markers of HSCs, such as OLFML3 and TPM2, and revealed novel genes encoding quiescence-defining features, including QSOX1, required for the incorporation of laminin into the extracellular matrix (ECM), and STEAP4, a gatekeeper of normal metabolic function in fat-storing cells (Figure 3B). HSCs have the capacity to undergo a process of activation during which they acquire a fibrogenic, myofibroblast-like phenotype. To assess whether HSC clustering was driven by a difference in activation status, the expression in both subpopulations of key activation marker genes (i.e. ACTA2, COL1A1, COL1A2, COL3A1, LOX, NOTCH3, PDGFRA, TIMP1, and VCL) and quiescence-associated genes (i.e. BAMBI, ETS1/2, GATA4/6, IRF1/2, LRAT, LHX2, NOTCH1, PLIN2, PPARG, RBP1, SPARCL1, and TCF21) was analysed (Figure 3C).,[47], [48], [49], [50], [51], [52], [53], [54], [55], [56], [57], [58] Although there were differences in expression, those differences were not uniformly distributed, suggesting that the activation status of the cells is not the primary determinant of clustering. Cells of both clusters also displayed similar cell cycle gene expression profiles (Figure S1C). Both subpopulations differentially expressed 340 genes (Figure 3D and Table S7), and were characterised by specific gene signatures differentiating them from each other as well as from the other analysed liver cell populations (Figure 3E, Figure S3, and Table S7). HSC1 specifically expressed the cell surface proteoglycan GPC3, the neurotrophic receptor NTRK2 as well as other genes previously reported to be expressed in human or murine HSCs, including EFEMP1, GEM, CCL2, and THBS1.,, By contrast, HSC2 expressed significantly higher levels of the dopamine norepinephrine-converting enzyme DBH, the hedgehog signalling modulator HHIP as well as several G-protein-coupled receptors targeted by vasorelaxation peptide hormones and previously shown to be abundantly expressed by both murine and human HSCs (i.e. VIPR1, PTH1R, RAMP1, EDNRB, and AGTR1A).

Fig. 3

Transcriptomic and spatial heterogeneity of human HSCs.

(A) UMAP clustering of 246 HSCs. (B) Expression levels of HSC-specific genes in LPs, MPs, ECs, CCs and HEPs, as well as in the HSC subpopulations HSC1 and HSC2. (C) Expression levels of key HSC activation and quiescence marker genes in HSC1 and HSC2. (D) Heatmap of the significantly differentially expressed genes between HSC1 and HSC2. (E) UMAP plots for a selection of 30 HSC1- and HSC2-specific genes. (F) Representative pictures of brightfield in situ hybridisation using GPC3- and DBH-specific probes in human liver tissue. (G) Representative pictures of fluorescent in situ hybridisation using GPC3-, DBH- and NGFR-specific probes in human liver tissue. CC, cholangiocyte; CV, central vein; EC, endothelial cell; HEP, hepatocyte; HSCs, hepatic stellate cells; LP, lymphoid cell; MP, macrophage; PC, pericentral; PP, periportal; PV, portal vein; UMAP, uniform manifold approximation and projection.

The HSC1- and HSC2-specific gene signatures were used to determine whether both subpopulations could be robustly identified in scRNA-seq data sets from additional donors. To do so, healthy liver mesenchyme scRNA-seq data published by Ramachandran et al. (Figure S4A) were imported and the HSC transcriptomes were extracted by opposing the expression of a panel of vascular smooth muscle cell- (i.e. CALD1, CNN1, HEXIM1, MYH11, MYLK, SMTN, and TAGLN)[61], [62], [63] and HSC-associated genes (i.e. LRAT, RGS5, NGFR, PTN, FBLN5, DPT, CXCL12, STEAP4, SPARC, QSOX1, OLFML3, IGFBP3, IGFBP7, HGF, DCN, CYGB, COLEC11, and BGN) (Figure S4B,C). Three HSC clusters were identified (i.e. clusters 2, 5, and 6) (Figure S4D), of which cluster 2 and 5 comprised cells originating from 5 different donors, whereas cluster 6 originated from 1 donor and had a lower library size compared with clusters 2 and 5 (Figure S4E,F). Comparative analysis of the 3 HSC clusters for HSC1- and HSC2-specific genes (Figure S4G,H) found clusters 2 and 5 to be enriched for HSC1- and HSC2-specific genes, respectively. When merged, cells from both data sets co-localised (i.e. HSC1 with cluster 2 and HSC2 with cluster 5) (Figure S4I). No overlap was found with cells from cluster 6.

Transcriptomic alterations in response to tissue dissociation can lead to faulty identification by scRNA-seq of cell populations that do not exist in vivo. To validate the existence of HSC1 and HSC2 and query their distribution in the native tissue, RNA in situ hybridisation (ISH) was performed using liver tissue from 5 different donors and GPC3- and DBH-specific probes, part of the HSC1- and HSC2-specific gene signatures, respectively (Figure 3F and Figure S5). GPC3 and DBH cells were found in all donor tissues. Their expression was mutually exclusive and displayed a different distribution pattern across the liver lobule. GPC3 cells were mostly restricted to the portal and central vein areas, whereas DBH cells were distributed more diffusely in the perisinusoidal space, but not in the portal and central vein areas. The perisinusoidal distribution pattern of DBH cells was confirmed using a HHIP-specific probe, a second HSC2-specific gene (Figure S6). To validate whether GPC3 and DBH cells are both HSC subpopulations in the native tissue, multiplex RNA-fluorescence ISH was performed using NGFR-specific probes (Figure 3G). GPC3 and DBH expression co-localised with NGFR+ cells, confirming their HSC identity. Given that a different distribution pattern was observed for HSC1 and HSC2 in the human liver, the expression of genes reported to discriminate murine portal vein-associated HSCs and central-vein associated HSCs was analysed (Figure S7). Although several genes were significantly differentially expressed between the HSC1 and HSC2 cell populations, including NGFR, IGFBP3, TNFRSF11B, TPM2, CCL2, GEM, CXCL14, and RSPO3, a consistent pattern was not found. Of note, NGFRGPC3DBH cells were also found, suggesting the existence of additional transcriptional HSC subpopulations in the native tissue.

To investigate the potential functional implications of these observations, a comparative GO analysis was conducted to identify specifically enriched gene networks in HSC1 and HSC2 (Table S9). This analysis revealed genes related to different physiological processes, including antigen presentation in HSC2 (HLA-DRB1, HLA-DRA, CD74, HLA-DQB1, HLA-DPA1, HLA-DPB1, and HLA-DRB5), and the metabolism of GAGs (PIM1, LUM, SDC2, GPC3, OGN, VCAN, UGDH, CD44, HAS1, and HAS2) and elastic fiber constituents in HSC1 (FBLN5, ELN, FBLN1, and MFAP4) (Figure 4A).

Fig. 4

Functional heterogeneity of human hepatic stellate cells.

(A) Distribution of the expression of genes associated with selected functions in HSCs. (B) Representative pictures of VIM and HLA-DR immunostaining in liver tissue. (C) Representative pictures of CD45 and HLA-DR immunostaining in liver tissue. (D) Distribution of the expression of CD80 and CD86 in HSCs. (E,F) Representative pictures of GAG and Elastin staining in liver tissue. (G) Expression levels of liver extracellular matrix glycoproteins, collagens and proteoglycans in LPs, MPs, ECs, CCs, and HEPs, as well as in the HSC subpopulations HSC1 and HSC2. (H) Expression levels of relevant growth factors and cytokines in HSC1 and HSC2. CC, cholangiocyte; CV, central vein; EC, endothelial cell; GAGs, glycosaminoglycans; HEP, hepatocyte; HLA, human leukocyte antigen; HSCs, hepatic stellate cells; LP, lymphoid cell; MP, macrophage; PC, pericentral; PP, periportal; PV, portal vein; VIM, vimentin.

To assess the relevance of these findings, co-immunostaining was performed for HLA-DR and the HSC marker VIM (Figure 4B)., VIM+HLA-DR+ cells were distributed throughout the liver lobule, without matching particularly to either HSC1 or HSC2 localisation. Although VIM potentially also marks MPs, the existence of CD45-HLA-DR+ cells was confirmed in the different areas of the liver lobule (Figure 4C), including in portal tracts, a finding that is in line with a previous report. However, no expression of the HLA class II co-stimulatory molecules CD80 and CD86 was found in either of the 2 HSC subpopulations (Figure 4D), calling into question their antigen-presenting capacity under physiological conditions.

By contrast, both GAGs and elastin were found to accumulate in PC and PP areas along big vessels (Figures 4E,F), which correlated with the location of HSC1. Given that HSCs are considered the main ECM-producing cell type of the liver, this prompted the evaluation of the expression of all glycoproteins, collagens and proteoglycans reported to be constituents of the human liver ECM.[68], [69], [70] Gene expression levels were detected for >90% of those proteins and ~45% were specifically expressed by HSCs compared with all other analysed cell types (Figure 4G). Interestingly, of those specific to HSCs, ~60% were specifically expressed by HSC1 or HSC2, indicating that both subpopulations collaborate in the synthesis and organisation of ECM components. Akin to this observation, HSC1 and HSC2 expressed different levels of key secreted cytokines (CCL2, CCL21, and IL32), chemokines (CXCL12 and CXCL14), angiopoietins (ANGPTL1, ANGTPTL2, and ANGPTL6), and mitogens (HGF, HDGF, VEGFC, and PGF) (Figure 4H). Altogether, these results show that HSCs in the human liver are heterogeneous, spatially zonated and characterised by unique gene expression signatures suggestive of important functional differences.

Discussion

The current study examined the near-native transcriptome of thousands of freshly isolated human liver cells in an unbiased way by scRNA-seq to unravel the heterogeneity of the main hepatic cell populations. In particular, these data highlight metabolic zonation as a key determinant of the HEP transcriptomic heterogeneity and, to our knowledge for the first time in human, transcriptomic, and functional heterogeneity of HSCs.

HEPs, through their organisation in concentric circles around the central vein, are exposed to various biochemical gradients that regulate the expression of genes and induce functional heterogeneity. The current study generated, to the best of our knowledge, the largest human HEP droplet-based scRNA-seq data set currently available. Based on the zone-specific HEP markers validated by Halpern et al. and the subsequently identified PC and PP gene signatures, the transcriptomes were classified according to their predicted position in PC, MZ, and PP areas. The relevance of this approach was validated at the protein level by immunohistochemistry. This classification confirmed the PP expression of HAL and SDS, as previously reported in mice. Both enzymes catalyse a reaction releasing ammonia, which can feed the urea cycle and, in the case of SDS, pyruvate, which can feed gluconeogenesis (i.e. 2 PP metabolic pathways)., Genes encoding secreted plasma proteins and, in particular, genes associated with iron transport (TF and HP) and its regulation (HAMP) showed increased expression in PP HEPs, whereas genes associated with iron storage (FTH1 and FTL) were upregulated in PC HEPs. PC HEPs also displayed higher levels of drug metabolism-related genes, including CYP1A2, CYP2E1, CYP3A4, and CYP3A5 (i.e. haemic iron-containing enzymes), with the notable exception of CYP2A6, CYP2A7, and CYP3A7, which were higher in the PP area. PC CYPs can also take part, with additional enzymes, in retinoid and lipid metabolism., Although bile production has been previously zonated in PC HEPs, the current study revealed the PC upregulation of genes associated with IgA transport from blood to bile (PIGR and FCGRT). In the current data set, 19% of the genes from the PC HEP signature (including AP-1 subunits FOS and JUN) and 31% of the genes from the PP HEP signature were reported as up- and downregulated by Wnt signalling in mice, respectively., Thus, the identified zonated signatures are in line with previous reports describing the Wnt centroportal gradient as a major regulator of HEP zonation. Additional studies are required to confirm the relevance of identified zonated pathways and their interplay from a functional standpoint.

Although scRNA-seq offers the advantage of in silico cell sorting, the optimal isolation of HSCs from liver tissue requires a pronase digestion step, which is incompatible with the isolation of HEPs because it induces their destruction.,, The absence of pronase in the current digestion protocol probably negatively affected the yield of HSCs in the total cell suspension. Nevertheless, this data set still presents 1 of the largest collection of individual HSC transcriptomes sourced from human liver tissue.

Similarly to previous reports in murine and porcine HSCs,, human HSCs were found to be transcriptionally heterogeneous and spatially zonated. To the best of our knowledge, this is the first report of heterogeneity in the expression of HLA class II genes in human HSCs, and the localisation HLA-DR+ HSCs in the portal area, the entry point of exogenous antigens. Previous observations in primary murine, rat, and human HSCs showed that basal expression of HLA class II decreases upon myofibroblastic activation in culture but could be restored by proinflammatory cytokines, including interferon γ.[79], [80], [81], [82] HSC-derived myofibroblasts in proinflammatory conditions were shown to be capable of exogenous antigen internalisation, to express co-stimulatory molecules,,, and to elicit a HLA class II-dependent T cell response.[79], [80], [81], [82] Additional studies further reported HLA class II-independent stimulatory or immunomodulatory functions of HSCs. Of note however, those experiments were performed on culture-activated HSCs and the physiological relevance of the antigen-presenting cell properties of HSCs in human liver physiology requires further investigation. In the current data set, expression of co-stimulatory CD80 and CD86 was not detected. Therefore, the role of HLA class II molecules in HSCs in physiological conditions requires more investigation.

The current work also suggests that the synthesis of a significant number of key liver ECM components is reserved to HSCs, and that this labour is asymmetrically distributed over subpopulations with differential localisation within the liver lobule. This finding is of particular interest, given the central role of HSCs in the development of fibrosis, and in light of the recent discovery that murine HSC subpopulations unequally contribute to the development of CCl4-induced liver fibrosis.

In conclusion, the current work contributes to recent efforts to achieve a more refined understanding of the heterogeneity and intralobular organisation of human liver cells, and calls for further research on the implications of HSC heterogeneity in health and disease.

Financial support

This research was funded by grants from La Wallonie (SPW Recherche DGO6 WALInnov program, convention no. 1710035) and the Brussels-Capital Region (Innoviris, RBC/2020-SPIN-206), together with the support of Promethera Biosciences SA. Promethera Biosciences SA had no influence on the design, the collection, analysis and interpretation of data, the writing of the report, or the decision to submit the article for publication.

Authors’ contributions

Acquired financial support for the project leading to this publication: W.C., B.C., E.M.S. and A.E.T. Conceptualised the study: W.C., B.C. and A.E.T. Designed the experiments: V.L.P., A.L. and A.E.T. Performed the experiments and/or collected the data: V.L.P., A.L., N.A.S., L.K., M.D., M.C., M.N. and A.E.T. Wrote the manuscript: V.L.P., A.L., N.A.S., and A.E.T. Reviewed and edited the manuscript: all authors.

Data availability statement

All raw sequencing data are available in GEO database under accession GSE158723.

Conflict of interest

The authors declare no conflicts of interest that pertain to this work.

Please refer to the accompanying ICMJE disclosure forms for further details.