# Single cell transcriptomics of Atlantic salmon (Salmo salar L.) liver reveals cellular heterogeneity and immunological responses to challenge by Aeromonas salmonicida

Richard Taylor, Rose Ruiz Daniels, Ross Dobie, Shahmir Naseer, Thomas Clark, Neil Henderson, Pierre Boudinot, Samuel A.M. Martin, Daniel Macqueen

<figure><img src="image_1.png" /></figure>

## To cite this version:

Richard Taylor, Rose Ruiz Daniels, Ross Dobie, Shahmir Naseer, Thomas Clark, et al. Single cell transcriptomics of Atlantic salmon (Salmo salar L.) liver reveals cellular heterogeneity and immunological responses to challenge by Aeromonas salmonicida. Frontiers in Immunology, 2022, 13, pp.984799. 10.3389/fimmu.2022.984799. hal-04011974

HAL Id: hal-04011974

https://hal.inrae.fr/hal-04011974

Submitted on 2 Mar 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

Check for updates

OPEN ACCESS

EDITED BY
Uwe Fischer,
Friedrich-Loeffler-Institute, Germany

REVIEWED BY
Alexander Rebl,
Research Institute for Farm Animal
Biology (FBN), Germany,
Charlotte Scott,
Flanders Institute for Biotechnology,
Belgium

*CORRESPONDENCE
Daniel J. Macqueen
daniel.macqueen@roslin.ed.ac.uk

†These authors have contributed
equally to this work and share
first authorship

SPECIALTY SECTION
This article was submitted to
Comparative Immunology,
a section of the journal
Frontiers in Immunology

RECEIVED 02 July 2022
ACCEPTED 03 August 2022
PUBLISHED 24 August 2022

CITATION
Taylor RS, Ruiz Daniels R, Dobie R, Naseer S, Clark TC, Henderson NC, Boudinot P, Martin SAM and Macqueen DJ (2022) Single cell transcriptomics of Atlantic salmon (*Salmo salar* L.) liver reveals cellular heterogeneity and immunological responses to challenge by *Aeromonas salmonicida*. *Front. Immunol.* 13:984799.
doi: 10.3389/fimmu.2022.984799

COPYRIGHT
© 2022 Taylor, Ruiz Daniels, Dobie, Naseer, Clark, Henderson, Boudinot, Martin and Macqueen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Single cell transcriptomics
of Atlantic salmon (*Salmo
salar* L.) liver reveals
cellular heterogeneity
and immunological
responses to challenge by
*Aeromonas salmonicida*

Richard S. Taylor¹†, Rose Ruiz Daniels¹†, Ross Dobie²,³,
Shahmir Naseer⁴, Thomas C. Clark⁵, Neil C. Henderson²,³,
Pierre Boudinot⁵, Samuel A.M. Martin⁴
and Daniel J. Macqueen¹*

†The Roslin Institute and Royal (Dick) School of Veterinary Studies, University of Edinburgh, Edinburgh, United Kingdom, ´Centre for Inflammation Research, The Queen's Medical Research Institute, Edinburgh BioQuarter, University of Edinburgh, Edinburgh, United Kingdom, ³MRC Human Genetics Unit, Institute of Genetics and Cancer, University of Edinburgh, Edinburgh, United Kingdom, ⁴School of Biological Sciences, University of Aberdeen, Aberdeen, United Kingdom, ⁵ Université Paris-Saclay, INRAE, UVSQ, VIM, Jouy-en-Josas, France

The liver is a multitasking organ with essential functions for vertebrate health spanning metabolism and immunity. In contrast to mammals, our understanding of liver cellular heterogeneity and its role in regulating immunological status remains poorly defined in fishes. Addressing this knowledge gap, we generated a transcriptomic atlas of 47,432 nuclei isolated from the liver of Atlantic salmon (*Salmo salar* L.) contrasting control fish with those challenged with a pathogenic strain of *Aeromonas salmonicida*, a problematic bacterial pathogen in global aquaculture. We identified the major liver cell types and their sub-populations, revealing poor conservation of many hepatic cell marker genes utilized in mammals, while identifying novel heterogeneity within the hepatocyte, lymphoid, and myeloid lineages. This included polyploid hepatocytes, multiple T cell populations including γδ T cells, and candidate populations of monocytes/macrophages and dendritic cells. A dominant hepatocyte population radically remodeled its transcriptome following infection to activate the acute phase response and other defense functions, while repressing routine functions such as metabolism. These defense-specialized hepatocytes showed strong activation of genes controlling protein synthesis and secretion, presumably to support the release of acute phase proteins into circulation. The infection response further involved up-regulation of numerous genes in an immune-cell specific manner, reflecting functions in pathogen recognition and killing, antigen presentation, phagocytosis, regulation of inflammation, B cell differentiation
and T cell activation. Overall, this study greatly enhances our understanding of the multifaceted role played by liver immune and non-immune cells in host defense and metabolic remodeling following infection and provides many novel cell-specific marker genes to empower future studies of this organ in fishes.

KEYWORDS

liver, single cell transcriptomics, bacterial infection, salmonid fish, immune-metabolism
cross talk

**Introduction**

The vertebrate liver is a multitasking organ with diverse physiological functions, including nutrient metabolism, transport and storage, growth signaling, endocrine regulation, and immunity (1). In mammals, these roles are performed by the cooperative actions of several distinct cell types including hepatocytes, cholangiocytes (epithelial cells of the bile duct), stellate cells, Kupffer cells (resident liver macrophages) and lymphocytes (1). Recent advances in single cell transcriptomics have revealed functional heterogeneity within the major hepatic cell types of mammals (2, 3), providing insights into liver spatial organization (4, 5) while revealing cellular and molecular drivers of disease and malignancy states (6, 7).

The adult liver of all vertebrates contains both immune and non-immune cells with important immunological functions (2, 7-9) that support immune homeostasis and tolerance (10) and the generation of inflammatory responses upon pathogen challenge, leading to secretion of acute phase proteins (APPs) into circulation by hepatocytes (7, 11). The liver is the major site of haematopoiesis in the mammalian fetus and hence an important organ for early immune cell development (12), though this feature is not conserved in fishes (13). The multifaceted functions of the liver must demand tight coordination of different cell types to achieve appropriate responses to prevailing physiological and environmental conditions, inclusive of immune system status following pathogen challenge. Immunological functions may need to be prioritized at the cost of investment into metabolic functions in such scenarios (14). In this regard, the role of the liver in co-regulating metabolism and immunity makes it an interesting organ to understand the coordinated responses of different cell types following pathogenic challenge.

Bulk RNA-Seq and proteomics have been widely applied to understand liver functions in commercially important fishes. For example, in Atlantic salmon (*Salmo salar* L.), which is among the most important aquaculture species globally (15), such work has uncovered the role of this organ in innate immune defense and the acute phase response (11, 16). In the biomedical model

zebrafish (*Danio rerio*), single cell transcriptomics has recently been used to reveal hepatic cellular heterogeneity, including within the hepatocyte, myeloid and lymphoid lineages (17, 18). However, zebrafish are distantly related to salmonids, leaving a gap in knowledge on the role of liver cellular heterogeneity in this key group of fishes, where such information can be applied to understand and manipulate health and immune traits relevant to sustainable aquaculture and food production.

The aim of this study was to reveal the major cell lineages within the liver of Atlantic salmon and uncover the role of hepatic cellular heterogeneity in the host response to bacterial infection. Using single-nuclei RNA-Seq (snRNA-Seq), we report a comprehensive single cell transcriptomic atlas of the Atlantic salmon liver, identifying novel heterogeneity across multiple hepatic cell types. By comparing cell-specific responses in control animals to those challenged with *Aeromonas salmonicida*, we uncover a dramatic transcriptomic remodeling of hepatocytes that underpins the acute phase response, alongside major changes in gene expression specific to distinct immune cell populations.

**Materials and methods**

**Disease challenge and sampling**

Animal work was carried out in compliance with the Animals (Scientific Procedures) Act 1986 under Home Office license PFF8CC5BE and was approved by the ethics committee of the University of Aberdeen. Atlantic salmon were kept in 250 L freshwater tanks at aquarium facilities of the University of Aberdeen. Water temperature was maintained at 14°C and fish were fed a commercial pellet diet at 2% body weight per day. After two weeks, twenty fish were anaesthetized using 2-phenoxyethanol (2.5 mL in 10L water/0.0025% v/v) and given an intraperitoneal injection of either PBS (0.5 mL) (n=10), or the pathogenic Hooke strain of *A. salmonicida* (2 × 10⁵ colony-forming units/mL in PBS; 0.5 mL/fish) (n=10). Sampling was performed 24 h post-injection (after ref 11). All twenty fish were
killed using a Schedule 1 method following anesthesia overdose using 2-phenoxyethanol (0.1% v/v) and destruction of the brain. Fish were immediately sampled, with liver samples (approx. 100mg) flash frozen on dry ice before storage at -80°C prior to snRNA-Seq library construction. Separate liver samples were placed in 1.5ml of Tri Reagent (Sigma-Aldrich) and used for quantitative PCR (qPCR) to validate the expected response in infected fish.

## Validation of immune response by qPCR

qPCR validation was performed on the twenty liver samples (comparing n=10 control vs. n=10 infected animals) to confirm an inflammatory response in the infected fish. RNA was extracted from 100 mg liver tissue using 1 mL Trizol reagent. The tissue homogenization was performed using two Tungsten Carbide Beads (Qiagen) (3 mm) on a Tissuelyser II (Qiagen) at a frequency of 30.0 I/sec for 2.5 min. 200 µl of chloroform (Sigma-Aldrich) was added and the mixture centrifuged at 12,000 g for 20 minutes at 4°C to separate the aqueous phase, which was retained. RNA precipitation was performed using 700 µl of isopropanol (Sigma-Aldrich) and centrifuged at 12,000 g for 20 min at 4°C and washed 3 times with 80% ethanol. The concentration and purity of total RNA was estimated using a NanoDrop 1000 Spectrophotometer (Thermo Scientific). A QuantiTect Reverse Transcription kit (Qiagen) was used to synthesize first-strand cDNA from 1 µg total RNA per sample, with a genomic DNA removal step included, in a total volume of 20 µl. The resulting cDNA was diluted 20-fold (working stock) with RNase/DNase free water. qPCR was performed using primers targeting the APP encoding genes *saa* and *hamp* (11). For normalization of gene expression, primers targeting *rps13* and *rps29* were used (19). qPCR was performed with 2x SYBR Green I (Invitrogen) master mix on a Mx3005P System (Agilent Technologies). Each assay was run in triplicate using 15 µl of reaction mix containing 7.5 µl Brilliant III Ultra-Fast SYBR Green (Agilent Technologies), 500 nM forward and reverse primers and 5 µl of cDNA (2.5 ng of total reverse-transcribed RNA). Assays were run with 1 cycle of 95°C for 3 min, followed by 40 cycles of 95°C for 20 s and 64°C for 20 s. A melting curve (thermal gradient 55°C - 95°C) was used to confirm single qPCR products. Each qPCR plate included two no-template (i.e. water) controls. LinRegPCR (20) was used to establish the efficiency of each assay. Gene expression data were analyzed using GeneX 5.4.3 (MultiD Analysis), correcting for differences in efficiency across genes, before target gene expression was normalized to *rps13* and *rps29* and placed on a scale of relative expression. Separate Kruskal-Wallis tests were performed in Minitab® 17.1.0 to compare gene expression values for *saa* and *hamp* between the control (n=10) and infected (n=10) liver samples.

## snRNA-Seq library construction

Based on the results of qPCR, n=4 samples were taken forward for snRNA-Seq library construction, representing n=2 control fish and n=2 infected fish (see Results; Supplementary Figure 1). A protocol adapted from ref (21) was used for nuclear extraction, employing a tween with salt and tris (TST) buffer. Approximately 45 mg of each frozen liver sample was placed in a 6-well tissue culture plate (Stem Cell Technologies) with 1 mL TST (2 mL of 2X ST buffer + 120 µL of 1% Tween-20 + 20 µL of 2% BSA brought up to 4 mL with nuclease-free water). The tissue was minced using Noyes spring scissors for 10 min on ice. The resulting homogenate was filtered through a 40 µm Falcon cell strainer, and a further 1 mL of TST was added to wash the well and filter. The volume was brought up to 5 mL using 3 mL of 1X ST buffer (diluted from 2xST buffer [292 µl of 146 mM NaCl, 100 µl of 10 mM Tris-HCl pH 7.5, 10 µl of 1 mM CaCl₂, 210 µl of 21 mM MgCl₂, brought up to 10 mL with nuclease-free water]). The sample was centrifuged at 4°C for 5 min at 500g before the resulting pellet was re-suspended in 1 mL 1X ST buffer and the recovered nuclei were filtered through a 40 µm Falcon cell strainer, Hoechst stained, visually inspected under a fluorescent microscope, and counted using a Bio-Rad TC20. Liver nuclei were processed through the 10X Chromium™ Single Cell Platform using the Chromium™ Single Cell 30 Library and Gel Bead Kit v3.1 and Chromium™ Single Cell A Chip Kit (both 10X Genomics) as per the manufacturer's protocol. For each sample, the nuclei were loaded into a channel of a Chromium 3' Chip and partitioned into droplets using the Chromium controller before the captured RNA for each nucleus was barcoded and reverse transcribed. The resulting cDNA was PCR amplified for 14 cycles, fragmented, and size selected before Illumina sequencing adaptor and sample indexes were attached. Libraries were sequenced on a NovaSeq 6000 by Novogene UK Ltd (2x150bp paired end reads).

## Generation of snRNA-Seq count matrix

Raw sequencing data were aligned to the unmasked ICSASG_v2 reference assembly (Ensembl release 104) of the Atlantic salmon genome (22). The analysis was restricted to protein coding genes. Mapping of reads to the genome, assignment of reads to cellular barcodes, and collapsing of unique molecular identifiers (UMIs) was performed with StarSolo v2.7.7a (23). The genome index was generated with standard settings and "sjdbOverhang" set to 149. The reads were then mapped with the "STAR" command and following settings: "soloType=CB_UMI_Simple", "outFilterMultimapNmax = 20", "outMultimapOrder(random", "soloUMIdedup=1MM_Directional", "soloFeatures = GeneFull", "soloBarcodeReadLength = 0", "outFilterMatchNminOverLread = 0", "soloCellFilter = TopCells 100000". The top 100,000 cell
barcodes ranked by UMI number were retained to ensure the
capture of transcriptionally quiet nuclei, lost when using the
automated StarSolo filtering algorithm. Mapping statistics for each
snRNA-Seq sample are provided in **Supplementary Table 1**.

Nuclei filtering and quality control

Nuclei filtering was performed manually on each of the four
samples using Seurat v3.1 (24). Ranked barcode plots of UMIs
and gene counts were used to identify the lower ‘elbow point’
and cell barcodes with UMI count or gene counts below the
elbow point removed as empty droplets. Further steps were used
to identify additional empty droplets: the “SCTransform”
function (25) was used to normalise the data, prior to
centering, scaling, principal component analysis (PCA) and
high-resolution clustering (using 30 PCs and resolution of 2).
Wilcoxon rank sum differential gene expression tests were used
to identify up-regulated genes in each cluster (see Methods
Section ‘Differential gene expression tests’). Cell clusters that
both lacked distinguishing markers and had a low median
UMI or gene counts (typically 2 median absolute deviations
lower than the median across all nuclei) were removed as likely
empty droplets or poor-quality nuclei. This process was repeated
iteratively on each sample until all such low-quality populations
were removed. Likely doublets were identified and removed later
in the analysis, after cell identity was established.

Assignment of cellular identity

Each cell cluster was assigned to one of the major liver cell
lineages using a *priori* marker genes (**Supplementary Table 2**).
Populations identified as erythrocytes and thrombocytes were
not included in downstream analyses due to their likely origin
from contaminating blood. Populations identified as
hepatocytes, cholangiocytes, mesenchymal cells, endothelial
cells, and immune cells were merged into five separate Seurat
objects (using the Seurat *merge* function) for separate analyses
of each cell lineage. Batch effects across samples were removed in
each merged object with Harmony (26). Gene annotations were
taken from the Ensembl annotation for Atlantic salmon, and in
cases where the gene was not assigned a name, the name of the
nearest annotated putative orthologue in rainbow trout
(*Oncorhynchus mykiss*), zebrafish or mouse (*Mus musculus*)
was used, or in some cases informed by BLASTp searches
against the NCBI non-redundant database.

For each of the five major liver cell lineages retained, data
was log normalised, then scaled and centered, before PCA was
performed and used as the input to graph-based clustering,

using the established Seurat pipeline. The appropriate number
of PCs to use for clustering of each sample to minimise
technical noise, was determined through visualisations
generated through the “ElbowPlot” and “DimHeatmap”
functions. The resolution parameter in the “FindClusters”
function was tuned to return biologically meaningful
heterogeneity with each cell type. At this stage, a further
quality control step was performed to remove doublets from
each lineage. Specifically, differential gene expression tests were
performed and clusters that exhibited canonical markers
(Supplementary Table 2) or lineage distinguishing markers
(Supplementary Table 3) from two distinct cell lineages yet
lacked any unique distinguishing marker genes of their own,
were removed as likely doublets, before the data was re-
clustered. This process was repeated until no doublet
populations remained. For the immune cells, the Seurat
object was further split into T cell, Myeloid, NK-like,
Neutrophil and B cell objects (see Results) and the same
process was repeated (using markers in Supplementary
Table 3) and cell sub-clusters identified, before each object
was merged back into a single “immune cell” Seurat object.
Finally, all five major lineages were merged into a final global
liver cell atlas object, with cell identities retained from the cell
lineage specific analyses. We used this strategy as opposed to
global clustering, as markedly more biologically meaningful
heterogeneity could be established in the cell lineage
specific analyses.

Differential gene expression tests

Throughout this study, differentially expressed genes were
defined using the Seurat function “FindAllMarkers”, applying
the Wilcoxon rank sum test with default cut-offs (multi-test
adjusted p-value < 0.05, log2-fold change > 0.25, with expression
of the gene in at least 20% of nuclei in the cluster tested). For the
global analyses of major cell types, differential gene expression
tests were performed for each defined cluster in turn versus all
other nuclei in the dataset as the background. For the
hepatocyte-specific analysis we compared each hepatocyte
cluster to all other hepatocyte nuclei. For the immune cell
nuclei, the comparisons made were between each T cell sub-
cluster and all other T cell nuclei, each myeloid sub-cluster and
all other myeloid nuclei, and between each of the other immune
cell clusters (i.e. NK-like, Neutrophil and B cell clusters) and all
other immune cell nuclei. To assess which genes were up-
regulated in infection in the defined immune sub-clusters, a
differential gene expression test was performed within each
cluster to compare nuclei that originated from control vs.
infected fish (Wilcoxon rank sum test, p-value < 0.05, log2-
fold change > 0.25, with expression of the gene in at least 20% of
nuclei from either the infected or control nuclei).
# Results

## Single-nuclei RNA-Seq atlas of the Atlantic salmon liver

Hepatocytes, the major epithelial cell type within the liver, have been underrepresented in mammalian scRNA-Seq analyses (2), which may be caused by damage occurring during the dissociation step. As Atlantic salmon are ectotherms, we were also concerned that enzymatic dissociation (typically done at >30°C) would activate cell stress and heat shock responses. We consequently decided to generate a snRNA-Seq atlas of salmon liver, using nuclei isolated from freshly flash frozen samples, an approach expected to provide an accurate representation of cell diversity (e.g. 27, 28). The profiled liver nuclei were from control fish (n=2) and animals infected by a pathogenic strain of *Aeromonas salmonicida* (n=2), the bacterial agent of furunculosis, a long-standing problem disease in salmonid aquaculture (29). The infected group were sampled 24 hours post *Aeromonas* challenge, previously shown to capture

the inflammatory and acute phase response (11). The fish used for sequencing were further selected based on gene expression data using marker genes for the acute phase response (11), which were robustly and significantly up-regulated in the infected group (Supplementary Figure 1).

Across all samples, we generated 47,432 nuclei transcriptomes with median UMI and gene counts per nucleus of 2,105 and 1,065, respectively. This was split across control fish as 11,679 and 11,433 nuclei, and *Aeromonas*-challenged fish as 19,148 and 5,172 nuclei, respectively (Supplementary Table 1). The major cell lineages were identified using a guided graph-based clustering strategy (24), with cell identity assigned using a priori defined marker genes (Figures 1A, B; markers in Supplementary Table 2). Hepatocytes comprised most of the profiled nuclei (88.1%), followed by cholangiocytes (4.3%), immune cells (3.5%), mesenchymal cells (2.6%) and endothelial cells (1.5%) (Supplementary Figure 2; Supplementary Table 1). Clustering and identification of these major cell lineages was repeatable across individual samples, with each sample contributing a large proportion of the nuclei (Supplementary Figure 2). Differential gene expression analysis

<figure><img src="image_4.png" /><figcaption>FIGURE 1</figcaption></figure>

Major cell types in the Atlantic salmon liver defined by 47,982 nuclear transcriptomes. (A) UMAP highlighting five main liver cell-type clusters according to a *priori* defined marker genes (see Supplementary Table 2). (B) Bubble plots showing the expression of a *priori* marker genes for the five main liver cell types, including mean expression (bubble intensity) and percentage of nuclei expressing gene (bubble size). (C) Heatmap of the top 20 differentially expressed genes per each liver cell type defined against the background of all other cell types. (D) Violin plots demonstrating the expression of the most specific marker genes per each of the main liver cell types (colours are matched to the colours of the 5 cell lineages defined in part A).
revealed markers for each major liver cell type (Figures 1C, D; Supplementary Table 3). Marker genes for candidate cell types and sub-populations are hereafter reported according primarily to annotations provided by Ensembl, or in cases where no Ensembl annotation was available, using supplementary BLAST homology searches against the NCBI database to support our inferences. While Ensembl annotation utilizes phylogenetic information to inform homology relationships, it may nonetheless fail to correctly capture orthology of Atlantic salmon genes to mammalian species, particularly for fast evolving and complex gene families. Adding to this challenge, genetic orthology is not a prerequisite for conservation of gene function or expression, which limits our ability to transfer knowledge about cell marker genes from mammalian studies to Atlantic salmon. While the reader must be aware of these caveats, they represent a general issue in functional genomics studies using non-model taxa like salmonids.

While it was possible to identify the five major liver cell lineages of Atlantic salmon using orthologues to marker genes defined in mammalian liver scRNA-Seq datasets (e.g. 6), many markers were notably absent or expressed at very low levels in our snRNA-Seq dataset (Supplementary Figure 3). For example, Atlantic salmon orthologues of the widely used epithelial marker *Epcam* did not show expression in the cholangiocyte cluster (Supplementary Figure 3). Likewise, salmon orthologues of *Pecam1* and *Pdgfrb*, which are excellent markers of mammalian endothelial and mesenchymal cells, were not detected at significant levels in these cell types in our dataset (Supplementary Figure 3). These results may be explained by differences in transcriptome composition between snRNA-Seq and scRNA-Seq datasets from liver, as identified in humans (30). In support of this idea, *epcam* and *pecam1* were respective markers of cholangiocytes and endothelial cells in a recently published zebrafish scRNA-Seq dataset (31), making it less likely that the lack of expression of these markers in our Atlantic salmon dataset represents a true evolutionary difference between mammals and teleosts. Differences in expression were also observed between predicted Atlantic salmon orthologues of mammalian marker genes for the major hepatic cell types, with, for example, only one of two *cdh5* co-orthologues marking the endothelial population, and only one of two *hnf4a* co-orthologues marking hepatocytes (Supplementary Figure 3).

Higher resolution clustering captured varying degrees of transcriptomic heterogeneity for each of the five major liver cell types (Figures 2A, B; Supplementary Table 4), which was consistent across the four samples (Supplementary Figure 4). We observed a split of hepatocytes into sub-populations explained largely by infection status, with 75.9% of ‘control-associated’ hepatocytes deriving from control fish and 73.3% of ‘infection-associated’ hepatocytes deriving from *Aeromonas*-challenged fish (expanded in next section). Significant markers for the latter were dominated by genes encoding APPs, many of which were observed to be expressed across all cell types (Figure 2B). This is likely a consequence of the numerically dominant hepatocytes (Figure 1A) leaking mRNA from highly

expressed genes into the ambient RNA. Similar examples of this type of leakage can be observed in many liver scRNA-Seq datasets with, for example, hepatocyte expressed genes *Alb* and *Hp* being widely ‘expressed’ across non-hepatocyte cells in past mammalian studies (e.g. 3, 32).

The immune compartment contained identifiable T and B cells, along with candidate populations of neutrophils, myeloid cells, and NK-like cells (Figure 2B). The transcriptome of cholangiocyte nuclei was homogeneous, while limited heterogeneity was identified in the endothelial and mesenchymal cells (Figure 1A; Supplementary Figure 5 and 6). The endothelia sub-clusters had a clear biological interpretation, with Atlantic salmon orthologues to marker genes from mammals distinguishing arterial and venous derived endothelial cells (33) (Supplementary Figure 5, marker genes in Supplementary Table 5). However, the mesenchymal sub-clusters were not readily biologically interpretable (Supplementary Figure 6; Supplementary Table 6).

## Hepatocyte remodelling dominates the liver response to bacterial infection

To explore how hepatocyte heterogeneity contributes to the response to *Aeromonas* infection, we analysed 41,792 available hepatocyte nuclei transcriptomes. Clustering using the most variable genes in this compartment identified nine sub-populations (H1-H9) (Figure 3A; marker genes in Supplementary Table 7), with several showing marked differences in abundance between control and infected fish (Figure 3B). H1-H4 comprised 90.2% of hepatocyte nuclei, with H1 and H2 deriving mainly from control fish and showing highly correlated transcriptomes (Figures 3B, C). H3 and H4 comprised 70.0% of nuclei from infected fish and showed closely related transcriptome profiles (Figures 3B, C).

Hepatocyte nuclei derived from infected fish increased from H1 (16.9%), to H2 (28.2%), to H3 (64.3%) to H4 (80.8%), with 2,842 genes differentially expressed on this gradient (Figures 3D, E; Supplementary Tables 8 and 9). 379 genes were up-regulated in infection-dominated H4 vs. control-dominated H1, showing overrepresented functions linked to host defense and the acute phase response (‘complement activation’, ‘defense response to other organism’, and ‘cellular iron homeostasis’), in addition to translational processes (e.g. ‘translational elongation’) (Figure 3F, Supplementary Table 10). This response was dominated by genes encoding APPs including hepcidin, haptoglobins, ferritin, transferrin, fibrinogens, ceruloplasmin, angiotensinogen, serum albumins, apolipoproteins, and c-reactive protein (Supplementary Table 8). One of the top up-regulated genes (3.7-fold up-regulated; ENSSSAG00000046715) encodes mechanistic target of rapamycin kinase (mTOR), a master regulator of translation (34). mTOR biasedly promotes translation of ribosomal protein genes (34), which is notable as
<figure><img src="image_5.png" /><figcaption>FIGURE 2</figcaption></figure>

Higher resolution atlas of Atlantic salmon liver cells defined by snRNA-Seq. (A) Unbiased graph-based clustering reveals varying heterogeneity levels across the major liver cell types, presented on a UMAP. Each cell population retains the gene signature of the parent lineage (Supplementary Figure 2), while also displaying its own distinct transcriptomic profile, presented here as a heatmap (B), inclusive of the top 10 marker genes based on differential gene expression against all other cell clusters. The colour bars above columns on the heatmap illustrate the cell types to which the genes shown were identified as markers (matched to part A).

28 such genes, encoding many proteins comprising the large and small ribosomal subunits of salmonid fishes (35), were up-regulated in H4 vs. H1, with only one downregulated (Supplementary Tables 8 and 9). We further observed up-regulation of *copb1* (ENSSSAG00000068484), *grp78* (ENSSSAG00000054661), *sec61a1* (ENSSSAG00000039008), and *srprb* (ENSSSAG00000046456), encoding major components of the COPI and translocon complexes, representing key protein secretion pathways (36).

H5 nuclei were mainly from infected fish (Figure 3B) and showed highly correlated transcriptomes to H3/H4 (Figure 3C), sharing many of the same key markers up-regulated in H4 vs. H1, but also specific markers (Supplementary Table 7) associated with NF-κB signaling, including *relb* (ENSSSAG00000052551), encoding a component of the NF-κB transcription factor complex. This is notable, as H5 also showed the highest expression among all hepatocyte sub-populations of *stat3* (ENSSSAG00000003657), encoding signal transducer and

activator of transcription 3, which acts downstream of NF-κB in mammalian hepatocytes to activate inflammation driving the acute phase response (37).

2,278 genes were downregulated in H4 vs. H1 (Supplementary Table 9), enriched in functions related to signaling (e.g. 'intracellular receptor signaling pathway'), metabolism (e.g. 'gluconeogenesis'), cell differentiation (e.g. 'stem cell differentiation') and transcription (e.g. 'regulation of transcription, DNA-templated') (Figure 3F) (Supplementary Table 11). The top downregulated gene was *prkag2a* (4.0-fold down-regulated; ENSSSAG00000079550), encoding a subunit of the AMPK complex - a master regulator of metabolism including gluconeogenesis (38). Also downregulated were master hepatic transcription factors for lipid and glucose metabolism pathways that interact with AMPK (38), including genes encoding hepatocyte nuclear factor 1 (HNF1) (ENSSSAG00000006158) and HNF4 (ENSSSAG00000047055) (39), carbohydrate-responsive element-binding protein (ENSSSAG00000039257) (40), and forkhead box proteins O1/O3
A

<figure><img src="image_6.png" /></figure>

B

<figure><img src="image_7.png" /></figure>

C

<figure><img src="image_8.png" /></figure>

D

<figure><img src="image_9.png" /></figure>

E

<figure><img src="image_10.png" /></figure>

F

<table><caption>FIGURE 3</caption><thead><tr><th>GO term</th><th>P-value</th></tr></thead><tbody><tr><td>Intracellular receptor signaling pathway.</td><td>3.1e-06</td></tr><tr><td>Small molecule catabolic process.</td><td>1.6e-04</td></tr><tr><td>Gluconeogenesis.</td><td>5.0e-03</td></tr><tr><td>Cellular carbohydrate catabolic process.</td><td>6.1e-03</td></tr><tr><td>Regulation of transcription, DNA-templated.</td><td>8.2e-03</td></tr><tr><td>Regulation of cell population proliferation.</td><td>7.8e-03</td></tr><tr><td>Signal transduction.</td><td>2.0e-02</td></tr><tr><td>Complement activation.</td><td>4.0e-07</td></tr><tr><td>Translational elongation.</td><td>5.3e-07</td></tr><tr><td>Cellular iron ion homeostasis.</td><td>7.5e-06</td></tr><tr><td>Negative regulation of fibrinolysis.</td><td>2.3e-04</td></tr><tr><td>Defense response to other organism.</td><td>3.5e-03</td></tr><tr><td>Response to bacterium.</td><td>5.0e-03</td></tr><tr><td>Ribosome biogenesis.</td><td>1.9e-02</td></tr></tbody></table>

Striking remodeling of the hepatocyte transcriptome in response to bacterial infection. (A) UMAP visualisation of 41,792 hepatocyte nuclei, with sub-clustering performed using the most variable genes restricted to this cell lineage. (B) Shows the same UMAP with nuclei coloured by infection state (top) and the proportion of nuclei originating from each sample after normalising for different nuclei numbers across samples (bottom). (C) Pearson correlation of the expression values for the top 2,000 most variable genes across the nine hepatocytes populations H1-H9. (D) Violin plots of mean expression for the 20 most down-regulated (left) and up-regulated (right) genes based on log-fold change in H4 vs. H1. (E) Heatmap of the top 20 most upregulated and top 20 most downregulated genes in H4 relative to H1, illustrating a gradient of expression from H1 — H2 — H3 — H4 (F) Example enriched GO terms in H4 based on all up-regulated and all down-regulated genes in H4 vs. H1 (full data provided in Supplementary Tables 10 and 11).

(ENSSSAG00000054669/ENSSSAG00000055241) (41) (Supplementary Table 9). Further evidence for repression of anabolism included downregulation of genes encoding growth hormone receptor (GHR) (ENSSSAG00000065355 and ENSSSAG00000081526) and Stat5a/b (ENSSSAG00000010616, ENSSSAG00000003584 and ENSSSAG00000048873), which act

downstream of GHR to activate growth and cell proliferation genes (42). In addition, a gene was downregulated encoding glucocorticoid receptor (GR) (ENSSSAG00000062169), which promotes expression of gluconeogenesis genes, while its interaction with Stat5 is required for transcriptional activation of growth genes *via* GHR signalling (42).
## Additional hepatocyte heterogeneity includes polyploid cells

The remaining hepatocyte subclusters were not linked to infection status (Figure 3B). H7 likely represents hepatocytes that have undergone polyploidization, which occurs progressively during aging in mammals, such that 4n-16n cells make up a large fraction of liver cells by adulthood (43). Polyploidy in H7 is indicated by a striking concordance of genes representing highly specific markers for H7 and those shown elsewhere to be up-regulated in 4n vs. 2n mammalian cells, encoding DNA primase subunit 2 (ENSSSAG00000001875), replication protein A 70 kDa DNA-binding subunit (ENSSSAG00000050927, ENSSSAG00000050209) and DNA polymerase (ENSSSAG00000078390), among others (44). The close relationship of H7 to H1 (Figures 3A, C) implies that polyploid hepatocytes derive from those supporting routine metabolism.

H9 showed many specific markers encoding mitosis proteins, e.g. cytoskeleton-associated protein 5 (ENSSSAG00000066206), centromere protein E (ENSSSAG00000073454), abnormal spindle-like microcephaly-associated protein homolog (ENSSSAG00000053226) and kinesin family member 23 (ENSSSAG00000044381), indicating these are cycling hepatocytes. H6 expressed a small number of highly specific markers, including genes encoding ligand of numb-protein X1 (lnx1) (ENSSSAG00000072728 and ENSSSAG00000070275), a E3 ubiquitin ligase that targets a wide range of proteins, including CD8 expressed on T-cells in mammals (45). H8 was biasedly represented by control fish and most correlated with H1/H2 in transcriptome profile (Figure 3C), expressing marker genes described as ‘novel’ in the Ensembl annotation (Supplementary Table 7).

## Immune cell heterogeneity in the Atlantic salmon liver

3.4% (n=1,620) of the liver nuclei were derived from immune cells (Figures 4A, B). Combinations of canonical marker expression was used to classify T cells (cd3e, tox2 and tcf7), B cells (ighm, cd37, cd79a), NK-like cells (prf1.3 and runx3 with absence of cd3e) and myeloid cells (mpeg1, cd63, csf1r, lyz2) (Figure 4B; marker genes provided in Supplementary Table 12). We also identified a candidate population of neutrophils based on marker genes that showed highest specificity or expression for neutrophils among different immune cells in the human protein atlas, namely itgax (ENSSSAG00000049715), ncf1 (ENSSSAG00000079828), and mmp9 (ENSSSAG00000069874) (Figure 4B). Mmp9 was an effective marker for neutrophils in other teleost species, and has a role in driving neutrophil migration in mammals (46, 47). The NK-like cells represent a tentative annotation owing to a lack of certain markers for NK cells, i.e. genes encoding granzymes. Each immune cluster had

substantial contributions from all samples, except the candidate neutrophils, which derived mainly from one Aeromonas-challenged fish (Supplementary Figure 4). Sub-clustering revealed heterogeneity in the T and myeloid cells, but not the NK-like cells, B cells (with no evidence of plasma cells) or candidate neutrophils.

Five T cell sub-clusters were compared using differential expression tests (Figure 4A; Supplementary Table 12). Only T1 and T2 expressed cd4 genes to low levels, specifically ENSSSAG00000076631, encoding CD4-1, which is also expressed by macrophages, and ENSSSAG00000076595, encoding CD4-2; shown elsewhere to be expressed by all CD4⁺ T cells (48). We did not identify any CD8 expressing T cells, likely due to low expression levels. T1 was the largest sub-cluster but showed few specific markers. T2 expressed genes encoding receptors involved in T cell activation. This included cd28 (ENSSSAG00000060163), encoding the main co-stimulatory T cell receptor (49) and cd44 (ENSSSAG00000076128), a widely used T cell adhesion, co-stimulation and activation marker (50). However, T2 cells did not specifically express any genes annotated as ctla-4, an IgSF member induced during T cell activation that regulates CD28 activity. Furthermore, T1 and T2 both expressed tcf7 (ENSSSAG00000006857) at a much higher level than T3-T5, encoding the master Wnt pathway transcription factor, which is most highly expressed on naïve mammalian T cells (51) and was a specific marker for resting CD4⁺ T cells in humans (52). T1 and T2 also expressed foxp1b (ENSSSAG00000077820) at a higher level than T3-5; a gene essential for quiescent naïve T cells in mammals (53). T1 and T2 therefore appear to be constituted mainly of resting T cells. T3 appears to be an activated T cell population based on specific expression of slamf1 (ENSSSAG00000043093) encoding CD150 (54). T3 expressed many highly specific markers, including the integrin coding gene itgal (ENSSSAG00000046537; encoding CD11a) and itgb2 (ENSSSAG00000022772; encoding CD18), whose products form lymphocyte function-associated antigen 1 (LFA-1), a molecule with key roles in T cell activation and migration, in addition to cytotoxic and memory responses (55). In mammals, LFA-1 is an established marker for the migration of liver T resident memory (T~RM~) cells into liver, which is also the case for CD103 (also known as ITGAE) (56), which was expressed more highly in T3 (ENSSSAG00000076346) than other T sub-clusters. The human orthologue of a highly expressed T3 specific marker, cacna2d2a (ENSSSAG00000071299), encoding a calcium voltage-gated channel, was not detected in any immune cell in the human protein atlas, implying a teleost-specific T cell marker. T4 expressed several activation markers including pou2f2a (aka oct2) (ENSSSAG00000071136) (57), CD226 (58), along with two distinct ctram genes (encoding cytotoxic and regulatory T-cell molecule), previously shown in mammals to be required for differentiation of cytotoxic CD4⁺ T cells (59). T5 specifically expressed two paralogues of sox13 (ENSSSAG00000077869,
<figure><img src="image_11.png" /><figcaption>FIGURE 4
Heterogeneity in Atlantic salmon immune cells. (A) UMAP visualization of 1,620 immune nuclei. (B) Heatmap of top 5 markers genes for each immune sub-cluster, sorted by log2 fold change. (C) Number of cell-specific genes up-regulated by infection in immune sub-clusters (D) Examples of genes showing cell-specific up-regulation in response to Aeromonas infection across the breath of immune cell heterogeneity identified.</figcaption></figure>

ENSSSAG00000058488), the defining vertebrate transcription factor for the γδ T lineage (60).

My1 and My2 markers are associated with monocytes and macrophages (Supplementary Table 12). My1 specifically expressed *cd4-1* (ENSSSAG00000076631), at a level higher than any T sub-cluster, likely representing a phagocytic CD4⁺ macrophage characterized in rainbow trout (48). My2 expressed specific monocyte marker genes including *csf3r* (ENSSSAG00000041566) and *timp2* (ENSSSAG00000042353 and ENSSSAG00000064056). Two *csf1r* copies were identified with reciprocal higher expression in My1 (ENSSSAG0000004088) and My2 (ENSSSAG00000061479). High *flt3* (ENSSSAG0000009390)

expression in My4 supports an annotation as dendritic cells (DCs) (61), with specific up-regulation of *cd9* (ENSSSAG00000059637) and *lamp2* (encoding CD107b) (ENSSSAG00000074801) genes, consistent with monocyte-derived DCs in mammals (62, 63). My3 expressed the second highest level of *flt3*, while the top My3 marker gene, *ptprsa* (ENSSSAG00000051752), is a specific marker for plasmacytoid DCs (pDCs) in mice and human (64). The relationship of My3 to pDCs is also supported by specific expression of *sbf1* (ENSSSAG00000071635) (65). However, another My4-specific marker gene, *tcf4* (ENSSSAG00000071044), encodes a transcription factor required for pDC development (66).
Immune cell-specific responses to
*Aeromonas* insult

To uncover the role of hepatic immune cell heterogeneity in the response to *Aeromonas* challenge, we performed differential expression tests comparing nuclei from control and infected fish within each immune sub-cluster (Figures 4C, D). 819 genes showed significant up-regulation in nuclei from infected fish (criteria: *P* < 0.05; Log2FC > 1), among which 274 (33%) and 72 (8.8%) were up-regulated by most (≥7 of the 12 sub-clusters) or all immune cell types, respectively, and 271 (33%) showed immune sub-cluster specific up-regulation (Figure 4C; Supplementary Table 13).

B cells specifically up-regulated two genes, namely *stat3* (2.4-fold up-regulated; ENSSSAG00000003657), which is essential for B cell differentiation (67) and a gene annotated zgc:174904 (2.1-fold up-regulated; ENSSSAG00000070511), encoding a 304 amino acid protein with a CD209/DC-Sign-like, C-type lectin-like domain (InterPro domain: IPR033989). C-type lectin/DC-SIGN is a broad-specificity PRR that detects bacteria by binding mannose or carbohydrate structures (68). In the candidate neutrophils, the top up-regulated gene was *ladderlectin* (38.9-fold up-regulated; ENSSSAG00000039613), encoding a soluble lectin that bound *Aeromonas* in salmonids, leading to bacterial killing actions (69). The NK-like cells up-regulated few genes specifically, one of which was *pglyrp2* (2.3-fold up-regulated; ENSSSAG00000054105), a peptidoglycan recognition protein with enzymatic activity targeting and limiting the inflammatory effects of bacterial peptidoglycan (70).

Among the T cell populations, T3, T4, and the γδ T cells (T5), showed the strongest responses to *Aeromonas* infection. The top up-regulated genes in T3 included *csl1* (4.3-fold up-regulated; ENSSSAG00000004327), encoding a L-rhamnose-binding lectin that binds bacteria and enhances phagocytosis in salmonids (71), and *c8b* (4.1-fold up-regulated; ENSSSAG00000073702), encoding a core component of the complement membrane attack complex. In T4, up-regulated genes included *tln1*, encoding Talin-1 (4.8-fold up-regulated; ENSSSAG00000063331), which is known to regulate the integrin LFA-1 complex (defining T3; see last section), and is required for sustained interactions between APCs and T cells, as well as T cell proliferation (72). T4 also up-regulated *itgb1* (encoding CD29, also called β1-integrin) (4.8-fold up-regulated; ENSSSAG00000007621), a signature marker for cytotoxic T cells in humans (73). In T5, among the top-up-regulated genes was *prex1* (4.8-fold up-regulated; ENSSSAG00000044871), a signaling molecule that promotes expression of key interleukin cytokines in activated human T cells, including IL-2 (74). T5 also up-regulated *catl1*, encoding cathepsin L (8.0-fold up-regulated; ENSSSAG000000077309), which regulates T cell cytotoxicity (75) and an unannotated gene encoding a protein with saposin-like domains, which is annotated as *prosaposin*-like in NCBI (74.4-fold up-regulated; ENSSSAG00000009411). In mammals,

*prosaposin* encodes the precursor to all saposin lysosomal proteins, which are known to have antibacterial activity and play a key role in presenting lipid antigens to CD1-restricted T cells (76).

My1 specifically up-regulated 17 genes, including a different *itgal* paralogue (encoding CD11A) to that noted as a marker for T3 (4.2-fold up-regulated; ENSSSAG00000046996), encoding a component of LFA-1 essential to the immune response of mice to *Mycobacterium tuberculosis*, supporting T cell-mediated activation and recognition of infected macrophages (77). My1 also specifically up-regulated *olfm4* (4.1-fold upregulated; ENSSSAG00000046003), encoding a glycoprotein induced in mice macrophages by *Helicobacter pylori* infection, which regulates inflammatory responses (78). My1 further specifically up-regulated *cats* (ENSSSAG00000070942) and *ctsd* (ENSSSAG00000027269), encoding cathepsin S and D, proteolytic enzymes with established macrophage roles in bacterial killing and antigen processing (79). My2 specifically up-regulated 32 genes, including *cath2* (6.6-fold upregulated; ENSSSAG00000049319), encoding an antimicrobial peptide that increased in abundance in response to *Aeromonas* infection in salmon plasma (80) and *abr* (6.6-fold upregulated; ENSSSAG00000080204), encoding a GTPase-activating protein that down-regulates the inflammatory actions of macrophages (81).

My3 specifically up-regulated 16 genes, including *il12p40b2* (5.7-fold up-regulated; ENSSSAG00000069633), best known as a component of IL-12 and/or IL-23 heterodimers, but that also has defined cytokine functions as a monomer protein, including promotion of DC migration in response to bacterial infection in mammals (82, 83). Another induced gene was *afp4* (3.7-fold up-regulated; ENSSSAG00000072959), encoding type IV ice-structuring protein LS-12, an apolipoprotein-like molecule that dramatically increases in abundance in salmon plasma following *Aeromonas* infection (84). My4 specifically up-regulated 64 genes, including *gas7* (8.56-fold up-regulated; ENSSSAG00000076110), which has a crucial role in phagocytosis (85), *dnase113* (8.4-fold up-regulated; ENSSSAG00000066441), which controls inflammasome-induced cytokine secretion (86), and *cd82* (6.9-fold up-regulated; ENSSSAG00000052206), previously shown to be up-regulated during DC activation, where is promotes stable interactions between DCs and T cells, and MHC-II maturation (87). My4 further upregulated *irf4* (6.28-fold up-regulated; ENSSSAG00000039730), a gene essential to the ability of DCs to promote Th2 differentiation and inflammation (88).

Discussion

This study greatly enhances our knowledge of liver function in a salmonid fish with global commercial and scientific importance. The major advancement compared to previous
work comes from our application of snRNA-Seq, which, in
contrast to previously past bulk transcriptomic or proteomic
studies, allowed us to identify multiple hepatic cell populations,
before dissecting the role of this heterogeneity in host defense
following bacterial infection. Furthermore, a plethora of novel
marker genes are reported for developmentally and functionally
diverse hepatic cell types, which will be useful for future studies
investigating traits relevant to salmonid health and
immunological status.

Our results demonstrate the essential contribution of
hepatocytes to the antibacterial and acute phase response in
Atlantic salmon. Transcriptomic heterogeneity in the dominant
hepatocyte population (i.e., H1-4; Figure 3) was inconsistent
with distinct hepatocyte populations. Instead, our data supports
a single hepatocyte population that can exist in radically distinct
transcriptional states dependent on infection status. At one
extreme are the hepatocytes that dominate the liver of healthy
fish (H1/2), which appear to be performing routine functions
controlled by master hepatic transcription factors and signaling
pathways. Conversely, the liver of bacterially infected Atlantic
salmon was dominated by hepatocytes (H3/4) that
downregulated master hepatic pathways (e.g. controlling
metabolic functions), and up-regulated a suite of genes
encoding APPs and innate immune molecules. This includes
many APPs routinely detected in Atlantic salmon plasma by
proteomics (89) implying extremely high abundance. As most
plasma proteins derive from liver, this aligns with our finding
that these ‘defense-specialized’ hepatocytes strongly up-
regulated mTOR, its target ribosomal protein-coding genes,
and genes from protein secretion pathways, presumably to
boost APP translation and secretion rates during the acute
phase response. This striking repurposing of hepatocyte
function upon infection illustrates the vital role these cells play
as a hub for cross-talk between metabolism and immunity,
presumably allowing energetic resources to be allocated
towards clearing a pathogen at the short-term cost of limiting
investment into routine hepatic functions e.g. supporting
growth (14).

Interestingly, a sub-population of defense-specialized
hepatocytes (i.e. H5; Figure 3), almost exclusively derived from
infected fish, specifically expressed genes associated with NF-κB
signaling (e.g. relb was ~10-fold more highly expressed in this
population compared to the average across the eight other
hepatocyte subpopulations). In addition, H5, among all
hepatocyte sub-populations, showed the highest expression of
stat3, which is indispensable for activation of APP and protein
secretion pathway gene expression during bacterial infection in
mice, acting downstream of NF-κB (37, 90). As APP and
secretory protein pathway genes were strongly up-regulated in
the dominant sub-populations of defense-specialized
hepatocytes (H3/4), which lacked significant relb expression,
H5 may represent an intermediate hepatocyte state, where the
activation of the APP response and associated secretory pathway

is first initiated. A previous scRNA-Seq study of zebrafish liver
failed to identify hepatocytes showing any equivalent
specialization towards host defense (18), while another
identified a minor hepatocyte population enriched for immune
functions (17), which may be analogous to H3-H5. Differences
with these past zebrafish studies may reflect the fact that both
studies utilized control zebrafish lacking immune stimulation.
However, the fact that H3-5 comprised a significant fraction of
hepatocytes in our control fish could also be explained by
differences in liver function, potentially indicative of a higher
baseline inflammatory state in the liver of the Atlantic salmon
population we studied.

While representing a small population of liver nuclei, we also
offer evidence for polyploid hepatocytes in Atlantic salmon. The
functional role of polyploidy in mammalian liver remains ill-
defined, despite extensive study over decades (43). Work done
over 40 years ago showed that the liver of several teleosts
contained polyploid hepatocytes (91), so our result is perhaps
not unexpected. Polyploidy increases with aging in mammals
(43), which may explain why this hepatocyte population was so
limited in the fish used in our study, which were juveniles.
However, polyploid hepatocytes were not reported in past
scRNA-Seq studies of zebrafish liver, which included an 18-
month-old adult population (17), representing half the adult
lifespan for this species. More work is required to understand the
role of hepatocyte polyploidy in teleost health and disease.

It is important to acknowledge that comparing our results
with liver scRNA-Seq studies in other species has limitations due
to fundamental differences in experimental execution, which, for
example, is linked to striking differences in the composition of
the cells captured. For example, a preprint by ref (18)
surprisingly identified T cells as the dominant liver population
in zebrafish, with hepatocytes reflecting a smaller proportion
than expected. While underrepresentation of hepatocytes has
been observed in several mammalian liver scRNA-Seq studies, a
separate zebrafish scRNA-Seq study identified hepatocytes as the
most abundant liver cell type (17), albeit at a smaller fraction
than for our snRNA-Seq atlas in Atlantic salmon. This perhaps
illustrates a recognized benefit of snRNA-Seq compared to
scRNA-Seq; a more accurate representation of the true tissue
cell diversity (e.g. 27, 28, 94). Considering our limited knowledge
of cellular diversity in most fish species, including salmonids,
careful comparisons of results from scRNA-Seq and snRNA-Seq
will be required to establish baseline expectations for
future studies.

This is the first single cell transcriptomic study to
characterize immune cell heterogeneity in the liver of a
salmonid fish, and the first single cell study in any teleost to
characterize transcriptomic responses of specific hepatic
immune cell subtypes to infection. Past scRNA-Seq studies of
zebrafish liver paint very distinct pictures of lymphocyte
diversity, with one reporting no B cells, and little T cell
heterogeneity (17). Conversely, a recent preprint reported a
small B cell population, and six T cell sub-clusters capturing distinct CD8 and CD4 subsets (18). While we also identified a single B cell population and multiple T cell sub-clusters, the identity of T cells was markedly less clear in our data, due partly to a general lack of *cd4* and *cd8* expression, which may reflect a limitation of snRNA-Seq, or the lower sequencing depth in our study compared to (18). However, unlike these previous studies, we identified a small population of *sox13*⁺ *γδ* T cells, an ancient vertebrate T cell subtype with roles bridging adaptive and innate immunity. In zebrafish, *γδ* T cells possess phagocytic ability and act as APCs that activate CD4⁺ T cells, inducing B cell proliferation (93). In mammals, *γδ* T cells increase dramatically in liver during inflammatory conditions (94) and produce IL-17 essential for the innate response to bacterial infection (95). While *γδ* T cells have not been reported among the plethora of immune cells reported to date in scRNA-Seq studies spanning different teleost species (96), we find them readily identifiable in Atlantic salmon by *sox13* expression in multiple tissues (not shown). Consistent with their known functions, salmon liver *γδ* T cells up-regulated genes with roles spanning innate and adaptive immunity during the early response to *Aeromonas* infection.

Our study also identified evidence for myeloid heterogeneity within the Atlantic salmon liver, including two candidate DC populations that showed a strong response to bacterial infection, up-regulating genes required for interactions with T cells, phagocytosis and inflammasome activation, suggesting conserved roles between mammals and salmonids, as shown elsewhere (97). DCs were reported in a recent scRNA-Seq study of zebrafish liver (31), but not in two other liver scRNA-Seq analyses from the same species (17, 18). DCs have also been reported in Atlantic cod spleen (98). We observed two distinct candidate macrophage populations, matching the level of heterogeneity recently reported in two previous zebrafish studies (17, 18). A natural question relates to the relationship of teleost macrophages and mammalian Kupffer cells (KCs). Past work largely agrees that KCs are amniote-specific and hence not present in teleost liver (99), including for Atlantic salmon (100). Consistent with this notion, a recent study in zebrafish identified that resident macrophages were not located in the sinusoidal space, and lacked phagocytic ability (101). However, a recent cross-vertebrate scRNA-Seq analysis of liver cell repertoires reported putative KCs in the liver of zebrafish based on a conserved transcriptomic signature with mammals (31). This study did not report any heterogeneity within the macrophages beyond the putative KCs, notably lacking any monocyte-derived macrophages known to make up a large proportion of macrophages in mammalian livers (2, 3), and which share marker genes with our Atlantic salmon macrophage clusters. We failed to identify conserved expression of the signature marker genes for KCs (in addition to transcription factors associated with KCs) defined in this past study (31) in any of our myeloid sub-clusters from Atlantic salmon (Supplementary

Figure 9). While KCs were not expected among our macrophage populations, it is worth noting the expression of the predicted Atlantic salmon orthologues of mammalian *Clec4f* (ENSSSAG00000040735 and ENSSSAG00000076658). *Clec4f* encodes a C-Type Lectin (also known as Kupffer cell receptor) expressed specifically in KCs to the exclusion of non-liver macrophages (102, 103). In our Atlantic salmon dataset, one predicted orthologue of *Clec4f* is expressed across all four myeloid populations (ENSSSAG000000076658), showing highest expression in My2, whilst the other (ENSSSAG00000040735) was expressed across all immune cell types (Supplementary Table 11). Furthermore, both genes were among the top up-regulated genes in hepatocytes following *Aeromonas* infection (Supplementary Table 8). These findings highlight the potential for major differences in marker gene cell-specific expression and presumably function for homologous genes shared by mammals and teleosts.

Having robust marker genes for different cell-types is clearly essential for accurate studies of cell biology in any species. However, single cell transcriptomics studies performed to date, including this study, demonstrate that marker genes can vary markedly across species, especially for immune cells (96, 104). This demands a broader uptake of single cell transcriptomics in more species to define conserved from non-conserved marker genes, and to separate true biological or evolutionary differences in cell heterogeneity and associated marker genes from differences introduced by technical reasons discussed earlier, including the use of snRNA-Seq vs. scRNA-Seq. In salmonid fishes, the presence of an ancestral whole genome duplication, and the associated retention of numerous duplicated genes that have diverged extensively in tissue expression (22, 105–107), further challenges the transfer of knowledge on marker genes from model species - more work is required in this area. In summary, our comprehensive dissection of the Atlantic salmon liver using single cell transcriptomics has generated many species-specific marker genes for a range of immune and non-immune cell types that contribute to health and immunological traits of relevance to sustainable aquaculture. Our results can further be used to extract cell-specific information from existing and future bulk gene expression studies.

## Data availability statement

All data generated in this study has been made available in the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Accession: GSE207655.

## Ethics statement

The animal study was reviewed and approved by Animal work was carried out in compliance with the Animals (Scientific
Procedures) Act 1986 under Home Office license PFF8CC5BE and was approved by the ethics committee of the University of Aberdeen.

## Author contributions

Designed study: DM, RRD, RT; coordinated sampling and disease challenge experiment: SM; performed fish sampling: SN, SM, DM; performed quantitative PCR: SN; optimized nuclear isolation: RRD; generated snRNA-Seq libraries: RD, RRD; provided infrastructure for snRNA-Seq: NH; performed bioinformatics: RT; interpreted immunological data: TC, PB; drafted manuscript: DM, RT, RRD; made figures and tables: RT, DM; contributed to data interpretation and finalization of manuscript: all authors.

## Funding

This study was funded by grants from the Scottish Universities Life Sciences Alliance (Technology Seed Funding Call), the University of Edinburgh's Data Driven Innovation Initiative (Scottish Funding Council Beacon 'Building Back Better' Call), and the Biotechnology and Biological Sciences Research Council, including the institutional strategic programme grants BBS/E/D/10002071 and BBS/E/D/20002174 and the responsive mode grant BB/W005859/1. NH is supported by a Wellcome Trust Senior Research Fellowship in Clinical Science (ref. 219542/Z/19/Z).

## References

1. Trefts E, Gannon M, Wasserman DH. The liver. *Curr Biol* (2017) 27(21): R1147-51. doi: 10.1016/j.cub.2017.09.019

2. MacParland SA, Liu JC, Ma XZ, Innes BT, Bartczak AM, Gage BK, et al. Single cell RNA sequencing of human liver reveals distinct intrahepatic macrophage populations. *Nat Commun* (2018) 9(1):4383. doi: 10.1038/s41467-018-06318-7

3. Ramachandran P, Dobie R, Wilson-Kanamori JR, Dora EF, Henderson BE, Luu NT, et al. Resolving the fibrotic niche of human liver cirrhosis at single-cell level. *Nature* (2019) 575(7783):512-8. doi: 10.1038/s41586-019-1631-3

4. Dobie R, Wilson-Kanamori JR, Henderson BE, Smith JR, Matchett KP, Portman JR, et al. Single-cell transcriptomics uncovers zonation of function in the mesenchyme during liver fibrosis. *Cell Rep* (2019) 29(7):1832-1847.e8. doi: 10.1016/j.celrep.2019.10.024

5. Ben-Moshe S, Itzkovitz S. Spatial heterogeneity in the mammalian liver. *Nat Rev Gastroenterol Hepatol* (2019) 16(7):395-410. doi: 10.1038/s41575-019-0134-x

6. Ramachandran P, Matchett KP, Dobie R, Wilson-Kanamori JR, Henderson NC. Single-cell technologies in hepatology: new insights into liver biology and disease pathogenesis. *Nat Rev Gastroenterol Hepatol* (2020) 17(8):457-72. doi: 10.1038/s41575-020-0304-x

7. Stamataki Z, Swadling L. The liver as an immunological barrier redefined by single-cell analysis. *Immunology* (2020) 160(2):157-70. doi: 10.1111/imm.13193

8. Robinson MW, Harmon C, O'Farrelly C. Liver immunology and its role in inflammation and homeostasis. *Cell Mol Immunol* (2016) 13(3):267-76. doi: 10.1038/cmi.2016.3

## Acknowledgments

For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) license to any Author Accepted Manuscript version arising from this submission.

## Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.984799/full#supplementary-material

9. Zhao J, Zhang S, Liu Y, He X, Qu M, Xu G, et al. Single-cell RNA sequencing reveals the heterogeneity of liver-resident immune cells in human. *Cell Discovery* (2020) 6:22. doi: 10.1038/s41421-020-0157-z

10. Heymann F, Tacke F. Immunology in the liver-from homeostasis to disease. *Nat Rev Gastroenterol Hepatol* (2016) 13(2):88-110. doi: 10.1038/nrgastro.2015.200

11. Causey DR, Pohl MAN, Stead DA, Martin SAM, Secombes CJ, Macqueen DJ. High-throughput proteomic profiling of the fish liver following bacterial infection. *BMC Genomics* (2018) 19(1):719. doi: 10.1186/s12864-018-5092-0

12. Popescu DM, Botting RA, Stephenson E, Green K, Webb S, Jardine L, et al. Decoding human fetal liver haematopoiesis. *Nature* (2019) 574(7778):365-71. doi: 10.1038/s41586-019-1652-y

13. Mahony CB, Bertrand JY. How HSCs colonize and expand in the fetal niche of the vertebrate embryo: An evolutionary perspective. *Front Cell Dev Biol* (2019) 7:34. doi: 10.3389/fcell.2019.00034

14. Rauw WM. Immune response from a resource allocation perspective. *Front Genet* (2012) 3:267. doi: 10.3389/fgene.2012.00267

15. Houston RD, Macqueen DJ. Atlantic Salmon (*Salmo salar* L.) genetics in the 21st century: taking leaps forward in aquaculture and biological understanding. *Anim Genet* (2019) 50(1):3-14. doi: 10.1111/age.12748

16. Moraleda CP, Robledo D, Gutiérrez AP, Del-Pozo J, Yáñez JM, Houston RD. Investigating mechanisms underlying genetic resistance to salmon rickettsial syndrome in Atlantic salmon using RNA sequencing. *BMC Genomics* (2021) 22(1):156. doi: 10.1186/s12864-021-07443-2
17. Morrison JK, DeRossi C, Alter IL, Nayar S, Giri M, Zhang C, et al. Single-cell transcriptomics reveals conserved cell identities and fibrogenic phenotypes in zebrafish and human liver. *Hepatol Commun* (2022) 6(7):1711-24. doi: 10.1002/hep4.1930

18. Huang Y, Liu X, Wang H-Y, Chen JY, Zhang X, Li Y, et al. Single-cell transcriptome landscape of zebrafish liver reveals hepatocytes and immune cell interactions in understanding nonalcoholic fatty liver disease. *bioRxiv* (2022). doi: 10.1101/2022.02.06.479276

19. Macqueen DJ, Kristjánsson BK, Johnston IA. Salmonid genomes have a remarkably expanded akirin family, coexpressed with genes from conserved pathways governing skeletal muscle growth and catabolism. *Physiol Genomics* (2010) 42(1):134-48. doi: 10.1152/physiogenomics.00045.2010

20. Ruijter JM, Ramakers C, Hoogaars WM, Karlen Y, Bakker O, van den Hoff MJ, et al. Amplification efficiency: linking baseline and bias in the analysis of quantitative PCR data. *Nucleic Acids Res* (2009) 37(6):e45. doi: 10.1093/nar/gkp045

21. Slyper M, Porter CB, Ashenberg O, Waldman J, Drokhlyansky E, Wakiro I, et al. A single-cell and single-nucleus RNA-seq toolbox for fresh and frozen human tumors. *Nat Med* (2020) 26(5):792-802. doi: 10.1038/s41591-020-0844-1

22. Lien S, Koop BF, Sandve SR, Miller JR, Kent MP, Nome T, et al. The Atlantic salmon genome provides insights into redilapidization. *Nature* (2016) 533 (7602):200-5. doi: 10.1038/nature17164

23. Kaminow B, Yunusov D, Dobin A. STARsolo: accurate, fast and versatile mapping/quantification of single-cell and single-nucleus RNA-seq data. *bioRxiv* (2021). doi: 10.1101/2021.05.05.442755

24. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM, et al. Comprehensive integration of single-cell data. *Cell* (2019) 177(7):1888-1902.e21. doi: 10.1016/j.cell.2019.05.031

25. Hafemeister C, Satija R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. *Genome Biol* (2019) 20(1):296. doi: 10.1186/s13059-019-1874-1

26. Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, et al. Fast, sensitive and accurate integration of single-cell data with harmony. *Nat Methods* (2019) 16(12):1289-96. doi: 10.1038/s41592-019-0619-0

27. Wu H, Kirita Y, Donnelly EL, Humphreys BD. Advantages of single-nucleus over single-cell RNA sequencing of adult kidney: Rare cell types and novel cell states revealed in fibrosis. *J Am Soc Nephrol* (2019) 30(1):23-32. doi: 10.1681/ASN.2018090912

28. Denisenko E, Guo BB, Jones M, Hou R, de Kock L, Lassmann T, et al. Systematic assessment of tissue dissociation and storage biases in single-cell and single-nucleus RNA-seq workflows. *Genome Biol* (2020) 21(1):130. doi: 10.1186/s13059-020-02048-6

29. Park SY, Han JE, Kwon H, Park SC, Kim JH. Recent insights into *Aeromonas salmonicida* and its bacteriophages in aquaculture: A comprehensive review. *J Microbiol Technol* (2020) 30(10):1443-57. doi: 10.4014/jmb.2005.05040

30. Andrews TS, Atif J, Liu JC, Perciani CT, Ma XZ, Thoeni C, et al. Single-cell, single-nucleus, and spatial RNA sequencing of the human liver identifies cholangiocyte and mesenchymal heterogeneity. *Hepatol Commun* (2022) 6 (4):821-40. doi: 10.1002/hep4.1854

31. Guilliams M, Bonnardel J, Haest B, Vanderborght B, Wagner C, Remmerie A, et al. Spatial proteogenomics reveals distinct and evolutionarily conserved hepatic macrophage niches. *Cell* (2022) 185(2):379-396.e38. doi: 10.1016/j.cell.2021.12.018

32. Aizarani N, Saviano A, Sagar, Mailly L, Durand S, Herman JS, et al. A human liver cell atlas reveals heterogeneity and epithelial progenitors. *Nature* (2019) 572(7768):199-204. doi: 10.1038/s41586-019-1373-2

33. Kalucka J, de Rooij LP, Goveia J, Rohlenova K, Dumas SJ, Meta E, et al. Single-cell transcriptome atlas of murine endothelial cells. *Cell* (2020) 180(4):764-779.e20. doi: 10.1016/j.cell.2020.01.015

34. Saxton RA, Sabatini DM. mTOR signaling in growth, metabolism, and disease. *Cell* (2017) 169(2):361-71. doi: 10.1016/j.cell.2017.03.035

35. Causey DR, Kim JH, Stead DA, Martin SA, Devlin RH, Macqueen DJ. Proteomic comparison of selective breeding and growth hormone transgenesis in fish: Unique pathways to enhanced growth. *J Proteomics* (2019) 192:114-24. doi: 10.1016/j.jprot.2018.08.013

36. Ahyi AN, Quinton LJ, Jones MR, Ferrari JD, Pepper-Cunningham ZA, Mella JR, et al. Roles of STAT3 in protein secretion pathways during the acute-phase response. *Infect Immun* (2013) 81(5):1644-53. doi: 10.1128/IAI.01332-12

37. He G, Karin M. NF-κB and STAT3 - key players in liver inflammation and cancer. *Cell Res* (2011) 21(1):159-68. doi: 10.1038/cr.2010.183

38. Herzig S, Shaw RJ. AMPK: guardian of metabolism and mitochondrial homeostasis. *Nat Rev Mol Cell Biol* (2018) 19(2):121-35. doi: 10.1038/nrm.2017.95

39. Lau HH, Ng NH, Loo LSW, Jasmen JB, Teo AK. The molecular functions of hepatocyte nuclear factors - in and beyond the liver. *J Hepatol* (2018) 68(5):1033-48. doi: 10.1016/j.jhep.2017.11.026

40. Kawaguchi T, Osatomi K, Yamashita H, Kabashima T, Uyeda K. Mechanism for fatty acid "sparing" effect on glucose-induced transcription: regulation of carbohydrate-responsive element-binding protein by AMP-activated protein kinase. *J Biol Chem* (2002) 277(6):3829-35. doi: 10.1074/jbc.M107895200

41. Tikhanovich I, Cox J, Weinman SA. Forkhead box class O transcription factors in liver function and disease. *J Gastroenterol Hepatol* (2013) 28 Suppl 1(0):125-31. doi: 10.1111/jgh.12021

42. Mueller KM, Themanns M, Friedbichler K, Kornfeld JW, Esterbauer H, Tuckermann JP, et al. Hepatic growth hormone and glucocorticoid receptor signaling in body growth, steatosis and metabolic liver cancer development. *Mol Cell Endocrinol* (2012) 361(1-2):1-11. doi: 10.1016/j.mce.2012.03.026

43. Donne R, Saroul-Aïnama M, Cordier P, Celton-Morizur S, Desdouets C. Polyploidy in liver development, homeostasis and disease. *Nat Rev Gastroenterol Hepatol* (2020) 17(7):391-405. doi: 10.1038/s41575-020-0284-x

44. Barrett MT, Pritchard D, Palanca-Wessels C, Anderson J, Reid BJ, Rabinovitch PS. Molecular phenotype of spontaneously arising 4N (G2-tetraploid) intermediates of neoplastic progression in barrett's esophagus. *Cancer Res* (2003) 63(14):4211-7.

45. D'Agostino M, Tornillo G, Caporaso MG, Barone MV, Ghigo E, Bonatti S, et al. Ligand of numb proteins LNX1P80 and LNX2 interact with the human glycoprotein CD8α and promote its ubiquitylation and endocytosis. *J Cell Sci* (2011) 124(Pt 21):3545-56. doi: 10.1242/jcs.081224

46. Parker J, Guslund NC, Jentoft S, Roth O. Characterization of pipefish immune cell populations through single-cell transcriptomics. *Front Immunol* (2022) 13:820152. doi: 10.3389/fimmu.2022.820152

47. Bradley LM, Douglass MF, Chatterjee D, Akira S, Baaten BJ. Matrix metalloprotease 9 mediates neutrophil migration into the airways in response to influenza virus-induced toll-like receptor signaling. *PloS Pathog* (2012) 8(4):e1002641. doi: 10.1371/journal.ppat.1002641

48. Takizawa F, Magadan S, Parra D, Xu Z, Korytář T, Boudinot P, et al. Novel teleost CD4-bearing cell populations provide insights into the evolutionary origins and primordial roles of CD4+ lymphocytes and CD4+ macrophages. *J Immunol* (2016) 196(11):4522-35. doi: 10.4049/jimmunol.1600222

49. Glinos DA, Soskic B, Williams C, Kennedy A, Jostins L, Sansom DM, et al. Genomic profiling of T-cell activation suggests increased sensitivity of memory T cells to CD28 costimulation. *Genes Immun* (2020) 21(6-8):390-408. doi: 10.1038/s41435-020-00118-0

50. Mou D, Espinosa J, Lo DJ, Kirk AD. CD28 negative T cells: is their loss our gain? *Am J Transplant* (2014) 14(11):2460-6. doi: 10.1111/ajt.12937

51. Zhao X, Shan Q, Xue HH. TCF1 in T cell immunity: a broadened frontier. *Nat Rev Immunol* (2022) 22(3):147-57. doi: 10.1038/s41577-021-00563-6

52. Szabo PA, Levitin HM, Miron M, Snyder ME, Senda T, Yuan J, et al. Single-cell transcriptomics of human T cells reveals tissue and activation signatures in health and disease. *Nat Commun* (2019) 10(1):4706. doi: 10.1038/s41467-019-12464-3

53. Feng X, Wang H, Takata H, Day TJ, Willen J, Hu H. Transcription factor Foxp1 exerts essential cell-intrinsic regulation of the quiescence of naive T cells. *Nat Immunol* (2011) 12(6):544-50. doi: 10.1038/nli.2034

54. Schwartzberg PL, Mueller KL, Qi H, Cannons JL. SLAM receptors and SAP influence lymphocyte interactions, development and function. *Nat Rev Immunol* (2009) 9(1):39-46. doi: 10.1038/nri2456

55. Walling BL, Kim M. LFA-1 in T cell migration and differentiation. *Front Immunol* (2018) 9:952. doi: 10.3389/fimmu.2018.00952

56. Bartsch LM, Damasio MPS, Subudhi S, Drescher HK. Tissue-resident memory T cells in the liver unique characteristics of local specialists. *Cells* (2020) 9(11):2457. doi: 10.3390/cells9112457

57. Kang SM, Tsang W, Doll S, Scherle P, Ko HS, Tran AC, et al. Induction of the POU domain transcription factor Oct-2 during T-cell activation by cognate antigen. *Mol Cell Biol* (1992) 12(7):3149-54. doi: 10.1128/mcb.12.7.3149-3154.1992

58. Burns GF, Triglia T, Werkmeister JA, Begley CG, Boyd AW. TLIISA1, a human T lineage-specific activation antigen involved in the differentiation of cytotoxic T lymphocytes and anomalous killer cells from their precursors. *J Exp Med* (1985) 161(5):1063-78. doi: 10.1084/jem.161.5.1063

59. Takeuchi A, Badr Mel S, Miyauchi K, Ishihara C, Onishi R, Guo Z, et al. CRTAM determines the CD4+ cytotoxic T lymphocyte lineage. *J Exp Med* (2016) 213(1):123-38. doi: 10.1084/jem.20150519

60. Melichar HJ, Narayan K, Der SD, Hiraoka Y, Gardiol N, Jeannet G, et al. Regulation of gammadelta versus alphabeta T lymphocyte differentiation by the transcription factor SOX13. *Science* (2007) 315(5809):230-3. doi: 10.1126/science.1135344

61. Karsunky H, Merad M, Cozzio A, Weissman IL, Manz MG. Flt3 ligand regulates dendritic cell development from Flt3+ lymphoid and myeloid-committed
progenitors to Flt3+ dendritic cells *in vivo*. *J Exp Med* (2003) 198(2):305-13. doi: 10.1084/jem.20030323

62. Leone DA, Peschel A, Brown M, Schachner H, Ball MJ, Gyuraszova M, et al. Surface LAMP-2 is an endocytic receptor that diverts antigen internalized by human dendritic cells into highly immunogenic exosomes. *J Immunol* (2017) 199 (2):531-46. doi: 10.4049/jimmunol.1601263

63. Rocha-Perugini V, Martínez Del Hoyo G, González-Granado JM, Ramírez-Huesca M, Zorita V, Rubinstein E, et al. CD9 regulates major histocompatibility complex class II trafficking in monocyte-derived dendritic cells. *Mol Cell Biol* (2017) 37(15):e00202-17. doi: 10.1128/MCB.00202-17

64. Bunin A, Sisirak V, Ghosh HS, Grajkowska LT, Hou ZE, Miron M, et al. Protein tyrosine phosphatase PTPRS is an inhibitory receptor on human and murine plasmacytoid dendritic cells. *Immunity* (2015) 43(2):277-88. doi: 10.1016/j.immuni.2015.07.009

65. Worah K, Mathan TS, Vu Manh TP, Keerthikumar S, Schreibelt G, Tel J, et al. Proteomics of human dendritic cell subsets reveals subset-specific surface markers and differential inflammasome function. *Cell Rep* (2016) 16(11):2953-66. doi: 10.1016/j.celrep.2016.08.023

66. Grajkowska LT, Ceribelli M, Lau CM, Warren ME, Tiniakou I, Nakandakari Higa S, et al. Isoform-specific expression and feedback regulation of e protein TCF4 control dendritic cell lineage specification. *Immunity* (2017) 46(1):65-77. doi: 10.1016/j.immuni.2016.11.006

67. Deenick EK, Pelham SJ, Kane A, Ma CS. Signal transducer and activator of transcription 3 control of human T and b cell responses. *Front Immunol* (2018) 9:168. doi: 10.3389/fimmu.2018.00168

68. van Kooyk Y, Geijenbeek TB. DC-SIGN: escape mechanism for pathogens. *Nat Rev Immunol* (2003) 3(9):697-709. doi: 10.1038/nri1182

69. Young KM, Russell S, Smith M, Huber P, Ostland VE, Brooks AS, et al. Bacterial-binding activity and plasma concentration of ladderlectin in rainbow trout (*Oncorhynchus mykiss*). *Fish Shellfish Immunol* (2007) 23(2):305-15. doi: 10.1016/j.fsi.2006.10.014

70. Royet J, Dziarski R. Peptidoglycan recognition proteins: pleiotropic sensors and effectors of antimicrobial defences. *Nat Rev Microbiol* (2007) 5(4):264-77. doi: 10.1038/nrmicro1620

71. Watanabe Y, Tateno H, Nakamura-Tsuruta S, Kominami J, Hirabayashi J, Nakamura O, et al. The function of rhamnose-binding lectin in innate immunity by restricted binding to Gb3. *Dev Comp Immunol* (2009) 33(2):187-97. doi: 10.1016/j.dci.2008.08.008

72. Wernimont SA, Wiemer AJ, Bennin DA, Monkley SJ, Ludwig T, Critchley DR, et al. Contact-dependent T cell activation and T cell stopping require talin1. *J Immunol* (2011) 187(12):6256-67. doi: 10.4049/jimmunol.1102028

73. Nicolet BP, Guislain A, van Alphen FP, Gomez-Eerland R, Schumacher TN, van den Biggelaar M, et al. CD29 identifies IFN-γ-producing human CD8+ T cells with an increased cytotoxic potential. *Proc Natl Acad Sci U.S.A.* (2020) 117 (12):6686-96. doi: 10.1073/pnas.1913940117

74. Kremer KN, Dinkel BA, Sterner RM, Osborne DG, Jevremovic D, Hedin KE. TCR-CXCR4 signaling stabilizes cytokine mRNA transcripts via a PREX1-Rac1 pathway: implications for CTCL. *Blood* (2017) 130(8):982-94. doi: 10.1182/blood-2017-03-770982

75. Perišić Nanut M, Sabotič J, Jewett A, Kos J. Cysteine cathepsins as regulators of the cytotoxicity of NK and T cells. *Front Immunol* (2014) 5:616. doi: 10.3389/fimmu.2014.00616

76. Darmoise A, Maschmeyer P, Winau F. The immunological functions of saposins. *Adv Immunol* (2010) 105:25-62. doi: 10.1016/S0065-2776(10)05002-9

77. Ghosh S, Chackerian AA, Parker CM, Ballantyne CM, Behar SM. The LFA-1 adhesion molecule is required for protective immunity during pulmonary mycobacterium tuberculosis infection. *J Immunol* (2006) 176(8):4914-22. doi: 10.4049/jimmunol.176.8.4914

78. Liu W, Yan M, Liu Y, Wang R, Li C, Deng C, et al. Olfactomedin 4 down-regulates innate immunity against helicobacter pylori infection. *Proc Natl Acad Sci USA* (2010) 107(24):11056-61. doi: 10.1073/pnas.1001269107

79. Szulc-Dąbrowska L, Bossowska-Nowicka M, Struzik J, Toka FN. Cathepsins in bacteria-macrophage interaction: Defenders or victims of circumstance? *Front Cell Infect Microbiol* (2020) 10:601072. doi: 10.3389/fimbci.2020.601072

80. Sun B, van Dissel D, Mo I, Boysen P, Haslene-Hox H, Lund H. Identification of novel biomarkers of inflammation in Atlantic salmon (*Salmo salar* l.) by a plasma proteomic approach. *Dev Comp Immunol* (2022) 127:104268. doi: 10.1016/j.dci.2021.104268

81. Cunnick JM, Schmidhuber S, Chen G, Yu M, Yi SJ, Cho YJ, et al. Bcr and abr cooperate in negatively regulating acute inflammatory responses. *Mol Cell Biol* (2009) 29(21):5742-50. doi: 10.1128/MCB.00357-09

82. Cooper AM, Khader SA. IL-12p40: an inherently agonistic cytokine. *Trends Immunol* (2007) 28(1):33-8. doi: 10.1016/j.it.2006.11.002

83. Khader SA, Partida-Sanchez S, Bell G, Jelley-Gibbs DM, Swain S, Pearl JE, et al. Interleukin 12p40 is required for dendritic cell migration and T cell priming after mycobacterium tuberculosis infection. *J Exp Med* (2006) 203(7):1805-15. doi: 10.1084/jem.20052545

84. Russell S, Hayes MA, Simko E, Lumsden JS. Plasma proteomic analysis of the acute phase response of rainbow trout (Oncorhynchus mykiss) to intraperitoneal inflammation and LPS injection. *Dev Comp Immunol* (2006) 30 (4):393-406. doi: 10.1016/j.dci.2005.06.002

85. Hanawa-Suetsugu K, Itoh Y, Ab Fatah M, Nishimura T, Takemura K, Takeshita K, et al. Phagocytosis is mediated by two-dimensional assemblies of the f-BAR protein GAS7. *Nat Commun* (2019) 10(1):4763. doi: 10.1038/s41467-019-12738-w

86. Shi G, Abbott KN, Wu W, Salter RD, Keyel PA. Dnase1L3 regulates inflammasome-dependent cytokine secretion. *Front Immunol* (2017) 8:522. doi: 10.3389/fimmu.2017.00522

87. Jones EL, Wee JL, Demaria MC, Blakeley J, Ho PK, Vega-Ramos J, et al. Dendritic cell migration and antigen presentation are coordinated by the opposing functions of the tetraspanins CD82 and CD37. *J Immunol* (2016) 196(3):978-87. doi: 10.4049/jimmunol.1500357

88. Williams JW, Tjota MY, Clay BS, Vander Lugt B, Bandukwala HS, Hrusch CL, et al. Transcription factor IRF4 drives dendritic cells to promote Th2 differentiation. *Nat Commun* (2013) 4:2990. doi: 10.1038/ncomms3990

89. Bakke FK, Monte MM, Stead DA, Causey DR, Douglas A, Macqueen DJ, et al. Plasma proteome responses in salmonid fish following immunization. *Front Immunol* (2020) 11:581070. doi: 10.3389/fimmu.2020.581070

90. Quinton LJ, Blahna MT, Jones MR, Allen E, Ferrari JD, Hilliard KL, et al. Hepatocyte-specific mutation of both NF-κB RelA and STAT3 abrogates the acute phase response in mice. *J Clin Invest* (2012) 122(5):1758-63. doi: 10.1172/JCI59408

91. Brasch K. Endopolyploidy in vertebrate liver: an evolutionary perspective. *Cell Biol Int Rep* (1980) 4(2):217-26. doi: 10.1016/0309-1651(80)90077-6

92. Bakken TE, Hodge RD, Miller JA, Yao Z, Nguyen TN, Aevermann B, et al. Single-nucleus and single-cell transcriptomes compared in matched cortical cell types. *PloS One* (2018) 13(12):e0209648. doi: 10.1371/journal.pone.0209648

93. Wan F, Hu CB, Ma JX, Gao K, Xiang LX, Shao JZ. Characterization of γδ T cells from zebrafish provides insights into their important role in adaptive humoral immunity. *Front Immunol* (2017) 7:675. doi: 10.3389/fimmu.2016.00675

94. Hammerich L, Tacke F. Role of gamma-delta T cells in liver inflammation and fibrosis. *World J Gastrointest Pathophysiol* (2014) 5(2):107-13. doi: 10.4291/wjgp.v5.i2.107

95. Hamada S, Umemura M, Shiono T, Tanaka K, Yahagi A, Begum MD, et al. IL-17A produced by gammadelta T cells plays a critical role in innate immunity against listeria monocytogenes infection in the liver. *J Immunol* (2008) 181 (5):3456-63. doi: 10.4049/jimmunol.181.5.3456

96. Chan JT, Kadri S, Köllner B, Rebl A, Korytář T. RNA-Seq of single fish cells - seeking out the leukocytes mediating immunity in teleost fishes. *Front Immunol* (2022) 13:798712. doi: 10.3389/fimmu.2022.798712

97. Bassity E, Clark TG. Functional identification of dendritic cells in the teleost model, rainbow trout (*Oncorhynchus mykiss*). *PloS One* (2012) 7(3):e33196. doi: 10.1371/journal.pone.0033196

98. Guslund NC, Solbakken MH, Brieuc MS, Jentoft S, Jakobsen KS, Qiao SW. Single-cell transcriptome profiling of immune cell repertoire of the Atlantic cod which naturally lacks the major histocompatibility class II system. *Front Immunol* (2020) 11:559555. doi: 10.3389/fimmu.2020.559555

99. Seternes T, Bøgwald J, Dalmo RA. Scavenger endothelial cells of fish, a review. *J Fish Dis* (2021) 44(9):1385-97. doi: 10.1111/jfd.13396

100. Speilberg L, Evensen O, Nafstad P. Liver of juvenile Atlantic salmon, *Salmo salar* l.: a light, transmission, and scanning electron microscopic study, with special reference to the sinusoid. *Anat Rec* (1994) 240(3):291-307. doi: 10.1002/ar.1092400302

101. Cheng D, Morsch M, Shami GJ, Chung RS, Braet F. Observation and characterisation of macrophages in zebrafish liver. *Micron* (2020) 132:102851. doi: 10.1016/j.micron.2020.102851

102. Haltiwanger RS, Lehrman MA, Eckhardt AE, Hill RL. The distribution and localization of the fucose-binding lectin in rat tissues and the identification of a high affinity form of the mannose/N-acetylglucosamine-binding lectin in rat liver. *J Biol Chem* (1986) 261(16):7433-9.

103. Yang CY, Chen JB, Tsai TF, Tsai YC, Tsai CY, Liang PH, et al. CLEC4F is an inducible c-type lectin in F4/80-positive cells and is involved in alpha-galactosylceramide presentation in liver. *PloS One* (2013) 8(6):e65070. doi: 10.1371/journal.pone.0065070

104. Perdiguero P, Morel E, Tafalla C. Diversity of rainbow trout blood b cells revealed by single cell RNA sequencing. *Biol (Basel)* (2021) 10(6):511. doi: 10.3390/biology10060511
105. Macqueen DJ, Johnston IA. A well-constrained estimate for the timing of the salmonid whole genome duplication reveals major decoupling from species diversification. *Proc Biol Sci* (2014) 281(1778):20132881. doi: 10.1098/rspb.2013.2881

106. Gundappa MK, To TH, Grønvold L, Martin SA, Lien S, Geist J, et al. Genome-wide reconstruction of rediploidization following autopolyploidization

across one hundred million years of salmonid evolution. *Mol Biol Evol* (2022) 39 (1):msab310. doi: 10.1093/molbev/msab310

107. Gillard GB, Grønvold L, Røsæg LL, Holen MM, Monsen Ø, Koop BF, et al. Comparative regulomics supports pervasive selection on gene dosage following whole genome duplication. *Genome Biol* (2021) 22(1):103. doi: 10.1186/s13059-021-02323-0