DOI: 10.5061/dryad.8cz8w9gmh
title: Genomic and immunogenic changes of Piscine novirhabdovirus (Viral hemorrhagic septicemia virus) over its evolutionary history in the Laurentian Great Lakes

description: A unique and highly virulent subgenogroup (VHSV-IVb) of Piscine
 novirhabdovirus, also known as Viral Hemorrhagic Septicemia Virus (VHSV),
 suddenly appeared in the Laurentian Great Lakes, causing large mortality
 outbreaks in 2005 and 2006, and affecting &gt;32 freshwater fish
 species to date. Periods of apparent dormancy have punctuated smaller and
 more geographically-restricted outbreaks in 2007, 2008, and 2017. In this
 study, we conduct the largest whole genome sequencing analysis of VHSV-IVb
 to date, evaluating its evolutionary changes from 48 isolates in relation
 to immunogenicity in cell culture. Our investigation compares genomic and
 genetic variation, selection, and rates of sequence changes in VHSV-IVb,
 in relation to other VHSV genogroups (VHSV-I, VHSV-II, VHSV-III, and
 VHSV-IVa) and with other Novirhabdoviruses. Results show that the VHSV-IVb
 isolates we sequenced contain 253 SNPs (2.3% of the total 11,158
 nucleotides) across their entire genomes, with 85 (33.6%) of them being
 non-synonymous. The most substitutions occurred in the non-coding region
 (NCDS; 4.3%), followed by the Nv- (3.8%), and M- (2.8%) genes.
 Proportionally more M-gene substitutions encoded amino acid changes
 (52.9%), followed by the Nv- (50.0%), G- (48.6%), N- (35.7%) and L-
 (23.1%) genes. Among VHSV genogroups and subgenogroups, VHSV-IVa from the
 northeastern Pacific Ocean has shown the fastest substitution rate
 (2.01x10-3), followed by VHSV-IVb (6.64x10-5) and by the VHSV-I, -II and
 -III genogroups from Europe (4.09x10-5). A 2016 gizzard shad (Dorosoma
 cepedianum) from Lake Erie possessed the most divergent VHSV-IVb sequence.
 Yet. the in vitro immunogenicity analysis of that sample displayed reduced
 virulence (as did the other samples from 2016), in comparison to the
 original VHSV-IVb isolate (which had been traced back to 2003, as an
 origin date). The 2016 isolates that we tested induced milder impacts on
 fish host cell innate antiviral responses, suggesting altered phenotypic
 effects. In conclusion, our overall findings indicate that VHSV-IVb has
 undergone continued sequence change and a trend to lower virulence over
 its evolutionary history (2003 through present-day), which may facilitate
 its long-term persistence in fish host populations.
description: <b>Sampling and
 nomenclature</b> All fishes were
 collected by federal, provincial, and state agency collaborators,
 following their respective permits and regulations, and according to the
 University of Toledo Institutional Animal Care and Use Committee (IACUC)
 protocol #106419. VHSV-IVb samples from 2015–2016 were collected by us and
 archived isolates were provided by G. Kurath (USGS, Seattle, WA) as frozen
 infected media from Bluegill Fry (BF-2) cell culture (ATCC: CCL-91).
 Additional complete VHSV (I–IVa; Sup. Table 1) and
 <i>Novirhabdovirus</i> genomes were downloaded from NCBI
 GenBank and aligned using the program MEGA X. VHSV-IVb isolates are named
 here with unique identifiers, using the first letter of their lake name,
 the last two digits of the sampling year, followed by the first two
 letters of the host species’ common name. Example: the original isolate,
 collected from a Lake St. Clair muskellunge (<i>Esox
 masquinongy</i>) in 2003 here is termed C03MU (it also is known as
 M103GL; GenBank: GQ385941). If more than one isolate shared the above
 identifiers, an additional lower-case letter was added. We defined a
 haplotype as “a unique gene sequence differing by one or more nucleotide
 substitutions” from C03MU. Population genetic relationships, based on this
 dataset, are evaluated separately. <b>Virus
 isolation</b> Frozen media was thawed on ice, and
 30 µl from each sample was diluted 1:5 with serum-free Minimal Essential
 Medium (MEM; Gibco), and added to individual wells of a 12 well plate with
 confluent BF-2 monolayers. Cells were incubated 1 h with the infected
 medium at 20°C, which then was replaced with complete MEM 10% (v/v) cosmic
 calf serum (GE Healthcare) and 1% penicillin/ streptomycin antibiotics
 (Invitrogen). Infected cells were incubated at 20°C for ≤one week and
 sampled when ≥80% cytopathic effects (CPE) were achieved. Media was
 collected and added to a 1.5 ml tube with 250 µl of versene (Gibco) for 10
 min, followed by 4 min 4,000 rpm centrifugation at 4°C. The supernatant
 was discarded, and the versene/cell mixture spun again, as before. The
 remaining supernatant again was discarded, and 250 µl Trizol® (Invitrogen)
 was added to the remaining pellet. <b>Viral
 genome sequencing</b> cDNA was synthesized from
 total RNA extracted from the tissue samples using SuperScript IV
 (Invitrogen), following manufacturer’s instructions. Genomic cDNA was
 amplified in four segments using primers from Schönherz, substituting
 VHSV_Frag1I_nt18_+s with a more specific primer (5’GAGAGCTGCAGCACTTCACCG
 C3’), and 1 μl cDNA in 25 µl polymerase chain reactions (PCRs) with One
 Taq DNA polymerase (New England Biolabs). Amplicons were examined under UV
 light on 1% agarose gels stained with ethidium bromide. Target PCR
 products were gel-excised and purified using QIAquick Gel Extraction kits
 (Qiagen). Genomic
 sequencing was outsourced to Ohio State University’s Molecular and
 Cellular Imaging Center (Wooster, OH). Sequences then were uploaded by us
 to the Galaxy web platform and analyzed with usegalaxy.org programs [61].
 Segments were aligned to the reference VHSV-IVb genome (C03MU, GenBank:
 GQ385941) using MAP WITH BWA-MEM [62]. For each of the 48 isolates we
 sequenced here, consensus sequences were generated followed by manual
 checking of each single nucleotide polymorphism (SNP) and read coverage
 was determined using the Integrative Genomics Viewer (IGV) [63,
 64]. Additional
 PCRs were conducted to amplify the front 700 nucleotides (NTs) and end 400
 NTs of the genome, with 45 s extension time, in order to complete the
 whole genomes. The front segment utilized VHSV_Frag1I_nt18_+s
 (5’GAGTTATGTTACARGGGACAGG3’) and anti-sense 5’TGACCGAGATGGCAGATC3’, and
 end primers were designed by us based on the VHSV-IVb original reference
 genome (GenBank: GQ385941) (End sense: 5’CCCAGATGCTATCACCGAGAA3’, End
 anti-sense: 5’ACAAAGAATCCGAGGCAGGAG3’). Those cleaned products were
 Sanger-sequenced at Cornell DNA Services (Ithaca, NY), and checked and
 aligned by us. The consensus sequences, front, and end segments then were
 concatenated, aligned, and trimmed using MEGA X, to complete the whole
 genome sequences. <b>Genetic
 analyses</b> The Basic Local Alignment Search
 Tool (BLAST) was used to align the whole genome sequences. JMODELTEST v3.7
 employed the Akaike Information Criterion (AIC) to determine the best-fit
 evolutionary models [66]. Phylogenetic analyses evaluated the most
 parsimonious evolutionary relationships with Maximum Likelihood (ML)
 (PHYML v3.0) and Bayesian (MRBAYES v3.1) algorithms. For the latter,
 Metropolis-coupled Markov Chain Monte Carlo (MCMCMC) analyses were run for
 five million generations, sampling every 100 to obtain posterior
 probability (pp) values. Burn-in was determined by plotting the log
 likelihood values to identify when stationarity was reached, discarding
 the first 25%. Branch support for ML was calculated with non-parametric
 bootstrapping replications, using 500 for the VHSV analysis and 1450 for
 the VHSV-IV analysis. Haplotype networks were analyzed by us using PopART
 (Population Analysis using Reticulate Trees; http://popart.otago.ac.nz)
 with TCS. <b>Evaluating evolution and
 selection</b> Comparative divergence times among
 VHSV taxa were estimated using BEAST v1.10.4, using JMODELTEST output and
 a relaxed molecular clock with lognormal distribution, sampled every
 50,000 of 500,000,000 generations. Outputs were assessed with TRACER v1.5
 (in BEAST) to ensure stationarity. Collection dates were used as
 calibration points and the tree branches were set, following the PHYML
 output. Numbers of nucleotide substitutions per site per year
 (<i>k</i> = substitutions site–1 yr–1) were determined from
 the pairwise (p) distances. Nucleotide (NT) and amino acid (AA)
 substitutions were evaluated for all VHSV-IVb isolates and compared to the
 reference C03MU sequence (GenBank: GQ385941). Results from the VHSV-IVa
 isolates were compared to KRRV9601 (GenBank: AB179621), and VHSV-I+III
 with the European VHSV-Ia isolate Hededam (GenBank: Z93412).
 Two codon-based methods
 examined the possibility of selection pressures for the entire genome and
 separately for each gene. Fast, unconstrained Bayesian approximations
 (FUBAR) were used to identify positive or purifying selection, which are
 the pressures that select for beneficial traits. However, FUBAR’s
 assumption of constant selection may not accurately represent VHSV-IVb
 since selection pressures may differ across hosts and with environmental
 factors. To remedy this, we additionally used MEME (mixed effects model of
 evolution), which can detect positive selection under strong purifying
 selection or with removal of detrimental variants. FUBAR and MEME were run
 with HyPHY on DataMonkey (www.datamonkey.org), with their significance
 evaluated using the posterior probability of &gt;0.95 for FUBAR and
 <i>p</i>&lt;0.05 for MEME. Amino acid sequences were submitted
 to the Phyre2 web portal (www.sbg.bio.ic.ac.uk/phyre2/) for protein
 modeling, prediction, and analyses. For open reading frame (ORF) analysis,
 full length cDNA nucleotide sequences per isolate were submitted to the
 NCBI ORF Finder (www.ncbi.nlm.nih.gov/orffinder/) to search for alternate
 products &gt;300 NTs in length. Any sequences that returned as coding
 in reverse were discounted because VHSV is single-stranded.
 <b>Virulence and immunogenicity
 assessments</b> <i>Epithelioma papulosum
 cyprinid</i> (EPC) (ATCC: CRL-2872) and BF-2 cells were cultured as
 previously described by Niner (2019). Three viable independent viral
 stocks from 2016 – isolates E16GSa (Cell16a) and E16LB (Cell16b and
 Cell16c) – were derived from pooled organs (kidney, liver, and spleen) of
 sampled fish. Cell culture amplified stocks are labelled “Cell” in our
 results. Due to low concentrations, the original E16LB sample was only
 partially sequenced, which matched E16GSa in those sequence regions [51,
 57]. The cell culture amplified reference control (CellC03) was derived
 from C03MU. VHSV-IVb isolates were amplified for subsequent purification
 by infecting confluent monolayers of BF-2 cells in 15 cm tissue culture
 dishes with a 1:1000 (v/v) dilution of un-purified virus stock in
 serum-free MEM (=EMEM). Viral adsorption occurred for 1 h, then the medium
 was replaced with complete EMEM. The plates were incubated at 20°C, until
 &gt;75% CPE occurred at ~72 hours post infection (hpi). The virus
 containing media and the attached cells were collected and subjected to
 one freeze-thaw cycle, followed by a 30 min 4,000 rpm centrifugal removal
 of the debris at 4°C. The supernatant from each isolate was clarified
 using 0.22 μm syringe-tip filters, and the viral particles were purified
 through a 25% (w/v) sucrose pad upon 3 h of 25,000 rpm
 ultra-centrifugation at 4°C. Virus-containing pellets were re-suspended
 overnight at 4°C in phosphate buffered saline (PBS). Virus stocks were
 titered by 1:10 serial dilution using confluent EPC cells, divided into
 100 μl aliquots, and stored at -80°C. Virus-induced CPE were quantified
 with a sulforhodamine B (SRB) assay [42]. Infected EPC cells were sampled
 respectively at 0, 18, 50, 72, or 96 hpi at a multiplicity of infection
 (MOI)=1.0 for C03MU or Cell16-a, -b, or -c, and stored at -80°C. A viral
 yield assay compared the replication ability of the three test isolates of
 E16GSa to C03MU virus, by titering the media from each time point in 1:10
 serial dilutions on BF-2 cells. Plaques were counted at 96 hpi, and the
 final viral concentration in plaque forming units per ml (pfu/ml) was
 calculated per each time point. Antiviral assays followed Ke et al. [42].
 UV-irradiated media from each harvested time point was added to EPC cells
 in 1:3 serial dilutions for 24 h. Cells then were challenged with
 sucrose-purified VHSV-IVb C03MU for 96 h, fixed, and stained with crystal
 violet (Sigma-Aldrich). CPE plaques were counted and normalized to counts
 obtained from untreated and uninfected wells. One unit of IFN was defined
 as the dilution that conferred 50% protection from viral CPE, whose value
 then was used to measure the IFN activity (ml) in the testing culture
 media per time point. Total RNA was extracted from the
 infected EPC cells, following a previously optimized TRIzol protocol [59].
 One μg of each RNA sample was reverse transcribed into cDNA using a 10 min
 incubation at 70°C with 100 ng of random hexamer primer (Thermo
 Scientific) and water, for 7 μl total volume. Reactions were cooled to 4°C
 before adding 13 µl of Moloney Murine Leukemia Virus Reverse Transcriptase
 (M-MLV-RT) mixture [10X First Stand buffer, 10 mM dNTPs, 0.05 mM random
 hexamers, 25 U/µl RNasin Plus (Promega), and 200 U/µl M-MLV (Promega)],
 which then was incubated at 42°C for 1 h. cDNA was diluted 10-fold in
 water and archived at -80°C. 1 μl of each cDNA was tested using RT-qPCR,
 with 5 μl of the Radiant Green Lo-ROX 2X qPCR kit (Alkali Scientific), 50
 ng of each oligonucleotide, and water to total 10 μl. Primers used were:
 VHSV-Nse/as, EPC IFN se/as, and β-actin se/as. Reactions and data
 collections were performed on a C1000 Real Time Thermocycler (Bio-Rad),
 with initial 3 min denaturation at 95°C, followed by 40 cycles of 15 s at
 95°C, and 30 s elongation at 60°C. Readings were recorded at the end of
 each elongation cycle, and the threshold values were obtained from an
 automated single point threshold within the log-linear range. Detections
 of VHSV-<i>N</i> and EPC IFN were normalized to EPC β-actin
 detection, and to the gene expression of uninfected cell samples. Relative
 gene expression used the 2<sup>-ΔΔCT</sup> method.
description: Review GenBank files (https://www.ncbi.nlm.nih.gov/GenBank) MK782981–MK783014