Single-cell transcriptomics reveals conserved cell identities and fibrogenic phenotypes in zebrafish and human liver.

The mechanisms underlying liver fibrosis are multifaceted and remain elusive with no approved antifibrotic treatments available. The adult zebrafish has been an underutilized tool to study liver fibrosis. We aimed to characterize the single-cell transcriptome of the adult zebrafish liver to determine its utility as a model for studying liver fibrosis. We used single-cell RNA sequencing (scRNA-seq) of adult zebrafish liver to study the molecular and cellular dynamics at a single-cell level. We performed a comparative analysis to scRNA-seq of human liver with a focus on hepatic stellate cells (HSCs), the driver cells in liver fibrosis. scRNA-seq reveals transcriptionally unique populations of hepatic cell types that comprise the zebrafish liver. Joint clustering with human liver scRNA-seq data demonstrates high conservation of transcriptional profiles and human marker genes in zebrafish. Human and zebrafish HSCs show conservation of transcriptional profiles, and we uncover collectin subfamily member 11 (colec11) as a novel, conserved marker for zebrafish HSCs. To demonstrate the power of scRNA-seq to study liver fibrosis using zebrafish, we performed scRNA-seq on our zebrafish model of a pediatric liver disease with mutation in mannose phosphate isomerase (MPI) and characteristic early liver fibrosis. We found fibrosis signaling pathways and upstream regulators conserved across MPI-depleted zebrafish and human HSCs. CellPhoneDB analysis of zebrafish transcriptome identified neuropilin 1 as a potential driver of liver fibrosis.
Conclusion: This study establishes the first scRNA-seq atlas of the adult zebrafish liver, highlights the high degree of similarity to human liver, and strengthens its value as a model to study liver fibrosis.

INTRODUCTION

Liver fibrosis is the excessive accumulation of extracellular matrix that is produced in response to chronic liver injury, which can progress to cirrhosis and liver failure if untreated.[ ] Despite the increasing prevalence and high morbidity associated with liver fibrosis, there are currently no approved antifibrotic treatments. Using animal models to study liver function and disease is essential to better understand liver cell biology, identify new therapeutic targets, and establish preclinical models to test potential therapies.

The zebrafish liver is remarkably similar to the human liver, with conserved cellular composition and functionality. The zebrafish liver is fully functional by 5 days post‐fertilization, making it a high‐throughput tool for disease modeling.[ ] The majority of studies have focused on larval zebrafish to model acute liver injury,[ , , ] whereas the adult zebrafish has been underutilized to study chronic injury, particularly in fibrotic liver diseases.[ ] We previously generated a zebrafish mutant line, mannose phosphate isomerase (mpi) with knockdown of Mpi, to recapitulate the characteristic early onset liver fibrosis seen in children with a mutation in MPI, resulting in a congenital disorder of glycosylation (MPI‐CDG).[ ] We used this model to discover that MPI depletion could directly activate hepatic stellate cells (HSCs), the main driver cell of liver fibrosis, but the effects of MPI loss in other liver cell types in the regulation of liver fibrosis is not known.[ ]

Single‐cell RNA sequencing (scRNA‐seq) has emerged as a powerful tool to identify key marker genes and pathways in various cell types and uncover unique cell populations in tissues. Recent high‐impact studies on human livers have used scRNA‐seq to reveal the existence and functions of unique subpopulations of various hepatic cell types, including macrophages, epithelial progenitor cells, and myofibroblasts.[ , , ] In contrast, there has been no reported scRNA‐seq data on the adult zebrafish liver in physiologic or fibrotic conditions, which limits our use of this model organism to study molecular and cellular dynamics in the mature liver at a single‐cell level.

Here, we present the first full characterization of the adult zebrafish liver transcriptome at the single‐cell level and demonstrate its use through comparison with the mpi zebrafish model of fibrosis to uncover potential mechanisms of fibrogenesis. We compare and contrast the transcriptional cell identities with adult human liver, specifically zebrafish HSC signatures to that of human HSCs in physiologic and fibrotic conditions. Lastly, we leverage our scRNA‐seq data to identify neuropilin 1 (Nrp1) as a potential mediator of altered vascular endothelial growth factor (Vegf) signaling in cholangiocytes to drive liver fibrosis. This adult zebrafish liver cell atlas can be a valuable resource to expand the use of adult zebrafish as a tool to study liver fibrosis.

MATERIALS AND METHODS

Zebrafish maintenance

Procedures were performed in accordance with the Mount Sinai Institutional Animal Care and Use Committee guidelines. Adult fish were maintained on a 14 hour:10 hour light/dark cycle at 28°C. AB wild type (WT), mpi , [ ] mpi, and Tg(hand2:EGFP) zebrafish strains were used.

scRNA‐seq on adult zebrafish livers

Single‐cell suspension

Whole livers were dissected from 18‐month‐old adult male zebrafish. Single‐cell suspensions were generated according to our published protocol.[ ]

Alignment, transcriptome assembly, and quality control

Alignment was performed according to a published protocol.[ ] The raw (unfiltered) output matrices were used for the clustering and downstream analysis in R package Seurat, version 4.0.[ ] Data were filtered to include cells with unique molecular identifiers (UMIs) >200, 200

Differential gene expression on scRNA‐seq data

To determine differential gene expression (DGE) between clusters and samples, the FindMarkers function in Seurat version 4.0 was used.[ ] The default Wilcoxon rank‐sum test was used to determine significance.

Human and zebrafish liver integrated clustering

A list of human‐to‐fish orthologs was generated using the ZFIN zebrafish–human orthologs database, which contains human orthologs for roughly 80% of annotated zebrafish genes (~20,000 unique genes; (August 20, 2020)).[ ] The orthologs were curated considering the following three factors: conserved genome location, amino acid sequence comparison, and the phylogenetics tree. Zebrafish gene names were converted to the orthologous human gene name. Joint clustering of filtered human liver scRNA‐seq data from MacParland et al.[ ] and zebrafish data was performed by normalizing and identifying variable features for each species dataset, and then integration of datasets using FindIntegrationAnchors. Data were scaled and dimensionality reduction was performed. The first 15 principal components were used to perform clustering at a resolution of 0.8. For joint clustering, HSCs and endothelial cells (EC) were identified from the human and zebrafish data sets, respectively, and subset. The same workflow was applied, using the first seven principal components and a resolution of 0.8, to generate clustering.

Identifying ligand–receptor interactions

CellPhoneDB, version 2.0[ ] was run using its statistical analysis method with normalized counts and cell annotation input files from each zebrafish scRNA‐seq data set. Dot plots of expression scores and p values for specific ligand–receptor interactions between cell type pairs were created using the R script for the dot_plot function obtained from https://github.com/Teichlab/cellphonedb.

Identifying modules of coregulated genes in EC populations

The graph_test function in Monocle 3[ , , ] was used to identify genes that are variable in expression across zebrafish WT ECs and HSCs. We used find_gene_modules to perform uniform manifold approximation and projection (UMAP) for dimension reduction analysis on significant genes (q ≤ 0.05) and group genes into modules by using Louvain community analysis. The aggregate_gene_expression function was used to aggregate the expression of all genes in each module to visualize the expression of each module in each cluster.

Cell culture

LX‐2[ ] and TWNT‐4[ ] human HSC lines were maintained in Dulbecco’s modified Eagle’s medium with 10% fetal bovine serum (FBS), L‐glutamine, and penicillin/streptomycin and routinely tested for mycoplasma using the Venor GeM Mycoplasma Detection Kit (Sigma‐Aldrich). Small interfering RNAs (siRNAs) targeting MPI were transfected using Lipofectamine RNAiMAX transfection reagent (ThermoFisher, Waltham, MA), as described.[ ] Cells were collected 48 hours after transfection for RNA or protein.

Immunofluorescence staining

Adult zebrafish livers were dissected and fixed in 4% paraformaldehyde in phosphate‐buffered saline (PBS) overnight at 4°C. Livers were incubated in 10%, 20%, and 30% sucrose in PBS for 24 hours each at 4°C. Samples were processed using a published protocol.[ ] LX‐2 cells were seeded at 200,000 cells/mL on chamber slides (Nunc Lab‐Tek; ThermoFisher), grown overnight, fixed in 4% paraformaldehyde for 15 minutes, permeabilized with PBS + 0.4% Triton X‐100, and blocked with PBS + 0.25% Triton X‐100 + 5% FBS + 2% bovine serum albumin for 1.5 hours at room temperature. Samples were stained with primary antibody (Table S1) overnight at 4°C and followed by secondary antibody for 1.5 hours in the dark at room temperature. Samples were mounted with ProLong Diamond Antifade Mountant with 4´,6‐diamidino‐2‐phenylindole (DAPI; Life Technologies) and imaged with a Leica SP5 DMi.

MPI enzyme activity assay

MPI activity assay was performed on cell protein extract according to our published protocol.[ , ]

Ingenuity pathway analysis on DGE

DGE was calculated for mpi versus WT HSCs as described in the methods above. Differentially expressed genes, log‐fold changes, and p values were imported into ingenuity pathway analysis (IPA) for all data sets (https://digitalinsights.qiagen.com/IPA)[ ] Core analysis was conducted on each data set using log‐fold change values. Figures were generated using Prism 9.1.1 (GraphPad Software, San Diego, CA).

Statistics

We used Prism 9.1.1 for data analysis. A cutoff of p < 0.05 was considered significant. Data represent mean ± SD.

RESULTS

Creation of a single‐cell atlas of adult zebrafish liver

To develop a single‐cell atlas of adult zebrafish liver, we performed scRNA‐seq on livers dissected from WT zebrafish (n = 3), which yielded 13,630 cells for further analysis (Table S2). Samples were integrated to perform clustering using the Seurat pipeline,[ ] and 19 transcriptionally distinct clusters were identified (Figure 1A). Each cluster contained cells from all three samples, supporting that these clusters are representative of biological populations (Figure 1B; Table S3).

FIGURE 1

Single‐cell transcriptome of the adult zebrafish liver. (A) scRNA‐seq of WT adult zebrafish liver cells (n = 3). UMAP shows 19 unique clusters. (B) Bar graph showing the percentage of cells contributed to each cluster from each sample. (C) Heatmap of gene expression for three of the top differentially expressed genes in each cluster. (D) UMAP visualization of marker gene expression for hepatocytes (tfa and cp), ECs (kdrl), HSCs (hand2), lymphocytes (il7r), and erythrocytes (hbaa1)

Hepatocytes were the most abundant cell type recovered from zebrafish livers and were identified by high expression of hepatocyte marker fatty acid binding protein 10a, liver basic (fabp10a)[ ] (Figure 1C). Other hepatocyte marker genes, transferrin‐a (tfa)[ ] and ceruloplasmin (cp),[ ] further discriminated hepatocytes from other cell populations (Figure 1D). Hepatocytes comprised nine unique clusters of the 19 total identified (zebrafish [z] hepatocyte [Hep]1‐9) (Figure 1A,C). Gene ontology (GO) enrichment analysis revealed both shared and distinct metabolic pathways among hepatocyte clusters. Steroid metabolic processes and cholesterol transport were common across clusters (Figure S1A). Interestingly, clusters were distinctly enriched for pathways involved in oxidative phosphorylation (zHep1, zHep7, and zHep8); cholesterol, steroid, and lipid metabolism and transport (zHep2, zHep3, zHep4, and zHep5); fatty acid transport and glucose metabolism (zHep6); and immune processes (zHep9) (Figure S1A). These data demonstrate that the adult zebrafish liver is comprised of unique hepatocyte populations with distinct functions, similar to the specialization of human hepatocytes.[ , ]

Cholangiocytes (Chol) made up three of the 19 total clusters (zChol1, zChol2, and anterior gradient 2‐positive (zAgr2+) (Figure 1A). Cholangiocytes were identified by expression of annexin A4 (anxa4)[ ] (Figure 1C). Other top DEGs in cholangiocyte clusters included lectin, galactoside‐binding, soluble, 2b (lgals2b), transmembrane 4 L six family member 4 (tm4sf4), GRAM domain containing 1Bb (gramd1bb), growth arrest and DNA‐damage‐inducible, beta b (gadd45bb), and Kruppel‐like factor 6a (klf6a) (Figure 1C). Cluster zChol2 uniquely showed a higher degree of overlap with hepatocyte gene expression patterns (Figure 1C). Interestingly, cluster zAgr2+ expressed agr2, ictacalcin (icn), and transcobalamin beta a (tcnl), which were absent in other cholangiocyte clusters (Figure 1C). Cholangiocyte populations were enriched for GO processes that included translation, apoptosis, and metabolism (Figure S1B).

ECs and HSCs grouped into a single cluster (zEC/HSC), identified by zebrafish markers kinase insert domain receptor like (kdrl)[ ] and  heart and neural crest derivatives expressed 2 (hand2),[ ] respectively (Figure 1D). Other enriched genes in this cluster included suppressor of cytokine signaling 3a (socs3a), insulin‐like growth factor binding protein 7 (igfbp7), heat shock protein alpha‐crystallin‐related, 1, (hspb1), fatty acid binding protein 4a (fabp11a), and natriuretic peptide B (nppb).

Macrophages comprised two unique clusters and were identified by marker genes, macrophage receptor with collagenous structure (marco) and granulin a, tandem duplicate 1 (grn1)[ ] (Figure 1C). An inflammatory macrophage cluster (zInf mac) was identified by expression of inflammatory cytokines interleukin 1b (il1b) and chemokine (C‐X‐C motif) ligand 19 (cxcl19) [ ] (Figure 1C) and were enriched for immune processes, including respiratory burst and regulation of neutrophil migration (Figure S1C). In contrast, a noninflammatory macrophage cluster (zNon‐inf mac) was identified with limited expression of inflammatory cytokines, and these cells were enriched for negative regulation of intrinsic apoptotic signaling in response to DNA damage and hemoglobin metabolic process (Figure S1C).

Populations of natural killer (NK) and T cells (zNK/T), erythrocytes (zEry), and neutrophils (zNeu) were also captured (Figure 1A). These populations were identified by known marker genes interleukin 7 receptor (il7r), and T‐cell receptor alpha constant (trac),[ ] hemoglobin, alpha adult 1 (hbaa1),[ ] and myeloid‐specific peroxidase (mpx)[ ] and matrix metallopeptidase 13a (mmp13a),[ ] respectively (Figure 1C,D). GO analysis for clusters zNK/T and zNeu revealed enrichment for immune response processes, while zEry was enriched for carbon dioxide transport (Figure 1C).

Finally, there was a small population that clustered separately from all other clusters (apelin‐positive [zApln+]). This cluster highly expressed apln, a gene that was absent in all other captured cell types. Other marker genes for this cluster included heparin‐binding epidermal growth factor‐like growth factor b (hbegfb), thrombospondin 1b (thbs1b), glycoprotein Ib platelet subunit beta (gp1bb), nuclear factor, erythroid 2 (nfe2), and glucosaminyl (N‐acetyl) transferase 4a (gcnt4a) (Figure 1C). apln is highly expressed in endothelial tip cells, a type of EC found at the leading ends of sprouting vasculature, and is involved in angiogenesis and vascular organization in zebrafish.[ ]

This scRNA‐seq data set provides an atlas for the distinct transcriptional profiles of the cell types that comprise the adult zebrafish liver and establishes an important foundation for using zebrafish as a tool to study liver function and disease.

Key cell types and transcriptional profiles in the human liver are conserved in zebrafish liver

As the first scRNA‐seq data set in adult zebrafish liver, we sought to determine the similarity of the adult zebrafish liver to the human liver single‐cell transcriptome. We coclustered our adult zebrafish liver single‐cell data set with a published human liver single‐cell atlas.[ ] Joint clustering of these data sets resulted in 21 transcriptionally distinct clusters (Figure 2A,B). Cells of like identities clustered together across species, as evidenced by the presence of both human and zebrafish cells in most clusters (Figure 2C,D; Table S4). Hepatocytes accounted for six unique clusters (zebrafish human [zh]Hep1‐6). Human cholangiocytes largely clustered with the agr2+ cholangiocyte (zhChol2) population from zebrafish, with few human cholangiocytes found in cluster zhChol1 (Figure 2C,D). EC populations from human and zebrafish largely clustered separately by species (zhEC1‐2) (Figure 2C,D). Human and zebrafish HSCs clustered together in a single cluster, but zebrafish EC markers were also expressed in this cluster (zhEC/HSC) (Figure 2B,D). Joint clustering identified two macrophage populations, one inflammatory (zhInf mac) and one noninflammatory (zhNon‐inf mac), each composed of both human and zebrafish cells (Figure 2C,D). Furthermore, zebrafish NK/T cells mostly clustered with human liver NK cells (zhNK), with equal proportions of the remaining cells clustering with either αβ (zhαβ T) or γδ T cells (zhγδ T) (Figure 2C,D). Finally, erythrocytes (zhEry) from both species were found in a single cluster (Figure 2C,D).

FIGURE 2

Joint clustering of adult human and zebrafish liver scRNA‐seq. (A) UMAP visualization of 21 clusters comprised of 8,444 adult human liver and 13,630 adult zebrafish liver cells. (B) Heatmap of gene expression for three of the top differentially expressed genes in each cluster. (C) UMAP of joint human and zebrafish clustering split by species. (D) Bar graph showing the percentage of cells contributed to each cluster from each species. (E) Dot plot showing differential expression of human cell type marker genes in select clusters split by species to show marker conservation

To further explore the similarity of human and zebrafish liver cell types involved in fibrosis, we compared expression of specific marker genes for human HSCs and other cell types that are known to regulate HSC activation.[ ] We found conservation of marker genes[ ] for noninflammatory macrophages (zhNon‐inf mac), NK cells (zhNK), ECs (zhEC1‐2), inflammatory macrophages (zhInf mac), HSCs (zhEC/HSC), cholangiocytes (zhChol2), and γδ T cells (zhγδ T) (Figure 2E) across species. Taken together, joint clustering of scRNA‐seq data from adult human and zebrafish livers revealed a high degree of similarity in the transcriptional profiles and specificity of marker gene expression across species.

Subclustering of zebrafish cluster zEC/HSC reveals distinct HSC and EC populations

HSCs are the key cell type in the development and progression of fibrosis.[ ] We sought to use our scRNA‐seq data to better understand stellate cell biology and their cellular interactions in vivo. In our WT zebrafish liver single‐cell transcriptome, ECs and HSCs clustered together, consistent with known EC–HSC interactions that regulate HSC activation.[ ] To identify other key cell types that interact with HSCs, we used CellPhoneDB[ ] to predict receptor–ligand interactions between zEC/HSC and other populations in the WT zebrafish liver. Immune cell clusters (zInf mac, zNon‐inf mac, zNK/T, zNeu) and zAgr2+ had the highest number of significant interactions with zEC/HSC (Figure S2); immune cells and cholangiocytes have known roles in regulating HSC activation in human liver fibrosis.[ ] Furthermore, zApln+ exclusively had interactions with cluster zEC/HSC (Figure S2). These predictions identified immune cells, zAgr2+ cholangiocytes, and the zApln+ population as key zebrafish HSC‐interacting partners.

To examine changes that occur in HSCs in vivo during liver fibrosis, we used an Mpi‐depleted zebrafish line (mpi) in which heterozygous Mpi‐mutant zebrafish livers demonstrate reduced Mpi enzymatic activity, increased fibrogenic gene expression, increased collagen deposition, and liver fibrosis similar to MPI‐CDG individuals.[ ] We performed scRNA‐seq on mpi adult zebrafish livers (n = 3; 15,286 cells) (Table S2) and clustered these samples with scRNA‐seq data from our WT zebrafish livers (Figure S3A). To investigate reciprocal interactions between HSCs and other cell types in physiological and fibrotic contexts, we subset and reclustered the zEC/HSC and identified HSC‐interacting clusters. This resulted in clusters of zInf mac, zNK/T1‐2, zNeu, zNon‐inf mac, zAgr2+, and zApln+ and also allowed us to identify a distinct population of zHSCs, zEC‐HSCs, and three distinct EC populations (zEC1, zEC2, and zEC3) (Figure 3A). All clusters contained both WT and mpi cells and were identified by marker genes determined by DGE (Figure 3B,C; Figure S3). The zHSC cluster was identified by hand2 expression,[ ] which was specifically expressed in the zHSC and zEC‐HSC clusters, although only captured in a small percentage of cells (Figure 3B). Top DEGs for cluster zHSC included collagen, type I, alpha 1b (col1a1b), interferon nduced transmembrane protein 1 (ifitm1), chemokine (C‐C motif) liigand 25b (ccl25b), secreted protein, acidic, cysteine‐rich (osteonectin) (sparc), STEAP family member 4 (steap4), and collectin sub‐family member 11 ( colec11) (Figure 3B). mpi livers had a higher percentage of captured HSCs on average, although this difference was not significant (Figure 3C; Table S5).

FIGURE 3

Characterization of zebrafish HSCs and ECs. (A) UMAP visualization of 13 clusters comprised of WT and mpi adult zebrafish liver cells subset from total liver cell clustering. (B) Dot plot of gene expression for top DEGs in clusters zEC1, zEC2, zEC3, zEC‐HSC, and zHSC. (C) Bar graph showing the percentage of cells contributed to each cluster from each sample. (D) Heatmap showing expression of modules of coregulated genes correlating with EC and HSC clusters

Three distinct EC clusters (zEC1, zEC2, and zEC3) were identified by expression of kdrl [ ] (Figure 3B). Top DEGs for EC clusters included oncoprotein induced transcript 3 (oit3), protein tyrosine phosphatase receptor tybe b(ptprb), and receptor (G protein‐coupled) activity modifying protein 2 (ramp2) for zEC1 and claudin 5b (cldn5b), selectin E (sele), and nppb for zEC2 and zEC3 (Figure 3B).

Liver sinusoidal endothelial cells (LSECs) are known to physically interact with HSCs.[ ] We identified a cluster present in all WT and mpi samples that demonstrated expression for genes that were enriched in both EC and HSC clusters (zEC‐HSC) (Figure 3B). Furthermore, the cells in this cluster had a higher number of detected UMIs on average than in the zHSC or zEC clusters, indicating that this population may be made up of EC and HSC doublets (Figure S3C).

The human liver contains unique EC subtypes,[ ] and LSECs are known to play a critical role in the regulation of HSC activation.[ ] We sought to determine whether unique EC subpopulations similarly exist in the zebrafish liver. We subset the WT zebrafish EC and HSC populations (Figure S4A) and used the pseudotime analysis tool Monocle[ , , ] as an additional method to determine DGE across ECs and HSCs. Rather than identifying DGE between already defined clusters, Monocle identifies genes that are expressed in focal regions of cells in UMAP to create modules of gene co‐expression that are agnostic of cluster identity. We identified DEGs that grouped into unique expression modules that specifically correlated with zEC and zHSC clusters, indicating that these populations are transcriptionally distinct from one another (Figure 3D). Furthermore, cluster zEC‐HSC showed moderate expression scores for modules specific to zHSC and zEC clusters, further suggesting that these cells are likely EC‐HSC doublets (Figure 3D) We confirmed that distinct modules had cluster‐specific expression through visualizing module expression scores in single cells in UMAP and validating increased expression of module genes (Figure S4B,C). These results demonstrate that the zebrafish EC clusters in our data set have significantly different transcriptional profiles, indicating that the zebrafish liver contains distinct subtypes of EC populations similar to human liver.[ ]

Demonstrating conservation of HSC gene expression with human and zebrafish joint clustering

Because the use of adult zebrafish to study liver fibrosis has been limited, we sought to determine the conservation of EC and HSC populations between zebrafish and human livers. WT zebrafish ECs and HSCs were jointly clustered with human HSCs, LSECs, and portal ECs. Recapitulating the populations identified in MacParland et al.,[ ] this resulted in three distinct human EC populations (zone 1 LSECs [zhLSEC1], zone 2/3 LSECs [zhLSEC2/3], and portal ECs [zhPEC]) (Figure 4A). There were three other EC populations; two were primarily comprised of zebrafish cells (zhEC1 and zhEC2) and one was equal proportions of human and zebrafish cells (zhEC3) (Figure 4B; Figure S5A; Table S6). Overall, zebrafish and human ECs did not cluster together due to poor conservation of human EC marker gene expression in zebrafish (Figure 4B; Figure S5A,B), an important consideration when studying the role of ECs in the zebrafish liver.

FIGURE 4

Joint clustering of adult human and zebrafish ECs and HSCs. (A) UMAP visualization of nine clusters comprised of adult human liver and adult zebrafish liver ECs and HSCs. (B) Bar graph showing the percentage of cells contributed to each cluster from each species. (C) Dot plot showing differential expression of HSC marker genes split by species to show marker conservation. (D) Immunofluorescent staining for COLEC11 and COL1A1 on WT adult zebrafish liver cryosections (imaged at magnification ×63) and LX‐2 cells (imaged at magnification ×20). (E) Immunofluorescent staining for GFP under the hand2 promoter and COLEC11 on WT adult zebrafish liver cryosections (imaged at magnification ×63). zhEC/Imm, endothelial cell/immune cell mix

In contrast to the ECs, the majority of human and zebrafish HSCs clustered together (zhHSC) (Figure 4B; Figure S5A). Expression of human HSC marker genes α‐actin 2 (ACTA2), collagen type I alpha 1 chain (COL1A1), collagen type I alpha 2 chain (COL1A2), SPARC, and COLEC11 [ , ] was observed in over 75% of all human HSCs (Figure 4C). Of these genes, only SPARC and COLEC11 were highly captured in zebrafish HSCs, while known zebrafish HSC marker, HAND2, and markers determined from our zebrafish scRNA‐seq data set, angiopoietin‐like 6 (ANGPTL6) and STEAP4, were poorly captured in human HSCs (Figure 4C). Commonly used protein markers for human HSCs, including glial fibrillary acidic protein (GFAP), platelet‐dervied growth factor recepor beta (PDGFRB), lecithin retinol acyltransferase (LRAT), desmin (DES),[ , ] were poorly captured by scRNA‐seq in both human and zebrafish HSCs (Figure 4C). These findings indicate that traditional markers of human and zebrafish HSCs at the transcriptional level show lower levels of conservation. However, using our scRNA‐seq data, we found COLEC11, a C‐type lectin containing collagen‐like and carbohydrate recognition domains, to be highly expressed in human and zebrafish HSCs (Figure 4C). We tested its expression and specificity for zebrafish HSCs by immunofluorescent staining in WT adult zebrafish liver tissue. COLEC11‐expressing cells colocalized with COL1A1 in both zebrafish liver and human HSCs by immunofluorescent staining (Figure 4D) and colocalized with an established zebrafish HSC marker, Hand2[ ] (Figure 4E). To our knowledge, this is the first use of Colec11 as a novel marker for zebrafish HSCs and demonstrates conservation of COLEC11 across human and zebrafish HSCs.

Taken together, while human and zebrafish ECs have distinct gene expression patterns, human and zebrafish HSCs show conservation in transcriptional profiles and marker gene expression, with Colec11 demonstrating specific expression in HSCs in zebrafish. The similarity of human and zebrafish HSCs further substantiates the use of the zebrafish to study HSC activation and liver fibrosis.

Key pathways in HSC activation, angiogenesis, and immune response are altered in MPI‐depleted HSCs

To identify key shared canonical pathways and upstream regulators that are altered in the context of fibrosis in both humans and zebrafish, we compared responses to MPI depletion in zebrafish and human HSCs. We first created a stable MPI mutant human HSC line (MPI MT LX‐2s) using clustered regularly interspaced short palindromic repeat (CRISPR)/Cas9 (Figure S6A,B). MPI MT LX‐2 cells had a residual MPI activity of 21.4% compared to WT controls (Figure 5A) and decreased amounts of MPI protein (Figure 5B). MPI MT LX‐2 cells also demonstrated elevated COL1A levels, indicative of an activated phenotype (Figure 5B). We also used HSCs with acute depletion of MPI using siRNA (siMPI) in TWNT‐4 human HSCs, which had a residual MPI activity of 44.4% compared to controls (Figure 5C). We then conducted bulk RNA‐seq on the stable MPI mutant LX‐2s and siMPI TWNT‐4s.

FIGURE 5

Activation phenotypes are conserved across human and zebrafish HSCs. (A) Bar graph showing MPI enzymatic activity in WT and MPI MT LX‐2 cells. (B) Western blot for MPI and COL1A in WT and MPI MT LX‐2 cells. (C) Bar graph showing MPI enzymatic activity in control and siMPI TWNT‐4s. (D,E) IPA conducted on DGE from the following three data sets: mpi zebrafish single‐cell HSCs (zHSCs), MPI mutant LX‐2 HSCs, and siMPI TWNT‐4 HSCs. Heatmaps display (D) −log10(p values) for canonical pathway activity alteration and (E) z score of activation/inhibition of upstream regulators. *p < 0.05, **p < 0.01. ILK, integrin‐linked kinase; LXR, liver X receptor; RXR, retinoid X receptor; STAT3, signal transducer and activator of transcription 3

To compare the effects of MPI depletion in human and zebrafish HSCs, we performed IPA on DEGs from bulk RNA‐seq of both human MPI‐depleted HSC lines and DEGs from mpi HSCs in our single‐cell data set. We compared significance scores for all altered pathways in IPA and found conserved significantly altered pathways in MPI‐depleted HSCs across zebrafish and human data sets. The hepatic fibrosis signaling pathway was significantly activated in all data sets, which further supports that loss of Mpi activates HSCs in vitro and in vivo. Other pathways included those involved in HSC activation, angiogenesis, and immune response, all of which are known to play key roles in the pathogenesis of liver fibrosis[ ] (Figure 6A). We also compared the activation z score of upstream regulators across all three data sets; significant predicted alterations included EGF, fibroblast growth factor (FGF), extracellular signal‐regulated kinase 1 and 2 (ERK1/2), and Pdgf (complex), which are known to be key signaling pathways in HSC activation[ ] (Figure 6B). Taken together, these data indicate that MPI depletion in in vivo zebrafish and in vitro human HSCs lead to similar significant alterations in fibrogenic pathways, demonstrating the power of zebrafish HSCs to recapitulate human HSC phenotypes.

FIGURE 6

Identifying key alterations in Mpi‐depleted cholangiocytes. (A) Bar graph showing Mpi enzymatic activity in WT and mpi adult zebrafish livers. (B) Immunofluorescent staining for 2f11 and Col1a1 on WT and mpi adult zebrafish liver cryosections imaged at magnification ×63. (C) Dot plot of receptor–ligand interaction scores determined by CellPhoneDB analysis of scRNA‐seq gene expression in WT and mpi zebrafish liver cells. Dot color indicates log2(mean expression score) with values closest to zero indicating the highest expression score. −Log10(p value) is represented by dot size, with increasing dot size indicating increasing significance. (D) Immunofluorescent staining for 2f11 and Nrp1 on WT and mpi adult zebrafish liver cryosections imaged at magnification ×63. White arrows indicate representative bile ducts for each sample. (E) Violin plot for Nrp1 expression in Sox9+ cholangiocytes from scRNA‐seq data on livers from control, BDL, and CCl4‐treated mice. **p < 0.01, ***p < 0.001

Expression of Nrp1 is up‐regulated in cholangiocytes in liver fibrosis

In order to validate our phenotypic findings from the mpi mutant line, we generated a second mpi mutant zebrafish line, mpi, harboring a deletion mutation in exon 3 (Figure S7A‐C). mpi adult zebrafish livers demonstrated a 44.2% reduction in Mpi enzymatic activity (Figure 6A) as well as notable Col1a1 deposition surrounding dilated bile ducts, marked by anti‐neurofilament monoclonal antibody (2f11),[ ] when compared to WT livers (Figure 6B); this phenocopied mpi adult zebrafish and MPI‐CDG patient histology with distorted and dilated bile ducts that are embedded in dense fibrosis.[ , ] Given the shared portal‐based fibrosis, we sought to investigate the impact of Mpi loss on cholangiocytes, the cell type that makes up intrahepatic bile ducts.[ ]

Our scRNA‐seq data demonstrated MPI to be expressed broadly across liver cell types in human and zebrafish (Figure S8A,B), and we sought to investigate the effects of Mpi depletion in other liver cell types, with a focus on cholangiocytes given the biliary‐type fibrosis. We performed CellPhoneDB[ ] analysis on our adult zebrafish liver scRNA‐seq data to identify altered receptor–ligand interactions between cholangiocytes and other cell types across WT and mpi livers. We found a predicted significant interaction for NRP1 (a VEGFA receptor) and VEGFA between zHSC and zChol1 in mpi livers but not WT (Figure 6C). Using our scRNA‐seq data set, we found nrp1a to be expressed by cholangiocytes, ECs, HSCs, and macrophages, while nrp1b appeared to be restricted to ECs and HSCs (Figure S9). To confirm these findings at the protein level, we performed immunofluorescent staining for Nrp1 and demonstrated a strong increase in Nrp1 in the cholangiocytes (marked by staining of 2f11 along the luminal membrane), comprising the dilated bile ducts surrounded by increased collagen deposition of mpi livers when compared to WT controls (Figure 6D).

To determine if induction of Nrp1 in cholangiocytes was only observed in Mpi depletion or was more broadly found in other models of biliary fibrosis, we queried published scRNA‐seq data on livers from mice with bile duct ligation (BDL),[ ] a known model of biliary fibrosis, compared to controls and CCl4‐induced liver fibrosis conditions. We analyzed the SRY‐box transcription factor 9 (Sox9)+ cholangiocyte population found in this data set and found that Nrp1 was significantly up‐regulated in livers from BDL mice when compared to control and CCl4 livers (Figure 6E). These data suggest that Nrp1 may be a mediator of Vegf signaling in cholangiocytes during biliary fibrosis and warrants further investigation.

DISCUSSION

The mechanisms regulating liver fibrogenesis are extremely complex and include multiple cell types and cell–cell interactions.[ ] We have created the first single‐cell atlas of the adult zebrafish liver, comprised of transcriptionally unique populations of hepatic cell types, including hepatocytes, biliary cells, ECs, immune cells, and HSCs. When choosing an animal model to study a human disease, the degree of conservation for specific cell types is critical. Although conventional studies have shown that the zebrafish liver is similar to the human liver,[ ] our comparative analysis using scRNA‐seq reveals highly conserved marker genes in like cell types demonstrating shared identity and functional roles in the zebrafish and human liver. While the zebrafish liver is structurally different from the human liver and not known to have architectural zonation,[ , ] our scRNA‐seq data set revealed three main but transcriptionally distinct groups of hepatocytes enriched for pathways involved in oxidative phosphorylation (zHep1, zHep7, and zHep8); cholesterol, steroid, and lipid metabolism (zHep2, zHep3, zHep4, and zHep5); and fatty acid transport and glucose metabolism (zHep6). Future studies to determine whether these differences are related to the spatial location of these hepatocytes in the zebrafish liver or indicate zonation of the zebrafish liver similar to that of the human liver can include spatially resolved transcriptomics[ ] and RNAscope.

Given that HSCs are the primary cell type in liver fibrosis,[ ] conservation of this cell type in zebrafish is critical for establishing the adult zebrafish as an effective model for studying liver fibrosis. We found that zebrafish HSCs demonstrate a high similarity to human HSCs. This is seen through both coclustering and conservation of marker genes across species. Furthermore, tools to identify and study HSCs in the zebrafish liver have been limited thus far. Only a transgenic zebrafish line with green fluorescent protein under control of the hand2 promoter (Tg[hand2:GFP])[ ] and GFAP immunofluorescent staining[ ] have been shown to identify HSCs in zebrafish. Our work expands this armamentarium by identifying Colec11 as a new zebrafish HSC marker that could be used to both visualize HSCs in vivo and isolate zebrafish HSCs.

ECs have established roles in both activation and reversion of HSCs in fibrogenesis.[ ] The human liver is known to have different EC types, including LSECs and portal ECs,[ ] but little is known about whether zebrafish liver ECs are a homogeneous or heterogeneous population. Our data revealed multiple EC populations characterized by unique modules of gene expression, demonstrating that discrete populations of ECs do exist within the zebrafish liver. However, conclusions are limited as human marker genes that discriminate between different EC populations do not show the same specificity in zebrafish. Further studies are needed to explore how these zebrafish EC subpopulations differ from one another in spatial location and function.

We have previously shown that loss of MPI induces HSC activation and liver fibrosis.[ ] Pathway analysis comparing three transcriptomic data sets across zebrafish and human HSCs demonstrates conserved fibrogenic responses to loss of MPI and highlights the role of MPI in the regulation of HSC activation. Furthermore, we leveraged our scRNA‐seq data set to uncover increased Nrp1 expression in zebrafish cholangiocytes involved in biliary fibrosis. We validated this finding in cholangiocytes from the livers of mice with BDL‐induced fibrosis, an alternative model of biliary fibrosis. VEGF has been implicated as important for cholangiocyte proliferation, which is associated with biliary fibrosis; however, only VEGF receptor 2 (VEGFR‐2) and VEGFR‐3 have been shown to be mediators of this increased signaling.[ , ] This induction suggests that NRP1 may also play a role in mediating VEGF signaling and promoting biliary fibrosis. Furthermore, NRP1 has been implicated in other liver cell types in fibrosis, including HSCs and LSECs,[ , ] making NRP1 a potentially interesting target for future study and development of antifibrotic therapies.

In summary, we demonstrate the feasibility and utility of scRNA‐seq in the adult zebrafish liver. Our findings highlight the similarities and differences to human liver, support its use as a valuable tool to study liver fibrosis, introduce new tools to study HSCs, and lend insights into important in vivo interactions as potential therapeutic candidates for further investigation.

CONFLICT OF INTEREST

Nothing to report.

Fig S1

Click here for additional data file.

Fig S2

Click here for additional data file.

Fig S3

Click here for additional data file.

Fig S4

Click here for additional data file.

Fig S5

Click here for additional data file.

Fig S6

Click here for additional data file.

Fig S7

Click here for additional data file.

Fig S8

Click here for additional data file.

Fig S9

Click here for additional data file.

Supplementary Material

Click here for additional data file.