Part 1: Impact of different freezing techniques on spatial transcriptomic analysis of human liver samples.
The constraints of clinical workflows, including the timing of biopsies and surgical procedures, can impact the ability to perform experiments using fresh liver samples. Also, many biorepositories bank samples of different disease states for future studies. In both instances, liver tissue is often snap frozen and stored at -80 C or in liquid nitrogen for future experiments.18 The Visium platform enables spatial transcriptomic analysis of frozen tissue specimens, and the manufacturer recommends samples be frozen using OCT at the time of collection. Based on our prior work with scRNAseq of human liver, we hypothesized that liver tissue that was snap frozen and stored would generate different transcriptomic profiles than liver tissue that was initially frozen using OCT.19 Fresh samples of normal liver tissue and BA with advanced fibrosis were collected and subject to two differing tissue preparation techniques at the time of collection: (i) snap freezing with storage in liquid nitrogen, followed by embedding in OCT at the time of retrieval from the freezer for transcriptomic analysis (SNAP), versus (ii) embedding in OCT up front with storage in liquid nitrogen (OCT) (Fig. 1a). A spatial map of gene expression for each 55 µm spot was generated for each sample, and clustering algorithms were utilized to identify groups of spots with similar gene expression patterns (Fig. 1b). The analysis resulted in a diverse set of 36 clusters of spots representative of the spatial complexity in the liver in both SNAP and OCT conditions. Next, DE analysis was performed to compare SNAP and OCT preservation techniques in normal and BA with advanced fibrosis samples. When comparing the two freezing techniques in normal liver, DE genes were largely upregulated, including activating transcription factor 3 (ATF3), DnaJ heat shock protein family member B1 (DNAJB1), and phosphoenolpyruvate carboxykinase 1 (PCK1) (fold change > 2.5) in SNAP relative to OCT (Fig. 1c). When comparing BA with advanced fibrosis between freezing conditions (SNAP vs. OCT), the subset of genes that were DE were primarily downregulated, including hepcidin antimicrobial peptide (HAMP) and metallothionein 1G (MT1G) (fold change < -2.2) in SNAP compared to OCT (Fig. 1c). We next performed pathway analysis of DE genes, highlighting effects on the proteasome, oxidative phosphorylation, and reactive oxygen species in SNAP relative to OCT in both normal and BA samples, indicating that preservation technique likely impacts downstream transcript quality in these pathways (Fig. 1d). In BA with advanced fibrosis, genes associated with oxidative phosphorylation, metabolic associated fatty liver disease (MAFLD)20, and thermogenesis had even larger fold changes in SNAP relative to OCT than normal livers, indicating that diseased tissue with advanced fibrosis may be more sensitive to tissue processing techniques (Fig. 1d). Given these findings, only the OCT samples were used for further analysis.
Part 2: Integration of snRNAseq and spatial transcriptomic data via deconvolution enables the creation of a spatially resolved single-cell atlas of the human liver.
Liver tissue is highly spatially organized with portal triads and central veins connected by sinusoids (Fig. 2a). The Visium spatial transcriptomic platform generates output ‘spots’ of 55µm diameter, which can represent 1–10 cell types per spot (Fig. 2b). When attempting to identify cell types based solely on Visium spots, multiple cell types are present in each spot, even when intentionally overclustering (Fig. 2c), as shown in the heatmap in Fig. 2d. Recent advances in informatics methods have been developed to map single cell transcriptomes from scRNAseq or snRNAseq outputs onto spatial transcriptomic datasets using a computational approach called ‘deconvolution’. There are two general approaches to deconvolution, which include identification of the percent of each cluster or cell type in a spatial spot (CARD21 or Seurat10) versus mapping a single cell type to a specific spot (Celltrek15). We implemented both deconvolution techniques, CARD and Celltrek (Supplemental Fig. 1). Given the complexity of liver tissue and the desire to understand spatial relationships between individual cell types and their transcriptomes, rather than just their location, we employed the second approach to map single cell types to Visium spots using Celltrek for downstream analysis since it improved output interpretability. To perform deconvolution, we first generated matched snRNAseq datasets from both normal and BA with advanced fibrosis samples. Unsupervised clustering was used to generate 24 clusters within the snRNASeq data, and cell lineage marker expression was used to annotate individual cell types (Fig. 3a,b). Relative shifts in cell proportions were evaluated between normal and BA with advanced fibrosis samples and revealed that BA tissue had increased cholangiocytes, liver sinusoidal endothelial cells (LSEC) zone 1, T/NK cells, B/Plasma cells, hepatic stellate cells (HSC)/fibroblasts, and portal endothelial cells (Fig. 3d). BA with advanced fibrosis had also decreased epithelial progenitor cells, LSEC zones 2&3, and hepatocyte progenitor cells. These annotated snRNAseq datasets were then used to deconvolute the spatial OCT samples (Fig. 3c). Each cell type was projected onto the spatial transcriptomic output and visualized relative to H&E staining to demonstrate density and spatial distribution in normal liver tissue (Fig. 4) and BA with advanced fibrosis tissue (Fig. 5). To increase our confidence in the deconvolution approach, we observed that the proportion of cells identified in each cell type cluster in deconvoluted data mirrored the snRNA seq data in both normal and BA liver tissues. In normal liver, the highest proportion of cells belonged to the LSEC zones 2&3 subcluster 1 cluster and the central-venous hepatocyte clusters (Fig. 4a). This is also seen in the visualization of the deconvoluted data on the tissue sample (Fig. 4c,d). Additionally, the annotated cells can be compared to pathologist-labeled periportal and central venous regions (Fig. 4b). While not spatially exclusive, cells identified in the portal endothelial cell gene expression cluster and the LSEC zone 1 fall in the periportal regions (Fig. 4d). In BA with advanced fibrosis, cells identified in the intrazonal 1 cluster were the largest proportion of cells (Fig. 5a), mirroring the distribution of cells in the deconvoluted data (Fig. 5c,d). When compared to microanatomy annotations (Fig. 5b), the immune cell types in BA were usually in periportal areas (Fig. 5d) whereas their spatial distribution was more diffuse in normal liver (Fig. 4d). Using these deconvoluted, spatially resolved single-cell data, colocalization analysis was completed to identify significant cellular interactions and avoidances (Fig. 3e). This revealed a loss of complexity in the cellular interactions in BA with advanced fibrosis, with immune populations interacting more with LSEC and hepatocyte populations in normal liver tissue.
Part 3: Spatial transcriptomic data augments Ligand-Receptor Analysis in human livers.
Next, L-R analysis was performed, first on the snRNAseq data and then on the spatial transcriptomic data. Overall, the L-R pairs among interacting cell types showed a decrease in signaling in BA with advanced fibrosis when compared to normal liver, in either the snRNAseq data, spatial transcriptomic data, or both datasets (Fig. 6a). Overall, 51% of the L-R interactions were downregulated in only the spatial transcriptomic dataset, whereas 14.5% were downregulated only in the snRNAseq dataset with the remaining 34.5% downregulated in both (Fig. 6b). This highlights that the snRNAseq data alone are not sufficient to understand potential L-R interactions in a liver tissue sample.
Comparison with spatial deconvolution adds additional information regarding communication pathways. For example, the Laminin Subunit Gamma-3 (LAMC3)-CD44 interaction present between HSC/fibroblasts and Kupffer cells is downregulated in BA with advanced fibrosis in both the snRNAseq and spatial transcriptomic data (Fig. 6c). In this instance, spatial deconvolution analysis adds confidence to the results from the snRNAseq L-R analysis that BA with advanced fibrosis samples show potential interference with the communication pathway between these two cell types; specifically between LAMC3, a major component of basement membranes and CD44, a cell surface receptor that has been implicated in fibrotic and wound healing processes.22–24 In contrast, the Fibroblast Growth Factor-23 (FGF23) – Fibroblast Growth Factor Receptor-1 (FGFR1) interaction from LSEC zones 2&3 subcluster 1 to subcluster 2 is only downregulated in BA spatial transcriptomic data (Fig. 6d). For this interaction, the spatial transcriptomic data added information not previously seen when analyzing snRNAseq data. This suggests that within LSEC zones 2&3 of BA tissue, there is a distinctively spatial dysregulation of the FGF23-FGFR1 pathway, which has been implicated in both liver metabolism and the development of fibrosis.25
Ultimately, when comparing L-R interactions downregulated using both datasets, the relative change in communication strength in interacting cell groups between snRNAseq and spatially deconvoluted data was low, with a Spearman correlation of 0.57, suggesting that predicted L-R interactions from snRNAseq should be verified using deconvoluted spatial transcriptomic data (Fig. 6e).