Identification of differentially expressed Atlantic salmon miRNAs responding to salmonid alphavirus (SAV) infection

Background MicroRNAs (miRNAs) control multiple biological processes including the innate immune responses by negative post-transcriptional regulation of gene expression. As there were no studies on the role(s) of miRNAs in viral diseases in Atlantic salmon, we aimed to identify miRNAs responding to salmonid alphavirus (SAV) infection. Their expression were studied at different time points post infection with SAV isolates associated with different mortalities. Furthermore, the genome sequences of the identified miRNAs were analysed to reveal putative cis-regulatory elements, and, finally, their putative target genes were predicted. Results Twenty differentially expressed miRNAs (DE miRNAs) were identified. The expression of the majority of these increased post infection with maximum levels reached after the viral load were stabilized or decreasing. On the other hand, some miRNAs (e.g. the miRNA-21 family) showed decreased expression at the early time points post infection. There were significant differences in the temporal expression of individual miRNA associated with different SAV isolates. Target gene prediction in SAV responsive immune network genes showed that seventeen of the DE miRNAs could target 24 genes (e.g. IRF3, IRF7). Applying the Atlantic salmon transcriptome as input 28 more immune network genes were revealed as putative targets (e.g. IRF5, IRF4). The majority of the predicted target genes promote inflammatory response. The upstream sequences of the miRNA genes revealed a high density of cis-regulatory sequences known as binding sites for immune network transcription factors (TFs). A high expression in the late phase could therefore be due to increased transcription promoted by immune response activated TFs. Based on the in silico target predictions, we discuss their putative roles as early promotors or late inhibitors of inflammation. We propose that the differences in expressions associated with different SAV isolates could contribute to their differences in mortality rates. Conclusions This study represents the first steps in exploring miRNAs important in viral-host interaction in Atlantic salmon. We identified several miRNAs responding to SAV infection. Some likely to prohibit harmful inflammation while other may promote an early immune response. Their predicted functions need to be validated and further studied in functional assays to fully understand their roles in immune homeostasis. Electronic supplementary material The online version of this article (doi:10.1186/s12864-017-3741-3) contains supplementary material, which is available to authorized users.


Background
MicroRNAs (miRNAs) are important regulators of gene expression at the post transcriptional level. Pre-miRNAs are transported from the nucleus and processed into two small mature miRNAs (5p and 3p) that usually are 20-24 nt in length. The mature miRNAs are incorporated into the miRNA-induced silencing complex (miRISC). Their function in miRISC are to recognize the target genes by partial base pairing to their transcripts [1,2]. Due to their role as regulators of gene networks, the majority of miRNAs are highly conserved in vertebrates, but there are also miRNAs that seem to be species specific [3][4][5].
MiRNAs are involved in regulation of almost all cellular processes including development, growth, maintenance of tissue-specific functions and apoptosis [6][7][8][9]. Some miRNAs are also important key regulators of normal immune function and immune responses [10,11]. Several studies in vertebrates have shown that some miRNAs play important roles in the virus-host interaction following virus infection. Viruses may themselves encode miRNAs (viral miRNAs) that may be involved in cellular reprogramming to promote a cellular environment favorable to viral life cycle. To achieve such an advantage to viral life cycle other viral factors may also repress or induce expression of host miRNAs (cellular miRNAs) [12]. On the other hand, altered expression of cellular miRNAs may also be a consequence of host responses to viral infection where the hosts utilize miR-NAs as part of its defense mechanisms [13,14]. Some of the affected cellular miRNAs have been shown to regulate genes in gene networks associated with innate immune response [15]. The modulation of cellular miR-NAs following virus infection could therefore be favorable for either the virus or the host depending on the function of the particular miRNA modulated and the direction of the change in expression [12,13,[16][17][18]. As demonstrated in some studies, detailed knowledge about miRNAs, their target genes and the disease mechanisms may potentially lead to the discovery of new diagnostic markers [19] and novel therapeutic approaches [20].
The role of miRNAs following RNA virus infection is less well studied in teleosts than in higher vertebrates. However, some recent studies of economically important teleost species have identified miRNAs responding to viral infection (e.g. [21][22][23][24]). We are particularly interested in miRNAs that affect disease development following viral infection in the economically important aquaculture species Atlantic salmon (Salmo salar L.). Recently, mature miRNAs and their corresponding miRNA genes were characterized in Atlantic salmon [25,26], and RT-qPCR methods to study the expression of single miRNAs have been validated [27]. So far, no studies of the putative role(s) of miRNAs in viral disease have been carried out in Atlantic salmon, but one study in rainbow trout (Onchorhyncus mykiss) identified miRNAs associated with Viral hemorrhagic septicaemia virus infection [28]. A better understanding of the particular miR-NAs showing a change in expression levels, the direction of change (reduced or increased expression) and prediction of their target genes could provide knowledge on the particular role the modulated miRNAs have and which gene regulatory networks are affected.
Salmonid alphavirus (SAV), a member of the genus Alphavirus, family Togaviridae, is a spherical, enveloped, single-stranded positive-sense RNA virus of approximately 60-70 nm in diameter with a 12 kb genome. Pancreas disease (PD) in Atlantic salmon and sleeping disease (SD) in Rainbow trout are caused by infection with SAV [29]. Six SAV subtypes (1-6) have been identified based on the nucleic acid sequences encoding proteins E2 and nsP3 [30]. Two of these, SAV2 and SAV3 have been reported in Norway [31]. Outbreaks of PD leads to increased mortality, reduced growth, affect feed conversion and eventually lead to low product quality. Histological examination usually reveals complete loss of exocrine pancreatic tissue, cardiac myocytic necrosis and inflammation, and degeneration and/or inflammation of skeletal muscle. PD outbreaks leads to large economic losses and is a major concern in the Atlantic salmon aquaculture. A recent challenge study using Norwegian SAV isolates has demonstrated that fish infected with SAV3 have higher mortality than those infected with SAV2. In addition, significant differences in mortality within the SAV3 genotype was observed. The cause for such differences in mortality is unknown [29,32].
The aim of this study was to identify cellular miR-NAs that are differentially expressed (DE miRNAs) following infection with SAV. The expression of DE miRNAs was studied in samples from fish infected by one SAV2 isolate (SAV2-i1) as well as samples from fish infected with two isolates of SAV3 with diverging mortality (SAV3-i4 and SAV3-i6) [32]. Expression of DE miRNAs was measured in each SAV group at three time points post virus challenge. This allowed us to identify miRNAs responding to SAV infection, to study their expression following viral infection by each of the three SAV isolates, and explore whether there were any associations between viral load and miRNA expression. The upstream genome sequences of the differentially expressed miRNA genes were studied in order to identify conserved cis-regulatory sequences associated with genes responding to innate immune responses, and, finally, the sequences of the mature DE miRNAs were used to predict their putative target genes. These in silico predictions could further contribute to reveal the particular role of the DE miRNAs in the host-virus interaction.

RT-qPCR analyses of SAV in cardiac tissue materials
The 98 cardiac tissue samples from controls and fish challenged with SAV2-i1, SAV3-i4 and SAV3-i6 collected at week 1-4 post onset of challenge experiment (poc) were successfully analysed for SAV by RT-qPCR. All eight control samples and the nine samples collected at week 1 poc were SAV negative (Ct > 45). Among the samples collected at week 2 poc, five of the nine samples from fish challenged with SAV2-i1 were SAV positive (55%), eight of the nine fish challenged with SAV3-i4 was SAV positive (89%) while all nine fish challenged with SAV3-i6 were SAV positive. All other samples collected at week 3 and 4 poc were SAV positive. Figure 1 illustrates the changes in viral load (assessed by measurements of Ct) in controls and all SAV groups at the four time points (including only the SAV positive samples at week 2 poc).
The results were in agreement with findings in Taksdal et al [32]. In Taksdal's study all controls were negative (Ct > 45), while at week 2 poc there were 63% SAV positive fish in the SAV2-i1 group, 88% SAV positive fish in the SAV3-i4 group and 99% SAV positive fish in the SAV3-i6 group while all fish from week 3 and 4 poc were SAV positive.
Comparison of Ct values in the SAV positive samples from each of the SAV groups indicated large increases in viral loads from week 2 to week 3 poc. The changes were less pronounced from week 3 to week 4 poc where the obtained Ct values indicated no change in virus load in the SAV3-i4 group, a small increase in the SAV2-i1 group and a small decrease in the SAV3-i6 group. The comparisons of Ct values in the SAV positive samples showed no significant differences at week 2 poc. However, a number of samples were SAV negative at week 2 poc, and most of these were from the SAV2-i1 group. This indicated that the time needed for infection by cohabitation was somewhat longer in the SAV2-i1 group. Ct measurements at week 3 poc did show significant differences (p < 0,05) between all three SAV groups (Fig. 1). These Ct-measurements at week 3 poc indicate an approximately ten times higher viral load in the SAV3-i6 group than in the SAV3-i4 group, while there was about seventeen times higher viral load in the SAV3-i4 group than in the SAV2-i1 group. Measurements at week 4 poc revealed smaller differences in Ct values, and only the difference between the SAV2-i1 group and the SAV3-i4 group was significant.

Identification of DE miRNAs following SAV infection Deep sequencing and identification of DE miRNAs by DESeq2 analysis
Total RNA was successfully extracted from 98 cardiac samples, and six of these were subsequently used for small RNA deep sequencing analysis. These six samples were from two groups; three healthy individuals (control group samples) and three individuals infected with SAV3-i4 (the SAV isolate associated with highest mortality rate in Taksdal et al [32]) collected at week 4 poc (see Materials and Methods). The samples were successfully Fig. 1 Shows the increase in viral loads assessed by RT-qPCR. The mean Ct-values are given for each of the three groups challenged with SAV2-i1 (blue), SAV3-i4 (dark orange) or SAV3-i6 (green) at week 1-4 poc. The controls are the SAV negative healthy fish sampled at the initiation of the challenge trial (Ct values > 45). At week 1 poc all samples were SAV negative. At week 2 poc the plot show results from only those samples that were SAV positive. All fish sampled at week 3 and 4 poc were SAV positive subjected to next generation sequencing using Illumina Genome Analyzer IIx sequencing platform as described in Materials and Methods. Quality control by FASTQC showed good sequence quality and similar size distributions in all samples. The results from deep sequencing including number of size filtered and adapter trimmed reads, and reads mapped as ssa-miRNAs in each sample are summarized in Additional file 1. Two out of the six samples had a number of reads recovered as miRNAs less than one million. We can not rule out that this may have affected the sensitivity to detect smaller fold changes in low expressed miRNAs. All results have been submitted to the NCBI SRA data base (Genbank SRP048613). Genbank accession numbers to the individual samples are SRR3928621 and from SRR3928623 to SRR3928627.

Validation of DESeq2 results by RT-qPCR analysis
DESeq2 is a robust and widely used statistical tool for quantitative analysis of deep sequencing data. However, it is still common practice to validate any findings by use of an RT-qPCR method that can target each of the differentially expressed mature miRNAs individually. Such independent validations also allow for measurements of their expression in a larger population sample. Thus, for the purpose of validating our results we developed RT-qPCR assays to measure expression of the different DE miRNAs identified by DESeq2 analysis. Two miRNAs were not included in this validation (ssa-miR-462b-5p and ssa-miR-146a-3-3p) as two other miRNAs (ssa-miR-462a-5p and ssa-miR-146a-3p) with very similar mature sequences from the same miRNA families were included. RT-qPCR assays were successfully developed for 17 DE miRNAs. The miRNA specific primer sequences along with measurements of assay performances (R-squared and efficiency) for all these RT-qPCR assays are given in Additional file 2. One miRNA assay (ssa-miR-92a-3-5p) did not pass the validation of the RT-qPCR assay due to non-specific amplification and was not analysed further.
The comparative Ct method (ΔΔCt-method) [35] were used to measure the relative difference in miRNA levels between the healthy controls (n = 8) and SAV3-i4 positive individuals (n = 9) at week 4 poc. The results are summarized in Table 1. The log2 fold changes (ΔΔCt) as well as adjusted p-values from each miRNA analysed by RT-qPCR is given. The results from RT-qPCR analysis agreed very well with the results from DESeq2 analysis for 15 of the 17 miRNAs tested. They all showed similar changes in expression which was significantly different compared to the controls (significant threshold p = 0.05). The two remaining miRNAs (ssa-miR-155-5p and ssa-miR-734-3p) did reveal differences in the same direction as those from deep sequencing, but they were less pronounced and not significant. Considering the good agreement between DESeq2 analysis and RT-qPCR for the other 15 miRNAs it may be possible that the lack of correlation between the two methods for these two miRNAs was due to suboptimal RT-qPCR assay performance. An efficiency in the lower range for these two RT-qPCR assays (Additional file 2) support that the RT-qPCR was the method providing the less accurate measurements of fold changes for these two miRNAs. Figure 2 shows the results (log2 fold changes) for the 17 miRNAs that were analysed by both deep sequencing (DESeq2) and RT-qPCR (ΔΔCt-method). The figure illustrates the good agreement between the relative expression changes detected by the two methods. Spearman's nonparametric test for correlation showed ρ = 0.83 (p = 0.001). The good correlation between the methods also demonstrates that the RT-qPCR assays were suitable for measurements of miRNA expression in all SAV positive samples collected (three SAV groups, different time points).
RT-qPCR analysis of miRNA expression in SAV2-i1, SAV3-i4 and SAV3-i6 groups The individuals in the three groups challenged with either SAV2-i1, SAV3-i4 or SAV3-i6 were successfully The results from all time points showed that there was only one miRNA (ssa-miR-734-3p) that could not be verified by RT-qPCR as significantly different from the controls in any of the SAV groups. A majority of the miRNAs showed increased expression across the time points with the largest increases between week 3 and 4 poc, and the higher level of expression at the latest time point (week 4 poc) (Fig. 3). Viral load, on the other hand, showed the largest Fig. 3 Shows the relative changes in miRNA expression at all time points for each of the SAV groups for the following miRNAs: ssa-miR-462a-5p, ssa-miR-731-5p, ssa-miR-146a-3p, ssa-miR-146b-5p, ssa-miR-223-5p, ssa-miR-181c-5p, ssa-miR-155-5p and ssa-miR-7132-3p. The corresponding mature miRNAs from the miRNA gene cluster miRNA-462/731 (ssa-miR-462a-3p and ssa-miR-731-3p) showed very similar results, and for simplicity reasons these are not given in the figure, but in Additional file 3 increases from week 2 to week 3 poc and modest changes from week 3 to week 4 poc (Fig. 1). Comparing the changes in miRNA expression with the changes in viral load, the expression of most of the miRNAs seem to mirror the change in viral load as if the large increases to the highest expression level at week 4 poc were a time delayed response to the large increase in viral load from week 2 to week 3 poc.
There were, however, differences in magnitude of the increases between groups. There were also differences in the direction of change between different miRNAs as some of the miRNAs revealed other changes than increases only (e.g. the miRNA-21 family and ssa-miR-1338-3p) (Fig. 4). These differences between groups are interesting as fish infected with either SAV2-i1, SAV3-i4 or SAV3-i6 have revealed differences in mortality rates [32].

miRNAs with decreased expression at early time points
The miRNAs from the miRNA 21 family (ssa-miR-21a-1-3p, ssa-miR-21b-3p and ssa-miR-21a-2-3p) and miR-1338-3p revealed a decreased expression across the early time points in the SAV2-i1 group. This was in contrast to the two other groups that showed unchanged or increased expression compared to controls at week 2 and 3 poc (Fig. 4). Interestingly, the SAV3-i4 group with the highest mortality rate also showed a higher expression of the three miRNAs from the miRNA 21 family than the SAV3-i6 group.
The changes in expression in ssa-miR-1338-3p were small in all SAV groups, but only the SAV2-i1 group showed a reduced expression, and the SAV2-i1 group was significantly different from SAV2-i6 at week 3 poc (Fig. 4). Ssa-miR-7132-5p and ssa-miR-2188-3p did both show a decreased expression at one or more time point in the SAV3-i4 group (Fig. 4). Ssa-miR-7132-5p showed a significant decrease at week 2 poc while ssa-miR-2188 revealed a decrease in expression across all time points in the SAV3-i4 group that was significantly different from both the controls and the SAV3-i6 group at week 4. Ssa-miR-734-3p showed a small decrease that was not significantly different to controls in any of the SAV groups (Additional file 3). However, this miRNA may have been more accurately measured in the DESeq2 analysis (4x times decrease, Table 1).

Prediction of miRNA target genes
The four miRNA target prediction tools RNAhybrid, TargetSpy, PITA and Miranda [36,37] were applied to predict the target genes. Two different approaches were used for selecting the genes used as input. First, we used a knowledge based approach where immune responsive genes known to respond to SAV infection [38,39] were used as input. A few immune relevant genes known to be targets of the identified DE miRNAs in other vertebrates [40][41][42][43][44][45] were also included. Thus, a total of 50 genes (Additional file 4) were analysed in this first approach. Then we performed a second target gene prediction without any prior selection based on function or response to SAV. Applying such a holistic approach all known Atlantic salmon transcripts in the Refseq database of Genbank were considered putative targets and used as input. The latter approach would predict a large number of false positives [46], but could also potentially identify immune response related genes not previously identified as differentially expressed following SAV infection that are true targets.
The results from analyzing genes known to respond to SAV infection showed that 24 of the genes were predicted as putative targets by at least three of the prediction tools applied. The genes along with the mature miRNAs predicted to target each of these genes are shown in Table 2. Thirteen of the predicted target genes revealed multiple target sites to several different DE miRNAs in their 3'UTR's. Three of these genes, neutrophil cytosolic factor 1, myeloperoxidase and IKKa, were predicted to have target sites to five of the DE miRNAs. One of the miRNAs (miR-223-5p) predicted to target IKKa has also been reported as targeting this gene in other vertebrates [43].
The target gene predictions applying all Atlantic salmon mRNAs (Refseq, Genbank) as input (holistic approach) generated matches to a total of 1335 transcripts. Each DE miRNA was predicted to target several transcripts with numbers ranging from 30 transcripts (miR-7132-3p) to 389 (miR-181c-5p). The complete results from this analysis, showing the predicted target transcripts along with the matching mature miRNAs are given in Additional file 6.
A subset of 28 genes predicted as targets applying the holistic approach were immune response network genes ( Table 3) that was not discovered using the SAV responsive set of genes as input in the target gene predictions. Although an association between these genes and SAV has not been reported, the in silico prediction analysis indicates that these genes may also be among the genes targeted by the DE miRNAs. The majority of these transcripts encode proteins with roles as activators (IRF1, IRF4, IRF5) and effectors (chemokine receptors) in the immune responses while a few of them are inhibitors (e.g. NFKappaB inhibitor) ( Table 3).

Genome sequence analysis of DE miRNA genes to identify putative cis-elements
The genome sequence of the miRNA precursors and their location within the Atlantic salmon genome was described in Andreassen et al. [25], and an updated version of the Atlantic salmon genome sequence was recently uploaded in Genbank (Genbank # GCF_000233375.1). Together these sources were utilized to identify and analyze the upstream sequences of the DE miRNA genes. The 20 short mature DE miRNAs could be processed from 25 miRNA genes as some of the mature miRNAs could originate from highly similar duplicated miRNA genes. While the upstream genome sequence was not available for one of the genes (miRNA-155-2), the genome sequences from the other 24 miRNA genes could be included in the genome sequence analysis. They were all analyzed regarding presence of the consensus sequences of eight conserved promoter motifs known to bind immune-responseactivated transcription factors (Additional file 7). The presence of such elements could indicate that their transcription was enhanced by transcription factors that are part of the innate immune response gene network.
The results from the search for response elements in the immediate upstream sequences of the DE miRNA genes are summarized in Table 4. Ten of the genes revealed TATA-box motifs, a finding in agreement with Long et al. [47]. More interestingly, two or more of the eight promoter motifs were observed in the immediate upstream sequences of all genes except miRNA-734 and miRNA-2188 that did not reveal any of the motifs in their upstream sequences. The miRNA-462a/731 gene cluster as well as the single miRNA-462b gene revealed the "classical" ISRE motif (ISRE-IRF1 in Table 4), the PU.1 motif and the GAS motif which have also been identified in the upstream sequence of the miRNA-462a/ 731 gene cluster in rainbow trout [28]. In addition, there were three more out of the eight promoter motifs present upstream of this gene cluster ( Table 4). The "classical" ISRE motif was also present in the miRNA-146b and miRNA-155-1 genes along with three or more of the other motifs. In agreement with findings that miRNA-155 transcription is regulated through binding to NF-kappaB and ETS2 motifs [48], these motifs were also present in the Atlantic salmon miRNA-155-1 gene. In summary, the genome analysis showed that there was a high density of motifs known to bind transcription factors activated by the innate immune system in the upstream genome sequences of the DE miRNA genes. This indicates that transcription of many of the DE miRNA genes may be enhanced by some of the key transcription factors of the innate immune system. The two miRNA genes, miRNA-734 and miRNA-2188, that did not have any matches to the conserved motifs in their upstream sequences were the ones that revealed a decrease of their mature miRNAs following infection (Table 1).

DE miRNAs may participate in negative feedback loops in the host innate immune system
The quantitative analysis of the results from the deep sequencing followed by RT-qPCR validations identified a subset of mature miRNAs that are differentially expressed following SAV infection in Atlantic salmon (DE miRNAs). The relative expression of the DE miR-NAs was subsequently measured across three time points post onset of challenge in fish infected with either SAV2-i1, SAV3-i4 or SAV3-i6 (Figs. 3 and 4).
The DE miRNAs were predicted to target a number of key genes of the innate immune system. Some are among the important activators of the innate immune responses like the ATP-dependent RNA helicase DHX58 as well as the transcription factors IRF3 and IRF7. Other innate immune response network genes that were predicted as targets were interferon gamma (IFN-γ) and interferon stimulated genes like viperine (RSAD2), C-C Table 3 Predicted target genes with immune response related function from in silico analysis using the Atlantic salmon transcriptome as input (holistic apporach)  (Table 2). These genes are all known to respond to SAV infection [38,39]. Extending the search for putative target genes to the Atlantic salmon transcriptome identified 28 additional immune network genes as putative targets ( Table 3). The majority of the putative target genes act as activators or effector molecules in the innate immune system pathway. If they are true target genes they are negatively regulated by DE miRNAs. The role of the DE miRNAs in the host-virus interaction would not only depend on what target genes they regulate, but would also depend on what time post viral infection their expression changes as well as the direction of change. Several recent reviews [11,[49][50][51] suggest that some miRNAs targeting activator and effector transcripts could have a role in promoting early immune response if their expression was decreased at the early stage of an immune response. In the late phase of inflammation, some miRNAs that target activator and effector transcripts have been suggested to have another role. In this phase increased expression of these miRNAs may help prevent pathological inflammatory processes. In this role the miRNAs stabilize and prevent further increase of activator and effector transcripts, and both the targets and miRNAs would have increased expression compared to the normal cellular state. The expression of some miRNAs believed to negatively regulate the immune response in the late phase of infection seem to be induced by immune network transcription factors. In such cases the miRNAs participate in negative feedback loops as their expression is enhanced by their targets [11,49].
The viral loads increased across the three first time points (week 1, 2 and 3 poc) in all the SAV groups, but seemed to decline at week 4 poc in SAV3-i6. The viral loads in the two other groups also seemed to have reached their maximum level at week 4 poc as the increases were much less pronounced from week 3 to week 4 poc (Fig. 1). These observations are in agreement with findings in Taksdal et al. [32] where the viral loads declined in all SAV groups after week 4 poc. Ten of the miRNAs (e.g. the eight miRNAs in Fig. 3) showed minor increases across the first time point (week 2 poc), but larger increases towards the last time point (week 4 poc). The findings from measurements of viral load showed that all fish had been infected for more than a week at the last time point (week 4 poc, Fig. 1). Comparing the increases in viral loads with the change in miRNA expression the latter appear as a time delayed response to viral infection. This could indicate that the purpose of the increases, resulting in a high expression of the DE miRNAs more than a week post maximum measurements of viral loads, were to negatively regulate the expression of inflammatory mediators to prevent excess inflammation [50,51]. Many DE miRNAs showed a time delayed increase of the DE miRNAs (Fig. 3) and they were predicted to target immune response activators. This could suggest that their role were to prevent a prolonged or elevated immune response that would be harmful to the host [11,49]. Table 3 Predicted target genes with immune response related function from in silico analysis using the Atlantic salmon transcriptome as input (holistic apporach) (Continued) C-X-C chemokine receptor type 3 Receptor for C-X-R chemokines, regulates biological processes such as inflammation, immunity, and would healing   The initial screening to identify DE miRNAs was carried out in materials from the latest time point post challenge (week 4 poc). It would therefore be expected that the majority of the miRNAs identified would be those involved in preventing pathological inflammation and tissue damage. Since we did not include earlier time points in the initial screening for DE miRNAs, we cannot rule out that other miRNAs are important regulators early in the infection. The role of such early responding miRNAs could be to promote an immune response. The decreased expression at the early time points of a few DE miRNAs identified in this study, e.g. the mature miRNAs from the miRNA 21 family (Fig. 4), could indicate that they may have such a role where they contribute to promote the early immune response. The downregulation at the early time points in the SAV2-i1 group of the mature miRNAs from the miRNA 21 family (Fig. 4) would be expected to contribute to an immediate increase in the immune response as downregulation of the miRNAs would lead to a proportional downregulation of the translational repression of the predicted target gene IRF3 ( Table 2). The downregulation of another miRNA, ssa-miR-734-3p, was detected in the DESeq2 analysis, but could not be validated by RT-qPCR in our study (Table 1). However, if this miRNA is downregulated in the very early phase of infection this would also lead to an increased immune response as it was predicted to target viperine and interferon gamma ( Table 2).
Analysis of the genomic sequences spanning the DE miRNA genes revealed a high density of response elements that may bind transcription factors associated with an innate immune response. Similar findings have been reported by Schyt et al. [28] for the teleost specific miRNA 462/731-gene cluster analysed in rainbow trout. Other conserved miRNA genes, like miRNA-146, miRNA-181 and miRNA-155 are also known to be INF-regulated in vertebrates with increased expression upon stimulation [11]. If the expression of the DE miRNAs identified in our study is enhanced by binding transcription factors activated by the immune response, this means that many of the DE miRNAs are the effector molecules participating in negative feedback loops. A recent study of miR-731-5p in Japanese flounder [44] demonstrated such a negative feedback where miR-731-5p was shown to target the transcription factor IRF7. The analyses of target genes in our study showed that, not only was ssa-miR-731-5p predicted to target the salmon ortholog of IRF7, but it was also predicted to target the ATP-dependent RNA helicase DHX58 (RIG-I-like 3). The RIG-I-like receptors are believed to be upstream activators of IRF7 [39]. Targeting not only IRF7, but also the upstream activator of IRF7, an increased expression of ssa-miR-731-5p would, as described in Zhang et al. [44], negatively regulate the immune response.
Together, the target gene predictions and the analysis of upstream genome sequences of the miRNA genes suggest that many of the DE miRNAs are important participants of the immune system gene pathways and some may act in negative feedback loops to maintain immune homeostasis. However, as demonstrated in our study, the expression of these DE miRNAs may differ depending on what kind of virus subtype that causes the immune response. Even virus infection caused by virus isolates belonging to the same subtype (SAV3-i4 and SAV3-i6) resulted in significant differences in expression of most of the DE miRNAs.

Different expression patterns of the DE miRNAs may help explain differences in mortality
The high expression at the late time point and the predicted target genes suggest that the ten miRNA (eight of these shown in Fig. 3) may act to prevent excess inflammation. The expression of these miRNAs (ssa-miR-462a-5p and 3p, ssa-miR-731-5p and 3p, ssa-miR-146a-3p, ssa-miR-146b-5p, ssa-miR223-5p, ssa-miR-155-5p, ssa-miR-181c-5p and ssa-miR7132-3p) were higher in the fish infected with SAV3-i6 (Fig. 3, Additional file 3) than in the other SAV3 isolate (SAV3-i4) at week 4 poc. This observation is interesting as Taksdal et al. [32] reported that there was a difference in mortality between the two isolates of SAV3 with a higher mortality in the SAV3-i4 group than in the SAV3-i6 group that could not be explained by their results. A higher expression of the ten DE miRNAs in the SAV3-i6 group than in the SAV3-i4 group could lead to a more pronounced down regulation of the inflammatory process in the SAV3-i6 group. The effect of this could be less pathological inflammation in the cardiac tissue in the following weeks that might result in a reduced mortality in the SAV3-i6 group compared to the SAV3-i4 group. The observation that ssa-miR-2188-3p has a significantly lower expression in the SAV3-i4 group than in the SAV3-i6 group (Fig. 4) could also contribute to a higher mortality in the SAV3-i4 group. The two predicted targets of this miRNA promote inflammatory response and cell death (Table 2). Thus, a prolonged high expression of these predicted target genes would not benefit the host. The observed differences in expression of DE miRNAs between the SAV3-i4 and SAV3-i6 group are therefore in agreement with, and offer an explanation to, the observed differences in mortality between fish infected with these isolates of SAV3 subtype.
Fish challenged with SAV2-i1 was reported with a significantly lower mortality than those challenged with the two SAV3 isolates [32]. The SAV2-i1 group did, like the SAV3-i6 group, show a significantly higher expression at the latest time point in some miRNAs compared to the SAV3-i4 group (ssa-miR-462a-5p, ssa-miR155-5p and ssa-miR-7132-3p, Fig. 3). However, in other miRNAs they revealed similar expression levels as those from the SAV3-i4 group, and comparing the SAV2-i1 group to the SAV3-i6 group all ten miRNAs (Fig. 3) showed a much lower expression in the SAV2-i1 group. The miRNA expression at the late time point in infection (week 4 poc) in these miRNAs could therefore not explain a mortality rate that was lower in the group infected by SAV2-i1 than in both of the groups infected with the SAV3 subtypes (SAV3-i6 and SAV3-i4). However, the expression of other miRNAs at the early time points may offer an explanation to these differences in mortality rate. There was a decreased expression in the miRNAs from the miRNA 21 family at week 2 and 3 poc in the SAV2-i1 group (Fig. 4). In contrast, there was an increase in expression across the same time points in the groups challenged with both isolates from the SAV3 subtype. These differences in expression, a reduction of miRNAs from the miRNA 21 family in the SAV2-i1 group versus an increase of these miRNAs in both of the SAV3 groups would lead to a faster and much more pronounced downregulation of the negative repression of the target gene in the SAV2-i1 group. The predicted target gene is IRF3, a constitutively expressed key activator of the innate immune system [52]. Thus, the differences in miRNA 21 family expressions could contribute to a faster and more effective early immune response in fish infected by SAV2-i1 and lead to a lower mortality in this group.
Although the observed differences in expression between DE miRNAs may offer some explanations to differences in virulence and mortality observed in fish infected by different SAV isolates they rely on the correctness of the target gene predictions. Thus, an obvious way forward would be to validate the target genes of the identified DE miRNAs and further explore the miRNAtarget gene interactions by experimental approaches.

Conclusion
This is the first miRNA profiling study in Atlantic salmon that has identified and characterized cellular miR-NAs with modulated expression following virus infection indicating that they are important participants in the virus-host interaction. The target gene predictions and analysis of upstream genome sequences suggests that many of the miRNAs may regulate innate immune responses as effector molecules in negative feedback loops. Our findings represent a first step in exploring and suggesting a role for Atlantic salmon miRNAs in virus-host interactions. Further functional studies will help to disclose the particular role of these miRNA genes, and may contribute more knowledge on how these miRNAs act in concert with other regulator molecules to fine-tune the innate immune system responses.

Materials
A challenge trial was carried out where healthy Atlantic salmon (smolt) were exposed to SAV by cohabitating with salmon shedders injected with virus (carrier fish). The cohabitant challenge experiment and sampling of tissue is described in detail in Taksdal et al. [32]. In short, the fish used for the challenge experiment was unvaccinated, clinically healthy Atlantic salmon from VESO Vikane hatchery. Health conditions was also confirmed by examining ten fish for common fish viruses. They all tested negative for piscine myocarditis virus (PMCV), piscine reovirus (PRV), infectious pancreatic necrosis virus (IPNV) and SAV. The fish was challenged with 3 SAV-isolates derived from independent outbreaks of Pancreas Disease in Norwegian fish farms, either marine SAV2 (SAV2-i1) or one of two isolates of SAV3; SAV3-i4 and SAV3-i6 [32].
Cardiac tissue samples were collected from the cohabitating (non-injected) fish. Co-habitant samples were collected at four time points; at the initiation of the experiment (controls), at week one, two, three and four post onset of the challenge experiment. Eight fish collected at the initiation of the experiment were used as healthy controls, while nine fish, three from each SAV group was collected at week 1 poc. Nine fish challenged with either SAV2-i1, SAV3-i4 or SAV3-i6 were sampled at week two, three and four poc. Together, cardiac tissue from 98 fish was sampled. All samples from controls, week one, two, three and four poc were analysed to measure viral load by RT-qPCR. All samples from week 1 poc, four samples from fish challenged with SAV2-i1 and one sample from fish challenged with SAV3-i4 at week 2 poc were SAV negative and were removed from further analysis of miRNA expression. The samples from the remaining individuals at week 2 poc, and all fish from week 3 and 4 poc, tested positive for SAV (n = 78). The eight control samples and the SAV positive samples from other time points (total of 84 samples) were all included in our study of miRNA expression (RT-qPCR of miRNA).
Three healthy controls and three samples collected at week 4 poc from the SAV3-i4 group were used for deep sequencing and DESeq2 analysis. Fish infected with SAV3-i4 were choosen as this was the SAV isolate associated with the highest mortality rate [32]. The late time point (week 4 poc) was choosen as the virus loads were high and 100% of fish from all three groups were infected at this time point. The experimental study was approved by the National Research Authority in Norway (NARA). All salmon used for sampling in the experiment were euthanized according to standard protocols approved by the Norwegian Food Safety Authorities prior to sampling.

Small RNA isolation
Sampling of cardiac tissue was carried out as described in Taksdal et al. [32]. Samples were fixed in RNA later (Life technologies, Carlsbad, USA) immediately after sampling. Total RNA was isolated from each of the samples by use of the mirVana miRNA isolation kit (Ambion, Life Technologies, Carlsbad, CA, USA) following the manufacturer's protocol. The RNA concentration and purity was determined using Nanodrop spectrophotometer (Thermo Fisher Scientific, WI, USA) following the manufacturer's protocol. The concentration of total RNA in samples ranged from 21-255 ng/μl (total volume 100 μl) and with A260/280 ratios above 1,9. The integrity of the RNA in samples used for deep sequencing was analyzed using the Agilent Bioanalyzer system (Agilent, CA, USA). The RIN values were from 6.8 to 10 (Additional file 1).

RT-qPCR detection of SAV
SAV in the cardiac samples was detected by RT-qPCR with a Qiagen Onestep RT-PCR kit using TaqMan assay with primers and probe targeting the nsP1 gene as described in Hodeland and Endresen [53]. Viral load (Ct values) was measured by RT-qPCR in all samples from controls, week one, two, three and four poc (n = 98).

Library construction, deep sequencing and quantitative analysis
The library construction was performed at the Norwegian Genomics Consortium's genomics core facility. The Illumina TruSeq Small RNA Library Preparation Kit (Illumina, San Diego, USA) was used in the preparation of the libraries as described by the manufacturer with 1 ug total RNA input. After adapter ligation and cDNA synthesis the products were purified on a gel and the fractions between 145 and 160 bp were used for sequencing. The small RNA libraries constructed from six samples (three from controls and three from the SAV3-i4 group sampled at week 4 post onset of challenge experiment) were successfully subjected to next generation sequencing using Illumina Genome Analyzer IIx sequencing platform as described in Andreassen et al. [25]. Reads quality were assessed using FastQC toolkit (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Cutadapt [54] was used for trimming of adapter sequences from raw sequence reads and removal of adapter-only sequences (5'TGGAATTCTCG GGTGCCAAGGAACTCCAGTCAC 3'). Finally, an additional filtering of reads outside of 18-25 nucleotide range was applied on all samples. All deep sequencing samples have been submitted to the NCBI Sequence Read Archive database (SRA) and accession numbers are given in Additional file 1.
Next, the reads from each of the six samples were aligned to a set of miRNAs that consisted of all known Salmo salar mature miRNAs [25,26] using Novoalign (http://www.novocraft.com/products/novoalign/). Custom made scripts were used to report number of successfully assigned reads to each miRNA. The reads mapping with edit distance 1 or less to the mature reference sequences were considered true mature miRNAs. Differentially expressed miRNAs were identified with DESeq2 package [33]. Rows with less than 2 reads for each condition were discarded from the analysis. A threshold of p-adjusted value <0.1 (adjusted according to Benjamini-Hochberg procedure) was applied to report differentially expressed miRNAs. The two groups compared were controls (three samples) and SAV3-i4 infected individuals (three samples) from week 4 poc.

cDNA synthesis and RT-qPCR of miRNAs
The miScript assays were used for cDNA synthesis (miScript II RT Kit) and qPCR (miScript SYBR Green PCR kit) as described by the manufacturer (Qiagen, Hilden, Germany). The reverse primer (universal primer) was provided with the miScript qPCR kit while the forward primer was custom designed and specific to each of the mature miRNA gene assays. The sequence of the forward primers were designed based on the mature sequence of each of the miRNA genes to be amplified [25]. All forward primers were purchased from Sigma Aldridge. They were purified by desalt only and provided as liquid solution of 100 μM. The primers were diluted to 10 μM for use in each of the qPCR assays. All forward primer sequences are given in Additional file 2. The qPCR analysis was run on an Mx3000p (Stratagene, Agilent Technologies, USA). The qPCR reaction mixture consisted of 12.5 μL 2x Quantitec Syber Green Master Mix, 2.5 μL 10x miScript Universal Primer, 2.5 μL of 10 μM forward miRNA gene specific primer, 5 μL Rnase free water, and 2.5 μL cDNA (template). Amplification was performed in 96-well plates. The following program was used: one thermal cycle at 95°C for 15 min followed by 40 cycles of 94°C for 15 sec, 55°C for 30 sec and 70°C for 30 sec. The cybergreen assay module was used for qPCR analysis. This module includes a final melting point analysis that follows the 40 cycles of quantitative PCR. A manual inspection of the plots from the melting point analysis of all miRNA gene assays tested was carried out to verify that forward primers were specific (data not shown).
In all cases of measurements of miRNA expression by RT-qPCR, the samples that were compared were run in the same set-up (within run measurements). Three miR-NAs ssa-miR-25-3p, ssa-miR-455-5p and ssa-miR-17-5p were used as reference genes [27]. Results (Ct-values) of target miRNAs were normalized by use of the geometric mean from the measurements of the three reference genes [55]. The relative increase in expression in target tissue of each miRNA was calculated using the principle of the ΔΔCt-method [35]. Applying this method the average relative difference in target tissue was calculated as the difference between mean normalized Ct from target tissue and mean normalized Ct from reference tissue. To test whether the relative differences in expression revealed by the ΔΔCt-method between controls and SAV infected groups was significant we performed Mann-Whitney tests using the statistical package SPSS. Likewise, Mann-Whitney tests were also used to test for significant difference in expression of miRNAs at the same time point between the three SAV groups. P-values were adjusted by use of the Benjamini-Hochberg procedure taking into account there were 153 comparisons between SAV groups at different time points (P-adjusted values equal to a 0.05 significance level). Test for correlation between estimates of relative differences by DESeq2 and by RT-qPCR of the 17 mature miRNAs analysed by both methods was carried out by use of Spearman's nonparametric test for correlation (SPSS package).

miRNA target gene predictions
The mature sequences from the DE miRNAs were used as input in the target gene prediction analysis. The 3'UTRs from two sets of genes were used as input for target genes. One set was genes reported as responding to SAV infection and associated with immune response or response to viral infection [38,39]. In addition, the Atlantic salmon orthologs to some genes that are key regulators of immune responses and reported in other species as targets of the same miRNAs as those identified in our study were included in this set of genes [40][41][42][43][44][45]. Applying these selection criteria based on prior knowledge about their relevance to SAV infection and immune response (knowledge based approach) a total of 50 genes were analysed to detect target site matches (Additional file 4). The other set of genes used as input were 3'UTRs from all Atlantic salmon mRNA transcripts in the Refseq database in Genbank (https:// www.ncbi.nlm.nih.gov/). This holistic approach could possibly identify target genes among immune network genes not previously known to respond to SAV infection.
The in silico analysis to identify target genes were carried out by use of four different target gene prediction tools; RNAhybrid, TargetSpy, PITA and Miranda. RNAhybrid analysis was performed with conditions helix constraint 2-8 and no G:U in seed allowing only target genes that had perfect "seed" matches to be detected [36]. Minimum free energy threshold for RNA hybrids was set to -18 kcal/mol to retrieve results (target site matches) from RNA hybrids that had a high stability [56]. RNA hybrid was used to generate the figure in Additional file 5. The target site analysis with the tools TargetSpy, PITA and Miranda was carried out applying default settings using the sRNAtoolbox (http://bioin-fo5.ugr.es/srnatoolbox/mirconstarget/) [37].
Functional motif analysis of the genome sequences spanning the DE miRNA genes The upstream genome sequences of all DE miRNA genes were retrieved by using precursor sequences [25] as input in BLASTn searches against the current version of the Salmo salar genome sequence in Genbank (Genbank # GCF_000233375.1). Several of the mature miR-NAs could originate from duplicated miRNA genes. Thus, the total number of miRNA genes analysed was 24. One thousand base-pairs of the upstream sequence and 200 base-pairs of the downstream sequence from the mature 5p sequences of each the DE miRNA genes were analysed regarding presence of conserved response elements known to bind key transcription factors from the innate immune system gene network. Consensus sequences of vertebrate response elements were retrieved from the MotifMap database (http://motifmap.ics.uci.edu/). In addition, Salmonid specific motifs were retrieved from Robertsen [57] and Schyth et al. [28]. The consensus sequences of the nine response elements used for upstream genome analysis are given in Additional file 7. The consensus sequences of conserved immune related response elements (ISRE-IRF, NF-kappaB and ETS2), the TATA-box consensus sequence as well as the response elements GAS and PU1 characterized in a study of the rainbow miRNA 462/731 gene cluster [28] are all detected with the general consensus sequences used in our genome sequence analysis. ISRE and NFkap-paB motifs in zebrafish are also detected using these motifs [58,59] indicating that the consensus sequence applied detect these motifs in teleost fish (including Atlantic salmon). Presence of any motif was detected by use of Sequencher software (Gene Codes Corporation, Michigan, USA).