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 experiments13. 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 OCT14. 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 using the Visium platform, 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, differential gene expression (DE) analysis was performed to compare SNAP and OCT preservation techniques in normal and BA with advanced fibrosis samples from the same patient. 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, Supplementary Fig. S1A). 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, Supplementary Fig. S1B). 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)15, 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 single-cell spatial mapping 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). Superimposing the spatial spots on the actual H&E stained tissue section confirms that a single Visium spot can represent multiple nuclei (Fig. 2C). When attempting to identify cell types based solely on Visium spatial transcriptomics data (Fig. 2D), we see multiple cell types present in each cluster, as shown in the heatmap in Fig. 2E, affirming that clustering of spatial spots cannot resolve individual cell populations. Recent advances in informatics methods have been developed to map single cell transcriptomes from scRNAseq or snRNAseq outputs onto spatial transcriptomic datasets. There are two general approaches: identification of the percent of each cluster or cell type in a spatial spot (CARD16 or Seurat10), also known as deconvolution, versus mapping a single cell type to a specific spot (Celltrek17). We implemented both techniques, CARD and Celltrek (Supplementary Fig. S2). 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 single-cell spatial mapping, 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 established cell lineage markers were used to annotate individual cell types (Fig. 3A,B)18. 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 endothelial progenitor cells, LSEC zones 2&3, and hepatocyte progenitor cells. These annotated snRNAseq datasets were then used to spatially map single cells onto the spatial OCT samples (Fig. 3C). Using these spatially resolved single-cell data, colocalization analysis was completed to identify spatial distribution of different cell types in the tissue (Fig. 3E). This revealed a loss of overall complexity in the cellular interactions in BA with advanced fibrosis. Looking closer at individual colocalization interactions, BA with advanced fibrosis is associated with an increase in the colocalization strength between different types of endothelial cells and hepatocytes as well as between HSC/fibroblasts, cholangiocytes, and immune cells. In normal liver tissue, immune populations and HSC interact more with LSEC and hepatocyte populations.
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 (Fig. 4) and BA with advanced fibrosis tissues (Fig. 5). To increase our confidence in the single-cell spatial mapping approach, we observed that the proportion of cells identified in each cell type cluster in the single-cell spatially mapped data mirrored the snRNAseq 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 single-cell spatially mapped 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 cluster fall in the periportal regions (Fig. 4D). In BA with advanced fibrosis, cells identified in the interzonal 2 cluster were the largest proportion of cells (Fig. 5A), mirroring the distribution of cells in the single-cell spatially mapped 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).
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 single-cell spatially mapped 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). No significantly upregulated L–R interactions were found in BA with advanced fibrosis when compared to normal liver when analyzing the snRNAseq samples. 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 may not be sufficient to understand potential L–R interactions in a liver tissue sample.
Comparison with single-cell spatially mapped data 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, single-cell spatial mapping 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 processes19,20,21. 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 fibrosis22. As FGF23 is known to be an antagonist of myofibroblast activity, a decrease in FGF23 signaling, as observed in BA tissue, leads to increased myofibroblast activity and therefore fibrosis progression23.
When comparing L–R interactions downregulated in both datasets, we see a low correlation in communication strength (Spearman correlation of 0.57). Looking at the relative change in communication strength in these interacting cell groups further corroborates that there is a difference in the L–R findings of snRNAseq and single-cell spatially mapped data (Supplementary Fig. S3A). This highlights that the snRNAseq data alone are not sufficient to understand potential L–R interactions in a liver tissue sample, especially given that these findings cannot be validated in snRNAseq data. Predicted L–R interactions from snRNAseq are enhanced by comparing single-cell spatially mapped transcriptomic data; this not only adds information on communication pathways, but the findings can be validated by projecting L–R pairs onto the actual tissue.
While SNAP spatial samples were not used in downstream analysis, we did compared L–R interaction results between the two different techniques for each patient (Supplementary Fig. S3B). Overall, 25% of the L–R interactions were downregulated in only the OCT dataset and 35% in only the SNAP dataset, while 40% were downregulated in both (Supplementary Fig. S3C). The difference in downregulated interactions between the two samples is not surprising given that we already saw a difference in transcriptional changes when individually comparing the two freezing techniques in normal and BA with advanced fibrosis and that the tissue blocks are different, resulting in differences in the number of portal triads, central veins, etc. (Fig. 1). Despite using different tissue blocks and the impact of freezing techniques on some gene expression pathways, the communication strength of the overlapping L–R signaling pairs between samples from the same patient was high when compared to the snRNAseq vs spatial comparison (Spearman correlation coefficient of 0.98 versus 0.57). This increase in the correlation of spatial interactions further supports comparing and confirming L–R analysis findings with spatial transcriptomics where possible.