Paired-cell sequencing enables spatial gene expression mapping of liver endothelial cells

Abstract

Spatially resolved single-cell RNA sequencing (scRNAseq) is a powerful approach for inferring connections between a cell's identity and its position in a tissue. We recently combined scRNAseq with spatially mapped landmark genes to infer the expression zonation of hepatocytes. However, determining zonation of small cells with low mRNA content, or without highly expressed landmark genes, remains challenging. Here we used paired-cell sequencing, in which mRNA from pairs of attached mouse cells were sequenced and gene expression from one cell type was used to infer the pairs' tissue coordinates. We applied this method to pairs of hepatocytes and liver endothelial cells (LECs). Using the spatial information from hepatocytes, we reconstructed LEC zonation and extracted a landmark gene panel that we used to spatially map LEC scRNAseq data. Our approach revealed the expression of both Wnt ligands and the Dkk3 Wnt antagonist in distinct pericentral LEC sub-populations. This approach can be used to reconstruct spatial expression maps of non-parenchymal cells in other tissues.

Figure 1. Strategy for paired-cell reconstruction of liver LEC zonation.

a) A diagram of the liver lobule. Blow up represents a typical porto-central sinusoidal unit. A typical unit consists of 10-15 hepatocytes and was coarse-grained into 8 or 4 concentric layers when analyzing paired cells and single cells respectively. b) Paired-cell RNA sequencing utilizes the hepatocyte zonation to determine tissue localization and single cell RNAseq of LECs to extract LEC-specific genes (1). LEC zonation is obtained by averaging expression of LEC genes in the spatially-localized pairs. This dataset is also used to extract LEC landmark genes for localizing single LECs (2).

Figure 2. Single cell RNAseq reveals the expression signatures of liver non-parenchymal cells.

a) tSNE map colored by the identified seven clusters, consisting of endothelial cells (b, 1,203 cells), T cells (c, 958 cells), plasmacytoid dendritic cells (pDCs, d, 119 cells), Kupffer cells (e, 340 cells), liver capsule macrophages (LCMs, f, 164 cells), B cells (g, 307 cells) and Neutrophils (h, 60 cells). Blue is high expression, gray is low expression in b)-h).

Figure 3. Sorting strategy to isolate pairs of attached hepatocytes and LECs.

a-c) FACS gates used in the isolation. a) FSC-A and SSC-A are used to select hepatocytes based on size. b) FSC-W is used to gate out clusters of hepatocytes. c) CD31 and CD45 fluorescense is used to select pairs that contain a LEC. Numbers in a-c) represent the percent of gated cells. d) Imagestream example of a sorted pair. BF – bright field. Hoechst staining in blue, CD31 in magenta. The experiment was repeated 3 independed times with similar results. e) Expression of hepatocyte pericentral genes (Cyp2e1, Oat, Gulo, Lect2, Cyp1a2) is enriched in the Rspo3+\Wnt9b+ pairs (p=3.3e-29, two-sided Wilcoxon rank-sum test) and depleted in the Efnb2+\Dll4+ pairs (p=6.45e-13). Expression of hepatocyte periportal genes (Cyp2f2, Sds, Hsd17b13, Aldh1b1, Alb, Cps1, Arg1, Pck1) is enriched in the Efnb2+\Dll4+ pairs (p=9.35e-9) and depleted in the Rspo3+\Wnt9b+ pairs (p=1.5e-29, n=4602 paired-cells). Expression is in units of fraction of total cellular UMIs. Box plot elements: center line, median; box limits, first to third quartile (Q1 to Q3); whiskers, extend to the most extreme data point within 1.5× the interquartile range (IQR) from the box; circles, data points. f) Reconstructed zonation profiles of hepatocyte genes based on the pcRNAseq data overlap profiles validated with smFISH . Patches are standard errors of the means. smFISH plots were based on n=10 lobules from 2 mice, pcRNAseq plots were based on n=4602 paired-cells.

Figure 4. Global zonation of LEC genes.

a) Reconstructed zonation profiles of the 475 significantly zonated LEC genes. Genes are sorted by the centers of mass of their zonation profiles and normalized by their maximum. Plots on the left show examples of pericentral (Wnt9b, Rspo3), non-monotonic (Pcdhgc5, Ecm1) and periportal (Ltbp4, Efnb2) profiles. Patches are standard errors of the means, plots were based on n=4602 paired-cells. b) smFISH image of the zonated gene Dll4. LEC marker Aqp1 in magenta, the periportal hepatocyte marker Acly in green and Dll4 in red. CV – central vein, PN – portal node. Blue – DAPI stained nuclei, gray – phalloidin stained membranes. Insets show a pericentral LEC (magenta outline), a mid-lobule LEC (yellow outline) and a periportal LEC (green outline). Red dots, transcripts of the Dll4, green dots, transcripts of Acly. Large orange blobs outside of LECs are lipofuscin – non-specific fluorescence common in liver tissue. Red outlines and white arrowheads mark LECs, green outlines mark hepatocytes. Scale bar is 20μm in the scan and 5μm in the insets. c) Validations of the zonation profiles using smFISH for 12 zonated LEC genes. Blue curves are the inferred zonation profile based on pcRNAseq, patches are s.e.m. Red curves are the average zonation profiles based on the smFISH of the sinusoidal LECs, patches are s.e.m. KruskalWallis p-values confirming the significance between the pericentral and periportal endothelial cells are shown at each plot title. Box plot elements: center line, median; box limits, first to third quartile (Q1 to Q3) of the smFISH expression; Gray patches in each plot mark the pericentral layers (CV, left) and periportal layers (PN, right). As opposed to the pcRNAseq, within these layers smFISH enabled distinguishing between LECs that line the central and portal vessels (the extreme boxes) and the sinusoidal LECs. N=30 cells for each layer from 2 mice were quantified for every gene.

Figure 5. Expression signature of pericentral LECs.

a) Zonation profiles of representative pericentral genes, selected out of the 60 zonated LEC genes with the highest ratio of expression between pericentral layer 1 and pericentral layer 8. Patches are s.e.m. plots were based on n=4602 paired-cells. b) Rspo3 (green dots) is highly expressed both in LECs that line the central vein (CV, arrows), and in sinusoidal pericentral LECs (arrowheads). c) Kit (green dots) is repressed in LECs that line the central vein (CV, arrows), and up-regulated in sinusoidal pericentral LECs (arrowheads). Red dots are mRNAs of the LEC marker Aqp1, yellow blobs are hepatocyte lipofuscins that fluoresce in both red and green channels. Scale bar is 10μm. d) Dkk3 (green dots) is expressed in a subset of LECs that line the central vein (green arrow), distinct from cells that express Wnt2 (red dots, cell marked by red arrow). Scale bar is 5μm. In (b-d) micrographs are representative of 10 lobules and two mice exhibiting similar results. e) Expression of Dkk3 and Wnt2 is anticorrelated in LECs that line the central vein, Spearman R=-0.37, p=3.9e-4 (n=92 cells).

Figure 6. Spatial reconstruction of single LECs using landmark genes obtained from pcRNAseq.

a-b) pcRNAseq-based zonation profiles of the panel of LEC landmark genes used to localize cells in the scRNAseq data. a) 70 pericentral LEC genes. b) 70 periportal LEC genes. Profiles were scaled between 0 and 1. c) Overlap of zonation profiles based on scRNAseq (blue) n=1163 single-cells and pcRNAseq (red) n=4602 paired-cells. Patches are s.e.m. Zonation profiles were normalized to their means across all layers (pcRNAseq profiles contain 8 layers whereas scRNAseq profiles contain 4). d) Ctsl (green dots) is highly expressed in hepatocytes (white dashed outline) and also expressed in the adjacent LECs (arrows). micrographs are representative of 10 lobules and two mice exhibiting similar results. e) scRNAseq spatial reconstruction reveals zonated expression of Ctsl in LECs, with a reduced expression level in pericentral LECs (1) compared to mid-lobule LECs (2). Box plot elements: center line, median; box limits, first to third quartile (Q1 to Q3) of the smFISH expression, horizonatal lines are medians. Quantification based on n= 30 cells from each layer from 2 mice. Gray patches mark the pericentral layers (CV, left) and periportal layers (PN, right). Blue is the scRNAseq-based zonation profile. Patches are s.e.m. Panels on the left show representative images. Blue is DAPI nuclear stain. Scale bar is 5μm in d)-e).