Understanding host response to infectious salmon anaemia virus in an Atlantic salmon cell line using single-cell RNA sequencing

Abstract

Background

Infectious Salmon Anaemia Virus (ISAV) is an Orthomixovirus that represents a large problem for salmonid aquaculture worldwide. Current prevention and treatment methods are only partially effective. Genetic selection and genome engineering have the potential to develop ISAV resistant salmon stocks. Both strategies can benefit from an improved understanding of the genomic regulation of ISAV pathogenesis. Here, we used single-cell RNA sequencing of an Atlantic salmon cell line to provide the first high dimensional insight into the transcriptional landscape that underpins host-virus interaction during early ISAV infection.

Results

Salmon head kidney (SHK-1) cells were single-cell RNA sequenced at 24, 48 and 96 h post-ISAV challenge. At 24 h post infection, cells showed expression signatures consistent with viral entry, with genes such as PI3K, FAK or JNK being upregulated relative to uninfected cells. At 48 and 96 h, infected cells showed a clear anti-viral response, characterised by the expression of IFNA2 or IRF2. Uninfected bystander cells at 48 and 96 h also showed clear transcriptional differences, potentially suggesting paracrine signalling from infected cells. These bystander cells expressed pathways such as mRNA sensing, RNA degradation, ubiquitination or proteasome; and up-regulation of mitochondrial ribosome genes also seemed to play a role in the host response to the infection. Correlation between viral and host genes revealed novel genes potentially key for this fish-virus interaction.

Conclusions

This study has increased our understanding of the cellular response of Atlantic salmon during ISAV infection and revealed host-virus interactions at the cellular level. Our results highlight various potential key genes in this host-virus interaction, which can be manipulated in future functional studies to increase the resistance of Atlantic salmon to ISAV.

Introduction

Aquaculture is currently the fastest growing food industry worldwide [1], and its products are a fundamental component of healthy and sustainable human diets. Atlantic salmon is the most important aquaculture fish species by value [1], and is a high-tech industry underpinned by large research and development programmes. Despite the high degree of technological progress, infectious diseases remain a major issue for salmon production. A notorious example is infectious salmon anaemia (ISA), a potentially fatal disease caused by virulent strains of the ISA virus (ISAV). ISAV is a fish Orthomyxovirus, and therefore, from the same family as Influenza viruses [2]. It is a segmented negative stranded RNA virus, with a genome that comprises 8 RNA segments coding for at least 10 proteins [3]. The clinical presentation of an ISAV infection is a severe anaemia, commonly accompanied by haemorrhage and necrosis of various organs [4]. ISA disease can cause up to 90% mortality in sea pens [5], and severe outbreaks have been responsible for the decimation of entire national aquaculture industries [4, 6].

Strategies to control and mitigate ISA varies across countries, although they mainly focus on disease containment [7]. As a notifiable pathogen, upon detection of a clinical ISAV outbreak, whole stocks have to be culled and the farm quarantined for a period of time [8]. Vaccines for ISA are available. Nevertheless, they do not offer complete protection and have proven insufficient to fully control the disease [9]. Encouragingly, resistance to ISA has been shown to be moderately heritable, which has been exploited by breeding companies to increase genetic resistance in salmon stocks based on regular experimental challenge testing of broodstock families and genomic selection [10,11,12,13,14,15]. While survival during such disease challenges is an important target trait for improvement, a detailed understanding of host-pathogen interactions holds great potential for fast tracking breeding efforts and the development of effective therapies against ISA. Firstly, it can assist in the identification of functional genetic variants to improve genomic selection accuracy and its persistency across distant relatives. Secondly, it can assist with the design of new vaccination, treatment, or management strategies. And thirdly, it can lead to targets for future genome editing studies to potentially develop fully ISA-resistant salmon strains.

Previous studies in Atlantic salmon have reported a notable up-regulation of the innate immune system in response to ISAV [15,16,17,18], albeit there is a large degree of tissue specificity in the responses [19]. This is expected since pathogen infections lead to complex and dynamic interactions with the host and its immune system. In vitro models provide simplified systems to study interactions between pathogens and the host cell machinery, which can help break down the multidimensional responses observed in vivo. The establishment of Atlantic salmon cell lines from various tissue types – e.g., gill, heart, and head kidney – has resulted in several in vitro models becoming available for the study of the main bacterial and viral diseases affecting salmon production [20,21,22,23]. In ISAV research, multiple in vitro studies have been performed, describing a rapid interferon response [24, 25] or the interplay between host and virus over the control of oxidative stress and apoptosis [26, 27]. In vitro studies using cell cultures have also been key to understand the role of the different ISAV proteins and their interaction with cellular mechanisms [28,29,30,31], further highlighting the importance of suitable model systems for the study of cell biology in Atlantic salmon.

While in vitro cell culture systems have led to key insights into cellular immune mechanisms and host-pathogen interactions, our understanding of ISA at the cellular and molecular level is still incomplete. One significant limitation of previous studies is that they were based on the analysis and interpretation of bulk RNA sequencing data, thus measured population-level cellular responses with limited resolution. In contrast, single-cell sequencing [32] allows us to obtain a more comprehensive picture of cellular heterogeneity, helping to unveil the complex dynamics of viral transcription and host molecular responses. Single cell sequencing has yet to be widely applied in Atlantic salmon, but has been used to gain insight into cellular transcription in gills [33]. Here, to better understand the molecular mechanisms of ISAV infection and the Atlantic salmon response to the virus, we have performed an in vitro ISAV infection in the salmon head kidney 1 (SHK-1) cell line. SHK-1 is a macrophage-like cell line that was established from adherent leukocytes isolated from the Atlantic salmon head kidney [34], the main hematopoietic organ in fish [35, 36]. Given SHK-1 is extensively used for ISAV research and viral diagnosis [26, 37,38,39], the high-resolution functional characterization of this cell line will contribute to build a rigorous model upon which to evaluate candidate targets or identify variants that may confer increased host resistance. The results of this study provide the first high-dimensional single-cell analysis of an ISAV infected fish cell line, revealing the heterogeneity of the infection process, and enabling the identification of early markers of disease and potential functional targets to impair ISAV infection.

Methods

Cell line

Atlantic salmon head kidney 1 (SHK-1) cells, an immortalised macrophage-like cell line derived from Atlantic salmon (Salmo salar), was obtained from the European Collection of Authenticated Cell Cultures (ECACC-97,111,106). The cells we obtained had been passaged 63 times, and the cells we used for the experiments were passaged 6 additional times for a total of 69 times. All cells were grown as a monolayer in L15 complete media (L15*), L15 (Sigma-Aldrich, St. Louis, USA) supplemented with 5% heat-inactivated foetal bovine serum (FBS) (Gibco, Waltham, USA), 40 µM β-mercaptoethanol (Gibco), 100 units/mL penicillin and 100 µg/mL streptomycin (Gibco). Cells were cultured in an incubator at 22 ± 1 °C without CO₂. SHK-1 was split 1:2 at 80% confluency with 1/3 conditioned media.

Viral stock

ISAV stock was obtained from Marine Research (Aberdeen, isolate V4782, derived from UK 2008/2009 outbreak). ISAV was passaged once on 1 x T175 SHK-1 cells (20 mL L15* but with reduced serum (2%) and the supernatant collected after 50% cell death, centrifuged, sterilised using a 0.45 μm filter) and stored in 1 mL aliquots at -80 ^C until use. The viral dose was estimated by flow cytometry. To do this, SHK-1 cells were infected with different doses of the virus for 48 h in L15* with reduced serum (2%) at 17ºC, followed by antibody staining (ADL anti-ISAv, AquaMab-P10 and anti-mouse-GFP, Invitrogen, CA, #A21202) and fluorescence quantification (BD Fortessa X20).

ISAV challenge

The SHK-1 cells were seeded at a density of 2.5 × 10⁵ SHK-1 in 6-well plates in L15 + 10% FCS and antibiotics (penicillin and streptomycin) 24 h prior to the challenge experiment. THE ISAV dose was determined using flow cytometry, and an ISAV infection dose of over 66% was used. The cells were left uninfected (control) or infected with 200 µL of ISAV in L15* with reduced serum (2%) and incubated at 17ºC for 24 h, 48 and 96 h. Cells (three infected time points and control) were collected by trypsin treatment, washed in phosphate buffered saline (PBS) and suspended to 10⁶ cells/mL in PBS + 0.05% BSA (Bio-Rad cell counter) (Fig. 1).

Single-cell RNA sequencing library preparation and sequencing

Cells were counted and checked for viability on a Cell Counting Slide (Bio-Rad Laboratories Ltd) and then appropriately diluted according to 10x Guidelines. Each individual sample was loaded separately into the 10x Chromium machine and 10x single cell 3’ GEM kit v3 was used to generate the libraries following the manufacturer’s protocol. The quality of the resulting libraries was assessed in a Bioanalyzer (Agilent) and sequenced in two lanes of a NovaSeq SP flowcell with a cycle setup of 28/8/91 at Edinburgh Genomics.

Single-cell RNA sequencing analysis

The raw single-cell RNA-seq samples control, [0 h (uninfected), 24 h, 48 and 96 h post-infection] were demultiplexed and mapped to a combined reference transcriptome of Atlantic salmon ICSASG_v2 (GCA_000233375.4) and ISAV (RefSeq, GCF_000854145.2) using Alevin/Salmon v1.4.0 [40]. The resulting raw count matrices were analysed using Seurat v3.1.5 [41] in R v3.6.3 [42]. Each library was loaded individually, discarding genes identified in fewer than 3 cells, and cells with fewer than 200 expressed genes. DoubletFinder V2.0.3 was used to assess the impact of doublets on cell clustering in each library, which was negligible. Thereafter all libraries were merged into a single Seurat object. For each cell, the percentage of mitochondrial and viral transcripts were estimated. After inspection of quality control parameters (Supplementary Fig. 1), cells with fewer than 1,000 or more than 10,000 genes or with percentage of mitochondrial transcripts above 25 were discarded. The infected/uninfected status of all cells was calculated as previously described [43]; based on the kernel density estimate of the distribution of percentage of viral counts on the log10 scale, which was used to find the first local minima. The percentage of viral transcripts in each cell was calculated using the “PercentageFeatureSet” function of Seurat using all viral genes. Cell cycle scores for each cell were estimated using the orthologues of Seurat´s (v3.1.5; Stuart at al. 2019) list of mammalian cell cycle markers. The gene counts of the filtered cells were normalised using ‘SCTransform’ with the percentage of mitochondrial RNA and cell cycle scores as regression variables. Dimensionality reduction was performed using the Uniform Manifold Approximation and Projection (UMAP) method using 30 principal components (determined using an Elbow plot; Supplementary Fig. 2), and clustered in groups using the ´FindNeighbors’ and ´FindClusters´ functions of Seurat (Stuart at al. 2019). Marker genes for each cell group were determined using the Wilcoxon Rank Sum test. Differentially expressed genes showing a log fold change > 0.25 and false discovery rate (FDR) corrected p-values < 0.05 were considered to be marker genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were carried out using KOBAS v3.0.3 [44]. Briefly, salmon genes were annotated against the KEGG protein database [45] to determine KEGG Orthology (KO). KEGG enrichment for gene lists was tested by comparison to the whole set of expressed genes (obtained from the Seurat object) using Fisher’s Exact Test. KEGG pathways with ≥ 5 DE genes assigned and showing a Benjamini-Hochberg FDR corrected p-value < 0.05 were considered enriched. Finally, the correlation between viral and host genes was calculated. To compute the correlation of each viral gene with the whole salmon transcriptome, only the cells where that viral gene was expressed were selected. Person correlation between viral genes and the host genes was computed using the “cor” function in R, considering only correlation with r > |0.3|.

Results and discussion

To assess the intracellular activity of ISAV and study the corresponding host response, single-cell RNA-Seq was performed on SHK cells infected with ISAV 24 h, 48 and 96 h post infection (Fig. 1). Uninfected SHK cells were used as controls. Raw reads were assigned to either Atlantic salmon genes or ISAV genes using a combined host-virus reference transcriptome. A total of 638, 512, 296 and 323 cells passed quality control filters for the control, 24 h, 48 and 96 h samples respectively. Clustering of the cells in these four libraries based on the host transcriptional response showed a clear separation between control and cells infected for 48 and 96 h. However, there was no transcriptional separation between control and cells infected for 24 h (Fig. 1).

Fig. 1

Experimental design and cell clustering. Diagrammatic figure showing the experimental design and workflow. Briefly, four SHK-1 cell plates were seeded at the same time, and three plates were infected 24 h, 48 and 96 h prior to harvesting. The four plates were collected 120 h post seeding, and immediately processed using the 10x workflow for single-cell RNA-seq library preparation. Libraries were sequenced using Illumina technology, and processed using Alevin and Seurat (see Methods) to cluster the cells of each sample according to their transcriptome

Viral transcription during infection

The percentage of infected cells (estimated following [43]) was 16% at 24 h post-infection, and 32–33% at 48 and 96 h (Fig. 2A). The infected cells formed two different groups, one with a small number of cells infected for 24 h, and the largest one with cells infected for 48 and 96 h where the percentage of viral transcripts was up to 22% of the transcriptome of the cells (Fig. 2B). Expectedly, the percentage of virus transcriptome, relative to the total cell transcriptome increased over the course of infection, mostly ranging between 2% and 15% (Fig. 2C). This is consistent with previous in vitro infections in other Orthomyxoviridae, although a larger heterogeneity has been described with viral transcript levels reaching up to 90% of the cell transcriptome [43, 46]. Viral transcript levels can vary substantially depending on the cell type [43], but also ISAV infection progression is typically slower than that of other Orthomyxoviridae, thus, later timepoints could show higher ISAV transcript levels.

The timing and level of expression of each viral gene also varies along the infection process (Fig. 2D). Of the 10 ISAV genes, 6 were found to be expressed in the cells. The viral polymerase basic protein 2 (PB2) gene is the most expressed throughout infection, followed by the polymerase acidic protein (PA) and the non-structural proteins 1 (NS1) and 2 (NS2), also known as nuclear exporting proteins or NEP (Fig. 2D). NS2 is a splicing product of NS1 [47] and is expressed at lower levels according to previous reports in ISAV [48] and Influenza A [49]. The expression of PB2, PA and NS2 showed a high positive correlation (0.92–0.96; Fig. 2E), suggesting their co-expression early during the infection, while NS1 showed a different expression pattern. Expression of the hemagglutinin (HA) is found in a small fraction of cells, and it shows the lowest correlation with the other (expressed) viral genes (Fig. 2E). Expression of the bicistronic mRNA encoding the matrix protein (M1) and a nuclear export protein (M2, also termed S8ORF2), of the polymerase basic 1 (PB1) and the viral membrane fusion protein (F0) was not detected in our dataset (Fig. 2D), which suggests they are either expressed late during the infection or expressed at low levels. In Influenza A, all the viral mRNAs show similar transcription patterns [50], however as mentioned above, ISAV infection typically progresses at a slower rate so it is possible that the viral mRNA dynamics are slightly different.

Fig. 2

ISAV infection and transcription. (A) UMAP dimensionality reduction plot showing the infected and uninfected cells, where the threshold for infection was the local minima of the kernel density estimate of the distribution of percentage of viral counts on the log10 scale [43]. (B) UMAP dimensionality reduction plot showing the viral load of each cell measured as percentage of viral transcripts, calculated using all viral genes. (C) Viral load in each sample, measured as percentage of viral transcripts, calculated using all viral genes. (D) Dotplot showing the expression level, calculated as the percentage of the transcriptome of the cell that the viral genes represents, and the percentage of host cells with transcripts of each viral gene in each sample. (E) Heatmap showing the Pearson correlation between the expression of the viral genes

Host response to ISAV

The cells clustered into 8 different groups according to host gene expression (Fig. 3A). Despite working with a relatively homogeneous cell line, major transcriptional differences between the cells were linked to infection status and host response. Clusters 1 to 4 are formed by cell from the control and 24 h samples; 1, 2 and 3 are predominantly composed by cells from the control sample, but also have large numbers of 24 h cells (27.5 to 44.8% of the total cells in the cluster). On the contrary, cluster 4 is composed mostly of 24 h cells (82.5%), and a high proportion of those cells are infected with ISAV. On the other hand, cluster 6 is formed by 48 and 96 h uninfected cells, while cluster 7 and 8 are 48 and 96 h infected cells, respectively. The clear separation of control-24 and 48 h-96 h uninfected cells suggests paracrine signalling by infected cells may have caused altered gene expression in uninfected cells at the two latest time points. Paracrine signalling has been described before in response to Influenza A [51] and other viruses [52, 53], eliciting an antiviral interferon-mediated response in uninfected cells. We further tested this hypothesis by comparing control cells to uninfected cells at the infected timepoints (Supplementary File 1). The number of differentially expressed genes is quite high (835, 2607 and 3402 for control vs. 24 h, 48 and 96 h, respectively). However, our experimental design presents two potential caveats. First, the viral stock may contain some cytosolic content of death cells, and this could also have an effect on the cells transcriptome, indistinguishable of paracrine signalling; however, this potential cytosolic content does not seem to have an effect on the transcriptome of the cells at 24 h, and therefore we consider its effect to be negligible. Secondly, while we aimed for a relatively high infection rate, the results suggest that only a small fraction of cells were infected. It is therefore plausible that there are some reinfection events in the later timepoints. However, the apparently slow replication of the virus in our study (Fig. 2) and the fact that early infection does not seem to have a large impact on the transcriptome of the cells suggest that the observed changes in the transcriptome of non-infected cells at 48 and 96 h post infection should be in a great extent due to other causes. In any case, it is important to keep in mind these caveats in the interpretation of the later timepoints. Finally, cluster 5 is the most difficult to interpret since it is formed by a low number of cells of each sample apparently not connected to the infection process (Fig. 3A); this cluster could represent a different cell type since this cell line was originally reported to be heterogeneous [54], however we do not believe that is the case since our cells have undergone a high number of passages and they are homogenous under the microscope and are clustered together by flow cytometry.

Each of these clusters is characterised by specific marker genes (Fig. 3B, Supplementary File 2). Clusters 1 and 2 are characterised by the high expression of certain ribosomal genes, which also show a high expression in half of the cells of group 3. This could be either (i) an “artefact” of the cellular response in infected / bystander cells (high expression of additional genes leads to a lower relative expression of ribosomal genes), or (ii) due to manipulation of the ribosomal machinery by ISAV, which can be connected to the hijacking of host translational machinery by the virus, a common strategy used by viruses to repress the cellular mRNA translation and allow the preferential translation of viral mRNAs [55]. On the other hand, a group of immune genes, including interferon alpha 2 (ifna2) or interferon regulatory factor 2 (irf2), are highly expressed in infected cells at 48 h (cluster 7), but moderately at 96 h (cluster 8), suggesting either viral repression or a shift in the immune response later during the infection. Further, this group of genes is not highly expressed in bystander cells at 48 and 96 h (cluster 6), suggesting that paracrine signalling by infected cells does not promote the activation of an interferon response in bystander cells.

To facilitate an unbiased interpretation of the marker genes of each cluster we performed a KEGG pathway enrichment analysis (Fig. 3C, Supplementary File 3). This confirmed the enrichment in ribosomal genes in the control / uninfected clusters (1 and 2). Clusters 3 and 4 (mostly 24 h cells, particularly in cluster 4) are enriched for “focal adhesion” and “ECM-receptor interactions”, which may reflect a consequence of the interaction of ISAV with the extracellular matrix and membrane receptors. The genes in these pathways include phosphoinositide 3-kinase (pi3k), focal adhesion kinase (fak) and c-Jun terminal kinase (jnk), which have been shown to be important for the entry of Influenza A into the cytoplasm [56,57,58]. This is consistent with the similarity in the entry mechanisms of Influenza A and ISAV [39, 59] and suggests that molecular interactions during viral entry are also conserved to some extent. These early genes are good candidates for functional studies aimed at disrupting ISAV entry into salmon cells and eventually to develop strategies to tackle ISAV in Atlantic salmon aquaculture, including the development of genome edited ISAV-resistant stocks.

Uninfected bystander cells at 48 and 96 h (cluster 6) showed enrichment for terms such as “proteasome”, “RNA degradation”, “mRNA surveillance pathway” or “cytosolic DNA-sensing pathway”, but they do not express typical immune genes (i.e. interferon or cytokines). This cluster shows expression of the seven lsm genes that compose the LSm1-7 ring involved in the degradation of messenger RNA in the cytoplasm [60]. The host RNA decay machinery can rapidly degrade viral mRNAs, and as such viruses have developed strategies to avoid or even use to their advantage parts of this pathway, including the LSm1-7 ring [61]. Many proteasome-specific genes are also expressed in this cluster, including two of the subunits of the immunoproteasome, involved in the degradation of intracellular proteins, including those of viral origin for presentation on major histocompatibility complex [62]. Again, many virus interact with the immunoproteasome as part of their infection strategy [63]. Connected to the proteasome, this cluster also expresses all the genes of the ubiquitination machinery: E1 ubiquitin-activating enzymes, E2 ubiquitin-conjugating enzymes and E3 ligases of all the types (HERC3, U-box and RING-finger and their adaptor proteins). Up-regulation of genes related to ubiquitination has been previously described in resistant fish in response to ISAV [64]. Moreover, proteins encoded by ISAV segment 8 (M1 and M2), which are known to interfere with interferon signalling, can be conjugated to ubiquitin and the ubiquitin-like interferon simulated gene 15 (ISG15) in Atlantic salmon cells [65]. Several viruses can hijack the host ubiquitination machinery, and in particular Influenza uses host ubiquitination as part of its cell entry and replication strategy [66,67,68]. Ubiquitination is also a fundamental part of the host immune response, controlling multiple cellular processes, and its role in the activation of the interferon response via TRIM25 - RIG-I interaction is well-documented [69]. Ubiquitination has been barely studied in fish; considering the large expansion of TRIM E3 ubiquitin ligases in fish [70], our recent results and the evidence in other species, ubiquitination deserves more attention to try to understand ISAV infection in Atlantic salmon.

Bystander cells at 48 and 96 h (cluster 6) also showed enrichment in ribosomal genes, as control and cells infected for 24 h (cluster 1 and 2), which was unexpected based on the expression pattern of ribosomal genes in the heatmap (Fig. 3B). However, we found that the underlying specific ribosomal genes are different, with only two genes in common between the two groups of cells (out of 55 in each list). The ribosomal genes enriched in Clusters 1 and 2 showed lower expression in the infected cells at 48 and 96 h (clusters 7 and 8), while the ribosomal genes enriched in Cluster 6 showed high expression also in infected cells at 96 h (cluster 8), but not so clearly in 48 h infected cells (cluster 7), although a wider expression level distribution was observed in these cells (Fig. 3D). A closer look at the gene lists revealed that most of the ribosomal genes more expressed in clusters 6 and 8 (48 h + 96 h bystander cells and 96 h infected cells) code for proteins of the mitochondrial ribosome (Supplementary file 2). Mitochondria activate anti-viral immune responses through mitochondrial antiviral signalling (MAVS) [71] and can also initiate apoptosis [72]. Influenza A NS1 protein has been observed in the mitochondria [73], and is able to alter mitochondria morphodynamics [74]. Our results suggests that paracrine regulation upon ISAV infection involves the mitochondria as a defence mechanism, and the fact that these mitochondria ribosomal genes are also overexpressed at 96 h in infected cells (Fig. 3D, cluster 8) suggests that they may play an active role in the immune response against ISAV.

Finally, clusters 7 and 8 (infected 48 and 96 h cells, respectively) show up-regulation of key immune response genes, such as interferon alpha 2 and 3 (ifna2, ifna3), interferon regulatory factor 2 (irf2) or C-X-C motif chemokine 10 (cxl10). These genes are slightly more up-regulated in the 48 h than in the 96 h cells (e.g. logFC 3.7 vs. 1.7 for ifna3, 3.8 vs. 1.6 for ifna2; Supplementary File 2), with pathways such as toll-like receptor signalling, RIG-I receptor signalling, TNF signalling or cytokine-cytokine interaction. However, these genes were not expressed in bystander cells at 48 and 96 h (cluster 6). Our results align with previous studies that have shown that ISAV triggers an early immune response characterised by the activation of interferon genes, both in vitro and in vivo [19, 75, 76]. IFNA2, IFNA3 and CXL10 are secreted proteins that act on other cells to promote an antiviral state and could explain the “activated” state of cluster 6 cells. However, surprisingly, they did not induce the up regulation of interferon genes in these cells, which is quite different to what has been reported during Influenza infections, where the paracrine signalling is important at 12 h post infection and increases the expression of interferon stimulated genes in bystander cells [51]. There are different potential explanations for this finding, for instance the virus may be able to inhibit the interferon-related paracrine signalling quite early during the infection, but it is also plausible that paracrine signalling functions differently in fish than in mammals. Further investigations are necessary to understand fish paracrine signalling, including studies with other fish species and viruses.

Fig. 3

Response to ISAV in the Atlantic salmon SHK-1 cell line. (A) UMAP dimensionality reduction plot showing the clustering of the cells according to the host cell transcriptome; (B) Heatmap showing the expression pattern of the top marker genes for each cell cluster; (C) Barplot with the fold enrichment of selected KEGG pathways in each cell cluster, ordered within each cluster according to their significance; (D) Ridgeplot showing the proportion of the transcriptome that represent the two sets of “Ribosome” genes (KEGG pathway) enriched in (i) cluster 1 and (ii) cluster 5

Interactions between host and viral genes

To evaluate whether specific viral proteins could have a direct impact on the expression of host genes, a correlation analysis was performed between expressed viral and host genes in infected cells (Fig. 4, Supplementary file 4). No strong correlations were observed (highest r = 0.51) and only those correlations with r > |0.3| were considered. PB2, PA and NS1 + NS2 expression patterns showed significant correlation with largely the same set of host genes. Considering these viral genes are the most expressed, these correlations do not necessarily represent regulation of host genes by viral proteins; they could also represent genes that respond proportionally to viral load. In fact, the host genes exhibiting the highest positive correlation include the immune genes mentioned in the previous section (ifna2, ifna3, irf2) and the regulator of antiviral signalling probable ATP-dependant RNA helicase DHX58 (dhx58). The highest correlation was observed with cholesterol 25-hydroxylaseI A (c25ha), an interferon-stimulated gene involved in the regulation of cholesterol biosynthesis, which has broad antiviral activity [77]. There were also host genes showing negative correlations with these viral genes, especially NS1 + NS2, for example with Phosphoinositide-3-Kinase Regulatory Subunit 1 (pik3r1), which may be exploited to tackle ISAV infections; Influenza A NS1 protein modulates the PIK3 / Akt pathway via interaction with PIK3, and disrupting this interaction led to reduced infectivity [78]. Several transcription factors also showed a negative correlation with NS1 + NS2, such as twist-related protein 2 (twst2) or zinc-finger protein SNAI1 (snai1).

The expression of over thirty cytoplasmic ribosomal genes is negatively correlated with the expression of the viral nucleoprotein NP (r = -0.3 to -0.4), suggesting a negative regulation of the host translation machinery by the virus. While NP is a structural protein, it has been shown that it interacts with host factors to facilitate viral replication [79, 80], although no interactions between NP and ribosomal genes have been previously described. The ribosomal machinery is an important part of viral replication since viral genomes do not usually harbour mRNA translation genes [81], and therefore they need to recruit host proteins including ribosomal proteins (RPs) to complete their replication cycle. Our results show that the viral protein NP may interact or at least drive the expression of specific mitochondrial RPs. While it is well known that viruses interact with ribosomal proteins as part of their infection strategy [82, 83], recent studies have also highlighted that some ribosomal proteins can have an antiviral function by either interacting with viral proteins to inhibit transcription/translation or by activating antiviral defence signalling pathways [82]. Non-infected cells in 48 and 96 h infected samples showed a clear up-regulation of mitochondrial RPs, which could suggest that they form an important part of a potential antiviral pathway in bystander cells, which may be repressed by the viral protein NP in infected cells considering the negative correlation.

There were also several chaperons showing positive correlation with viral genes, such as host heat shock protein 30 (hsp30) and clusterin (clu), which has antiapoptotic activity and is directly targeted by the Influenza A virus nucleoprotein [84]. Chaperones are used by viruses to facilitate the folding and assembly of the viral proteins [85]. Another protein of interest is sequestosome-1 (sqstm1), a receptor required for selected autophagy that inhibits Seneca valley virus and avian influenza replication [86, 87]; previous studies have reported that mammalian cells can use autophagy to restrict the replication of avian Influenza [88], but it has also been reported that it can promote the replication of influenza A [89, 90]. Sqstm1 is also known as ubiquitin-binding protein p62, being able to bind ubiquitin and providing the link between ubiquitination (relevant in bystander and infected cells; Fig. 3C and Supplementary File 2) and autophagy [73].

While it is not possible to discriminate whether the observed correlations between viral and host proteins are caused by (i) proportional response of the host to viral load or (ii) regulation of the host transcription by the virus, the genes and processes highlighted help us better understand the infection process and represent good targets for further functional studies aiming to develop ISAV-resistant Atlantic salmon. We have also identified genes without previous known associations to viral processes, such as the transcription factors twst2 or snai1. In vitro functional analyses are required to confirm the potential of these genes to sway the salmon-ISAV interaction in the fish favour (e.g. via disruption of their interactions with viral proteins), paving the way for their application in aquaculture to limit the impact of this disease and improving fish welfare and food security.

Fig. 4

Correlation between viral gene expression and host gene expression. All Pearson correlation values are between |0.3| and |0.5| (Supplementary File 4)

Conclusion

Here we applied single-cell RNA sequencing to study ISAV infection in Atlantic salmon SHK-1 cells. This simplified system, in combination with this novel genomics technology, has allowed us to finely characterise the transcriptomic changes that occur during early ISAV infection both in the virus and the host, and suggest potential interactions. Early ISAV infection of Atlantic salmon SHK cells is characterised by high expression of the PB2, NS1 + NS2 and PA viral genes, which show highly correlated expression patterns. After 24 h of infection infected SHK cells showed similar transcriptomic profiles to control cells, however genes and pathways connected to viral entry were identified, which could be good targets for functional studies to impair viral infection. At 48 and 96 h post-infection there was a clear host transcriptomic response to the virus which also affected uninfected cells, therefore suggesting paracrine signalling from infected cells. This putative paracrine signalling seemed to produce an unspecific alert state in the uninfected cells, involving mRNA surveillance and RNA degradation, ubiquitination and the proteasome, or mitochondrial activation. These pathways were also up-regulated in 48 and 96 h infected cells, but they also showed a clear antiviral interferon response, including genes such as ifna or irf2. Correlation between the expression of viral and host genes suggested potential negative regulation of certain host transcription factors and ribosomal genes by viral proteins. This study increases our understanding of the molecular interactions between ISAV and Atlantic salmon cells, and provides new targets for functional studies aiming to increase the resistance of Atlantic salmon stocks, leading to increased food security and fish welfare.