Transcriptomic characterization of transitioning cell types in the skin of Atlantic salmon

Abstract

Background

The skin maintains the body’s integrity and serves as the first line of defence against pathogens, stressors and mechanical injuries. Despite the global significance of salmon in aquaculture, how the transcriptomic profile of cells varies during wound healing remains unexplored. Teleost’s skin contains adult pluripotent cells that differentiate into various tissues, including bone, cartilage, tendon, ligament, adipose, dermis, muscle and connective tissue within the skin. These cells are pivotal for preserving the integrity of skin tissue throughout an organism’s lifespan and actively participate in the wound healing processes. In this study, we characterize the transcriptomic profiles of putative mesenchymal stromal cells (fibroblast-like adult stem cells) in healthy Atlantic salmon tissue and during the wound healing process.

Results

Single-nucleus sequencing and spatial transcriptomics were used to detect transcriptomic changes occurring during wound healing that are commonly associated with mesenchymal stromal cells. We followed the transcriptomic activity of these cells during an in vivo wound healing time course study showing that these cells become more transcriptionally active during the remodelling stage of wound healing. The changes detected give insights into the potential differentiation pathways leading to osteogenic and fibroblast lineages in the skin of Atlantic salmon.

Conclusions

We chart the transcriptomic activity of subclusters of putative differentiating stromal cells during the process of wound healing for the first time, revealing different spatial niches of the various putative MSC subclusters, and setting the stage for further investigation of the manipulation of transitioning cell types to improve fish health.

Background

Atlantic salmon is one of the most important aquaculture species worldwide, accounting for 32.6% of the global cultured marine finfish trade in 2020 [1], and with demand steadily increasing. However, sustainable expansion of the industry has been hindered by a range of disease and welfare issues. In particular, skin health issues, such as injuries caused by delousing operations, sea louse infections and bacterial winter wounds, are on the rise and represent one of the main current challenges in Atlantic salmon aquaculture [2]. These issues not only have major impacts on economic performance but also raise concerns about animal welfare.

The skin is a vital tissue for overall animal health, welfare and robustness. In Atlantic salmon (Salmo salar, Linnaeus 1758), the skin is a complex organ that serves various functions, including protection against pathogens, mechanical stress and homeostasis regulation. A key characteristic of skin is its remarkable ability to heal, which is enabled by specialized cells working in concert to repair and restore tissue integrity. Wound healing is a constant, resilient response that allows for the maintenance of homeostasis and scarless repair in adult fish [3, 4]. Among the diverse cell types in fish skin, including different populations of keratinocytes, immune cells, fibroblasts, iridocytes and others, mesenchymal stromal cells (MSCs) are essential for maintaining various tissue structures [5]. These cell types coordinate activities during the process of wound healing to maintain tissue integrity and homeostasis. Understanding cellular composition and the role of differentiating cell types in this barrier tissue is essential for understanding its function.

Of all the cell types in the skin of Atlantic salmon, MSCs are one of the least characterized. These are adult pluripotent stem cells with fibroblast-like morphology and phenotype [6,7,8]. They exhibit immunomodulatory properties: influencing the immune response at the wound site by regulating immune cell activity, reducing excessive inflammation and fostering an environment conducive to tissue repair [9]. These heterogeneous cells can not only self-renew but are also a highly transitioning cell type, being able to differentiate into adipocytes, fibroblasts, osteoblasts and myocytes to maintain tissue homeostasis and repair (wound healing) [10,11,12]. Recent single-cell sequencing data has identified clusters of cells with stem cell properties that fit this description, highlighting the fact that they are present in fish [13], and can be identified using transcriptomics [14, 15]. In salmonids, their differentiation into osteoblasts and adipocytes has also been demonstrated through primary cell culture of visceral fat [16]. Additionally, immunohistological analysis of wound healing has hinted at this process, with spindle-shaped proliferating cells observed to be particularly active on day 14 post-wounding (DPW) [17]. This indicates a role of these cells in homeostasis, with heightened activity during the tissue remodelling phase in the skin of Atlantic salmon. A better understanding of the ability of these cells to differentiate into a myriad of cell types and their role in maintaining tissue physiology and homeostasis can provide insights into the mechanisms underlying skin repair and pave the way for therapies aimed at improving wound healing and skin barrier function [6, 9, 18].

One caveat of this work that must be considered is that when using single-cell technologies to characterize cell types, we do so based on a priori markers for those cell types [19, 20]. It is important to note that the classification of cell types is inherently based on model species such as humans [21] and, to a more limited extent, zebrafish [22, 23]. Most functional annotations have been derived from mammals, which can lead to the assumption that mammalian data are directly applicable to fish. Clusters are typically assigned putative cell types based on available markers, which currently serve as the primary criterion for defining cell types in non-human, non-mouse systems such as salmon. However, these criteria are still evolving [24,25,26,27]. This challenge becomes more pronounced when dealing with nuanced cell types, particularly in non-homologous tissues like skin. Unlike amniote skin, such as human skin, teleost skin is non-cornified, with living cells in its outer layers that actively interact with the environment [28]. It is often assumed that function is conserved, but these assumptions have not been well tested because single-cell technologies are relatively new. As such, some clusters may be more clearly defined than others, and the classifications may evolve as further research clarifies the functional and molecular characteristics of these cells in salmonids.

Even with this caveat, the development of single-cell technologies, such as single-cell RNA sequencing (scRNA-seq), has revolutionized our ability to study individual cells. While the extraction of individual cells can be challenging in certain tissues, such as skin and fins [29], leading to biases in cell-type composition, these challenges can be addressed through single-nuclei RNA sequencing (snRNA-seq), allowing the dissociation of nuclei from frozen tissue without losing cellular diversity [24, 30]. While single-cell/nuclear technologies reveal the transcriptome of individual cells, the spatial location of the cells is lost during dissociation. Combining snRNA-seq with technologies that capture the positional context of gene expression, such as spatial transcriptomics [31], allows the contextualization of single-cell transcriptomic data. The power of combining these technologies to understand skin biology has already been exploited in humans, for example, to unravel the role of different skin cell types during wound healing processes [32]. Although both single-nucleus RNA-seq and spatial transcriptomics have been used separately to study skin in Atlantic salmon [26, 33], cell expression patterns have not yet been investigated in their tissue context, and this is crucial to understanding of biology of processes such as wound healing.

This study aimed to advance our understanding of cell types in the skin of Atlantic salmon and chart the transcriptomic activity of transitioning cell types during wound healing. Using snRNA-seq of nuclei isolated from flank skin (containing scales) and fins (lacking scales), we aimed to (i) identify differences in cellular composition between flank skin and fins, (ii) identify putative transitioning cells in Atlantic salmon skin, (iii) infer their heterogeneity and trajectory dynamics and provide marker genes for these cells and, finally, (iv) employ spatial transcriptomics to chart their transcriptomic activity during wound healing.

Results

Six snRNA-Seq Atlantic salmon flank skin and fin libraries were sequenced (two flank skin samples, two pelvic fins, one pectoral fin and one dorsal fin). Summary statistics for all 6 libraries can be found in Additional file 1: Table S1. After quality control, 27,989 nuclei remained, averaging 4665 nuclei per sample, with an average of 2813 (± 791) unique transcripts (unique molecular identifiers—UMIs) and 1566 (± 334) genes per nucleus. The number of nuclei detected per sample after filtering stages can be found in Additional file 1: Table S2.

Unsupervised clustering of 27,989 nuclei based on their transcriptomes revealed 18 distinct cell clusters (Fig. 1a), potentially representing 18 different cell types. To assign putative cell types to the identified clusters, we referred to lists of marker genes from previously published works on the same tissues [26], where it was also shown that cell types are relatively stable in salmon skin across species. A combination of transcriptomic clustering was used to classify cells. Differential gene expression analysis revealed markers for each major skin cell type Atlantic salmon skin cell types and their top 20 markers per cell type are presented in Additional file 1: Table S3a, with a review of these markers provided in Additional file 2: Table S3b.

Fig. 1

Comparison of cell types between flank skin and fins. a, b Uniform manifold approximation and projection (UMAP) plot showing the potential 18 cell types present across all flank skin and fin samples, respectively. Cells are coloured according to their cluster, and the putative cell types are provided in the legend. c Proportion of each cell type in the fin and skin normalized by dividing the number of cells of that cell type by the total number of cells. d Pearson rank correlation between the transcriptomes of each cell type in flank skin and fins

Clustering and identification of skin cell lineages were repeatable across individual samples and across flank skin and fin samples, with details found in Additional file 1: Table S4, which shows the amount of cells/percentage of cells per sample and per type of sample for all cell types, while Additional file 1: Table S5 provides the same information for putative mesenchymal stromal cell subtypes. The cell clusters broadly matched expectations based on previous work [17, 26]. Marker genes for candidate cell types and sub-populations are reported according to annotations provided by Ensembl when possible. The marker genes are referenced when they are known as a priori markers of that cell type.

Comparison of cellular composition between flank skin and fins

Two different sizes of fish were used in this study: 486 and 25 g. This size variation was introduced to capture cells from different life stages, as smaller fish are likely to have a higher turnover rate of skin cells. To ensure consistency, we sampled both flank skin and fin from fish of different sizes for this analysis. Cell types are generally preserved across tissues and taxa [34] including salmonids. While the flank skin (Fig. 1a) and fin (Fig. 1b) have broadly similar cell types, the proportion of the different cell types differ. Eighteen cell types were found in the flank skin, this being consistent with previous work [26]. The main three markers for each cell type are detailed in Additional file 2: Table S3b, with further markers shown in Additional file 1: Table S3a. The flank skin has a higher proportion of fibroblasts (cluster 7), skeletal muscle cells (cluster 12) and neural crest cells (cluster 13). In contrast, the fin has a higher proportion of intermediate keratinocytes (cluster 0), basal keratinocytes (cluster 3), two populations of fibroblasts (clusters 1 and 11) and neural cells (cluster 17) (Fig. 1c). Further details of cell type identity and proportion across sample and tissue type can be found in Additional file 1: Table S4.

Putative mesenchymal stromal cells

A particular cell type, previously characterized as fibroblasts in the skin of Atlantic salmon [26], emerged as a promising candidate for the mesenchymal stromal cell type. Markers associated with this cell type (Additional file 1: Table S3a) included classic fibroblast markers, such as fndc1.1 (ENSSSAG00000008223), antxr1.1 (ENSSSAG00000078628), itgbl1, (ENSSSAG00000044364), col6a3 (ENSSSAG00000042605) and aff2.1 (ENSSSAG00000057351) [21]. In addition, this cell type expressed markers associated with MSC activities: htra1b (ENSSSAG00000051286) which promotes fibroblast differentiation [35]; fbn2 (ENSSSAG00000057875) which regulates skeletal stem cell differentiation [36]; antxr1.1 (ENSSSAG00000078628) which is linked to stromal chondrogenesis [37]; and col6a3 (ENSSSAG00000042605) which is expressed during osteogenic differentiation [38]. Finally, this population showed an upregulation of wound-healing-like transcripts during infections with sea lice [26].

As another test of the potential for these cells to be MSCs, we looked at the expression of markers characterizing MSCs that have been previously identified in other species, specifically two paralogues of cd34 (ENSSSAG00000073342 (Fig. 2a), ENSSSAG00000074590) [39] (Fig. 2a), itga4 (ENSSSAG0000006416) [40] and itga5 (ENSSSAG00000098283) [41] (Additional file 2: Fig. S1a). These genes were also found to be expressed in fibroblast population 1 (cluster 1, Fig. 2a). Given its fibroblast like identity, infection response and expression of amniote stem cell markers, this population was identified as a putative MSC containing population. The cell population was subsetted out from the main data set and was reanalysed separately to explore potential heterogeneity within this cell type. Five subclusters of cells were identified within this putative MSC cell population (Fig. 2b), details of amount of cells/percentage of cells per sample found in Additional file 1: Table S5. Subcluster 2 included putative pure mesenchymal stromal cells with high expression of itga4, itga5 and cd34 (Fig. 2c).

Fig. 2

Identification and characterization of putative mesenchymal stromal cells. a Uniform manifold approximation and projection (UMAP) plot showing the expression of mammalian MSC markers, such as cd34, in a salmon skin and fin cell atlas. This panel displays expression levels for one cd34 paralogue; for other paralogues and stem markers refer to Additional file 2: Fig. S1a. b UMAP analysis illustrating distinct cell subclusters within the putative MSC-containing population. Each subcluster is colour-coded to represent different sub cell types. c Violin plots showing the expression levels of MSC markers across the various subclusters of the putative MSC populations. Each colour corresponds to one of the six samples analysed, ordered from left to right: fin1, fin2, skin1, fin3, fin4 and skin2. d Heatmap of the top 5 differentially expressed genes for each cellular subcluster, compared against the background of all other cell types. See Additional file 1: Table. S6 for the full list of marker genes

The remaining clusters were identified based on their marker genes (Fig. 2d, provided in more detail in Additional file 1: Table S6) as fibroblasts (subcluster 0; expression of collagen and keratin genes), undefined cells (subcluster 1; expression of nonspecific keratin genes), and two groups of differentiating MSCs (subclusters 3 and 4; based on the expression of fgfr4 (ENSSSAG00000017777, ENSSSAG00000020292) and wnt7 (ENSSSAG00000009090, ENSSSAG00000086070). Based on the expression of genes pth1r (ENSSSAG00000097651), ptprd (ENSSSAG00000096916, ENSSSAG00000111433), runx3 (ENSSSAG00000003444, ENSSSAG00000064246), fgfr4 (ENSSSAG00000017777, ENSSSAG00000020292), we classified cluster 3 as putative bone/muscle precursors. Based on the expression of genes wnt7bb (ENSSSAG00000009090, ENSSSAG00000086070), ebf1 (ENSSSAG00000070298, ENSSSAG00000079780), tafa5 (ENSSSAG00000088442), and myh10 (ENSSSAG00000101624), we classified 4 as putative adipocyte precursors. The presence of these subsets, which also expressed markers for cellular sub-types usually associated with MSCs, was taken as further proof of their putative stromal identity.

The dynamics of differentiating putative salmon mesenchymal stromal cells

UMAPs are useful visualization tools but are known to have limitations in faithfully capturing the dynamics and structure of high-dimensional data [42]. To further explore the transitions of cell types, we used PHATE as an alternative dimensional reduction method for the analysis of the putative MSC cluster. PHATE identified 7 different cell subclusters (Fig. 3a) with differential expression of marker genes (Fig. 3b), 2 more than that observed in UMAP. This extra clustering reflects the transition probabilities between cell subclusters or states, thus capturing more variability in the data in low-dimensional space while removing noise and retaining both the global and local structure [43]. The correspondence between the UMAP and PHATE clusters are shown in Additional file 2: Fig. S1b, where PHATE was run with UMAP parameters, with population 1 representing the putative pure mesenchymal stromal cells (Fig. 3b, Table 1). The top 20 marker genes for the PHATE subclustering were visualized using a UMAP to show the correspondence between the clusters across the methods (Additional file 2: Fig. S2, S3, S4, S5, S6 and S7). PHATE also revealed different transitioning states, represented by the spindly arms observed (populations 2, 4, 6 and, to a lesser extent, 3), with the length of the arm indicating the extent of differentiation/change in a particular cluster/subcluster of cells (Fig. 3a). These transitioning states were characterized based on their gene markers (Fig. 3b, Table 1 and were identified as putative fibroblasts (cluster 2 and 3), bone cell precursor cells (cluster 4), adipose precursor cells, muscle/bone precursor cells [6] and an undefined cell cluster (0) given its lack of biologically informative markers. This cluster was therefore characterized as “undefined” (Fig. 3b, see Table 1 for a list of key marker genes for each cluster Additional file 1: Table S7 a full list).

Table 1 The PHATE subclusters of the putative MSC-containing subpopulations, along with their characteristic marker genes and putative functions as reported in the literature. These observations provide insights into their transcriptomic identity. Collectively, the subclusters exhibit transcripts and characteristics typically associated with mesenchymal stromal cells

Fig. 3

Identification and characterization of mesenchymal stromal cells using PHATE embedding. a PHATE analysis showing the various states of the putative MSC containing cluster. b Heatmap showing the top marker genes for the seven identified sub-cellular states in putative MSC-containing cluster 1. More details of these markers are provided in Table 1

PHATE cluster 1 included, undifferentiated putative mesenchymal stromal cells (MSCs) based on their expression of the MSC markers cd34 [39] and itga5 [41]. This cluster also expressed klf2 (ENSSSAG00000006501), a known marker of human mesenchymal stem cells [49]. Consistent with their role as progenitors or bone cells and adipocytes [16], this cluster was characterized by multiple marker genes associated with osteogenesis and adipogenesis. For example, fbln2 (ENSSSAG00000008760) [51] is expressed by MSCs differentiating into bone cells, as is loxl1 (ENSSSAG00000004651) [52], which is involved in bone and adipocyte differentiation of MSCs (Table 1). This cluster, along with cluster 4 (putative adipocyte precursor cells), also highly expressed ebf1 (ENSSSAG00000070298, ENSSSAG00000079780), which is involved in adipocyte differentiation [46]. An additional marker gene for cluster 1 was cobll1 (ENSSSAG00000080500) [53], which is highly expressed in adipose-derived MSCs, and to a lesser extent in fibroblasts [50]. When these transcripts were mapped on to Seurat groups (Additional file 2: Fig. S2), most are expressed in cluster 2 “pure” MSC in Seurat.

PHATE cluster 2 was identified as fibroblasts by the expression of col12a1 (ENSSSAG00000076478, ENSSSAG00000070858), which was previously found to be a marker of fibroblasts in several salmonid species [26], as well as flrt2 (ENSSSAG00000118246, ENSSSAG00000116745) [54], a fibroblast growth factor. When these transcripts were mapped on to Seurat groups (Additional file 2: Fig. S3), most of them are expressed in cluster 0, the fibroblast progenitor in Seurat, with some expression in cluster 3, the bone/muscle progenitor. Cluster 3 in PHATE was marked as fibroblast with markers such as tenascinb (ENSSSAG00000078781) [57] and htra1b (ENSSSAG00000051286) [35]. When PHATE cluster 3’s top 20 makers were mapped on to Seurat clustering (Additional file 2: Fig. S4), most were expressed in cluster 0, the fibroblast progenitor in Seurat, and some in cluster 4, the adipocyte population in Seurat.

The identity of cluster 4 as a bone precursor cell type is suggested by its expression of many genes associated with bone development, including pth1r (ENSSSAG00000097651) [59], ptprd (ENSSSAG00000096916, ENSSSAG00000111433) [60], runx3 (ENSSSAG00000003444, ENSSSAG00000064246) [61, 70] and cbfb (ENSSSAG00000064481) [61]. Other marker genes supporting this identity include palld (ENSSSAG00000068098), a gene associated with osteoblast differentiation from MSCs [12], and fgfr4 (ENSSSAG00000017777), a marker for bone cells [58], particularly for osteoblasts in Atlantic salmon skin and fin samples [26]. When PHATE cluster 4’s top 20 makers were mapped on to Seurat clustering (Additional file 2: Fig. S5), most were expressed in cluster 3, the bone muscle progenitor in Seurat.

The identity of PHATE cluster 5 as an adipocyte precursor cell is suggested by its expression of multiple genes associated with adipogenesis, including wnt7bb (ENSSSAG00000009090, ENSSSAG00000086070), [62] (Fig. 4c), ebf1 (ENSSSAG00000079780, ENSSSAG00000070298) [46], myh10 (ENSSSAG00000101624) [63] and lpl (ENSSSAG00000072070) [65]. This cluster also expresses tafa5a (ENSSSAG00000109872, ENSSSAG00000088442), which is highly expressed in adipocytes [64]. This cluster is unlikely to contain fully differentiated adipocytes, however, given the expression of ttn (ENSSSAG00000095939), a key marker for muscle cells in salmonid skin samples [26]. Interestingly, sulf1 (ENSSSAG00000008807), which has been observed to be downregulated during adipogenesis [71], likely due to its role in driving osteogenesis [72], was also a marker for this cluster. This potentially reflects that this cluster is an early adipocyte precursor. When its top 20 makers were plotted on Seurat (Additional file 2: Fig. S6), most were expressed in cluster 4, the adipocyte population in Seurat.

Fig. 4

Differences in mesenchymal stromal cells between the flank skin and fin. a PHATE plot showing the skin cells. b PHATE plot showing the fin cells. c Violin plots of expression levels of three Wnt genes in the six PHATE subclusters (0 to 6)

Cluster 6 expressed markers associated with both bone and muscle differentiation. This cluster shared many marker genes with cluster 4 (bone precursor cells), including palld (ENSSSAG00000068098), fgfr4 (ENSSSAG00000017777, ENSSSAG00000020292), cbfb (ENSSSAG00000064481), runx3, ptprd (ENSSSAG00000096916, ENSSSAG00000111433), pth1r (ENSSSAG00000097651) and wnt5b (ENSSSAG00000113059) [58]. However, this cluster was also uniquely defined by multiple marker genes associated with muscle development: slit1 (ENSSSAG00000001168) [67], sox6 (ENSSSAG00000111310) [68], mmp23b (ENSSSAG00000077749) [69] and ttn (ENSSSAG00000095939) [26]. When its top 20 makers were plotted on Seurat UMAP (Additional file 2: Fig. S7). Most were expressed in Population 3, the bone/muscle progenitor in Seurat. Given that cluster 6 was not apparent in fin samples (Fig. 4b), we hypothesize about its identity in the section below. Collectively, the transitioning subclusters exhibit transcripts and characteristics typically associated with mesenchymal stromal cells. For a complete description of these transcripts and their putative functions, please refer to Table 1. Note that most of the functions associated with these cells have been validated in other organisms, primarily mammals but also in zebrafish; relevant references can be found in the table.

Differences in mesenchymal stromal cells in the fin and flank skin

The separate plotting of the putative MSC containing sub-populations in flank skin and fins indicated differences between the two tissues. The flank skin exhibited a more complex structure compared to the fin, which is biologically consistent given the skin’s need to regenerate a wider variety of cell types, such as scales that are absent in fins (Fig. 4a and b). The fin showed more quiescent cell subpopulations, particularly in clusters associated with bone (cluster 4) and fibroblasts (clusters 2 and 3). Bone/muscle precursors (cluster 6) were absent from fin.

The data also highlight the cluster-specific expression of genes. For instance, wnt7bb (ENSSSAG00000009090) and wnt7bb.1 (ENSSSAG00000086070), both key regulators in the pleotropic Wnt signalling pathway [22] that governs cell fate determination, are predominantly expressed in cluster 5, which is shared between flank skin and fin. Conversely, wnt5 (ENSSSAG00000113059), known to play a role in MSC differentiation and chondrocyte proliferation during bone development [73], is exclusively expressed in cluster 6, which is more prominent in the flank skin (Fig. 4c). This population expresses a mixture of bone and muscle transcripts (Table 1). Wnt5 is a pleiotropic gene involved in many developmental processes, including potential scale formation [74]. It could be hypothesized that this is a cell subtype involved in scale formation, but further investigation will be needed to confirm this.

Additionally, fibroblast trajectories appear more active in the skin compared to the fin. This is consistent with Seurat clustering on these two tissues separately (Additional file 2: Fig. S8 and Additional file 1: Tables S8 and S9). This suggests that MSCs in the flank skin actively differentiate and play a crucial role in maintaining tissue homeostasis and regeneration under normal conditions. In contrast, the MSC population in the fins shows lower levels of proliferation and differentiation but remains responsive to external injuries or factors compromising fin integrity [75, 76].

“Pure” mesenchymal stromal cells during wound healing

To investigate the activity of putative MSCs during wound healing, we used tissue sections of mechanically induced wounds at 2 and 14 days post wounding (DPW) (Fig. 5a and b, respectively). These time points represent an early phase of wound healing (2 DPW), characterized by the formation of new epidermis, inflammation and tissue degeneration (Fig. 5a), and a late phase of wound healing (14 DPW) characterized by the replacement of damaged fibres with granulation tissue, a mix of new connective tissue and small blood vessels (Fig. 5b) [17]. Using spatial transcriptomics, the expression of genes for the putative pure MSC cluster, as defined by snRNA-seq for Seurat subcluster 2 (Fig. 2), was assessed during wound healing as a proxy for MSC activity and proliferation. We chose the top 20 markers expressed by this population (Additional file 1: Table S6) and charted their transcriptomic activity during wound healing (Fig. 5). At 2 DPW (c), these markers exhibited a diffuse distribution throughout the entire tissue section, with mean expression value of 0.085 (replicate 1 at 0.09 and replicate 2 at 0.08) across all markers. By 14 DPW (d), MSC activity and expression were distinctly concentrated in the centre of the granulation tissue, displaying higher levels of expression, with average expression value of 0.28 (replicate 1 at 0.34 and replicate 2 at 0.22). For direct comparison with PHATE, the average expression of the top 20 most highly expressed genes (Additional file 1: Table S7) was used for plotting on the spatial slides in Fig. 5. The same trend was observed. At 2 DPW (e), with diffuse distribution throughout the entire tissue section, with average relative expression value of 0.11 (replicate 1 at 0.07 and replicate 2 at 0.15). At 14 DPW (f), MSC activity and expression were distinctly concentrated in the centre of the granulation tissue, displaying higher levels of expression, with average expression value of 0.32 (replicate 1 at 0.4 and replicate 2 at 0.23). Generally, all the “pure” MSC precursor subclusters exhibited lower average expression at 2 dpw (a, c, e and g) than at 14 DPW (b, d, f and h). The contribution of each marker identified by both methods to the overall expression levels is provided in Additional file 1: Table S10.

Fig. 5

Wound healing in Atlantic salmon. a Histological micrographs depicting incisional wound at 2 days post-wounding (DPW), with Alcian blue and periodic acid-Schiff (PAS) staining to facilitate viewing of the incisional wound during the inflammation stage b Histological micrographs depicting incisional wound at 14 DPW, with Movat staining to facilitate viewing of the granulation tissue during the remodelling stage. c-f 10x Visium Spatial Transcriptomics Slides. Expression of putative pure MSC population 2 top 20 transcript from Seurat snRNA-seq data (Additional file 1: Table S10) at 2 DPW and 14 DPW (c, d). Expression of putative pure MSC population 1 top 20 transcripts from PHATE snRNA-seq analysis (Additional file 1: Table S10) at 2 DPW and 14 DPW (e, f). Wound bed (Wb), epidermis (Epi), dense connective tissue (Dct), skeletal muscle fibres (Mu), damaged white muscle fibres (Mu*), myosepta (Myo) and newly formed epithelial tissue (“Neo Epi”). Scales (Sc), adipose tissue (Adi), polymorphonucleated inflammatory cells (InF), granulation tissue (Gt), blood vessel formation (Bv) and fibril formation (Ff). Scale bars: 500 µm (a–f). The colour of the scale in c-f indicates the expression of transcripts mapped on the slide from low (blue) to high (red)

Transcript activity of putative mesenchymal stromal cell-associated subclusters during wound healing

To better explore the transcriptomic dynamics of putative MSCs associated with subclusters at a finer scale during wound healing, the MSC subgroups identified during the PHATE embedding analysis (Fig. 3) were plotted individually (Fig. 6). For bone precursors, a comparable trend was observed to both “pure” MSC markers taken from Seurat and PHATE (Fig. 5): low relative average expression at 2 DPW (a), with mean of 0.38 (replicate 1 at 0.26 and replicate 2 at 0.51), and higher relative average expression at 14 DPW (b), with mean of 0.51 (replicate 1 at 0.64 and replicate 2 at 0.39). Visually, there was less transcriptomic activity in the granulation tissue, but higher expression in the dense connective tissue at this time point. Adipocyte precursors showed random distribution across the section at 2 DPW (c), with relative average expressions of 0.14 (replicate 1 at 0.08 and replicate 2 at 0.20). At 14 DPW (d), the average expression increased to 0.22 (replicate 1 at 0.27 and replicate 2 at 0.16). Visually, there was higher expression in dense connective tissue, scale pockets, the apoptogenic layer and the wound bed on day 14 DPW compared to 2 DPW. For bone/muscle precursors, aside from one dot in the scale pockets at 2 DPW (e), there was low relative average expression of 0.016 (replicate 1 at 0.008 and replicate 2 at 0.02). In contrast, at 14 DPW (f), the relative average expression increased to 0.045 (replicate 1 at 0.05 and replicate 2 at 0.04), with expression being quite widespread in the epidermis, dense connective tissue and granulation tissue. For the two fibroblast populations, the trend is similar to that observed for MSCs. Fibroblast 2 (g) at 2 DPW had relative average expressions of 0.18 (replicate 1 at 0.09 and replicate 2 at 0.19) and showed a dramatic visual increase in expression in the wound bed at 14 DPW (h), where the relative expression values were 0.51 (replicate 1 at 0.63 and replicate 2 at 0.38). A similar trend was observed for fibroblast 3, with relative average expressions at 2 DPW (i) of 0.63 (replicate 1 at 0.41 and replicate 2 at 0.85), compared to 14 DPW (j), where the values were 1.3 (replicate 1 at 1.59 and replicate 2 at 1.02) (marker contributions for each subset shown in Additional file 1: Table S11).

Fig. 6

Charting the transcriptomic activity of different putative MSC-associated subclusters (10x Visium Spatial Transcriptomic slides). a-h Expression of MSC subtype-specific transcripts, taken from the Phate analysis (average expression of the top 20 markers of each subtype, Supplementary Table 11). Expression of bone precursors on day 2 and 14 DPW (a, b). Expression of adipocyte precursors on day 2 and 14 DPW (c, d). Expression of bone/muscle precursors on day 2 and 14 DPW (e, f). Expression of fibroblast 2 on day 2 and 14 DPW (g, h). Expression of fibroblast 3 on day 2 and 14 DPW (i, j). 2 DPW represents inflammation stage (right column) and 14 DPW represents remodelling stage (left column) of wound healing. Scale bars: 500 um (a-j). The colour of the scales in a-j indicates the expression of transcripts mapped to the slide from low (blue) to high (red)

Discussion

This study provides a detailed transcriptomic analysis of the composition of Atlantic salmon skin from the body and fins, focusing on characterizing potential transitioning cell types that exhibit various characteristics indicative of a potential MSC identity. Although some caution is warranted as the number of examined samples is low due to the high costs of the technologies used, the implications of the results are highly relevant to Atlantic salmon aquaculture, providing a basis for further exploration and validation of the role of these putative MSCs in wound healing. Putative MSCs were identified among a previously identified population of fibroblasts that react to disease in Atlantic salmon [26]. It was found that they expressed MSC markers discovered in other species [41, 51, 77]. The subclustering revealed cell subclusters that expressed markers consistent with cell types associated with MSCs. Further trajectory inferences shed light on these putative MSCs transitioning to cell states/subtypes associated with MSCs. Among these MSC-associated subclusters, static “pure” MSCs (non-differentiated) were found in both methods and showed a similar trend when mapped on a spatial slide, with increased expression on day 14 DPW, compared to day 2 DPW, with this increase of expression being concentrated in the wound bed (Fig. 5). Trajectory state inference showed putative MSC-associated cell subclusters in a dynamic state expressing markers of fibroblast, osteocyte, muscle and adipocyte cell subtypes, cell subtypes associated with MSC across taxa [52, 78, 79]. The integration of sn RNA-seq data and spatial transcriptomics suggests that these cells are potentially present throughout the skin and are transcriptionally active, during wound healing, especially during the tissue remodelling phase.

Variation in the cellular composition of flank skin and fin

Fish flank skin and fins are generally very similar tissues. One of the main differences between the tissues is that the cell types annotated as “skeletal muscle”, “progenitor pigment cells” and “fibroblasts” are practically only present in the flank skin, likely reflecting the different anatomical architecture of the flank skin. Erythrocytes also show a relatively large transcriptomic divergence between the flank skin and fins. Erythrocytes in fish unlike in mammals are nucleated and are also involved in immune responses [80], and perhaps the relatively low correlation between transcriptomes is simply a signature of greater transcriptional variation in this cell type.

The largest differences were observed in the osteo lineage, with osteocytes being more numerous in the flank skin and showing the greatest transcriptomic differences between the two tissues. Putative MSCs in the flank skin appear to differentiate into two distinct osteo-like cell subtypes, tentatively identified as bone precursors and bone-muscle precursors, based on their gene expression and transition states. The bone precursors do not seem to be in as transcriptionally active a state in the fin (Fig. 4b) compared to the flank skin (Fig. 4a), and the bone-muscle precursors are not present in the fins. These differences could be due to the distinct processes of bone and scale regeneration in these tissues. The regeneration of scales in the flank skin differs from that of bones in the fin, as it is coordinated by the activity of osteoblasts (secretion and mineralization of the bone matrix) and osteoclasts (resorption of the bone matrix) to form an exoskeletal appendage. Bone cells also participate in the remodelling and vascularization of nerves and blood vessels within the flank skin [81]. The regeneration of bone in the fin is regulated by osteoblast migration and differentiation [82] into the area where the bone is regenerated. However, the Wnt signalling pathway is common to both processes of ossification [74, 81], as it plays a role in many cellular regenerative processes [73, 81, 83]. One hypothesis is that wnt5 appears to be expressed in population 6, which is found exclusively in the skin, possibly indicating a scale-specific population, this would require further validation. This work paves the way to further understanding how these cells are involved in the regeneration of scales and bones by providing markers to better understand the regeneration of bone in these two distinct bony structures.

Putative mesenchymal stromal cells in fish

In this study, putative MSCs were characterized in Atlantic salmon. These putative MSCs are embedded in cell populations previously classified as fibroblasts [26] and show high expression of the two paralogues of cd34, itga4 and itga5, which are known markers of MSCs in mammals [39,40,41]. MSCs have been shown to be very similar to fibroblasts in morphology, proliferation dynamics and immunomodulatory capacities, as well as in gene expression patterns [79]. Some authors have even suggested that they are fibroblasts [84]. Further exploration of this population of putative MSC-like cells revealed signatures not only of stem cells but also of osteocytes, adipocytes, fibroblast and muscle cells. This suggests that this cell type is, in fact, a heterogeneous mixture of cell lineages, composed of putative stem cells differentiating into different cell subtypes consistent with an MSC-like identity.

In this study, PHATE was used to identify and characterize potential MSCs subcluster. This method provided insights into cell transitional dynamics, including their differentiation into osteogenic cell subtypes, which has also been described in humans [85]. Atlantic salmon, MSC-like cells isolated from fat tissue were described as highly plastic and undergoing transitions from osteogenic to adipogenic states [16]. These findings suggest that MSCs in salmon likely undergo transition from a “pure” stromal cell state to various lineages, including bone and fat cells. This aligns with other studies on cell-type-specific responses to challenge in Atlantic salmon. For example, single-nuclei sequencing has shown that in response to bacterial infection, hepatocytes in salmon liver rapidly change their transcriptomic profile while still retaining their cell-type identity [19]. What defines a cell type and how to characterize this is a much debated topic in this emerging research area [27]. Using approaches such as PHATE to capture differentiation dynamics within cell populations provides a valuable alternative to methods that describe static cell states, as it enables the identification of potential transitioning MSC-associated subclusters, such as fibroblast, bone, adipocyte and bone muscle precursor cell states—hinting that fish MSCs have the potential to differentiate into other cell lineages when influenced by different spatial niche environments during wound healing. Entering a certain path (e.g. fibroblasts vs. osteocytes) is usually irreversible [7, 16], although our results cannot confirm whether these putative Atlantic salmon stem cells are transitioning from one state to another. Further studies will be necessary to confirm these dynamics in vitro, but this study has provided new salmon putative MSC markers as well as markers for different subclusters of putative cell types associated with MSCs with different transition identities.

Wound healing

During the wound healing process, tissue damage necessitates replacement, and fish skin exhibits remarkable regenerative ability. Zebrafish demonstrate complete scarless recovery from cutaneous wounds, even in adults [4]. Fish, including salmon, employ mechanisms of wound closure similar to those in mammals, but they possess the unique capability of scarless regeneration throughout their lives, akin to embryonic mammals [3]. Central to this phenomenon are MSCs, which play a crucial role in both wound healing and tissue regeneration [86]. These versatile cells differentiate into various essential cell types critical for regenerating tissues such as connective tissue, bone and cartilage. This study detected putative MSC-specific markers at 2 DPW, although their distribution within the wound area and the skin was scattered. In contrast, at 14 DPW, these markers were notably concentrated within the wound bed. These results align with the putative pure MSC markers identified using both subclustering in Seurat and PHATE analysis (Fig. 5). This concentration suggests cell proliferation and possibly cell migration following the clearance of damaged fibres. The analysis indicates that MSCs may differentiate into various tissue types. For instance, bone, adipocytes and bone muscle share marker genes, as shown in the heatmap of the top five expressed genes within the PHATE cluster (Fig. 3c). An example is the gene zeb2b, involved in MSC differentiation and formation [22], which is expressed in clusters identifies as adipocytes (cluster 5) and bone (cluster 4). Other genes, such as Wnt genes, are unique to specific lineages (Fig. 4c). This finding aligns with previous studies showing that osteogenic and adipogenic cells share common transcriptomic features but diverge when committing to a specific cell type [16].

These putative MSC-associated subclusters were transcriptomically charted during in vivo wound healing. Different MSC subclusters showed distinct transcriptomic activity during various phases of wound healing (Fig. 6). For example, on day 2 of wound healing, markers associated with MSCs, along with bone and adipocyte precursors, were randomly expressed throughout the tissue section, with bone muscle precursors expressed at low levels. However, during the remodelling phase, these cells were found in specific locations. Bone precursors were expressed in the scale beds around the wound granulation site, while adipocyte precursors concentrated on the outer adipogenic layer and were most highly expressed at the edge of the granulation site. Bone muscle precursors followed a similar pattern but were more highly expressed within the granulation site, resembling the concentration of MSC precursors previously observed at the same site (Fig. 5). By day 14, the rebuilding of connective tissue had begun [17], requiring high MSC activity. A notable trend was observed in the fibroblast 2 and 3 populations, where increased relative expression and localization to the wound bed were evident. This finding is consistent with current literature highlighting the significant role of fibroblasts in wound healing [87]. This study is the first to transcriptomically classify putative MSCs at the cellular level, characterize their differentiating niches as MSC-associated cellular subclusters, and chart their transcriptomic activity during wound healing. However, these interpretations are limited by the small sample size (n = 2), which restricts the conclusions. Nevertheless, this work lays the foundation for future studies by providing unique markers for putative MSCs in salmon and distinct transcriptomic profiles for MSC-associated subclusters (Additional file 1: Tables S6, S7, S8, S9, S10 and S11). Further in vitro and in vivo studies will be needed to confirm these findings.

Potential applications of MSC knowledge in aquaculture

In humans, it has been shown that diet influences the differentiation and action of MSCs, leading to a variety of different health issues in multiple tissue types [88], and that MSC differentiation can shift toward adipogenic or osteogenic lineages by changing gene expression patterns [89]. This opens an avenue for the exploration of nutritional modulation in salmon, which could have significant ramifications for fish health. Non-healing skin wounds, including those caused by sea lice, are a constant issue in Atlantic salmon aquaculture, and diseases such as winter ulcer diseases caused by Moritella viscosa are becoming increasingly prevalent [2]. Enhancing or manipulating MSC function may contribute to the rapid and effective repair of the skin barrier, potentially mitigating the impact of skin damage and infections, and even reducing the use of antibiotics [90]. These manipulations may also target other multifaceted contributions of MSCs, such as their role in promoting angiogenesis by stimulating the growth of capillaries, ensuring adequate blood supply, or their ability to secrete growth factors and cytokines, orchestrating the complex wound healing process and facilitating tissue regeneration and recovery.

Conclusions

This study advances our understanding of the cellular composition of the flank skin and fins of Atlantic salmon and charts the transcriptomic activity of putative mesenchymal stromal cells during wound healing. The study describes the transcriptomic dynamics of putative MSC differentiation in fish. The results also reveal the distinct transcriptomic activity of different putative MSC-associated cellular subclusters during the process of wound healing in fish skin. Although these findings are limited to transcriptomic activity in a small number of replicates, they provide a foundation and the necessary tools for further in vivo exploration and validation to chart the dynamics of adult pluripotent transitioning cells in the barrier tissues of Atlantic salmon.

Methods

Skin and fin sampling

The skin and fins of post-smolt Atlantic salmon (Salmo salar L.) were obtained from two different sources for single-nuclei extraction. Six Atlantic salmon, with a mean weight of 486 g, were randomly netted from a 2000 L stock tank at the University of Stirling’s Marine Environmental Research Laboratory (MERL) (Machrihanish, Scotland). The fish were maintained at ambient sea temperature (14 °C) in full-strength seawater (33%) in a flow-through system, with dissolved oxygen (DO) concentrations ranging from 8.6 to 8.8 ppm. They were fed commercial salmon pellets (Inicio Plus, BioMar, UK) at 1% of their body weight per day. The fish were euthanized by anaesthetic overdose using MS-222 tricaine methanesulfonate (100 mg/L), and the brain was destroyed according to UK Home Office Schedule 1 guidelines. Fish lengths and weights were recorded immediately following euthanasia. Approximately 45 mg (1 cm³) samples of pectoral fin, dorsal fin, pelvic fin and flank skin were dissected from the fish, snap-frozen in liquid nitrogen and stored at − 70 °C until further processing. Two additional Atlantic salmon, with a mean weight of 25 ± X g, were sampled from a recirculating system at the Center for Aquaculture Technologies (Prince Edward Island, Canada). The fish were housed in 135-L tanks with recirculating water at approximately 12 °C. They were sedated using an anaesthetic dose of MS-222 according to SOP/CATC/2085, followed by a lethal blow to the head. Skin and fin samples were collected, snap-frozen on dry ice and stored at − 70°C until further processing. One skin sample, one pelvic fin sample and one pectoral fin sample from the first sampling, and one skin sample, one pelvic fin sample and one dorsal fin sample from the second sampling were selected for snRNA-seq.

Nuclear isolation, library construction and sequencing

A tween with salt and tris buffer (TST) protocol adapted from [91] was used for nuclear extraction; this protocol had been previously optimized for Atlantic salmon fins and skin and is described in [30]. Approximately 45 mg of flash-frozen skin or/fin was placed in a 6-well tissue culture plate (Stem Cell Technologies, cat. no. 38015) with 1 mL of TST buffer composed of 2 mL of 2X ST buffer, 120 µL of 1% Tween-20 (Sigma‒Aldrich, catalog no. P-7949), 20 µL of 2% BSA (New England Biolabs, catalog no. B9000S) and 1.86 mL of nuclease-free water. The tissue was minced on ice for 10 min using Noyes Spring Scissors (Fine Science Tools, catalog no. 15514–12). The resulting tissue homogenate was filtered through a 40-µm Falcon™ cell strainer (Thermo Fisher Scientific, catalog no. 08–771-2), and a further 1 mL of TST was used to wash the well and passed through the filter. The volume was increased to 5 mL with 3 mL of 1X ST buffer diluted from 2X ST buffer, which was composed of 146 µL of mM NaCl (Thermo Fisher Scientific, catalog no. AM9759), 292 µL of 10 mM Tris–HCl, pH 7.5 (Thermo Fisher Scientific, catalog no. 15567027), 10 µL of mM CaCl₂ (Vwr, E506-100 mL), 210 µL of mM MgCl₂ (Sigma–Aldrich, catalog no. M1028) and 9388 mL of nuclease-free water. The sample was then centrifuged at 4 °C for 5 min at 500 g and 4 °C in a swinging bucket centrifuge. The resulting pellet was resuspended in 1 mL PBS and 0.02% BSA buffer. The nucleus solution was then filtered through a 40-µm cell strainer (Falcon™). The TST and PBS buffers contained 200 U mL⁻¹ of the ultrapure protector RNase inhibitor (Sigma‒Aldrich, catalog no. 3335399001), and the 1xST contained 20 U mL⁻¹.

Salmon skin and fin nuclei were encapsulated into tiny droplets, with each droplet containing one cell and a unique barcode, allowing for the individual profiling of each one through the Chromium Single Cell Platform using the Chromium Single Cell 3' Library and Gel Bead Kit. v3.1 (10X Genomics, Pn-10001221) and the Chromium Single Cell A Chip Kit (10X Genomics, PN-120236) according to the manufacturer’s protocol. Briefly, after isolation, the nuclei were counted using a disposable hemocytometer (Neubauer C-Chip DHC-N01). The chromium was loaded with the aim of recovering 7000 nuclei in a Chromium 3’Chip. The nuclei were then partitioned into droplets in the Chromium controller. Once the emulsions were formed, RNA was barcoded and reverse transcribed. The resulting cDNA was then amplified and fragmented, and adaptors and sample indices were added. Libraries were sequenced on a NovaSeq 6000 platform by Azenta, resulting in approximately 220 million paired-end 150-bp reads per sample.

Genome indexing and read alignment with STAR

Genome indexing and library mapping were performed with STAR (version 2.7.10a) [92]. We appended the mitochondrial genome from the ENSEMBL V2 Atlantic salmon genome (Salmo_salar.ICSASG_v2.dna_rm.toplevel.fa.gz, v2, release 105, masked genome, assembly ID: GCA_000233375.4) to the ENSEMBL V3 Atlantic salmon genome (Salmo_salar.Ssal_v3.1.dna_rm.toplevel.fa.gz, v3.1, release 106, masked genome, assembly ID: GCA_905237065.2) for both the. gff and.fna files prior to indexing. The genome annotation.gff files were converted to.gtf files using gffread (v0.10.1) [93]. The Atlantic salmon genome was indexed using STAR (–runMode genomeGenerate). Each library was then mapped against its genome with the 10X V3 cell barcode whitelist (3 M-February-2018.txt) using standard parameters for single-cell libraries (–soloMultiMappers Unique –soloBarcodeReadLength 150 –soloType CB_UMI_Simple –soloUMIlen 12 –soloCBwhitelist 3 M-February-2018.txt –soloFeatures GeneFull –clipAdapterType CellRanger4 –outFilterScoreMin 30 –soloCBmatchWLtype 1MM_multi_Nbase_pseudocounts –soloUMIfiltering MultiGeneUMI_CR –soloUMIdedup 1MM_CR –readFilesCommand zcat –outSAMtype BAM Unsorted). The raw (unfiltered) files (genes.tsv, barcodes.tsv and matrix.mtx) generated for each sample were then used for downstream analysis.

Quality control, clustering, integration

The raw STAR files were then analysed in the R (v4) environment using Seurat (v4.1), [94]. Seurat objects for each library were created after removing nuclei with fewer than 200 features and features occurring in fewer than three nuclei. Nuclei where mtDNA features represented more than 10% of their total UMIs were removed, and then, mtDNA features were removed from the dataset (leaving 46,832 features). Upper and lower thresholds for UMI and feature counts per nucleus were then determined for each sample based on knee plot visualization. For all samples, only nuclei with more than 500 UMIs but less than 6000 UMIs and more than 500 features and less than 3500 features were retained.

The samples were then merged into a single Seurat object. The “v2” SCTransform version with the glmGamPoi method (v 1.8.0, [95]) was used to normalize RNA counts for each sample prior to calculating cell cycle scores using the “CellCycleScoring” function based on the “SCT” assay (See Additional File 1: Table S12, for the list of genes used), regressing out scores for the S and G2M cell cycle stages. The SCTransformation was then repeated as described above, but additionally, the “S. Score” and “G2M. Score” variables.

Linear dimension reduction was conducted for each sample using the “RunPCA” function with 50 PCs. After consulting the elbow plots for each sample, a UMAP using 20 PCs was run for each sample, and the “FindNeighbours” function was applied using 20 PCs before using the “FindClusters” function with a resolution of 0.2. DoubletFinder (v 2.0.3, [96]) was then applied independently to each sample by selecting the pK values with the highest associated BCmvn values. We assumed a 4% doublet formation rate (based on the Chromium instrument specifications) and adjusted for homotypic doublets (See Additional File 1: Table S2) for the remaining cells per sample after doublet removal.

The samples were integrated using 10,000 features, and the anchors were identified via the “rpca” reduction method and the “FindIntegrationAnchors” function. PCA was performed on the integrated dataset using 50 PCs, and 30 PCs were used for subsequent UMAP generation and clustering with a resolution of 0.2 (See Additional File 2: Fig. S9).

Assignments of clusters and markers

Markers for each cluster were assessed using the logistic regression method. Differentially expressed genes were calculated using the Seurat FindAllMarkers function, which applies the Wilcoxon rank sum test with default cut-offs (multi-test adjusted p-value < 0.01, log2-fold change > 0.25, with expression of the gene in at least 25% of nuclei in the cluster tested). Sample ID was included as a latent variable to help reduce batch effects among samples. 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, using sample ID as the background.

Four clusters were removed from the dataset due to a small number of distinguishing marker genes. Clusters 1, 3 and 7 also had very low average feature/UMI counts (See Additional File 2: Fig. S9). The SCTransform was then repeated for each sample based on the RNA assay as described above, and the integration of samples for each species was conducted as described, using 30 PCs for UMAP generation and a resolution of 0.2 for clustering. Additional file 2: Fig. S10 shows the distribution of UMIs and features per sample and per cell type after all filtering.

A full list of the markers for 18 skin clusters, assigned with the criteria described above, can be found in Additional file 1: Table S3a, which includes the top 20 markers for each Atlantic salmon skin cell type. This list includes both a priori markers and links to the citations that describe these markers as associated with specific cell types.

Comparison between the epidermis of skin and fins

To account for differences in sample sizes and sequencing, cell counts in each cluster were normalized by dividing the total number of cells of each cluster in that tissue type by the total number of cells captured in that tissue type. The transcriptomic similarity between tissues for each cell population was determined by Pearson rank correlation, using the mean expression of the set of all marker genes. Additional file 1: Table S3a includes the top 20 markers, while Additional file 1: Table S3b contains the main representative markers. Furthermore, Salisbury et al. [26] give a further clarification of these clusters.

Inference and dynamics of mesenchymal stromal cells

In previous work, we were increasingly interested in the previously annotated fibroblast population 1, which was also shown to be reactive to infection during sea lice infection [26]. Known stem cell markers, specifically two paralogues of cd34 [39], integrin α5 (ITGA5) [41] and integrin alpha 6 (ITGA6) [97, 98], were plotted across the skin atlas and found to be expressed in this fibroblast population (Fig. 2 and Additional file 2: Fig. S1b). This cluster was further analysed to infer the expression of specific putative MSC markers within it. Fibroblast population 1 was subset into a new Seurat object and normalized with SCTransform to generate updated residuals. The clusters were re-clustered using the number of principal components (PCs) informing the clustering, with a resolution of 0.3 to provide biologically meaningful clusters (20 PCs). Differentially expressed genes for each cluster were then calculated as described, using the Seurat FindAllMarkers function (with the same parameters as for the global Assignments of Clusters and Markers detailed above). This new clustering identified five cell subsets of MSCs, each showing unique markers. These populations were then examined for stem markers, and population 2 was identified as containing putative “pure” stem markers (Additional file 1: Table S6). PHATE [43] was used to generate two-dimensional embeddings to identify differentiation processes between cell types. PHATE embeddings were used to recluster the populations using the FindNeighbors and FindClusters Seurat commands. Seven putative subclusters of cell subtypes/states were identified. The putative pure MSC subset homologous to that population 2 was taken to compare to the Seurat cluster in spatial analysis in Fig. 5 and equivalence of these two clustering: Additional file 2: Fig. S1b.

Wound healing samples and spatial transcriptomics

An Atlantic salmon wounding experiment was carried out as described by [17] with minor modifications. Briefly, the fish were fully anaesthetized with MS-222 (Sigma‒Aldrich). Incisional skin ulcers were introduced on the flank of the fish with a sterile scalpel blade, which was cut through all skin layers, resulting in deep cutaneous wounds with damaged muscle fibres. Samples (N = 2 fish per sampling point and group) were taken on days 2 and 14 post wounding at the wound site. Tissue samples were frozen immediately in liquid nitrogen and stored at − 80 °C. The fish were euthanized with an overdose of anaesthetic (MS-222) prior to tissue sampling. The fish were in good health prior to incisional wounding, and there were no mortalities. Spatial transcriptomic libraries were prepared as described previously [33]. Summary statistics for spatial transcriptomics libraries: Additional file 2: Table S1b.

To map the homologous putative pure stem markers (population 2 in Seurat and 1 in PHATE analysis), the top unique markers were selected (Additional file 1: Table S6) for the putative pure MSC subclusters from the Seurat population 2 (Fig. 2) and PHATE analysis population 1 (Fig. 3), as well as the average expression and the contribution of each individual marker to these: Additional file 1: Table S10. To map the different MSC-associated subclusters onto the spatial slide and considering the co-expression of these cell subtypes among the MSC subclusters, the top 20 expressed genes for each cell subtype were taken from the PHATE output (Additional file 1: Table S11). The average expression levels of 20 marker genes for each putative MSC cell subtype were visualized on a spatial transcriptomic slide (Fig. 6). These averages were computed by averaging the expression values of individual markers and the contribution of each marker to the overall average was determined by normalizing individual marker averages against the total average expression. However, not all marker genes were detected in the spatial transcriptomic data, resulting in variability in the number of markers included in the analysis for each subtype.

Data were normalized using the “v2” SCTransform which adjusts for sequencing depth and technical variation. Markers were selected as listed above. For each marker, the average expression across all spatial spots was calculated using the SCTransform-normalized expression matrix. The overall average expression of the selected markers was then computed by averaging these individual marker values. The contribution of each marker to the overall average expression was determined by normalizing individual marker averages against the total average expression. This calculation assessed the spatial distribution and relative contributions of putative MSC cell markers within the tissue samples. The results can be found in Additional file 1: Table S10 for “pure” MSC populations and Additional file 1: Table S11 for MSC-associated subtypes, where a higher relative abundance value indicates greater overall marker expression within the tissue and comparison. The mean average expression can be compared across different datasets or experimental conditions to assess the relative expression of these markers.