The virus isolate Glesvær/2/90 (Dannevig et al., 1995; Markussen et al., 2008; Merour et al., 2011) was used to verify the sequencing accuracy of the amplicon sequencing method. This virus isolate has previously been sequenced twice by two different research groups using Sanger sequencing (Markussen et al., 2008; Merour et al., 2011). The genome of the Glesvær/2/90 isolate was sequenced from both tissue samples and cell culture supernatants and compared to the two published genomes.

For this study, 12 virus samples from the northern part of Norway in 2013 and the subsequent two years were chosen for sequencing. Tissue homogenates were inoculated and cultivated in Atlantic Salmon Kidney (ASK, ATCC CRL-2747) cells, as described in the WOAH Manual (World Organisation for Animal Health, 2021).

A set of candidate ISAV-specific primers was generated with the program Primal Scheme specifying 400 bp as the target amplicon size (Quick et al., 2017). The specificity of each primer was assessed manually by aligning the primers to a multiple sequence alignment of publically available European ISAV sequences in the NCBI database (Benson et al., 2013). Some of the primers were degenerated or redesigned to match the variety of virus strains. The degenerated bases were written in the IUPAC ambiguity code (Cornish-Bowden, 1985; Johnson, 2010). Primers for the indexing PCR (PCR-2) were described by de Muinck et al. (2017). A series of pilot experiments were performed to assess primer performance relating to uniform sequencing depth for each amplicon. Some primers were redesigned based on the analysis of relative sequencing depth. Four multiplex PCR reactions with 10 primer pairs were employed to amplify the whole ISAV genome. Two short ISAV-specific primers to be used in the cDNA synthesis step were designed in an attempt to better capture the segment ends (Supplementary Table 1).

A MagNA Pure 96 system (Roche Diagnostics, Oslo, Norway) was used for nucleic acid extraction as described by the supplier. Briefly, 400 μl of MagNA Pure LC RNA Tissue Lysis Buffer (cat. no. 03604721001, Roche Diagnostics) was added to 150 μl of clarified cell culture supernatant or tissue homogenate before total RNA was extracted using the RNA Tissue FF Standard LV protocol with the MagNA Pure 96 Cellular RNA Large Volume Kit. Samples were eluted using 100 μl elution buffer and stored at −80°C.

Reverse transcription of total RNA was performed using SuperScript IV Reverse Transcriptase (cat. no. 18090050, Fisher Scientific AS, Oslo, Norway). The RT-reaction was primed with 2,500 nM random hexamers (cat no N8080127, Fisher Scientific AS) and 50 nM ISAV-specific primers, ISAV_RT_S5, ISAV_RT_C (TAG Copenhagen A/S, Denmark, Supplementary Table 1) and 5 μl of total RNA extract at 90°C for 5 min before Superscript IV reverse transcriptase was added and the incubation continued at 55°C for 30 min. The PCR1 amplification with gene-specific primers (Supplementary Table 1) was performed in four separate multiplex PCR reactions with multiplex primer pools 1–4 using the following reaction mixture: 11.25 μl of nuclease-free water, 5 μl of Q5 reaction buffer, 0.5 μl of 10 mM dNTPs (New England Biolabs, MA, USA), 0.25 μl of Q5 DNA Polymerase (New England Biolabs), primers as specified in Supplementary Table 1, and 5 μl of cDNA template. The following program was used for amplification: 98°C for 30 s, 32 cycles of 98°C for 15 s, 63°C for 5 min, and a final extension at 65°C for 10 min. The four PCR1 products were combined by equal volume and used as input for PCR2 (indexing PCR). The PCR2 reaction was set up using the following reaction mixture: 9.25 μl of nuclease-free water, 5 μl of Q5 reaction buffer, 0.5 μl of 10 mM dNTPs (New England Biolabs), 0.25 μl of Q5 DNA Polymerase (New England Biolabs), 2.5 μl of 1 μM forward primer, 2.5 μl of 1 μM reverse primer, and 5 μl of PCR1 template, as described by de Muinck et al. (2017). The following program was used for PCR2 amplification: 98°C for 30 s, 8 cycles of 98°C for 15 s, 63°C for 5 min, and a final extension at 65°C for 10 min. Individual PCR2-amplified libraries were pooled based on molarity measured using a Tapestation 4200 system (Agilent Technologies, CA, USA). The resulting library pool was cleaned using a double-sided size selection protocol with AMPure XP reagent (cat no A63881, Beckman Coulter INC, CA, USA), aimed at a target size range between 400 bp and 1,000 bp. The final library pool DNA concentration was quantified using the Tapestation 4200 system (Agilent Technologies, CA, USA) and with Qubit Fluorometric Quantification (Fisher Scientific, CA, USA).

Libraries were sequenced on a MiSeq (Illumina) using a Reagent Kit v3 (600 cycles) with 301 cycles paired-end reads and 7 cycles for i7/i5 index reads according to the manufacturer's recommendations. Illumina PhiX Control v3 (Illumina) was spiked in at a concentration of 5%.

Sequence quality and read counts were analyzed with FastQC version 0.11.9 [Andrews, 2010] and MultiQC version 1.14 (Ewels et al., 2016). Illumina adapter sequences were removed from the 3′ end of each read (if present) and quality trimmed with Trimmomatic version 0.39 (Bolger et al., 2014). Trimmomatic was run in paired-end mode with the following parameters: first, the last base was removed, followed by the removal of adapters, before low-quality bases were removed with a sliding window algorithm (4 bases and Q > 20). Unpaired reads were discarded. ISAV-specific primer sequences were removed from the 5′ end of each read with the BBDuk program from the BBMap package version 38.93. BBDuk was run with the following parameters: Ambiguous bases were expanded, the search was restricted to the first 30 bases, a kmer size of k=17 was used, and a single mismatch was allowed. Finally, the forward and reverse reads in each pair were merged into a “single amplicon read” with NGmerge version 1.4.2 (Gaspar, 2018). Dovetail alignments were allowed during the merging process. More than 90,000 reads per sample survived QC trimming and merging.

A consensus sequence for each segment was reconstructed by mapping reads from each sample to the Glesvær/2/90 genome (Merour et al., 2011) with Bowtie2 version 2.4.5 (Langmead and Salzberg, 2012) using the very sensitive option. The Glesvær/2/90 genome consisted of the concatenated sequences HQ259671.1, HQ259672.1, HQ259673.1, HQ259674.1, HQ259675.1, HQ259676.1, HQ259677.1, and HQ259678.1. The resulting mapping files (.sam) were sorted and compressed with Samtools version 1.10 with htslib version 1.10.2-3 (Li et al., 2009) and were used to generate a consensus sequence for each sample with iVar version 1.2.1 (Grubaugh et al., 2019a; Castellano et al., 2021), using all positions flag and requiring a depth of 10 at each position to call consensus. To be able to compare whole genomes efficiently were the segments concatenated with 50 Ns between each segment before multiple sequence alignments (MSA) were generated with MAFFT version 7.520 (Katoh and Standley, 2013) using the local alignment with an iterative refinement algorithm (L-INS-i). The HPR region in segment 6 was manually curated to represent the deletion as a single gap compared to the HP0 sequence. To the best of our knowledge, there is no alignment program that can correctly align the deletion automatically.

IQ-TREE 2 version 2.2.7 was used to generate maximum likelihood (ML) phylogenetic trees (Minh et al., 2020). First, the ModelFinder application in IQ-TREE was used to identify the model that best explained the data (Kalyaanamoorthy et al., 2017). A transversion model with equal transition and transversion rates, unequal base frequencies, and allowing for invariant sites (TVM+F+I) was selected based on the Bayesian information criterion (BIC) score. Ultrafast bootstrapping with 10,000 replicates was performed (Hoang et al., 2018). Patristic distances were extracted from the phylogenetic trees with the cophenetic function from the stats package in R version 4.1.2 and converted to percentage. The numbers of SNPs were counted with snp-dists version 0.8.2 (indels excluded, https://github.com/tseemann/snp-dists). Classification of HPR deletions was made by visual comparison (Devold et al., 2001; Nylund et al., 2003). A multiple sequence alignment was visualized with GGMSA version 1.0.0 under R version 4.1.2 (Zhou et al., 2022).

To sequence the ISAV genome, 40 primer pairs were designed by identifying conserved regions in an MSA created from European isolates of ISAV, with tiling amplicon lengths of 375–438 bp (Supplementary Table 1). The method can be used with both 300 bp PE (V3 chemistry) and 250 bp PE (V2 chemistry) read length. Pilot experiments were performed to assess sequencing depth across the set of amplicons. For amplicons with low sequencing depth, the primers were redesigned, and a new pilot experiment was performed (results not shown). To decrease the number of PCR reactions, PCR primer pairs were pooled in four multiplexes, each containing 10 PCR primer pairs.

The genome of the Glesvær strain (Glesvær/2/90) has previously been sequenced twice with Sanger sequencing (Markussen et al., 2008; Merour et al., 2011). The Markussen genome was generated by the Norwegian Veterinary Institute, and the Merour genome was generated by the National Research Institute for Agriculture, Food and the Environment, France (INRA) (2011). Both genomes were based on cultivated viruses used as templates. The two genomes were not identical and the Markussen genome was 12,141 bp long, whereas the Merour genome was 13,227 bp long. In addition, the two genomes differ by 9 SNPs in the pairwise alignment. These two genomes were used as reference sequences to assess sequencing accuracy for the amplicon method (Table 1). RNA was extracted from the same sample of the Glesvær strain as was used for both the Markussen and Merour genomes. RNA isolated from cell culture supernatant (CCS) was sequenced three times (Table 1, Glesvær-CCS-1,−2 and−3), and RNA extracted from tissue directly was sequenced twice with the amplicon method (Table 1, Glesvær-Tissue-1 and−2). The five resulting consensus sequences were aligned and compared to the Markussen and Merour genomes. The resulting MSA had an average length of 12,317 base pairs. The three CCS samples were identical to each other and had two SNPs compared to the Merour genome and 11 SNPs when compared to the Markussen genomes. The two tissue samples differed by one SNP to each other. They showed four and five SNPs compared to the Merour genome and 13 and 14 SNPs compared to the Markussen genome, respectively (Table 1). For the validation of sequencing accuracy, we defined the distance between the genomes generated by amplicon sequencing, and the comparators should not be >0.1%. The average distance was 0.062%, and the sequence accuracy of the amplicon sequencing method was concluded to be acceptable.

A local ISA epidemic including outbreaks at 12 farms occurred in the northern part of Norway (production area 9, Vestfjorden and Vesterålen) (Figure 1), starting in the summer of 2013 and continuing for 2 years. In total, 12 farms with ISA diagnosis were sampled and sequenced with the amplicon sequencing protocol (Table 2), and approximately 100,000 to 200,000 PE reads were generated per sample as specified in Table 2. Near complete ISAV genomes were generated for all 12 samples (Table 2) with accession numbers as listed in Table 3. The concatenated genomes were aligned with eight concatenated references (Supplementary Table 2), resulting in an MSA with an average length of 12,298 and a range of 11,973 to 12,341 bp. Based on the MSA, an ML tree was inferred, as described in the Materials and Methods section (Figure 2).

The genomes from the six farms in the Lofoten area clustered together in a single monophyletic group (Figure 2, cluster 1). These six farms, represented by samples 1, 2, 7, 8, 9, and 10, are located in the southern part of the production area 9 (Figure 1). The average distance between the six genomes was 0.065% SNPs/site (7.9 SNPs, range: 0 to 17 SNPs). In addition to cluster 1, samples 11 and 12 clustered together (Figure 2, cluster 2) with a distance between the two genomes of 0.017% SNPs/site (2 SNPs), while samples 3 and 6 clustered together (Figure 2, cluster 3) and were separated by 0.025% SNPs/site (3 SNPs). All three clusters had bootstrap support of ≥ 95%. Samples 4 and 5 have bootstrap support of <95% and were not considered a cluster. The distance between these two samples was 0.48% SNPs/site (54 SNPs).

At the time of the outbreaks, segments 5 (fusion protein) and segment 6 (hemagglutinin esterase) were sequenced as a part of the outbreak characterization. One phylogenetic tree was constructed based on each segment. This is still the most common strategy for the characterization of ISAV outbreaks. When segments 5 and segment 6 were analyzed separately with the amplicon data (Figures 3A, B), it became apparent that the two segments were not reflecting the same evolutionary history. In the tree for segment 6 (Figure 3A), sample 4 was closer to cluster 3, with an average genetic distance of 0.16% SNPs/site (2 SNPs) between sample 4 and cluster 3. In the tree for segment 5 (Figure 3B), sample 4 did not belong to cluster 3. The whole genome tree, based on the concatenated segments, supported that sample 4 did not belong in cluster 3 (Figure 2). Moreover, the deletion in the HPR region of segment 6 of sample 4 was classified according to the HPR classification system (Devold et al., 2001; Nylund et al., 2003) as HPR2, while the two samples (sample 6 and 3) of cluster 3 were classified as HPR3, which demonstrate that two independent deletion events had occurred. Clusters 1 and 2 merge in the segment 5 phylogenetic tree as one monophyletic group with bootstrap support of ≥ 95% (Figure 3B). The average distance between the isolates in the two clusters was 0.18% SNPs/site (average 2 SNPs, range 1–4). In the tree representing segment 6 (Figure 3A), clusters 1 and 2 are clearly distinct, as they are in the whole genome tree (Figure 2). Moreover, all isolates of cluster 1 were classified as HPR3, while the isolates in cluster 2 were classified as HPR5, which also indicates that clusters 1 and 2 are two distinct clusters, as indicated by the WGS tree.

Virulent ISAV strains have, in addition to a deletion in segment 6, a point mutation or insertion in a putative protease cleavage site on segment 5 (Markussen et al., 2008). Most of the samples had the expected glutamine-to-leucine (Q266L) mutation. The samples in cluster 1 had a 21 bp in-frame insertion (Figure 4), but not the Q266L mutation, supporting that the isolates in cluster 1 have the same origin and that cluster 1 is distinct from the two other clusters.