Genomic and immunogenic changes of Piscine novirhabdovirus (Viral hemorrhagic septicemia virus) over its evolutionary history in the Laurentian Great Lakes

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 >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.Review GenBank files (https://www.ncbi.nlm.nih.gov/GenBank) MK782981–MK783014 Funding provided by: National Science FoundationCrossref Funder Registry ID: http://dx.doi.org/10.13039/100000001Award Number: DBI-1354806Funding provided by: Agricultural Research ServiceCrossref Funder Registry ID: http://dx.doi.org/10.13039/100007917Award Number: 3655-31320-0000DSampling and nomenclature 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 Novirhabdovirus 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 (Esox masquinongy) 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. Virus isolation 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. Viral genome sequencing 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. Genetic analyses 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. Evaluating evolution and selection 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 (k = 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 >0.95 for FUBAR and p<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 >300 NTs in length. Any sequences that returned as coding in reverse were discounted because VHSV is single-stranded. Virulence and immunogenicity assessments Epithelioma papulosum cyprinid (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 >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-N and EPC IFN were normalized to EPC β-actin detection, and to the gene expression of uninfected cell samples. Relative gene expression used the 2-ΔΔCT method