A protocol for human liver dissociation for scRNA-seq

A central problem in liver tissue dissociation is that hepatocytes and cholangiocytes are sensitive to cell death, and other cells activate, in response to cell manipulation15,16,17. We developed a cell isolation approach (Fig. 1a, b) without density gradient or column purification, or flow cytometry. We found that cell sorting approaches led to loss of fragile cells, particularly hepatocytes and cholangiocytes, and the relative enrichment of more robust NPCs, such as endothelial cells/KCs. Preliminary scRNA-seq experiments showed an under-representation of NPC cell populations in the sorted NPC-enriched fraction compared to total liver homogenate (TLH) (Supplementary Fig. 1). Specifically, the NPC fraction contained two endothelial cell clusters compared to three in the TLH. Furthermore, when we targeted 6000 viable cells for scRNA-seq from TLH and NPC fractions, the NPC fraction yielded far fewer viable cells than the TLH (Supplementary Fig. 1) indicating liver cell sensitivity to flow sorting.

Fig. 1 ScRNA-seq profiling of parenchymal and non-parenchymal cells from the human liver. a Overview of the single-cell isolation and analysis workflow. Workflows for b the dissociation of human caudates to single-cell homogenates, c the generation of scRNA-seq cDNA expression libraries using the 10x Genomics genomics platform and, d data analysis Full size image

Using our liver dissociation protocol, cell viability of the TLHs obtained from five caudate lobes ranged from 49 to 90% viable by trypan blue exclusion (Supplementary Fig. 2a). However, the actual number of cells profiled per caudate after filtering for library size and mitochondrial transcript percentage (Fig. 2a) ranged from 1073 to 3255 cells per sample (between 17.9 and 54.3% passing quality control). The mean library size (number of UMIs detected/cell) for TLHs was 5227 (range 3122–6043 UMIs) and the mean number of genes detected per cell was 1313 (range 906–1537 genes) (Supplementary Figs. 2, 3). The variation in numbers of cells profiled may be attributed to differences in sample viability and technical differences in cell capture rates for each scRNA-seq run. Standard scRNA-seq filtering excludes cells with a high ratio of reads from mitochondrial genome transcripts, indicating potential plasma membrane rupture and dissociation-based damage18. This filter is often set at 10%, but hepatocytes can have very high mitochondrial content19, thus we chose a threshold of 50% to optimize keeping hepatocytes and removing dead and dying cells. Evaluation of six mitochondrial cut-offs ranging from 10–60% (Supplementary Figs. 4, 5) showed that all clusters (except cluster #6 at the 10% cutoff) are identified as unique populations in t-SNE plots at all cut-offs, indicating that our map is robust to this threshold. As expected, hepatocytes were most susceptible to removal by this filter (Fig. 3), while endothelial cells, macrophages were consistently represented in all livers profiled. Cell doublets are normally detected in scRNA-seq experiments as cells with abnormally high library size. However, we did not detect a natural threshold in library size per cell that we could select to remove doublets. Further, there are naturally occurring binucleated hepatocytes that are expected to have a high library size, which would make it difficult to distinguish between doublets and binucleated cells. A doublet filter set to remove cells with the top 99.9% library size (following standard protocols20), mostly removed hepatocyte (cluster #14) and plasma cell (cluster #7) populations. As both binucleated hepatocytes21 and tissue trafficking plasmablasts22 have been previously described, we suspect that the cells identified as “doublets” may be single biological cells containing two nuclei and thus did not apply doublet filtering to our map (Supplementary Fig. 6).

Fig. 2 20 distinct cell populations were revealed in healthy human livers. a Viable cells were identified from the single-cell libraries having a minimum library size of 1500 transcripts and a maximum of 50% mitochondrial transcript proportion. b t-SNE projection of 8444 liver cells (each point represents a single cell). Cells are colored by library size, with darker colors indicating larger libraries. c t-SNE projection where cells that share similar transcriptome profiles are grouped by colors representing unsupervised clustering results. d Heat map analysis using known gene expression profiles of hepatocytes/immune cells. The identity of each cluster was assigned by matching the cluster expression profile with established cell-specific marker gene expression for hepatocytes, endothelial cells, cholangiocytes, and immune cells. e Cell-cycle phase prediction showed that hepatocyte clusters were less proliferative than immune cell clusters. f Cluster map showing the assigned identity for each cluster defined in c. The cluster number of each potential cell population is indicated in parentheses. DE: differentially expressed, MACs: macrophages, PCA: principal component analysis, t-SNE: t-distributed stochastic neighbor embedding, PCs: principal components Full size image

Fig. 3 Contribution of cells to each scRNA-seq cluster by sample and subpopulation analysis. a The proportion of cells that contributed to each cluster by liver sample. b t-SNE projection of all cells, colored by the source donor, and labeled with cluster number. Most cell-type associated clusters are made up of multiple donors. Hepatocyte clusters, on the other hand, appear to segregate by donor. c Proportions of mature vs. plasma B cells across five liver samples as a percentage of total B cells. Similar subpopulation analyses were carried out for d αβ & γδ T cells e Macrophages (Macs), f Endothelial cells, and g Natural killer (NK)-like cells. LSECs liver sinusoidal endothelial cells Full size image

The landscape of cells in the healthy human liver

We applied our liver single-cell dissociation protocol and scRNA-seq technology to identify resident liver cells in TLHs from NDD liver transplant donors (donor characteristics found in Supplementary Table 1). Importantly, NDD grafts are subjected to the systemic inflammation that accompanies brain death and are thus themselves mildly inflamed23. However, we confirmed normal histological patterns in these livers (Supplementary Fig. 7). Pooling results from all donors, we captured 8444 cells (after filtering out low viability cells) that clustered into 20 discrete cell populations (Fig. 2c, f, Supplementary Data 1, 2) that are described below.

Hepatocytes

Hepatocytes are the building blocks of the liver and play a role in detoxification, lipolysis, and gluconeogenesis24. In mice and human, hepatocytes have been described as having distinct functions based on their location in the hepatic acinus, known as their zonation12,24. We found six distinct hepatocyte populations, comprising Clusters 1, 3, 5, 6, 14, and 15 (Figs. 2f, 4ci), that were generally less proliferative cells (Figs. 2e, 4a), and showed enriched ALB (Albumin) expression, a hallmark of hepatocytes.

Fig. 4 ScRNA-seq analysis of hepatocyte populations. a Distribution of hepatocytes by cell-cycle phase (G1, G2/M, S) and hepatocyte cluster (1, 3, 5, 6, 14, 15). b Box plot of library size for each hepatocyte cluster with median library size (top) and graphically denoted median (dark horizontal line). Outliers (black dots) and interquartile range (black box) are indicated. c t-SNE plots showing the expression of general hepatocyte markers based on PCA clustering of 8444 cells. c (i) ALB, c (ii) HAMP, c (iii) ARG1, c (iv) PCK1, c (v) AFP, c (vi) BCHE. Legend for relative expression of each marker from lowest expression (yellow dots) to highest expression (purple dots) (top left). c (vii) t-SNE projection showing a reference map of all six hepatocyte clusters. d Pathway enrichment analysis examining active cellular pathways in clusters 1, 3, 5, 6, 14 & 15. The size of the nodes represents the number of genes in a particular pathway. Highly related pathways are grouped into a theme (black circle) and labeled in Cytoscape (Version 3.6.1). Intra- and inter-pathway relationships are shown (green lines) and represent the number of genes shared between each pathway. Periportal and pericentral zones are assigned in relation to correlation analysis between mouse and human liver in Supplementary Fig. 8. Statistical significance of the correlation between mouse liver layers and human liver clusters (denoted under each pathway analysis) calculated using Pearson correlation. ∗∗∗P < 0.001, ∗∗P < 0.01, ∗P < 0.05. t-SNE t-distributed stochastic neighbor embedding Full size image

To directly infer hepatocyte cluster zonation patterns, single-molecule FISH (smFISH) and laser capture microdissection studies are required, these techniques were not a part of this study. However, a systematic comparison was carried out examining gene expression patterns in the identified human hepatocyte clusters with respect to the zonated gene expression patterns previously shown in mouse12 (Supplementary Fig. 8, Fig. 4d). This analysis revealed that gene expression patterns in four human hepatocyte clusters correlated significantly with zonated gene expression patterns identified in the mouse sinusoid (Supplementary Fig. 8, Supplementary Data 3). Two clusters showed weak correlation, possibly due to differences in the genes that define mouse and human zonated liver expression patterns. Specifically, human Cluster 5 correlated with the most periportal mouse liver layer, while human Cluster 3 correlated with the most central venous mouse liver layer. Human Cluster 1 was correlated with mouse layers 2 and 3 (more periportal) and Cluster 15 correlated with an interzonal mouse layer (Layer 4). Clusters 6 and 14 did not correlate significantly with any mouse layers, possibly because the top DE genes defining these clusters were not found in the top 94 genes which defined zonation in mouse livers.

To study the function of the human clusters, we applied pathway analysis to identify active cellular pathways in each cluster (Fig. 4d, Supplementary Data 4). Characteristic functional pathways in sinusoidal zones have been well described12,24. Using these approaches, we describe the possible zonation for human hepatocyte clusters with the expectation that natural hepatocyte heterogeneity will be better addressed as more human liver samples are profiled across multiple study sites.

Zone 1 (periportal) hepatocytes have a known role in gluconeogenesis and β-oxidation24. Cluster 5 showed enriched expression of genes involved in lipid and cholesterol synthesis such as SCD (Supplementary Fig. 9), HMGCS1, and ACSS2 with enhanced expression of PCK1 (Fig. 4c.iv) and CPS1 (Supplementary Data 1, Supplementary Fig. 9), genes that are enriched in the periportal hepatocytes of mice12, suggesting that this cluster may represent Zone 1 hepatocytes. A further indication of a periportal nature is enriched expression of urea cycle gene ARG1 (Fig. 4c.iii). (Top DE genes cluster 5: SCD, HMGCS1, ACSS2, TM7SF2, TMEM97, CP, CRP, SLPI, C2orf82, ACAT2, TM4SF5, MSMO1, LEPR). Pathway analysis revealed that Cluster 5 was enriched for pathways characteristic of liver periportal function including cholesterol and sterol biosynthesis along with numerous active immune pathways (Fig. 4d).

In mice, Zone 3 or central venous hepatocytes play a role in drug metabolism and detoxification, constitutively expressing high levels of cytochrome P450 enzymes12. Clusters 3 and 1 were enriched for the expression of BCHE (butyrylcholinesterase), involved in drug metabolism (Cluster 3 top DE genes: BCHE, G6PC, GHR, ALDH6A1, RCAN1, AR, RP4-710M16.2, LINC00261, PLIN1, RP11-390F4.3) (Fig. 4c.vi.). Cluster 1 also showed enriched expression of HAMP (Fig. 4c.ii), a gene that is enriched in mouse midzonal hepatocytes12. (Cluster 1 top DE genes: RPP25L, HSD11B1, HAMP, GHR, APOM, APOC4-APOC2, TKFC, G6PC, G0S2, PON3, C1orf53, TTC36, GOLT1A, RCAN1, RP4-710M16.2, FST, MCC, AQP9, PLIN1). In the comparison with zonated mouse liver data, Cluster 3 correlated most significantly with the most central venous mouse sinusoid layer (Supplementary Fig 8). Cluster 1 showed a significant correlation with mouse layers 2 and 3, suggesting that Clusters 3 and 1 may be more central venous. Pathway analysis revealed that the cells making up these two clusters were active in cellular pathways characteristic of zone 3 functions in mice and human including P450 pathways, drug metabolism, Wnt activation, hypoxia, amino acid biosynthesis, and glycolysis12,24 (Fig. 4d), supporting the notion that these cells might have a central venous origin.

Cluster 15 had the smallest library size of the hepatocyte clusters (2707 UMIs per cell) (Fig. 4b) (Top DE genes: HPR, GSTA2, AKR1C1, MASP2, NNT, SAA4, MRPS18C, OCIAD1, APOA5, TTR). Gene expression patterns in this cluster correlated significantly with mouse sinusoid layer 4, suggesting that these cells may be interzonal hepatocytes (Supplementary Fig 8). Future work will be needed to clarify the origin and role of these cells.

Activated cellular pathways in Cluster 14 included complement activation and immune activation, known periportal functions in mice12,24 (Top DE genes in Cluster 14: CYP2A7, ENTPD5, CYP3A7, CYP2A6, C4B, EID2, TP53INP2, SULT1A1, ATIC, SERPINH1, SAMD5, GRB14) (Supplementary Data 1&2). While this cluster is most transcriptionally similar to periportal mouse layers (Supplementary Fig. 8), there was no significant correlation to known layers, likely because the top DE genes (CYP2A genes) in this cluster have different roles in mouse and human25,26 and did not map well with the 94 genes that defined the mouse sinusoidal layers.

Cluster 6 hepatocytes (DE genes: SEC16B, SLBP, RND3, ABCD3, RHOB, EPB41L4B, GPAT4, SPTBN1, SDC2, PHLDA1, WTAP, ACADM) have activated cellular pathways in immune cell activation, fibrinolysis, and triglyceride biosynthesis, suggesting a periportal origin. However, since this cluster did not correlate with the mouse sinusoid regions (Supplementary Fig. 8), the origin of this cluster requires further examination.

Our analysis may suggest new aspects of hepatic stem cell biology. The liver displays a high regenerative capacity involving resident stem cells. There is evidence that hepatic stem cells can originate from niches in the canal of Herring, the terminal branch of the intrahepatic biliary system27. These stem cells are called hepatic progenitor cells in human and oval cells in rodents. Alpha-fetoprotein (AFP) is described as a marker for hepatic stem cells in human28. In our current understanding of the zonated rodent liver, the perivenous compartment is resistant to proliferative signals29, while the periportal region contains oval cells30. Our study challenges this paradigm, since AFP is expressed in all clusters of hepatocytes except Cluster 14 (Fig. 4c.v)—we were unable to identify a discrete cluster specific to hepatic progenitor cells. We selected hepatocytes with detectable AFP expression, and generated a ranked list of genes DE in AFP− vs AFP+ hepatocytes (Supplementary Data 5), which we used to examine enriched cellular pathways cells in the AFP-expressing cells. Enriched pathways in AFP+ cells included those for cellular division and IL-6/7 signaling (Supplementary Fig. 10, Supplementary Data 5). This supports the notion that these cells may be hepatic stem cells since IL-6 is a key cytokine for the proliferation of hepatocytes31. Our findings raise the possibility of a less localized and more heterogeneous model of hepatic progenitor cells in the human liver.

Liver endothelial cells

Our current knowledge of human LSECs is limited to immunofluorescence, flow cytometric studies32,33, bulk RNA-seq studies on sorted LSECs34, and functional properties assessed during brief in vitro cultures35,36. The liver vascular endothelium, made up of LSECs and the endothelium of blood vessels, provides a dynamic barrier between the blood and the liver microenvironment. LSECs are defined as scavenger endothelia that use clathrin-mediated endocytosis for the clearance of macromolecules from the blood36. Recently, immunofluorescent staining was used to describe the zonation of human LSECs and a population of CD36hiCD32B−CD14− LYVE-1− LSECs in Zone 1 of the hepatic acinus (the periportal area). A separate LYVE-1+, CD32Bhi, CD14+, CD54+, CD36mid-lo LSEC population was located in Zones 2 and 332 (central venous). While this work validates the notion of zone-specific heterogeneity within LSECs, it remains limited by the number of tested antigens and availability of validated immunological reagents.

We observed three endothelial cell populations (clusters 11–13; Figs. 2f, 3a, b, f, 5; Supplementary Fig. 11), which were less proliferative than immune cells (Figs. 2e & 5a), and expressed CALCRL (Fig. 5c.i) and RAMP2, suggesting sensitivity to adrenomedullin signaling37. The most abundant endothelial cell cluster (Cluster 11) displayed enriched expression of F8, PECAM1, with little expression of CD32B, LYVE-1, STAB2, and CD14 (Fig. 5c, Supplementary Fig. 11). (Top DE genes: MGP, SPARCL1, TM4SF1, CLEC14A, ID1, IGFBP7, ADIRF, CTGF, VWF, CD9, C7, SRPX, ID3, CAV1, GNG11, AQP1, HSPG2, EMP1, SOX18, CLDN5). In line with previous work32, we propose that these endothelial cells are likely periportal LSECs (Zone 1). The second most abundant endothelial cell population (Cluster 12) was characterized by the enriched expression of CD32B, LYVE1, STAB2, with little expression of VWF (Top DE genes: CCL14, CLEC1B, FCN2, S100A13, FCN3, CRHBP, STAB1, GNG11, IFI27, CLEC4G, CLDN5, CCL23, OIT3, RAMP3, SGK1, DNASE1L3, LIFR, SPARC, ADGRL4, EGFL7, PCAT19, CDKN1C) (Fig. 5c.ii–iii). Based on prior histological examinations of LSEC zonation32, we propose that this is an LSEC cluster that is central venous in origin (Zone 2/3). The least abundant hepatic endothelial cell population (Cluster 13) was characterized by low or no expression of LSEC markers (LYVE1, STAB2, CD32B). These cells are likely non-LSEC endothelial cells including central vein and portal arterial and venous endothelial cells based on the expression of ENG (protein alias CD105) and PECAM1 (protein alias CD31) as has been described in the human liver via immunohistochemistry32. (Top DE genes: RAMP3, INMT, DNASE1L3, LIFR, PTGDS, C7, CTGF, TIMP3, RNASE1, ID3, ENG, MGP, PCAT19, HSPG2, GPM6A, PTPRB, VWF, FAM167B, SRPX, LTC4S, IFI27). To confirm our DE findings at the protein level and to strengthen our assertions regarding LSEC zonation, we examined the expression of CD32B by immunohistochemistry and found that CD32B protein expression was concentrated in the central venous (Zones 2/3) areas of the sinusoid, including the central vein, with limited periportal staining (Fig. 5e.i–vi), confirming previous findings32. A pairwise pathway analysis comparing all endothelial cell clusters to each other revealed that clusters 11 and 13, periportal LSECs and portal endothelial cells, respectively, were functionally very similar (Supplementary Fig. 12), with periportal LSECs having significantly enriched pathways in vessel development, cell-cycle arrest, cell junction, and apoptosis. Both periportal LSECs and portal endothelial cells showed enriched pathways in translation, targeting ER, and TNF activation in comparison to central venous LSECs (Fig. 5d, Supplementary Fig. 13). Central venous LSECs (Cluster 12) displayed highly enriched immune pathways, including innate immunity, phagocytosis, leukocyte activation, bacterium defense, in comparison to Clusters 11 and 13 (Fig. 5d, Supplementary Fig. 13). These results present a rich marker description and functional profile of human liver endothelial cells. However, the characterization of endothelial cell clusters as sinusoidal, rather than arterial or vascular endothelium, requires visualization of fenestrations and specific scavenger receptor functions38 that are uniquely characteristic of sinusoidal endothelial cells. These descriptions will require optimization of LSEC culture conditions, identification of surface markers that distinguish these populations, and confirmation of protein expression of these markers by methods such as CITE-seq39, followed by flow sorting and in vitro culture to visualize fenestrations.

Fig. 5 ScRNA-seq analysis of LSEC populations. a Distribution of LSECs by cell-cycle phase (G1, G2/M, S) and LSEC cluster (11, 12, 13). b Box plot of library size for each LSEC cluster with median library size (top) and graphically denoted median (dark horizontal line). Outliers (black dots) and interquartile range (black box) are indicated. c t-SNE projection of the expression of established LSEC markers in the three identified clusters. c (i) CALCRL c (ii) CD32B, and c (iii) VWF. d Pairwise pathway enrichment analysis of genes DE between clusters 11 and 12 defined in Fig 2f. Pathways enriched in periportal LSECs (Cluster 11) are labeled in red and pathways enriched in central venous LSECs (Cluster 12) are indicated in blue. Colored circles (nodes) represent pathways, sized by number of genes they contain. Green lines depict intra- and inter-pathway relationships according to the number of genes shared between each pathway. Black circles group related pathways into themes that are labeled. e (i) Immunofluorescence of CD32B distribution in liver zone 1 (portal vein - PV), and 2/3 (central vein - CV). CD32B stains mainly LSEC cluster 12 and CK19 stains periportal ductal cells. e (ii, iii, iv, v) are magnified sections corresponding to the indicated roman numerals (white) in (e (i)). e (i) Scale bar represents 200 μm, e (v) scale bar represents 20 μm. Staining was performed on HIER (10 mM Citrate pH 6.0, 95 °C,15 min) treated slides visualized using the matching donkey anti-host antibody and counterstained with DAPI. Slides were scanned and lobules defined as in Supplementary Fig. 10. e (vi) Quantification of percent CD32B positive cells in liver zones 1–3. Error bars show the standard error of the mean for at least 10 replicates. Statistical significance evaluated using a one-way analysis of variance (ANOVA) with a Bonferroni post-test ∗∗∗P < 0.001, ∗∗P < 0.01, ∗P < 0.05. t-SNE: t-distributed stochastic neighbor embedding Full size image

Cholangiocytes

Cholangiocytes are epithelial cells which line bile ducts that comprise 3–5% of the cells in the human liver and generate 40% of the total bile volume40. Intrahepatic cholangiocytes originate from hepatoblasts during embryonic development. Based on the duct diameter size, cholangiocytes are categorized as small and large, each with different secretory functions and sensitivity to injury20. The circumference of small bile ducts is formed by 4–5 cholangiocytes and by 10–12 cholangiocytes in larger ducts41. Recent studies have outlined single-cell transcriptomic profiles from rodent cholangiocytes42,43, however the transcriptional profile of human cholangiocytes has not been described due to challenges in isolating these cells44. Using our dissociation protocol, we were able to identify a cholangiocyte cell cluster in all five livers profiled (Figs. 2f, 3a, b). Differentially expressed genes in this cluster (Cluster 17) included SOX9, EPCAM, and KRT19 (Fig. 6a, b, f), in accordance with genes that were previously described as enriched in primary human cholangiocytes45. Importantly, we found a population of cells which displayed upregulated expression of genes encoding secretory and inflammatory pathways (Top DE genes: TFF2, SCGB3A1, FXYD2, KRT7, DEFB1, CD24, TFF3, LCN2, KRT19, CXCL1, PIGR, TFF1, CXCL6, LGALS2, TACSTD2, ELF3, SPP1, MUC5B, LGALS4). This transcriptional profile of cholangiocytes can serve as a reference point for assessing the physiological nature of human stem cell-derived cholangiocytes.

Fig. 6 ScRNA-seq analysis of the cholangiocyte population. t-SNE plots showing relative distribution of commonly expressed cholangiocyte genes in the healthy NDD liver. Protein alias described in parentheses if different from gene name. Expression of a KRT19 (CK19), b EPCAM, c FXDY2, d CLDN4, e CLDN10, f SOX9, g MMP7, h CXCL1, i CFTR, j TFF2, k KRT7 (CK7), l CD24. Legend for relative expression of each marker from lowest expression (yellow dots) to highest expression (purple dots) (top left). t-SNE: t-distributed stochastic neighbor embedding Full size image

Hepatic stellate cells

Hepatic stellate cells (HSCs) are found in the subendothelial space of the liver sinusoid, known as the space of Disse46. HSCs are the main storage site for Vitamin A and are major contributors to tissue fibrosis. Upon activation, human HSCs express α-smooth muscle actin (ACTA2) and begin to lay down extracellular matrix, which is composed of collagen (e.g., COL1A1, COL1A2). Cluster 20 was identified as HSCs based on the expression of ACTA2, COL1A1, TAGLN, COL1A2, COL3A1, SPARC, and the expression of retinol binding protein 1 (RBP1), a vitamin A-associated transcript (Figs. 2f, 7). These genes have been previously described as being upregulated during hepatic stellate cell activation in human liver46. Additional genes specific to this cluster include DCN, MYL9 (Fig. 7h–i), TPM2, MEG3, BGN, IGFBP7, IGFBP3, CYR61, OLFML3, IGFBP6, CCL2, COLEC11, CTGF, HGF (Top DE genes). Thus, our data confirms and extends the human HSC transcriptional signature.

Fig. 7 ScRNA-seq analysis of hepatic stellate cell population. t-SNE plots showing the relative distribution of commonly expressed hepatic stellate cell genes in the healthy liver. Protein alias described in parentheses if different from gene name. Expression of a ACTA2, b COL1A1, c TAGLN, d COL1A2, e COL3A1, f SPARC, g RBP1, h DCN, and i MYL9. Legend for relative expression of each marker from lowest expression (yellow dots) to highest expression (purple dots) (top left). t-SNE: t-distributed stochastic neighbor embedding Full size image

Zonation of human hepatic tissue: the cellular basis for zone-specific gene signatures

Recently, the zonation of human hepatic tissue was examined by laser capture microdissection and RNA profiling of uninvolved liver tissue, obtained during hepatectomy for metastatic tumors to the liver24. In that study, gene expression patterns were described for each zone of the liver lobule. The limitation of such a study is in not knowing the cellular origin of the genes identified: our study is thus complementary and determines the cellular origins of the zonal gene expression signatures. Among the genes previously found to be upregulated in zone 1 (periportal area)24, we show enriched expression of CRP, SLPI, LEPR, A2M, CHI3L1, HAL in Cluster 5 (which correlates with mouse periportal zones); C7, MGP, ID1, LDB2, CD9, AQP1, SOX18 in portal endothelial cells and zone 1 LSECs (Clusters 11 &13); and KRT7, KRT19, CXCL6, SFRP5, CLDN10, BICC1, AQP1, ERICH5 in cholangiocytes (known to be enriched in the periportal areas). Among the zone 3 genes previously identified24, RELN was enriched in zone 2/3 LSECs and CYP3A4 and ADIRF were enriched in Cluster 15 (Supplementary Data 1, 2).

Intrahepatic immune cells

Increasing evidence shows that the liver possesses complex immunological properties1,47,48,49. However, the identity, frequency, and phenotype of hepatic immune cells is largely unknown due to the challenges in extracting immune cells from human liver tissue in a non-biased and viable manner48. In mice, a population of KCs has been identified by scRNA-seq12, but additional immune cell populations were not described. Using our gentle fractionation methods and scRNA-seq analysis, one key goal of our study was to overcome prior challenges and expand the current knowledge of intrahepatic immune cell populations.

Intrahepatic monocytes/macrophages

The liver is the solid organ in the body with the largest population of tissue resident macrophages1,50. Tissue resident macrophages have been described as the immunological sentinels of the liver51, but because they are difficult to isolate from humans and have a complex ontogeny52 they are poorly understood47. Traditionally, macrophages are classified as having inflammatory or immunoregulatory properties, as determined by their surface marker phenotype or cellular functions53. We previously employed flow cytometry to quantify expression of traditional inflammatory or regulatory macrophage surface markers and found a spectrum of expression in freshly isolated human KCs13. Surprisingly, our scRNA-seq analysis consistently revealed the presence of two distinct populations of intrahepatic CD68+ macrophages (Figs. 2f, 3a, b, Supplementary Fig. 14) in all livers studied. CD68+ macrophage population 1 (Cluster 4; Fig. 2f), was characterized by enriched expression of LYZ, CSTA, CD7454, suggesting that this cluster represents inflammatory macrophages (Top DE genes: S100A8, LYZ, S100A9, HLA-DPB1, S100A12, RP11-1143G9.4, EVI2A, HLA-DPA1, VCAN, S100A6, CXCL8, HLA-DRA, MNDA, TYROBP, HLA-DRB1, FCN1, HLA-DQA1, IL18, C1QC, CD74, HLA-DRB5).

CD68+ macrophage population 2 (Cluster 10; Fig. 2f) was characterized by enriched expression of CD5L, MARCO, VSIG4, CPVL, CD163, CCDC88A, C5AR1, LIPA, LILRB5, MAF, CTSB, MS4A7, VMO1, RAB31, SLC31A2, TTYH3, VCAM1, KLF4, HMOX1, AIF1l, TMIGD3 (top DE genes). The DE genes in this cluster suggest that these macrophages have a tolerogenic function. For example, VSIG4 is a co-inhibitory ligand and, in mice, is required to maintain an intrahepatic tolerogenic milieu55. Similarly, HMOX1 (hemoxygenase) knockdown in mice leads to hepatic inflammation56. Confirming the unique functionality of the KC populations, pathway analysis revealed that KCs in cluster 10 were enriched for pathways related to tolerance while KCs in cluster 4 were enriched for inflammatory pathways (Fig. 8c).

Fig. 8 ScRNA-seq identifies two distinct populations of human liver resident macrophages/monocytes. a, b t-SNE projection of 8444 liver cells, with each cell colored based on expression of a CD68 and b MARCO. c Pairwise pathway enrichment analysis comparing gene expression in the two CD68+ macrophage clusters defined in Fig 2f. Pathways enriched in non-inflammatory KCs are labeled in blue and pathways enriched in inflammatory KCs are indicated in red. Colored circles (nodes) represent pathways, sized by number of genes they contain. Green lines depict intra- and inter-pathway relationships according to the number of genes shared between each pathway. Black circles group related pathways into themes that are labeled. d Flow cytometry data showing the response of monocytes/macrophages in total liver homogenate cell suspensions to stimulation with 1 μg/ml LPS and 25 ng/ml IFN-γ. Cells were stained with anti-human CD45 (clone: HI30), anti-CD68 (clone: Y1/82 A), anti-MARCO (polyclonal; Invitrogen, PA5-26888, goat anti-rabbit secondary antibody), and anti-TNF-α antibodies. Full gating strategy and controls shown in Supplementary Fig. 15. e, f Distribution of MARCO-positive cells in liver zones 1–3. Scale bar represents 500 μm. Staining was performed on 5–7 μM slices cut from formalin-fixed, paraffin-embedded resected liver tissue. Using anti-MARCO (clone: Invitrogen, PA5-26888) and anti-CD68 (clone: PG-M1) at ×40 magnification. f Quantification of percent MARCO-positive cells in liver zones 1 to 3. Error bars show the standard error of the mean for at least seven replicates. Statistical significance evaluated using a one-way analysis of variance (ANOVA) with a Bonferroni post-test ∗∗∗P < 0.001, ∗∗P < 0.01, ∗P < 0.05 Full size image

A key finding of our study is the presence of two distinct populations of human liver macrophages, seeming to segregate into pro-inflammatory and immunoregulatory phenotypes. Our gene list identifies several markers that are unique to one or the other population (for example, MARCO (MAcrophage Receptor with COllagenous structure) is only expressed in non-inflammatory KCs), a finding that can be exploited in tissue-based studies. Using flow cytometry, we observed a subpopulation of macrophages which expressed MARCO on the cell surface (Fig. 8d, Supplementary Fig. 15). By immunohistochemistry staining for MARCO, MARCO+ cells are concentrated in the periportal areas (Fig. 8e, f; Supplementary Figs. 16, 17). We then examined the human macrophage expression profiles that characterize both CD68+ populations in the context of the mouse macrophage ontogeny literature. Human CD68+ MARCO+ cells appear to be transcriptionally similar to long-lived, sessile, liver resident KCs identified in mouse57,58,59, with these cells expressing high levels of VCAM160, CD5L60, HMOX1, MRC161, CD16360, M4SA7, and VSIG458,60. CD68+ MARCO−macrophages have a similar transcriptional profile to inflammatory, recently recruited macrophages including decreased expression of CD16358 and increased expression of PLBD157. To further examine this point using functional assays, we stimulated both populations of macrophages in vitro and examined cytokine secretion via intracellular cytokine staining. We found that MARCO-positive macrophages secreted less TNF-α in response to LPS/ IFN-γ stimulation than MARCO-negative CD68+ macrophages, suggesting that CD68+MARCO- cells are more pro-inflammatory (Fig. 8d, Supplementary Fig. 15). The expression of MARCO in the tumor microenvironment has been linked to poorer outcomes in human breast cancer62. MARCO has also been examined in preclinical mouse colon cancer models with the observation that MARCO expression defined a subtype of suppressive tumor-associated macrophages (TAMs). These TAMs could be polarized to an inflammatory phenotype by anti-MARCO antibody which promoted tumor immunogenicity63. These findings provide a point of reference for examining the role of intrahepatic monocyte/macrophage subsets in the establishment and progression of liver disease.

Liver resident T cells

The human intrahepatic T cell phenotype has been examined via flow cytometry evaluations of enzymatically and mechanically dispersed biopsies taken after reperfusion of donor organs during transplantation64. In this study, the frequency and phenotype of intrahepatic T cells differed from that of peripheral circulating T cells. In our study, we expand on these findings by examining the gene expression patterns that characterize T cell populations in the flushed human liver using cells obtained from a segment of liver tissue obtained prior to reperfusion.

The T cell repertoire within the human liver can be categorized into two broad groups: conventional and unconventional. Conventional T cells consist of CD8+ and CD4+ cells expressing αβ-chain TCRs that recognize antigen in the context of MHC class I and class II molecules respectively. Conventional T cells comprise up to 86% of all CD3+ T cells within the human liver64. Unconventional T cells can express either αβ or γδ TCRs but are not restricted by MHC and do not recognize classical peptide antigens65. In mice, the microbial antigens that enter the liver by the portal vein following digestion have been found to sustain γδ T cell homeostasis and activation66. The relative abundance of subsets of unconventional T cells within the human liver remains to be defined.

In our five healthy NDD livers, we identified three clusters of CD3+ T cells. The most abundant population of T cells (Cluster 2) was characterized by the expression of CD2, CD3D, TRAC, GZMK, CCL5, CCL4L2, PYHIN1, TRBC1, TRBC2, GZMA, CD3E, JUNB, CD69, IL7R, DUSP2, IFNG, LTB, IL32, CD52 (Top DE genes) (Fig. 9a, Supplementary Data 1, 2). These appear to be αβ T cells due to their expression of CD3 and αβ-chain TCRs. Within this population, cells show enriched expression of CD69 and CD8A (Supplementary Data 1, 2). CD69 is a marker for recently activated T cells, and is enriched in tissue resident memory T cells when compared to circulating memory T cells67,68. However, multiple studies have shown that CD69+ CD45RA− tissue resident memory T cells do not express activation markers such as CD25 (IL2RA) or CD38, indicating that the predominant cell in this cluster likely represents tissue resident memory T cells instead of newly activated T cells67. Furthermore, CD69 antagonizes the function and downregulates S1PR1, which is needed for T cell egress from tissues47,48. Since we observed no enrichment of expression of S1PR1, CD38, or IL2RA, all genes characteristic of newly activated non-memory T cells, this further supports our identification of these cells as tissue resident memory T cells. The second cluster of CD3+ cells (Cluster 9), was a population of CD3D expressing cells with enriched expression of TBX21 (aka T-bet), KLRB1 (protein alias CD161), FCGR3A (CD16), NKG7, and GNLY (NKG5) with enriched expression of the TCR delta chain TRDC and the TCR gamma chain TRGC1. The expression profile of these cells suggests that the predominant population in this cluster may be unconventional γδ T cells65 (top DE genes: GNLY, PTGDS, GZMB, S100B, FGFBP2, NKG7, PRF1, KLRF1, HOPX, CST7, KLRD1, CTSW, SPON2, IFITM1, GZMA, CD247, CLIC3, CD7, ADGRG1, CCL5, TRDC). A third cluster of T cells (Cluster 18) was identified as having enriched expression of CD3, little CD4 and CD8 expression with enriched expression of TRDC (Supplementary Data 1&2). The most highly DE gene in this cluster was STMN1, which is expressed on proliferating cells of various lineages69. Cell-cycle analysis revealed that 79% of the cells in this cluster were in G2M (Fig. 2e, Supplementary Fig. 18). MKI67 was also upregulated as a further indication that the cells in this cluster are dividing. Other highly expressed genes in this cluster were GNLY, NKG2A, TYMS, and TOP2A, suggesting that these cells might be phosphoantigen-reactive γδ T cells due to shared enriched genes previously identified in phosphoantigen-reactive γδ T cells isolated from human PBMCs70 (top DE genes STMN1, HMGB2, TYMS, KIAA0101, MKI67, UBE2C, TUBA1B, TRDC, ASPM, CENPA, TOP2A, GNLY, PCNA, AURKB, BIRC5, NUSAP1, TROAP, TUBB, H2AFX, CENPF, CCNB1, H2AFZ).

Fig. 9 The distribution of commonly expressed lymphocyte genes in the healthy liver. a–c i–iv t-SNE plots showing the relative distribution of common αβ T cell, γ∂ T cell, and NK cell markers. d i–iv Common markers of antibody-secreting B cells (plasma cells). e i–iv Common markers of mature B cells. Legend for relative expression of each marker from lowest expression (yellow dots) to highest expression (purple dots) (top left) Full size image

NK-like cells

Recently, three populations of innate-like lymphocytes, ILC1, ILC2, and ILC3 were described71. ILC1 cells include NK cells and other ILCs characterized by expression of T-bet and production of IFN-γ. ILC 2 include cells that express GATA3, BCL11B, and GFI1 and produce IL-4, IL-5, IL-13, and amphiregulin (AREG). ILC3 are thought to include lymphoid tissue inducer (LTi) cells and cells that are positive or negative for the natural cytotoxicity receptors (NCRs) NKp44; they are defined by the expression of RORγt and the production of IL-17A, IL-17F, and IL-22. In studies of human liver cells obtained by flow-based cell sorting of enzymatically and mechanically dissociated hepatic resection tissue, perforin, and granzyme staining revealed that human liver NK cells rapidly respond to antigens and malignant cells by quickly releasing lytic granules36. By flow cytometry, CD16− NK cells make up to 50% of NK cells in the human liver36. We found a population of NK-like cells (Cluster 8) characterized by enriched expression of CD7, C–type lectin receptor KLRD1 (CD94), GZMK (Granzyme K), NCR1(NKp46), and NCAM1 (CD56), without upregulated expression of FCGR3A (CD16) or ITGA1 (CD49a) (Supplementary Data 1&2). This population also showed enriched expression of EOMES, suggesting that the predominant population in this cluster may correlate with the long-lived NK cell population previously described in human liver biopsies by flow cytometry72. The top DE genes of Cluster 8 are: CD7, CMC1, XCL2, KLRB1, XCL1, KLRC1, KLRF1, IL2RB, CD160, CCL3, KLRD1, NKG7, TXK, ALOX5AP, TRDC, CD69, TMIGD2, CLIC3, GZMK, DUSP2, MATK, IFITM1, CCL4, CD247. Along with conventional T cells, γδ T cells, and NK-like cells, it is known that NKT cells, MAIT cells and atypical T cells populations are found in the liver73,74. In our initial analysis, we identified clusters that corresponded to T and NK-like cells, but our data does not currently have the resolution to characterize the less-frequent immune cell populations from the TLH analysis and we suggest that these clusters are likely multiple cell populations. MAIT cells in particular are a subset of T cells with an αβ TCR characterized by a semi-invariant TCR alpha (TCRα) chain and CD161 expression. These cells likely fall within Cluster 2, for example CD161 (KLRB1) is expressed on 52% of cells in Cluster 2 (Supplementary Data 1, 2). In the future, we will employ TCR clonotyping to examine the transcriptional signature of CD161+ αβ TCR+ cells expressing the known semi-invariant TCRs characteristic of MAIT cells (TRAV1-2/TRAJ12/20/3). As well, in cluster 8 (NK-like), 23% of the cells express NKp46 (NCR1) and 18% express CD56 (NCAM1)(Supplementary Data 1&2). This suggests that additional cell populations are clustering with CD56+ NKp46+ cells due to similarities in gene expression. To address this issue, we focused on the T and NK cell populations and carried out a sub-analysis comparing clusters 2, 8, 9, and 18 to one another. This analysis yielded nine populations of conventional and unconventional T cells and NK-like cells that comprise cells with inflammatory or immunoregulatory properties (Supplementary Data 2, 6 and Supplementary Fig. 19). These results provide a description of baseline T and NK cell transcriptomes in the human liver. Further characterization of these cell populations will require a more cells and additional single-cell analysis tools, including surface protein labeling (CITE-seq)39 and simultaneous characterization of the T cell receptor clonotype using single-cell V(d)J sequencing [https://www.10xgenomics.com/solutions/vdj/]. These approaches will improve our ability to detect and identify these important cell populations in future analyses.

B cells

The frequency and role of hepatic B cells remains a topic of debate48. Flow cytometric analyses of human liver tissue indicated that CD19+ B cells account for up to 8% of all lymphocytes64, while B cells have been found to make up as much as 37.5% of the mouse intrahepatic immune cell population75. Information on the phenotype and function of these B cells is limited due to their small number and the difficulty in isolating and analyzing hepatic B cells48. Recently, IgA-producing plasma cells of GALT (gut-associated lymphoid tissue) origin were found to be enriched in human liver periportal areas by histology76. Using scRNA-seq we were able to identify distinct populations of B cells with different stages of development. We found two populations of liver resident B cells which appeared to be plasma cells and antigen inexperienced B cells (Figs. 2f, 3a, b, 9d i–iv, e i–iv). Cluster 16 identified a subset of cells that were enriched for the expression of IGHD, CD19, MS4A1 (protein alias CD20), CD22, CD52 without expression of either CD27 or CD138 which suggest that this cluster is comprised of mature antigen inexperienced B cells77,78 (top DE genes: MS4A1, LTB, CD37, CD79B, CD52, HLA-DQB1, TNFRSF13C, TCL1A, LINC00926, STAG3, IGHD, BANK1, IRF8, BIRC3, P2RX5, RP11-693J15.5, RP5-887A10.1, VPREB3, CD22, CD74, SELL).

Plasma cells are terminally differentiated B cells that reside in tissue and continuously synthesize and secrete antibodies directed against specific antigens79. The cells of Cluster 7 were enriched for the expression of immunoglobulin heavy and light chains, CD27, CD3878, and had reduced expression of MS4A1(CD20)80, suggesting that this cluster is comprised of plasma cells (top DE genes: IGLC2, IGHG1, IGKC, IGHG2, IGHG3, IGHGP, IGLC3, JCHAIN, IGHA1, IGHG4, IGHA2, IGHM, IGLV3-1, IGLC7, MZB1, CD79A, SSR4, IL16). Our baseline description of B cell transcriptional profiles in the human liver will be useful in any examination of liver diseases modulated by B cells, such as alcohol-induced hepatitis76.