Increasing virulence, but not infectivity, associated with serially emergent virus strains of a fish rhabdovirus.

Surveillance and genetic typing of field isolates of a fish rhabdovirus, infectious hematopoietic necrosis virus (IHNV), has identified four dominant viral genotypes that were involved in serial viral emergence and displacement events in steelhead trout (Oncorhynchus mykiss) in western North America. To investigate drivers of these landscape-scale events, IHNV isolates designated 007, 111, 110, and 139, representing the four relevant genotypes, were compared for virulence and infectivity in controlled laboratory challenge studies in five relevant steelhead trout populations. Viral virulence was assessed as mortality using lethal dose estimates (LD50), survival kinetics, and proportional hazards analysis. A pattern of increasing virulence for isolates 007, 111, and 110 was consistent in all five host populations tested, and correlated with serial emergence and displacements in the virus-endemic lower Columbia River source region during 1980-2013. The fourth isolate, 139, did not have higher virulence than the previous isolate 110. However, the mG139M genotype displayed a conditional displacement phenotype in that it displaced type mG110M in coastal Washington, but not in the lower Columbia River region, indicating that factors other than evolution of higher viral virulence were involved in some displacement events. Viral infectivity, measured as infectious dose (ID50), did not correlate consistently with virulence or with viral emergence, and showed a narrow range of variation relative to the variation observed in virulence. Comparison among the five steelhead trout populations confirmed variation in resistance to IHNV, but correlations with previous history of virus exposure or with sites of viral emergence varied between IHNV source and sink regions. Overall, this study indicated increasing viral virulence over time as a potential driver for emergence and displacement events in the endemic Lower Columbia River source region where these IHNV genotypes originated, but not in adjacent sink regions.

1 Introduction

Infectivity and virulence are important viral traits that impact how a virus interacts with its host, and may ultimately influence how it succeeds in spreading and persisting in host populations at the landscape level. Variations in virulence and infectivity are thus hypothesized to be potential drivers of epidemiological success among different genotypes of a virus that co-occur in the same geographic region. The virus studied here, infectious hematopoietic necrosis virus (IHNV), is an important rhabdoviral pathogen of Pacific Salmon and trout and is endemic throughout the Pacific Northwest of North America, where it has caused disease epidemics in cultured fish in hatcheries and fish farms since the 1950s. Here, we have investigated infectivity and virulence as potential drivers of three natural viral emergence and displacement events observed in steelhead trout, Oncorhynchus mykiss. These studies were designed after previous work demonstrated that in-host replicative fitness, competitive co-infection fitness, superinfection fitness, induction of innate immunity, and virion stability in the external environment did not correlate with emergence and displacement events, but that virulence might (Breyta, Jones, and Kurath 2014; Kell, Wargo, and Kurath 2014).

Genetic typing and phylogenetic analysis of IHNV field isolates have revealed three major IHNV genogroups in North America, designated U, M, and L (Kurath et al. 2003). Both U and M genogroup viruses occur in the landscape studied here, and relevant to this work, the M genogroup has host-specificity demonstrated as high incidence and high virulence in steelhead and rainbow trout (anadromous and freshwater forms, respectively, of the fish species O. mykiss) (Garver, Conway, and Kurath 2006a; Peñaranda, Purcell, and Kurath 2009).

Long established fish health surveillance programs throughout the Pacific Northwest continually generate IHNV field isolates that are analyzed by genetic typing based on the sequence of the 303 nucleotide (nt) ‘midG’ region in the viral glycoprotein gene (Emmenegger and Kurath 2002; Kurath et al. 2003; Breyta et al. 2013). Full glycoprotein gene and genome sequences are available for selected IHNV isolates, providing analyses that are consistent with midG phylogenetic analysis (Breyta and Kurath, unpublished). However, midG typing has been done on thousands of isolates, allowing far more thorough near real-time monitoring of epidemiological patterns. In the course of ongoing genetic typing, four M genogroup IHNV genotypes were observed as dominant genotypes in steelhead trout, based on detection over multiple years and spread from the initial site of detection to multiple sites. These genotypes were each detected in both epidemic juvenile and asymptomatic adult fish, and they were involved in major virus emergence and displacement events in different regions of the Pacific Northwest (Fig. 1, Table 1) (Kell, Wargo, and Kurath 2014).

Figure 1.

Geographic regions of emergence and displacement events for M genogroup IHNV in steelhead trout, and origin of steelhead trout populations used in experimental studies. Three regions of the North American Pacific Northwest have different historical levels of M genogroup IHNV. The Lower Columbia River region (orange) is endemic for IHNV in steelhead trout, as is the Lower Snake River (blue), but the Washington Coast (green) is not endemic but has had periodic emergences if IHNV in steelhead trout (A). Four genetic types of virus, mG007M (007), mG111M (111), mG110M (110), and mG139M (139), have been observed in a serial pattern of emergence and displacement. The temporal incidence of these dominant genetic types is shown for each region (B). Five populations of steelhead trout were tested with representative isolates of the four virus genotypes, and their site of origin are indicated by two letter fish population abbreviation: Lake Quinault (LQ), Quinault National (QN), Skamania (SK), wild Elwha (EW), and Dworshak (DW).

Table 1.

Virus isolates used in controlled animal studies

Virus strain (genotype, subgroup)a	Geographic originb	Origin and emergence/displacement phenotype by region
007 (mG007M, MN)	Lower Columbia, WA	Emerged in lower Columbia steelhead during 1980–1994; detected in Lower Snake region in 1983; never detected in coastal WA
111 (mG111M, MD)	Lower Columbia, WA	Emerged in lower Columbia steelhead trout during 1994–1999 where it displaced type mG007M; never detected in Lower Snake; emerged in WA Coast 1997c
110 (mG110M, MD)	Coastal WA	Emerged in lower Columba steelhead trout in 2002–2013 where it displaced type mG111M; never detected in Lower Snake region; emerged in WA Coast in 2007–2009c
139 (mG139M, MD)	Lower Snake, ID	Emerged in lower Columbia steelhead during 2003–2008, where it did not displace type mG110M; emerged in Lower Snake region in 2008–2013; emerged in WA Coast in 2009–2011, where it displaced type mG110M

aIHNV midG genotype using the USD nomenclature (Emmenegger and Kurath 2002), which specifies unique sequences within the 303 nt midG region of the viral glycoprotein gene. MN indicates a sequence type in the M genogroup that does not fall into any defined subgroup; MD is a phylogenetically defined subgroup within the M genogroup.

bWA, Washington; Columbia, Columbia River Basin; ID, Idaho.

cIn the coastal WA region mG111M occurred as a single epidemic in juvenile steelhead trout that were culled and the genotype was not detected again. Therefore, the coastal emergence of mG110M in 2007 is not considered to have displaced mG111M due to the 10-year gap in M virus detection in the region).

In the virus-endemic lower Columbia River region, IHNV genotype mG007M emerged in 1980 and persisted through 1994, type mG111M emerged in 1994 and persisted through 1999, and type mG110M emerged in 2002 and persisted through at least 2013 (Garver, Troyer, and Kurath 2003; Breyta et al. 2013). This serial emergence of three genotypes comprises two displacement events within the lower Columbia River region in which genotype mG007M was displaced by mG111M, and mG111M was then displaced by mG110M. We use the term ‘displacement’ here to indicate an event in which a previously dominant viral genotype is no longer detected and a new dominant genotype arises within the same general geographic region. We do not assume that direct competition between displaced and displacing genotypes occurs in all displacements, although that is one possible mechanism for such events. Subsequent to emergence in the lower Columbia River region, genotype mG110M spread to infect steelhead trout in the coastal Washington region (Fig. 1) and was responsible for a sustained emergence involving numerous epidemics from 2007 to 2009 (Breyta et al. 2013). A fourth genotype, mG139M, emerged on a smaller scale in the lower Columbia River region during 2003–2008, but it did not displace mG110M there. Instead it emerged on a larger scale in the upstream Snake River region in 2008–2010 (Breyta et al. 2016), and then in the coastal Washington region during 2009–2011, where it displaced mG110M (Fig. 1, Table 1) (Breyta et al. 2013). Therefore, types mG110M and mG139M comprised a third ‘conditional’ displacement event, because displacement of mG110M occurred only in coastal Washington, not in the lower Columbia River region.

Four IHNV isolates were selected to represent the four dominant genotypes, mG007M, mG111M, mG110M, and mG139M, in controlled laboratory infection studies. These four virus isolates are referred to hereafter as 007, 111, 110, and 139, respectively. In Kell, Wargo, and Kurath (2014), these viruses were tested for variation in fitness. In this study, they were tested for virulence and infectivity, measured as mortality or infection after immersion exposure to four different doses of each virus. Differences in virulence were assessed for viruses within each displacement pair by estimating the dose of virus required to cause mortality of 50 percent of exposed fish (LD50), comparing survivorship kinetics, and conducting proportional hazard analysis. To test for possible contributions of host genetic factors, these studies were all conducted in five different steelhead trout populations selected to represent diverse host genotypes from the three geographic regions where the emergence and displacement events occurred, and one population naive to M genogroup IHNV (Fig. 1). The host populations tested included the Lake Quinault (LQ) population used in Kell, Wargo, and Kurath (2014), and the Quinault National (QN) population shown previously to be significantly more resistant to IHNV than the LQ population (Breyta, Jones, and Kurath 2014). Challenge studies in the LQ population were repeated in two different years as a test of year-to-year variation in host resistance levels. Infectivity of the viruses was compared in the LQ and QN populations by estimating the dose required to infect 50 percent of exposed fish (ID50).

The results demonstrated a consistent pattern, evident in all five host populations tested, of increasing virulence for the three virus strains that serially emerged and displaced one another in the endemic Lower Columbia region. This indicated evolution of increased virulence as a potential driver in the first two displacement events. This is a valuable observation because studies documenting evolution of viral virulence in the field are relatively rare (Kurath and Wargo 2015). In contrast, infectivity did not vary significantly among viruses. Fish host populations varied significantly in resistance, suggesting that host genetics play some role in emergence of IHNV. Observations presented here suggest that drivers for landscape virus emergence and displacement events may differ according to regional and ecological factors.

2 Methods

2.1 Virus isolates

The strains used to represent the emergence and displacement events in Fig. 1 are as in Kell, Wargo, and Kurath (2014) and have different temporal and geographic origins (Table 1). The 007 strain (type mG007M, previously referred to as LR80) was provided by J. Winton and has been described (Nichol, Rowe, and Winton 1995). The 111 strain (type mG111M, previously referred to as Mer95) was isolated from steelhead trout in the Columbia River basin in 1995 and was provided by J. Thomas (Washington Department of Fish and Wildlife). The 110 strain (type mG110M, previously referred to as Qts07) was isolated from epidemic juvenile steelhead trout at a hatchery on the Queets River in 2007, and was provided by J. Bertolini (Northwest Indian Fisheries Commission). These first three strains were used in studies described previously (Breyta, Jones, and Kurath 2014). The 139 strain (type mG139M, previously referred to as Dw09) was isolated at the Dworshak National Fish Hatchery on the Clearwater River in Idaho in 2007 from epidemic juvenile steelhead trout, and was provided by C. Samson (US Fish and Wildlife Service). As in Kell, Wargo, and Kurath (2014), these four viruses will be referred to hereafter by genotype as 007, 111, 110, and 139, respectively. Viruses were propagated in the Epithelioma papulosum cyprini fish cell line (Fijan et al. 1983) using a multiplicity of infection of 0.01–0.001 at 15 °C, and titered by high-accuracy plaque assay as previously described (Wargo, Garver, and Kurath 2010).

2.2 Steelhead trout (O. mykiss) populations

The five populations used in challenges were selected from geographic regions that differ in historical IHNV burden (Table 2). Steelhead trout populations were obtained as eggs as follows: Lake Quinault (LQ) from T. Jurasin (Quinault Nation, Lake Quinault Tribal Fish Hatchery); Quinault National (QN) from B. Gilles (US Fish and Wildlife Service, Quinault National Fish Hatchery); Elwha-wild (EW) from L. Ward (Lower Elwha Klallam Tribe, Elwha Tribal Fish Hatchery); Skamania (SK) eggs from J. Allen (Washington Department of Fish and Wildlife, Skamania State Fish Hatchery); and Dworshak (DW) from C. Samson (US Fish and Wildlife Service, Dworshak National Fish Hatchery). Eggs from each population were pooled from no fewer than 12 pairs of male/female parents to ensure a broad genetic representation of the population and to minimize potential for a founder resistance effect. Eyed eggs were transported, disinfected, and reared in the wet laboratory facilities of the Western Fisheries Research Center (WFRC). Eggs and fish were reared in single-pass flow through 10 °C filtered and UV-treated lake water and fed a moist pellet diet from Skretting at a feed rate of 1 to 1.5 percent average body mass. All fish populations were constantly monitored for consistent health. All challenge studies were performed with a constant water temperature of 10 °C.

Table 2.

Steelhead trout populations used in controlled animal studies

Fish population, abbreviation	Region, river of origina	Rationale for including in study; historic incidence with each virus type
Lake Quinault, LQ	Coastal WA, Quinault R	Region with newly emerged M group IHNV during 2007–2011, not endemic; recent history of mG110M and mG139M
Quinault National, QN	Coastal WA, Quinault R	Region with newly emerged M group IHNV during 2007–2011, not endemic; recent history of mG139M only
Elwha-Wild, EW	Salish Sea WA, Elwha R	Region with no history of M genogroup IHNV
Dworshak National, DW	Lower Snake ID, Clearwater R	Region with endemic history of M group IHNV since 1983, history of mG007M and mG139M
Skamania, SK	Lower Columbia WA, Washougal R	Region with endemic history of M group IHNV since 1981, history of mG007M, mG111M, mG110M, and mG139M

aCoastal indicates Washington State Pacific coastal watersheds

R, river; Columbia, Columbia River Basin; Snake, Snake River

2.3 Exposure experiments in steelhead trout

A total of nine different controlled laboratory challenges were performed in 2012 and 2014. In 2012, two experiments were conducted to compare the lethal dose (LD50) and infectious dose (ID50) estimates in the LQ steelhead trout population only. Subsequently in 2014, four experiments were conducted to determine the LQ and QN steelhead trout LD50 and ID50 estimates, and an additional three experiments determined LD50 estimates in other steelhead trout populations (SK, EW, DW) from a wider geographic distribution.

2.4 Immersion challenges for LD50 in juvenile fish

In 2012 and 2014, fish from each population were reared for 3 months after yolk sac absorption, to an average size of ∼1.5 g. Each fish population was tested independently, as the fish reached this target size. Batch challenge was conducted by the immersion method as previously reported (Garver, Batts, and Kurath 2006b) using triplicate tanks of randomly selected twenty-five fish for each virus dose, or mock treatment. In brief, fish were counted into 1 l of static, temperature controlled water with aeration before virus or mock solutions were added. Exposures consisted of four virus doses (2e+02, 2e+03, 2e+04, and 2e+05 PFU/ml) of each virus stock, or cell culture media alone, for 1 h before water flow was restored to the full volume of 4 l. All tanks were monitored daily for mortality for 30 days. From the fish that died during the challenge, 20 percent were randomly selected to be analyzed by plaque assay to confirm IHNV as cause of death; plaque assay lower limit of detection is 100 plaque-forming units per gram of whole fish (PFU/g).

2.5 Immersion challenges for ID50 in juvenile fish

In 2012 and 2014, fish from the LQ and QN populations were reared for 3 months after yolk sac absorption, to an average size of ∼1 g. Each fish population was tested independently, as the fish reached this target size. Batch challenge was conducted by the immersion method as above, using single tanks of randomly selected sixteen to twenty fish for each virus dose, or mock treatment. Exposures consisted of four or five virus doses, or cell culture media alone, for 1 h before water flow was restored to the full volume of 4 l for a 1-h rinse. Individual fish were then transferred to isolation beakers of 400 ml of static, temperature controlled water. Beakers were covered with a mesh barrier that prevented fish escape but allowed observation. All isolation beakers were monitored daily for 3–7 days. After 3 or 7 days, fish were euthanized within each beaker by the addition of 1 ml of 96 mg/ml MS-222 and 1 ml of 36 mg/ml sodium bicarbonate to reach a final concentration of 240 mg/l of buffered MS-222. All fish were then aseptically removed, weighed, and stored at  −80 °C until processing for viral RNA quantification.

2.6 RNA extraction, cDNA synthesis, and quantification via qPCR

Total RNA was extracted from whole fish as previously described (Wargo et al. 2010). In brief, 4 ml/g fish of guanidinium thiocyanate-based denaturing solution was added to each fish, and the fish was homogenized using a Seward Stomacher® 80 (Biomaster, UK). RNA was extracted from 1 ml of the homogenate and cDNA was prepared using M-MLV reverse transcriptase with random heximer primers (Wargo et al. 2010). Viral load in each fish was quantified using a general IHNV G-gene qPCR assay (Purcell et al. 2006) that detects all four IHNV genotypes used here. In brief, each cDNA sample was combined with forward and reverse primer and Taqman probe and then amplified on a 7900HT ABI Prism machine. This assay detects both genomic and messenger RNA (Purcell et al. 2006), and this combined quantity will be referred to as viral load per gram of host tissue. Lower limit of detection for this assay is 100 genome copies per gram of fish.

2.7 Statistical analysis

Final mortality in replicate exposure groups was tested for non-independence using the Chi-squares test in the Mass package of the R statistical program. LD50 and ID50 estimates were calculated by Generalized Linear Model as implemented in R using the dose.p function of the Mass package, with a quasi-binomial distribution prior (Venables and Ripley 2002). Significant differences between specific LD50 or ID50 estimates were calculated by the Welch-Satterthwaite (W-S) two-tailed t-test correction in R, as previously described (Breyta, Jones, and Kurath 2014). Survival analysis was conducted using the Survival package of R (Therneau and Grambsch 2000), with a type 1 right-censor prior. Survival analysis used the Kaplan–Meier (K–M) survivorship model for individual doses of virus in each fish population and used both Wilcoxon and log-rank tests, both of which produced similar results. The Cox proportional hazard model was used to test the significance of the covariates of fish population and virus strain on risk of death after exposure to virus for all doses and fish populations together and was assessed using the likelihood ratio test for significance. Testing the proportionality assumption of the Cox model was performed using the Cox–Aalen semi-parametric multiplicative model using the timereg package of R (Martinussen and Scheike 2006).

3 Results

3.1 Virulence measured as LD50 among four strains of IHNV

The four virus strains chosen to represent the serial emergence events observed during 1980–2013 (Fig. 1) were tested for ability to cause mortality in steelhead trout. In all experiments, there was a stable dose response between virus dose and cumulative mortality, and tight clustering of triplicate groups as measured by Chi-squared analysis. Calculating an LD50 value to explore relative virulence means that the lower the LD50 estimate, the higher the virulence of the virus. Overall the observed LD50 values varied more than seven orders of magnitude, from 3.51e+02 for the 110 virus in LQ 2014 fish to 7.31e+09 for the 007 virus in EW fish (Fig. 2). Within fish populations (i.e. within experiments) the LD50s for the four virus strains varied by as little as one (DW fish) or as much as four (EW fish) orders of magnitude. In all five fish populations, the first three virus strains showed a consistent pattern of decreasing LD50 value, indicating an increase in virulence with emergence time. The fourth emergent virus, 139, did not follow this pattern and was not consistent in all fish populations. In four of the five fish populations, the 139 virus had a higher LD50 estimate than the previous virus 110, but in EW steelhead trout, the 139 virus had a lower LD50 than 110 (Fig. 2E).

Figure 2.

Estimates of how much virus is needed to kill 50 percent (LD50) of immersion-exposed fish for four virus isolates in each of five fish populations. LD50 estimates (PFU/ml, 95% CI error bars) for each of four virus strains (bottom, arranged in order of emergence time from left to right) in Lake Quinault steelhead trout in two different years (A, B, 2012 and 2014, respectively), and Quinault NFH steelhead trout (C), Skamania SFH steelhead trout (D), wild Elwha River steelhead trout (E), and Dworshak NFH steelhead trout (F) in 2014 only. Significant differences between virus pairs indicated as follows: ε indicate that the LD50 estimates of the two viruses is significantly different (P < 0.05), and α# indicates that the K–M survivorship is significantly different (P < 0.05) at the exposure dose(s) indicated by #, where 2 = 2e + 02, 3 = 2e + 03, 4 = 2e + 04, and 5 = 2e + 05 PFU/ml.

Since each serially emergent virus type is thought to have displaced the previously emerged virus type, tests for significant differences between LD50 estimates of virulence focused on the comparison of results for viruses in three displacement pairs where the first virus strain is the previously established type and the second is the displacing type. Thus, the first displacement pair is 007:111, the second is 111:110, and the third is 110:139. While the overall trend of LD50 estimates for the first three viruses decreasing with viral emergence time was consistent in all experiments, the LD50 estimates were not significantly different in every displacement pair, or in every fish population (Table 3). In the 2012 LQ experiment, the LD50 estimates were only significantly different for the first pair, but in the 2014 LQ experiment both the first and second pairs were significantly different (Fig. 2A and B, ε, P = 0.012, P = 3.4e−4, P = 8.9e−4, respectively). Among experiments with the other four fish populations, only the first displacement pair in SK fish and the third pair in DW fish differed significantly by LD50 values (Fig. 2D and F, ε, P = 0.024, P = 0.049, respectively). In total, the LD50 analysis found significant differences for the first displacement pair in three of six experiments, the second pair in one of six experiments, and the third pair in one of six experiments.

Table 3.

Significance of virulence differences for each of three displacement pairs of IHNV, from LD50 and K–M survivorship analyses

Experiment	Analysis method	007:111a	111:110a	110:139a
2012 LQ	LD50	< , P = 0.012	<, ns	>, ns
	K–M	<, 2e3	<, 2e4	>, 2e5
2014 LQ	LD50	<, P = 0.000	<, P = 0.000	>, ns
	K–M	<, 2e3, 2e4, 2e5	<, 2e2, 2e4	>, ns
2014 QN	LD50	<, ns	<, ns	>, ns
	K–M	<, 2e3	<, ns	>, ns
2014 SK	LD50	<, P = 0.024	<, ns	>, ns
	K–M	<, 2e2, 2e4, 2e5	<, ns	>, 2e4
2014 EW	LD50	<, ns	<, ns	<, ns
	K–M	<, 2e2	<, 2e2	<, 2e2
2014 DW	LD50	<, ns	<, ns	>, P = 0.049
	K–M	<, ns	<, 2e5	>, ns

aQualitative trend of virulence difference is shown first in each entry indicating if the earlier virus had a lesser < or greater > virulence than the later emerging virus. For LD50 estimates significant WS t-test P-values < 0.05 are shown. For survivorship analyses the specific virus dose(s) at which K–M survival had P < 0.05 are shown; ns, not significant. For LD50 estimates, the lower the estimate, the higher the virulence; for K–M, the lower the survivorship, the higher the virulence.

Although not appropriate for statistical analysis, it was interesting to note that LD50 estimates varied among fish populations as well. Overall, the LQ steelhead trout had a broader range of LD50 values in 2014 than in 2012, with a higher LD50 for the 007 virus strain, but lower LD50s for the other three strains. All viruses had higher LD50s in QN fish than in LQ fish, with the exception of 139, which had a similar LD50 estimate in QN2014 and LQ2012 fish. Comparing LD50 values for all four viruses within each host population, the overall IHNV resistance ranking of the five fish populations would be EW > QN > SK > DW > LQ (in order of highest to lowest average LD50).

3.2 Virulence measured through estimates of dose-specific survivorship after infection.

K–M survival analysis was performed with the mortality data to determine if the instantaneous probability of survival (survivorship) differed significantly for fish exposed to displacement pair viruses at each specific virus dose. Using the K–M analysis to explore relative virulence means that the lower the survivorship, the higher the virulence of the virus. An example of this analysis is shown for the 2012 LQ experiment in Fig. 3, where survival curves for each displacement pair at three exposure doses of virus are shown. In this experiment, the survivorship of fish exposed to the displacing virus was significantly lower than fish exposed to the displaced virus for three of the twelve comparisons: the first displacement pair at the 2e+03 PFU/ml dose, the second displacement pair at the 2e+04 PFU/ml dose, and the third displacement pair at the 2e+05 PFU/ml dose (Fig. 2A–F, α, all P values presented in Supplementary Table S1). None of the virus pairs in 2012 LQ fish were significantly different at the lowest dose of 2e+02 PFU/ml (Supplementary Table S1). This K–M analysis was performed for each displacement pair at every dose in every experiment, and significant differences were found for at least one virus pair in each experiment (Supplementary Table S1). The first pair 007:111 was significantly different in nine of twenty-four comparisons, the 111:110 pair was different in five of twenty-four comparisons, and the 110:139 pair was different in three of twenty-four comparisons. Collectively, these K–M results revealed significant differences in virulence at some challenge dose(s) for the first displacement pair in five of six experiments, the second pair in three of six experiments, and the third pair in three of six experiments. In the combined results of LD50 and K–M analysis each displacement pair had significant differences in virulence, and for the first two displacement pairs the significant differences always indicated that the displacing virus had higher virulence (Table 3).

Figure 3.

Comparison of K–M survivorship of immersion-exposed fish for three virus pairs in 2012 Lake Quinault steelhead trout. At each of three exposure doses (2e + 03 top row, 2e + 04 middle row, and 2e + 05 bottom row) the survivorship of fish exposed to either the previously dominant (black) or newly emerging (red) virus strain is shown for each of the three virus pairs 007:111 (left column), 111:110 (middle column), or 110:139 (right column). Comparisons for the 2e + 02 dose, and all doses for all remaining experiments, are listed in Supplementary Table 1. Dashed lines are 95% CI around solid line estimate of instantaneous probability of survival; + indicates end of study censoring; P-values from Wilcoxon tests, see Section 2.

3.3 Virulence measured through estimates of the proportional hazard of not surviving virus exposure

To determine how much relative risk is associated with either specific virus strains or fish populations in the survival of steelhead trout upon virus exposure, the Cox proportional hazard model was fitted to the mortality data. Using the proportional hazards model to explore relative virulence means the higher the relative risk value, the higher the virulence of the virus. Virus exposure was a significant covariate of relative risk, relative to mock infection (P < 0.001). Each virus strain conferred a significantly different level of increased risk (relative to mock) of death upon exposure, regardless of which fish population was infected (Fig. 4A). The relative risk increased with emergence time of the first three viruses. Thus, 007 conferred a 6.2-fold increased risk over mock, 111 conferred a 9.1-fold increased risk, and 110 was still more virulent, conferring a 11.7-fold increased risk of death. The relative risk associated with infection with 139 was not greater than the 110 virus, however, with a relative risk of 9.8-fold. When this proportional hazard model was used to infer a survivorship function (to incorporate time), the difference between the least virulent virus strain 007 and the other viruses was apparent early in the course of disease, as soon as 10 days post-exposure (e.g. at 10 days the 95% CI for the 007 proportional hazard survivorship curve no longer overlapped those of the other viruses—Fig. 4B). The most virulent virus 110 took longer to be distinguishable from the two intermediate viruses 111 and 139, at around 15 days post-exposure. The kinetics of the relative risk indicated that all virus strains except 111 and 139 differed significantly in virulence from one another (Fig. 4B).

Figure 4.

Both virus strain and fish population are significant covariates for the probability of survival after exposure to IHNV. The relative risk of death after exposure to each of the four virus strains (A, arranged in order of emergence time from left to right) or for each of the exposed fish populations (C) was calculated from the Cox proportional hazard model (error bars are standard error). Each proportional hazard model was then used to fit a survivorship curve for either infecting virus (B) or exposed fish population (D). Dashed lines are 95% CI around solid line estimate of instantaneous probability of survival;+ indicates end of study censoring.

Fish population, of any origin, was also identified as a significant cofactor on the proportional hazard function (relative to pooled mock exposed fish of all populations) when infected with any virus (P < 0.001). The LQ steelhead trout population (combined results from 2012 to 2014) was the least resistant, suffering a 12.4-fold greater risk of death relative to mock exposure (Fig. 4C). The next lowest resistance fish population was DW, with a relative risk of 10.5-fold. The QN and SK fish were similar, with 8.1- and 8.6-fold risks, respectively. The EW fish were the most resistant to mortality caused by any virus, with only a 4.6-fold increased risk. When the proportional hazard model was fit to a survivorship curve to assess the impact of the relative risk over time, the most resistant EW fish had a shallower survival curve apparent as early as 8 days post-exposure (Fig. 4D). The least resistant LQ fish also resolved from the other fish populations early, around 11 days post-exposure, and had a much steeper survival curve. The fish populations with intermediate relative risk (QN and SK) resolved from the least resistant fish (DW and LQ) somewhat later, around 15 days post-exposure. The kinetics of the relative risk indicated that all fish populations except SK and QN differed significantly in resistance from one another (Fig. 4D). The overall host population resistance ranking by proportional hazards analysis was EW > QN ≈ SK > DW > LQ. This matched the host resistance ranking by average LD50 values.

An underlying assumption in the Cox proportional hazards model is that the effect of covariates is constant over time. To test if this assumption is violated in our data, we applied the Cox–Aalen multiplicative model to test for time-varying signal in each of the three covariates of virus type, virus dose, and fish population. Virus type was the only covariate that was always proportional with respect to time (Table 4, constant effects). Virus dose and fish population, in at least some tests, did show a significant degree of time-varying effect on relative risk of death (Table 4, test for time invariant effects). However, since fish in these experiments are artificially exposed and infected at the same time (synchronized) and infection leads to disease and death or recovery within a time frame of 7–21 days, it may be that the disease process itself is driving the signal, since the longer a fish is infected/diseased the greater its risk of dying. It is likely there are additional processes acting in the fish population covariate, but the implications of this effect are not clear.

Table 4.

Cox–Aalen semi-parametric multiplicative hazards model: simulation based tests for time-invariant effects and estimates of constant effects for the three covariates of virus dose, virus strain, and fish population

Test for time invariant effectsa			Constant effectsb
Covariate	Test statistic	P-value	Covariate	ϒ	Robust SE
200 PFU/ml	1.680	0.0000	111 virus	0.495	0.118
2,000 PFU/ml	1.990	0.0554	110 virus	0.764	0.115
20,000 PFU/ml	2.510	0.0070	139 virus	0.570	0.119
200,000 PFU/ml	2.860	0.0000
EW fish	−0.841	0.2270
LQ fish	0.137	0.0002
QN fish	−0.392	0.0844
SK fish	−0.311	0.0010

aTime invariance is the null hypothesis, so where the P-value < 0.05, that covariate demonstrates a time-variant effect.

a,bCox–Aalen model tests performed in the timereg package of R, see Section 2.

3.4 Four virus strains did not differ in ability to infect exposed fish

Since another possible factor impacting ability to displace an older dominant virus is infectivity, we conducted studies to determine how much of each virus was needed to infect 50 percent of exposed animals (ID50). The hypothesis was that a displacing virus would have a lower ID50 estimate, requiring fewer infectious virions to infect 50 percent of exposed fish. These studies were conducted only in the LQ and QN fish populations, and the LQ fish were tested in two different years. In 2012, LQ fish were exposed to four doses of each of the four virus strains, and assessed by qRT-PCR for infection status at 3 and 7 days post-infection. At 3 days post-exposure, the percent of infected fish followed a dose–response pattern for each virus, though the slope of the dose–response curve varied by virus strain (Fig. 5A). By 7 days post-exposure, the percent of fish infected at each dose of each virus was lower than at 3 days post-exposure, indicating that virus was already being cleared (Fig. 5A, dotted lines). Therefore, all subsequent ID50 analyses were based on data collected 3 days post-exposure.

Figure 5.

Percent infected and estimate of the dose necessary to infect 50 percent (ID50) of immersion exposed fish. Tests were conducted in Lake Quinault steelhead trout in two different years (A–B in 2012 and C–D in 2014), and Quinault NFH steelhead trout (E–F) in 2014 only, by exposure to one of four virus strains: 007/mG007M (black), 111/mG111M (orange), 110/mG110M (red), or 139/mG139M (purple). The percent of exposed fish found to be infected 3 days (A, solid lines) or 7 days (A, dashed lines) after exposure to indicated dose of virus (PFU/ml). The ID50 estimates were calculated from infection status at day 3 post-exposure (right column, PFU/ml), see Section 2.

In the 2012 LQ experiment, ID50 estimates varied by 5-fold, between the highest at 6.32e+02 PFU/ml for 110 and the lowest at 1.2e+02 PFU/ml for 139 (Fig. 5B). This was a much smaller range than observed in the LD50 estimates in the same fish in the same year, which differed between highest and lowest by 26-fold (Fig. 2A). There was a trend of higher ID50 estimates with emergence time for the first three viruses, while the 139 virus did not have a higher ID50 estimate than 110. Despite this trend, there were no significant differences among the ID50 estimates of viruses in any of the three displacement pairs. When the LQ steelhead trout were tested again in 2014, the percent of fish infected by each dose of each virus again showed a strong dose–response with each virus having a different shaped curve (Fig. 5C). Compared with the curves from the 2012 experiment, the shapes of the curves in 2014 were slightly different for 007, 111, and 110, but very different for 139 (compare Fig. 5A and C). The ID50 estimates for each virus in LQ fish in 2014 were similar to those obtained in 2012, ranging between 2.0e+02 and 6.35e+02 PFU/ml. This was a 3-fold difference between highest and lowest ID50 estimates, which was again a much narrower range than LD50 estimates in the same year in the same fish, which differed by 44-fold between highest and lowest (Fig. 2B). Unlike the results from 2012, the LQ2014 ID50s showed no apparent trend of decreasing infectivity with time of viral emergence, and the relative ranking of ID50 values for the four virus strains differed between the two experiments. Again, there were no displacement pairs that differed significantly by ID50 estimates.

The QN steelhead trout were tested in 2014 with exposure to five doses of each virus. The percent of fish infected by each virus showed a dose–response pattern, and as in LQ, the shapes of the curves for each virus differed (Fig. 5E). The ID50 estimates ranged from 7.9e+01 to 3.21e+02 PFU/ml, a 4-fold difference between highest and lowest. As with the LQ fish, the range of ID50 estimates was much narrower than LD50 estimates in the same fish in the same year, where a range of 34-fold was observed between the highest and lowest estimates (Fig. 2C). There was no apparent trend between ID50 estimates and virus emergence time, and none of the displacement pairs were significantly different from one another (Fig. 5F).

To determine if, upon infection, any virus might appear to have a replicative advantage over another, the mean viral loads within virus-positive fish exposed to each dose of each virus were calculated from the 2014 LQ and QN experiments (Fig. 6). For higher exposure doses where more than a few fish were infected, no significant differences were found between viral loads of any viruses. At the lowest exposure dose of each experiment, there were some significant differences in the viral loads of the different viruses, but only a few fish were infected in these cases. Within each fish population, there was a general trend of increasing mean viral load with increasing virus exposure dose (Fig. 6). Also, viral loads were different between the two fish populations, with ∼10-fold more virus in LQ fish than in QN for all viruses tested.

Figure 6.

Mean viral load of virus-positive fish 3 days after exposure to different doses of four virus strains. Tests were conducted in Lake Quinault (A), and Quinault NFH steelhead trout (B) in 2014, by exposure to one of four virus strains: 007 (black), 111 (orange), 110 (red), or 139 (purple). Quantitative RT-PCR was used to determine the mean viral load of positive fish, see Section 2; error bars are standard error; dashed line indicates the lower limit of detection of the assay.

4 Discussion

4.1 IHNV traits that correlate with observed emergence and displacement

Thorough genetic surveillance of IHNV incidence in steelhead trout of the Pacific Northwest reveals a clear pattern of serial emergence and displacement of four distinct M genogroup types, but does not indicate what forces might be driving this pattern. Previous work using the same four IHNV strains used here demonstrated that various components of viral fitness tested in vivo did not correlate with the observed ability to displace a previously established virus type (Kell, Wargo, and Kurath 2014). Here, we tested whether increased infectivity or virulence correlated with emergence and displacement. The infectivity of the four virus strains revealed no significant differences or patterns consistent with emergence time. However, virulence varied widely among the four viruses and showed a consistent pattern of increased virulence with emergence time for the first three virus strains. This suggests that the four viruses have a similar capacity to infect fish but a widely different capacity to kill fish after infection, as indicated by the 26–44-fold differences in LD50 values. This observation suggested that the ability to cause mortality after infection is a viral trait that varies between IHNV strains. This trait is one definition of the term ‘virulence’, although operational definitions of virulence often vary based on the trait measured experimentally. Our operational definition of virulence was mortality due to exposure to virus, which actually measured the combined result of both infection and the ability to cause mortality once infected. Interestingly, this higher virulence was not associated with greater virus replication fitness, as indicated by the lack of differences in viral loads for the four viruses in our infectivity studies, and the previously reported viral fitness studies (Kell, Wargo, and Kurath 2014). This rules out many components of overall viral fitness as drivers of IHNV displacement (Wargo and Kurath 2012) but other components such as differences in transmission fitness and viral persistence are currently being investigated.

The most consistent pattern so far observed between quantitative viral traits and naturally occurring IHNV emergence/displacement was the pattern of variation of virulence in the four virus strains. Year-to-year tests within the same population showed variation in absolute LD50 values, but the relative virulence rankings of the four virus strains were the same. This is consistent with a previous study (Breyta, Jones, and Kurath 2014), and demonstrates expected levels of natural variation. Overall, increased virulence was associated with the first two displacement events. In contrast, increased virulence was not apparent in the conditional displacement pair 110:139. Instead, a general trend of decreasing virulence was found for this displacement pair, but this trend was poorly supported statistically. Since genotype 139 did not displace 110 in the lower Columbia region, the relative virulence of all four virus isolates correlated completely with the occurrence of displacement or non-displacement in the endemic lower Columbia River region.

4.2 Steelhead trout populations vary in resistance to IHNV mortality

Another significant finding of this study was the variation of IHNV resistance among steelhead trout host populations. The relative rankings of the host populations based on magnitude of mean LD50 estimates and proportional hazard analyses were similar, indicating that overall the host resistance rankings in this study can be summarized as EW > QN ≈ SK > DW > LQ. This ranking confirmed the previous report that QN steelhead trout are more resistant than LQ steelhead, despite having a common ancestry ∼40 years ago, and inhabiting the same coastal watershed, presumably with similar virus exposure history (Breyta, Jones, and Kurath 2014). Among the other host populations however, the ranking was not completely consistent with a previously published challenge study using only the 110 virus strain in four of the same host populations tested here: the EW population differed dramatically and the QN and SK populations differed slightly (Breyta, Jones, and Kurath 2014). While different year classes of a host population can be expected to vary, the changes in relative resistance rankings make interpretation of the long-term impact of host resistance difficult.

The host resistance variation presented here suggests two hypotheses for how resistance might play a role in emergence and displacement. The first is that the virus is acting as a selection pressure in the field sufficient to drive increased host resistance. If this is true, then relative resistance of different host populations should correlate with varying levels of historic exposure to IHNV. This appeared to be true in the previous study (Breyta, Jones, and Kurath 2014), but it is less clearly supported here due to the shift of the IHNV-naive EW population from least resistant to most resistant in host population ranking. In addition, the DW population from the virus endemic lower Snake River region was the second least resistant host populations in the ranking. The lack of correlation between infectivity and displacement indicates that if co-evolution between virus and host has driven virus evolution, it has not had an effect on this component of transmission. Overall with the present data we are less confident that there is a general pattern indicating virus selection on host resistance in the field, but it may occur in the endemic lower Columbia River region. At present, we have tested only one host population from that region (SK) and it had an intermediate level of resistance, so testing of additional populations would be needed to clarify whether there is evidence for co-evolution with host resistance.

The second hypothesis relevant to variation in host resistance is that low resistance host populations may have a higher probability of facilitating viral emergence events. This is not supported in the endemic lower Columbia River SK population that has experienced emergence of all four virus types tested here despite having intermediate resistance. However, outside this endemic region, viral emergence events have occurred repeatedly in both the DW and LQ populations, which are both of relatively low resistance. Thus, we would speculate that the role of host genetic resistance in viral emergence events may vary among different regions.

4.3 Source and sink dynamics in the natural ecology and geography of IHNV

Within the landscape considered here, the lower Columbia River region has been endemic with M group IHNV since the early 1980s (Garver, Troyer, and Kurath 2003; Breyta et al. 2013). Similarly, M group IHNV has been detected in the lower Snake River region since 1983. In contrast, the coastal Washington region is not endemic: M group IHNV has only rarely been detected, mostly during a recent transient emergence 2007–2011 (Emmenegger and Kurath 2002; Breyta et al. 2013). Endemicity may influence viral emergence via niche availability. The three regions are also different epidemiologically. Transmission inference based on midG typing has indicated that the lower Columbia River, where all four of the dominant M genotypes described here were first detected (Garver, Troyer, and Kurath 2003; Breyta et al. 2013), is a source region because the virus that continually circulates and evolves there routinely transmits to neighboring regions. The lower Snake and coastal Washington regions, where M group viruses are repeatedly introduced from the source region, are therefore sink regions (Breyta et al. 2013, 2016). In our results, virulence correlates completely with displacement or non-displacement in the endemic source region, but not in the sink regions. Interestingly, an epidemiological model of human vector-borne diseases recently noted that selection on pathogen virulence is predicted to be weak in sink vector populations (Rascalou et al. 2012).

In both the lower Snake River and coastal Washington sink regions, factors other than increased viral virulence must function in emergence and displacement events. As noted above, it is possible that presence of low resistance host populations (DW and LQ, respectively) may have facilitated viral emergence. In general, the expression of virulence depends strongly on the host environment, and wide variation in host resistance such as observed for the four host populations from sink regions (EW, LQ, QN, and DW) could overwhelm the input of selection for higher virulence. Ecological factors including virus founder effects, genetic drift, and other stochastic factors could also be strong contributors in some cases. For example, in the lower Snake River sink region M virus incidence is less constant than in the lower Columbia River region, and a void of all M virus detection occurred for some years before emergence of mG139M in 2008–2009 (Breyta et al. 2016). This suggests availability of an open niche and a stochastic founder event as possible factors contributing to this emergence. In the non-endemic coastal Washington sink region, the incidence of the emergent mG110M was waning by 2009, and M virus was maintained by recurrent introduction instead of an endemic pattern of local transmission. In this situation, it is likely that the large-scale emergence of mG139M in the Snake River region in 2008–2009 created an increased source of virus carried by migratory fish that contributed to the displacement of mG110M by mG139M in the coastal region in 2010 (Breyta et al. 2013). In other words, it appears that type mG139M may not have been capable of displacing mG110M where that type was endemic (source region) but only emerged and became dominant in sink regions where type mG110M was absent (Lower Snake) or non-endemic and waning (coastal Washington).

4.4 Challenges in representing complex ecology in a controlled laboratory setting

In designing experiments to investigate drivers of natural virus emergence and displacement events within a virus:host landscape, there are many variables that require consideration. In this study, we have tested four virus strains representing three different major displacement events, the viruses were tested in five different host populations, and virulence was assessed by exposure to four doses of each virus. However, we have tested only one virus strain from each of the four serially emergent virus genotypes, so a caveat to the interpretations presented here is that we cannot be sure the phenotypes of these strains are completely representative of the dominant genotypes responsible for the events observed on the landscape. Based on the intriguing patterns of virulence variation seen here, additional virus strains from each of the four genotypes have been selected for testing in the near future to determine if the phenotypes are consistent. Another caveat inherent in ‘common garden’ infection studies such as those presented here is that, by necessity, viruses isolated in the past (as much as 35 years ago) are being used to challenge fish populations from today. Although this temporal asynchrony may influence the observed levels of viral virulence or host resistance, this was not evident in the relative rankings by individual virus strain or fish host population.

4.5 What would cause IHNV virulence to increase?

If we accept that increased virulence was a driver for the two IHNV displacement events that occurred in the lower Columbia River region, the question then becomes, what was the mechanism that led to evolution of increased virulence? Current virulence evolution theory suggests several factors that are predicted to lead to evolution of increased pathogen virulence (Kurath and Wargo 2015). Among these factors, co-infection by more than one virus genotype, and co-evolution with increasing host resistance are the two that seem most likely to apply to the IHNV:steelhead trout system. We know from controlled laboratory infection studies that co-infection is easily established when fish are exposed to two different IHNV genotypes at the same time, even if they differ in fitness and virulence (Troyer et al. 2008; Wargo et al. 2010; Peñaranda, Wargo, and Kurath 2011). Superinfection also occurs upon sequential experimental exposures, albeit with lower frequency (Kell, Wargo, and Kurath 2014). Detection of co- infections in field samples is relatively infrequent, but it is found more commonly in situations with active viral infections in high densities of fish, as might occur in the endemic lower Columbia River region. The potential contribution of co- evolution of IHNV virulence and steelhead trout resistance is uncertain at present since we used only one host population from the endemic lower Columbia River region as noted above, but in theory conditions conducive to co-evolution have been in place in the endemic region for several decades.

4.6 Mechanisms of viral displacement

In terms of viral displacement events there are numerous virus:host systems where molecular epidemiology has identified displacements of dominant genotypes similar to the events presented here. In Dengue virus, displacement of major genotypes is associated in some cases with higher virulence (Rico-Hesse 2007). Experimental studies of other arboviruses have shown that in some cases displacing genotypes had higher viral fitness associated with enhanced replication or dissemination in the insect vector (Moudy et al. 2007; Hanley et al. 2008; Heise, Tsetsarkin, and Weaver 2011). However, there were also some studies that did not find increased fitness associated with displacement, similar to our results with IHNV (Myat Thu et al. 2005; Lourenco and Recker 2010; OhAinle et al. 2011). Beyond arboviruses, displacement associated with higher virulence in the poultry herpesvirus infectious laryngotracheitis virus was based on increased replicative fitness and transmission potential (Stewart et al. 2015). The IHNV work presented here contributes to the efforts to understand viral displacements as an example of increasing virulence without the often-assumed aspect of increasing replicative fitness.

4.7 Relevance to viral virulence evolution

There are relatively few examples where long-term viral virulence evolution has been clearly documented using controlled laboratory infection studies in ‘common garden’ hosts (Kurath and Wargo 2015). Both myxomavirus and rabbit hemorrhagic disease virus have evolved long-term increases in virulence (for myxomavirus after a period of virulence decrease) as a result of co-evolution in response to increasing resistance in their rabbit hosts (Kerr 2012; Elsworth et al. 2014). Avian influenza virus in chickens and infectious salmon anemia virus in Atlantic salmon have each repeatedly evolved by mutation from asymptomatic reservoir genotypes to virulent genotypes selected by favorable conditions in intensive farming (Alexander 2007; Markussen et al. 2008). Although selection pressure associated with conditions in agricultural settings may also be relevant to IHNV in cultured fish, IHNV virulence evolution presented here is longer-term, sequential, incremental increases in virulence, starting from genotypes that were already quite virulent. This pattern is more similar to that of Marek’s disease virus in chickens, where long-term sequential virulence increases from a virulent progenitor have been observed in response to conditions in poultry farming, including vaccination (Witter 1997). However, there is no IHNV vaccine in use for steelhead trout. Thus, the overall the driving factors suggested here for IHNV emergence and displacement differ in various ways from those proposed for other viral:host systems, making IHNV a useful system for broadening our understanding of virulence evolution.

Overall, our results suggest evolution of increased viral virulence as a driver of IHNV displacement in the endemic lower Columbia River source region, and other factors including variation in host resistance and ecological or stochastic factors contributing to emergence and displacement in the Lower Snake River and coastal Washington sink regions. Viral infectivity and numerous features of viral fitness do not appear to be drivers of IHNV displacement. Studies to quantitate other viral traits that could be responsible, like transmission potential and long-term persistence or recovery, are currently underway.