An epidemiological model for proliferative kidney disease in salmonid populations

Abstract

Background

Proliferative kidney disease (PKD) affects salmonid populations in European and North-American rivers. It is caused by the endoparasitic myxozoan Tetracapsuloides bryosalmonae, which exploits freshwater bryozoans and salmonids as hosts. Incidence and severity of PKD in brown trout populations have recently increased rapidly, causing a decline in fish catches and local extinctions in many river systems. PKD incidence and fish mortality are known to be enhanced by warmer water temperatures. Therefore, environmental change is feared to increase the severity of PKD outbreaks and extend the disease range to higher latitude and altitude regions. We present the first mathematical model regarding the epidemiology of PKD, including the complex life-cycle of its causative agent across multiple hosts.

Methods

A dynamical model of PKD epidemiology in riverine host populations is developed. The model accounts for local demographic and epidemiological dynamics of bryozoans and fish, explicitly incorporates the role of temperature, and couples intra-seasonal and inter-seasonal dynamics. The former are described in a continuous-time domain, the latter in a discrete-time domain. Stability and sensitivity analyses are performed to investigate the key processes controlling parasite invasion and persistence.

Results

Stability analysis shows that, for realistic parameter ranges, a disease-free system is highly invasible, which implies that the introduction of the parasite in a susceptible community is very likely to trigger a disease outbreak. Sensitivity analysis shows that, when the disease is endemic, the impact of PKD outbreaks is mostly controlled by the rates of disease development in the fish population.

Conclusions

The developed mathematical model helps further our understanding of the modes of transmission of PKD in wild salmonid populations, and provides the basis for the design of interventions or mitigation strategies. It can also be used to project changes in disease severity and prevalence because of temperature regime shifts, and to guide field and laboratory experiments.

Background

Epidemics of emerging diseases may have large economic and ecological impacts, threaten livelihoods, elicit biodiversity losses and affect key ecosystem services. Understanding causes and patterns of disease emergence becomes even more topical as causative and correlative links with global climate and environmental change are established [1–3]. Proliferative Kidney Disease (PKD) of salmonid fish, caused by Tetracapsuloides bryosalmonae (phylum Cnidaria, class Malacosporea) [4], is highly problematic for hatcheries and fish farms, potentially reaching 100 % infection prevalence and high mortalities (up to 95 %) in affected fish [5, 6]. Although the impacts of PKD in wild fish populations are still poorly known, PKD is considered a key factor contributing to the decline of wild salmonid populations in Switzerland and Northern Europe [7–9].

The parasite life-cycle alternates between freshwater bryozoans, as primary hosts [10–12], and salmonid fish, where infection can cause PKD [4, 13]. The transmission between hosts occurs through spores that are released into water. Development and pathology of PKD are dependent on water temperature, as clinical signs and disease-related mortality increase with increasing water temperatures [4, 5, 9, 14–16]. Water temperature is also key to the development of T. bryosalmonae in its bryozoan host. Increased temperatures have been shown to promote the production of spores infective to fish [17, 18]. The connections between temperature, disease development and fish mortality suggest that, due to climate change, the prevalence, severity and distribution of PKD are likely to increase further [19]. Hence, PKD must be considered as an emerging disease having a sizeable and growing impact on the health of salmonid populations and, consequently, on the economy of fish farming industry and hatcheries. In the perspective of understanding drivers and controls of the disease, and in the search for mitigation strategies, the development of a dynamical model of PKD transmission becomes crucial.

Epidemiological models incorporating parasites with simple life-cycles have a long history [20, 21] and have been successfully applied to several human diseases and zoonoses. Conversely, models incorporating parasites and pathogens with complex life-cycles and multiple hosts are less frequently developed, although many major human and wildlife diseases are caused by parasites with such life-cycles. Successful application of epidemiological models for these diseases fundamentally relies on the incorporation of the ecological dynamics of the different host populations. In the context of fish diseases, epidemiological models have been exploited to address various issues, such as sea lice infection on salmon farms [22]; Koi herpes virus of the common carp (Cyprinus carpio) [23]; amyloodinosis, a disease of warm water mariculture [24]; fish-borne zoonotic trematodes in agriculture-aquaculture farms [25]; Ceratomyxa shasta, a myxozoan parasite endemic to river systems throughout the Pacific northwest region of North America [26]; and whirling disease, a myxozoan disease of farmed and wild salmonids [27]. As for PKD, the only modeling attempt used a Bayesian probability network to assess the decline of brown trout, citing PKD as a driver of increased mortality [28].

The present work proposes an epidemiological model capable of describing the intra- and inter-annual dynamics of PKD. The model aims to translate the current knowledge of the modes of transmission of the disease and of the life-cycle of its causative agent into a mathematical form, making the connection between temperature and epidemiological parameters explicit. In this way, it is possible to state under which conditions PKD can establish in a fish population. The development and the analysis of such a model is also instrumental to guiding further experimental and field studies on PKD.

Methods

In this section, we first describe the transmission cycle of PKD, and subsequently present the mathematical model that stems from such biological and epidemiological evidence.

PKD transmission cycle

Tetracapsuloides bryosalmonae has a complex life-cycle that exploits freshwater bryozoans and salmonids as hosts (for review, see [29]). Within bryozoans, the parasite can express either covert or overt infection stages. During covert infections, the parasite exists as non-virulent single-cell stages [30, 31]. The transition to overt infection implies increase in virulence, with the formation of multicellular sacs from which thousands of T. bryosalmonae spores (approximately 20 μm in diameter) are released into water. Spores are characterized by two amoeboid cells and four polar capsules [32, 33]. Parasite transmission from bryozoans to fish can thus take place only during the overt infection stage. Overt infection hampers bryozoan growth, whereas covert stages pose low energetic cost to bryozoans [17]. Peaks of overt infection have been observed in late spring and autumn [34]. In particular, overt stages develop when bryozoans are undergoing enhanced growth as a result of optimal temperatures or food levels [18, 35]. Overt infection also elicits temporary castration [36], with a severe reduction in the production of statoblasts (asexually produced dormant propagules). Covertly infected bryozoans can produce infected statoblasts thus allowing vertical transmission of T. bryosalmonae [37]. Recovery mechanisms for bryozoan colonies are poorly observed and understood [37].

Parasite spores released into water by overtly infected bryozoans infect fish through skin and gills [38, 39]. Tetracapsuloides bryosalmonae subsequently enters the kidney of fish hosts, where it undergoes multiplication and differentiation from extrasporogonic to sporogonic stages [40]. Spores developed in the lumen of kidney tubules contain one amoebid cell and two polar capsules [41], and are eventually excreted via urine [42]. Although PKD seems to develop in all salmonids to a varying degree, life-cycle completion may be highly species-specific. For example, in Europe, parasite spores infective to bryozoans develop in the brown trout (Salmo trutta) and the brook trout (Salvelinus fontinalis), but not in the rainbow trout (Oncorhynchus mykiss) [43–45]. Fish infected with PKD often die owing to secondary infections [13]; however, PKD alone has been shown to cause mortality [14, 15]. Fish that do not die during the acute phase of the infection may become long-term carriers of the parasite. Such carriers are reportedly able to infect Fredericella sultana, one of the most common bryozoan hosts of T. bryosalmonae, for a period up to 2 years after exposure [43].

Model

The proposed model couples PKD transmission and population dynamics of fish and bryozoan populations. When reproduction processes are concentrated in time (as in our case, in which spawning and hatching for fish, and statoblast hatch for bryozoans occur mainly during the cold season), population dynamics are traditionally analyzed with discrete-time models (see the widely known Ricker model [46], introduced to study stock and recruitment in fisheries). On the other hand, PKD transmission occurs continuously throughout the warm season, thus calling for a continuous modelling approach (e.g. the classical Susceptible-Infected-Recovered model [20]). We therefore propose a discrete-continuous hybrid framework that couples intra-seasonal and inter-seasonal dynamics. Seasons can be thought of as the periods of bryozoan proliferation (e.g. from April to November, depending on climate). Within-season dynamics are described by a set of coupled ordinary differential equations expressing how bryozoan biomass, fish density, and the abundances of infective spores and statoblasts change during the warm season. The transition between the end of a warm season and the beginning of the next one is modelled as a discrete-time update. Therefore, between-season dynamics are described by a set of difference equations. An interesting framework to formalize discrete-continuous hybrid models is the so-called time scale calculus [47]; however, the complexity of the model presented in the following paragraphs prevents an effective use of such formalism. In this work, we focus on a local-scale model, able to mimic PKD dynamics in an isolated water body, where fish and bryozoans are well mixed, and there are no additional inputs or outputs. An outline of the model is shown in Fig. 1. The variables of the model are listed in Table 1; all state variables are referred to as concentrations.

Fig. 1

Schematic representation of PKD transmission dynamics. a Within-season dynamics. The red-shaded area highlights the bryozoan sub-model. The circular arrows illustrate bryozoan growth. The dotted line indicates that the growth of overtly infected bryozoans is impaired by the parasite. The blue-shaded area highlights the fish sub-model. The dead-end arrow refers to excess PKD-induced fish mortality. In both sub-models, straight solid lines indicate fluxes among epidemiological classes, dashed arrows represent disease transmission between the two sub-systems. Natural mortality is not displayed. Parameters driving the transition between classes are displayed in gray. Letters T and τ indicate fluxes that are expected to be dependent on water temperature T and time τ. b Schematic representation of between-season (overwintering) dynamics relevant to PKD transmission. Colors as in (a). Continuous lines represent survival over winter. Dashed-dotted lines stand for reproduction processes. c Graphic timeline of the discrete-continuous hybrid model

Table 1 List of model variables

Within-season model

The dynamics of the state variables during the warm season (from time τ ₀ to τ ₁) of year y (see Fig. 1c) is described by the following system of first-order differential equations:

$$ \frac{d{B}_S}{d\tau }={g}_S\left({B}_S,{B}_C,{B}_O;T\right){B}_S-{\beta}_B{Z}_F{B}_S+\psi\;\left({B}_C+{B}_O\right); $$

(1a)

$$ \frac{d{B}_C}{d\tau }={g}_C\left({B}_S,{B}_C,{B}_O;T\right){B}_C+{\beta}_B{Z}_F{B}_S-\left[{d}_{CO}(T)+\psi \right]{B}_C+{d}_{OC}(T){B}_O; $$

(1b)

$$ \frac{d{B}_O}{d\tau }={g}_O\left({B}_S,{B}_C,{B}_O;T\right){B}_O+{d}_{CO}(T){B}_C-\left[{d}_{OC}(T)+\psi \right]{B}_O; $$

(1c)

$$ \frac{d{S}_S}{d\tau }={f}_B\left(\tau, T\right){B}_S+\varphi\;{f}_B\left(\tau, T\right){B}_C; $$

(1d)

$$ \frac{d{S}_I}{d\tau }=\left(1-\varphi \right){f}_B\left(\tau, T\right){B}_C; $$

(1e)

$$ \frac{d{F}_S}{d\tau }=-{\mu}_F{F}_S-{\beta}_F{Z}_B{F}_S+\zeta\;{F}_C; $$

(1f)

$$ \frac{d{F}_E}{d\tau }={\beta}_F{Z}_B{F}_S-\left[{\mu}_F+h\;(T)\right]{F}_E; $$

(1g)

$$ \frac{d{F}_I}{d\tau }=\left(1-\varepsilon \right)h\;(T){F}_E-\left[{\mu}_F+a\;(T)+\gamma \right]{F}_I; $$

(1h)

$$ \frac{d{F}_C}{d\tau }=\varepsilon\;h(T){F}_E+\gamma\;{F}_I-\left({\mu}_F+\zeta \right){F}_C; $$

(1i)

$$ \frac{d{Z}_B}{d\tau }={\pi}_B{B}_O-{\mu}_Z{Z}_B; $$

(1j)

$$ \frac{d{Z}_F}{d\tau }={\pi}_F\left({F}_I+\kappa\;{F}_C\right)-{\mu}_Z{Z}_F, $$

(1k)

where the dependence of the parameters on temperature T and time τ has been explicitly expressed.

The first terms in the right-hand sides of Eqs. (1a), (1b) and (1c) express the growth of the bryozoan biomass, which is assumed to be logistic-like. Thus we set

$$ {g}_X\left({B}_S,{B}_C,{B}_O;T\right)={r}_X(T)\left[1-\rho \left({B}_S+{B}_C+{B}_O\right)\right], $$

(2)

where r _X (T) and ρ are the baseline instantaneous growth rate of class X = {S, C, O} and the inverse of the carrying capacity of the bryozoan population, respectively. We assume that susceptible and covertly infected bryozoans have the same baseline growth rate r _S  = r _C  = r, whereas r _O  << r since the growth of overtly infected bryozoans is strongly impaired [36]. Both r and r _O are assumed to be monotonically increasing functions of temperature, as experimental evidence suggests [17]. The term β _B B _S Z _F in Eqs. (1a) and (1b) represents the flux of bryozoan biomass from the susceptible to the covertly infected class, with β _B being the exposure rate of bryozoans to fish-released spores. The transition from covert to overt infection is assumed to occur at a rate d _CO [see Eqs. (1b) and (1c)], which is defined as the inverse of the mean time necessary for the development of the overt stage of infection in a previously covertly infected bryozoan unit. The rate of transition from overt to covert infection is expressed by d _OC . Both parameters are temperature-dependent: in particular d _CO increases with increasing temperature, while d _OC decreases [18]. As previously stated, the development of an overt infection, which poses a high energetic cost to the bryozoan colony, also depends on host conditions and food availability. However, as a first approximation, temperature is considered as the only determinant of change of the two transition rates. Infected bryozoans can possibly clear the infection and become again susceptible at a rate ψ [Eqs. (1a), (1b) and (1c)]. The production of statoblasts is only achieved by susceptible and covertly infected bryozoans, whereas overtly infected colonies do not produce statoblasts. In particular, susceptible bryozoans produce uninfected statoblasts S _S [Eq. (1d)] at a rate f _B that is assumed to depend on both temperature and time; specifically, the release of statoblasts is enhanced towards the end of the season. Seemingly, uninfected bryozoans have been observed while producing infected statoblasts as well. However, Abd-Elfattah et al. [37] argued that their observation could be interpreted as a sign of a recent loss of infection. Instead, covertly infected bryozoans can produce both infected S _I [Eq. (1e)] and uninfected statoblasts [37]. We term φ the probability that infected bryozoans produce uninfected statoblasts. The total statoblast production rate of infected bryozoans is assumed to be equal to that of susceptible ones. T. bryosalmonae spores [Eq. (1j)] are produced by overtly infected bryozoans at a constant rate π _B . A constant mortality rate for spores (μ _Z , i.e. the inverse of the time span during which spores are viable) is also accounted for.

As for the fish population, no recruitment or immigration are considered to take place during the warm season. Therefore, the abundance of susceptible fish [Eq. (1f)] decays monotonically throughout the season, due to natural mortality (at a rate μ _F , defined as the inverse of the average lifespan of a fish) and infection from bryozoan-released spores Z _B . The term β _F F _S Z _B in Eqs. (1f) and (1g) represents the flux of fish from the susceptible to the exposed compartment, with β _F being the exposure rate of fish to bryozoan-released spores. Exposed fish have already contracted the disease but are not yet infective. The decrease in the abundance of these fish, besides natural mortality, is ruled by the temperature-dependent rate h [Eq. (1g)]; namely, 1/h is the average time for the development of the disease in fish. The parameter h increases with increasing temperature [48]. We assume that a fraction ε of fish exiting from the exposed class does not show clinical symptoms and is not subject to PKD-related mortality, thus directly entering the carrier class. The remaining part (1 − ε) becomes infected. The abundance of infected fish [Eq. (1h)] decreases because of both natural and PKD-caused mortality. The latter is assumed to occur at a rate a that is positively correlated with temperature [14, 15]. In addition, infected fish can enter the carrier compartment at a rate γ [Eq. (1i)], which corresponds to the inverse of the average duration of the acute phase of the infection. Fish infected by the parasite that did not die in the first year reportedly continue to release infective spores for a period up to two years [43] and no evidence of infection clearing has been observed. However, field and experimental observations on the duration of this carrier stage are still scarce and further experiments are underway (Strepparava, personal communication). Therefore we explore the possibility that carrier fish may recover by introducing a recovery rate ζ. Evidence of possible immunity is scant (e.g. Foott & Hedrick observed immunity to second infection in rainbow trout [49]), thus, as a safe assumption, recovered fish are assumed to enter the susceptible compartment again [Eq. (1f)]. Spores [Eq. (1k)] are released by both infected and carrier fish at rates π _F and κ π _F , respectively, with κ being an appropriate coefficient ranging from 0 to 1. Spore decay is accounted for through the parameter μ _Z .

Since the dynamics of disease transmission between bryozoan and fish are driven by the product between the concentration of spores (Z _B and Z _F ) and the rates (β _F and β _B ) at which susceptible organisms are actually exposed to the infectious agents, we can introduce two new state variables Z ^* _B  = β _F Z _B and Z ^* _F  = β _B Z _F . These new quantities are termed equivalent spores: namely, Z ^* _F (Z ^* _B ) [T⁻¹] is the concentration of spores needed to infect a unit concentration of susceptible bryozoan biomass (a single susceptible fish in a unit water volume) per unit time. With this definition, the exposure rates β _F and β _B can be discarded and two synthetic contamination rates π ^* _B  = β _F π _B [L³ M^-1 T^-2] and π ^* _F  = β _B π _F [L³ T^-2] are introduced. Note that both exposure and contamination rates as defined in model (1) are hardly measurable and would most likely need to be calibrated by contrasting model simulations against experimental or field data. The new parameter definitions thus reduce the number of parameters of the model and make the comparison between data and model predictions easier and more robust. The new set of equations accounting for the rescaled state variables therefore reads:

$$ \frac{d{B}_S}{d\tau }={g}_S\left({B}_S,{B}_C,{B}_O;T\right){B}_S-{Z}_F^{*}{B}_S+\psi \left({B}_C+{B}_O\right); $$

(3a)

$$ \frac{d{B}_C}{d\tau }={g}_C\left({B}_S,{B}_C,{B}_O;T\right){B}_C+{Z}_F^{*}{B}_S-\left[{d}_{CO}(T)+\psi \right]{B}_C+{d}_{OC}(T){B}_O; $$

(3b)

$$ \frac{d{F}_S}{d\tau }=-{\mu}_F{F}_S-{Z}_B^{*}{F}_S+\zeta {F}_C; $$

(3c)

$$ \frac{d{F}_E}{d\tau }={Z}_B^{*}{F}_S-\left[{\mu}_F+h\;(T)\right]{F}_E; $$

(3d)

$$ \frac{d{Z}_B^{*}}{d\tau }={\pi}_B^{*}{B}_O-{\mu}_Z{Z}_B^{*}; $$

(3e)

$$ \frac{d{Z}_F^{*}}{d\tau }={\pi}_F^{*}\left({F}_I+\kappa\;{F}_C\right)-{\mu}_Z{Z}_F^{*}. $$

(3f)

These equations, coupled with Eqs. (1c), (1d), (1e) (1h) and (1i) constitute the within-season model hereafter applied.

Between-season model

The following difference equation system relates the state of the model variables at the end of a season (y; τ ₁) with that at the beginning of the following season (y + 1; τ ₀).

$$ {B}_S\left(y+1;{\tau}_0\right)={\sigma}_B{B}_S\left(y;{\tau}_1\right)+\nu\;{S}_S\left(y;{\tau}_1\right); $$

(4a)

$$ {B}_C\left(y+1;{\tau}_0\right)={\sigma}_B{B}_C\left(y;{\tau}_1\right)+{\sigma}_O{B}_O\left(y;{\tau}_1\right)+\nu\;{S}_I\left(y;{\tau}_1\right); $$

(4b)

$$ {B}_O\left(y+1;{\tau}_0\right)=0; $$

(4c)

$$ {S}_S\left(y+1;{\tau}_0\right)=0; $$

(4d)

$$ {S}_I\left(y+1;{\tau}_0\right)=0; $$

(4e)

$$ {F}_S\left(y+1;{\tau}_0\right)={\sigma}_F{F}_S\left(y;{\tau}_1\right)+{f}_F\left({F}_S\left(y;{\tau}_1\right),{F}_E\left(y;{\tau}_1\right),{F}_I\left(y;{\tau}_1\right),{F}_C\left(y;{\tau}_1\right)\right); $$

(4f)

$$ {F}_E\left(y+1;{\tau}_0\right)=0; $$

(4g)

$$ {F}_I\left(y+1;{\tau}_0\right)=0; $$

(4h)

$$ {F}_C\left(y+1;{\tau}_0\right)={\sigma}_F\left[{p}_E{F}_E\left(y;{\tau}_1\right)+{p}_I{F}_I\left(y;{\tau}_1\right)+{F}_C\left(y;{\tau}_1\right)\right]; $$

(4i)

$$ {Z}_B\left(y+1;{\tau}_0\right)=0; $$

(4j)

$$ {Z}_F\left(y+1;{\tau}_0\right)=0. $$

(4k)

At the beginning of a new season (y + 1; τ ₀), no overtly infected bryozoans are present [Eq. (4c)]; the same condition applies for statoblasts [Eqs. (4d) and (4e)] and spores [Eqs. (4j) and (4k)]. Statoblasts surviving for more than one year are here neglected for the sake of simplicity. The biomass of susceptible bryozoans [Eq. (4a)] is composed of a fraction σ _B of the susceptible bryozoan biomass at the end of the previous season (y; τ ₁) that managed to survive over winter, and of newly hatched colonies from the uninfected statoblasts released during the previous season. The parameter ν is defined as the mean amount of bryozoan biomass produced by a single statoblast.

The population of covertly infected bryozoan at the beginning of a new season [Eq. (4b)] is given by the sum of three contributions: the fraction of the covertly infected biomass at the end of the previous season that survives over winter (with probability σ _B ), the fraction of overtly infected bryozoans at time (y; τ ₁) with survival probability σ _O , the newly established colonies generated by infected statoblasts, according to the parameter ν.

As for fish, we assume that a fraction σ _F of susceptible and carrier organisms at the end of the previous season survives over winter [Eqs. (4f) and (4i)]. Exposed and infected fish either die or enter the carrier class. This is modelled by computing additional coefficients p _E , p _I accounting for natural and PKD-induced deaths:

$$ {p}_I=\frac{\gamma }{\mu_F+\widehat{a}+\gamma };\kern2em {p}_E=\frac{\widehat{h}}{\mu_F+\widehat{h}}\;\left[\varepsilon +\left(1-\varepsilon \right)\frac{\gamma }{\mu_F+\widehat{a}+\gamma}\right]. $$

(5)

Specifically, p _E (p _I ) is the probability that an exposed (infected) fish survives and enters the carrier class in the first period of the winter season. In Eq. (5) â and ĥ are the rates of PKD-caused mortality and disease development averaged over the first period of the winter season, respectively. For the sake of simplicity, we assume $\widehat{a}=a\;\left(\widehat{T}\right)$ and $\widehat{h}=h\;\left(\widehat{T}\right)$, with $\widehat{T}$ being a representative value of temperature over the considered period.

Fish reproduction occurs during winter; we assume that the amount of newborn uninfected fish (termed f _F ) depends on the total fish population at (y; τ ₁):

$$ {f}_F\left(\tilde{F}\right)=\eta\;\tilde{F} \exp \left(-\xi\;\tilde{F}\right), $$

(6)

where $\tilde{F}={F}_S\left(y;{\tau}_1\right)+{p}_E{F}_E\left(y;{\tau}_1\right)+{p}_I{F}_I\left(y;{\tau}_1\right)+{F}_C\left(y;{\tau}_1\right)$. Eq. (6) assumes a Ricker growth model [46]. The parameter η is the baseline fertility rate, i.e. the average number of offspring produced by a single fish when the total fish population size is low, while ξ determines the strength of density dependence. For the aforementioned reasons, at the beginning of a season there are neither exposed nor infected fish [Eqs. (4g) and (4h)]. Table 2 summarizes all parameters of systems (1) and (4).

Table 2 List of parameters

Model simulation

In order to understand possible patterns of disease dynamics, we ran model simulations. The feasible range of several model parameters was estimated based on literature values and experts’ knowledge. Reasonable values were assumed for the remaining ones. Reference parameter values and feasible ranges are reported in Table 3. In the absence of detailed information about recovery mechanisms for fish and bryozoans, we assumed slow recovery rates as reference values (average recovery time equal to half of the lifetime for fish and half of the yearly proliferation period for bryozoan) and explored a wide range of possible values, including the absence of infection-clearing mechanisms (ψ = ζ = 0). As for temperature-dependent parameters, linear (r, r _O ) or parabolic (d _CO , d _OC , h, a) relationships have been assumed; a superlinear dependence on T for the latter parameters was deduced from literature review (see references in Table 3). The functional forms used are illustrated in Fig. 2. Note that, for these parameters, the corresponding reference values are obtained from the relationships of Fig. 2 by assuming T = 15 °C. For the numerical simulations, we use a time series of stream water temperature measured in River Langeten, Switzerland, in the period 2002–2013 (data provided by the Swiss Federal Office for the Environment - FOEN).

Table 3 Parameter ranges and reference values

Fig. 2

Functional forms for temperature-dependent parameters. The maximum values of these parameters are set to match the upper limits proposed in Table 3. Furthermore, we assume ${f}_B=0.005\kern0.22em T\kern0.28em \left[{}^{\mathrm{o}}\mathrm{C}\right]+0.00025\kern0.22em \tau \kern0.28em \left[\mathrm{d}\right]$, where τ is the time elapsed since January 1

Analysis of the within-season model

Parasite invasion in a previously PKD-free system most likely occurs during the warm season. For this reason, we first focus our analysis on the continuous, within-season model and investigate the conditions under which the introduction of the parasite in a fully susceptible population of fish and bryozoans (say through the introduction of spores, infected fish or infected bryozoans) leads to an outbreak of PKD. To this end, we analyze the linear stability of the disease-free equilibrium (DFE) in a simplified model that disregards population dynamics of both fish and bryozoans (i.e. both population sizes are treated as model parameters). Technical details are provided in Additional file 1. The instability of the within-season DFE would indicate long-term persistence of the parasite only if model parameters were kept constant at the value used for the computation of the stability criterion. However, in our system, population dynamics and water temperature variations constantly affect parameters values. Also, the alternation of warm and cold seasons, each endowed with different eco-epidemiological processes, de facto prevents a stable, nontrivial steady-state from being established in the system. Therefore, the instability of the within-season DFE can be seen as a useful indicator only of the short-term invasibility of the system. To investigate long-term parasite persistence and the establishment of endemic transmission conditions, a more complex analysis that accounts for both seasonal forcing and intra- and inter-seasonal population and epidemiological dynamics is required, as illustrated in the following section.

Stability analysis of the full model

As the seasonal cycle of water temperature critically controls PKD transmission and bryozoan population dynamics, we study the stability of the disease-free trajectory (DFT), a succession of system states in which the disease is not present and endemic transmission is not possible. To achieve this aim, we define a Poincaré map coupling between-season update with the within-season dynamics, with the latter being described through an integral operator derived from Floquet theory [50–52], a mathematical framework that allows the analysis of systems of periodically forced differential equations. As this theory assumes periodic forcing, the time evolution of water temperature is approximated via a sinusoidal signal with yearly period.

Disease-free trajectory

The DFT is defined as a discrete-continuous trajectory describing the yearly periodic evolution of biomass or abundances for the uninfected classes (B _S , S _S , F _S ) in the absence of PKD. Therefore, we first focus on a reduced system where only susceptible compartments are taken into account. Let x _S(y; τ) = {B _S (y; τ); S _S (y; τ); F _S (y; τ)} be the state of this system at season y and time τ. The time-evolution of this system can be represented as follows:

$$ {\mathbf{x}}_{\mathrm{S}}\left(y;{\tau}_0\right)\overset{v\left(\cdot \right)}{\to }{\mathbf{x}}_{\mathrm{S}}\left(y;{\tau}_1\right)\overset{w\left(\cdot \right)}{\to }{\mathbf{x}}_{\mathrm{S}}\left(y+1;{\tau}_0\right)\overset{v\left(\cdot \right)}{\to }{\mathbf{x}}_{\mathrm{S}}\left(y+1;{\tau}_1\right) $$

where v(⋅) is an operator describing the continuous-time dynamics of the system from τ ₀ to τ ₁, while w(⋅) is an operator describing the discrete-time update of the system between the end of season y and the beginning of season y + 1. Therefore, we obtain:

$$ {\mathbf{x}}_{\mathbf{S}}\left(y+1;{\tau}_0\right)=w\left(v\left({\mathbf{x}}_{\mathbf{S}}\left(y;{\tau}_0\right)\right)\right)=\left(w\circ v\right)\left({\mathbf{x}}_{\mathbf{S}}\left(y;{\tau}_0\right)\right), $$

where w ∘ v defines a so-called Poincaré map [53]. A fixed point of the Poincaré map is defined as an equilibrium $\overline{{\mathbf{x}}_{\mathbf{S}}}\left(y;\tau \right)$ such that $\overline{{\mathbf{x}}_{\mathbf{S}}}\left(y;{\tau}_0\right)=\overline{{\mathbf{x}}_{\mathbf{S}}}\left(y+1;{\tau}_0\right)=\left(w\circ v\right)\left(\overline{{\mathbf{x}}_{\mathbf{S}}}\left(y;{\tau}_0\right)\right)$.

The disease-free trajectory $\overline{{\mathbf{x}}_{\mathbf{df}}}\left(y;\tau \right)$ for systems (1) and (4) is defined under the hypothesis that water temperature can be approximated by a sinusoidal function of time, with period equal to one year. We therefore obtain

$$ \overline{{\mathbf{x}}_{\mathbf{S}}}\left(y;\tau \right)=\left\{\frac{B_{S,0}}{\rho {B}_{S,0}+\left(1-\rho {B}_{S,0}\right) \exp \left(-R\left(\tau \right)\right)},0,{F}_{S,0} \exp \left(-{\mu}_F\left(\tau -{\tau}_0\right)\right)\right\}, $$

(7)

where

$$ R\left(\tau \right)={\displaystyle {\int}_{\tau_0}^{\tau }}rdt $$

and B _S,0 (F _S,0) is the susceptible bryozoan (fish) population at τ ₀ along the DFT. Note in fact that $\overline{{\mathbf{x}}_{\mathbf{S}}}\left({\tau}_0\right)=\left\{{B}_{S,0},,,0,{F}_{S,0}\right\}$. B _S,0 is the solution of the following trascendental equation:

$$ 1-\frac{\sigma_B}{\rho {B}_{S,0}+\left(1-\rho {B}_{S,0}\right) \exp \left(-R\left({\tau}_1\right)\right)}-{\displaystyle {\int}_{\tau_0}^{\tau_1}}\frac{\nu {f}_B}{\rho {B}_{S,0}+\left(1-\rho {B}_{S,0}\right) \exp \left(-R\left(\tau \right)\right)}d\tau =0, $$

The value of B _S,0 can be obtained numerically for a given set of parameters. Instead, an analytical expression can be found for F _S,0:

$$ {F}_{S,0}={\xi}^{-1} \ln \left[\frac{\eta }{ \exp \left({\mu}_F\left({\tau}_1-{\tau}_0\right)\right)-{\sigma}_F}\right] \exp \left({\mu}_F\left({\tau}_1-{\tau}_0\right)\right). $$

(8)

Therefore, $\overline{{\mathbf{x}}_{\mathbf{df}}}\left(y;\tau \right)$ is the trajectory such that its susceptible components are equal to (7), while its infected components x _I  = {B _C ; B _O ; S _I ; F _E ; F _I ; F _C ; Z ^* _B ; Z ^* _F } are null. Note that Eq. (8) requires η > η _min  = exp(μ _F (τ ₁ − τ ₀)) − σ _F in order to avoid extinction of fish. For the reference parameter set of Table 3, one gathers η _min  ≈ 0.205.

Stability of the disease-free trajectory

In order to study the stability of the DFT $\overline{{\mathbf{x}}_{\mathbf{df}}}\left(y;\tau \right)$ of model (1) and (4) with respect to small perturbations, we refer to the reduced system where only infected compartments are considered. Hence, x _I (y; τ) denotes the state of such a system, whose time evolution is represented as follows:

$$ {\mathbf{x}}_{\mathbf{I}}\left(y;{\tau}_0\right)\overset{\mathbf{V}\left(\cdot \right)}{\to }{\mathbf{x}}_{\mathbf{I}}\left(y;{\tau}_1\right)\overset{\mathbf{W}\left(\cdot \right)}{\to }{\mathbf{x}}_{\mathbf{I}}\left(y+1;{\tau}_0\right)\overset{\mathbf{V}\left(\cdot \right)}{\to }{\mathbf{x}}_{\mathbf{I}}\left(y+1;{\tau}_1\right). $$

(9)

V and W are operators describing the evolution of the linearized system either within or between seasons. Note that the fixed point for this system is the null vector.

The propagation matrix V is given by the solution of the following matrix ODE system:

$$ \frac{d\;\mathbf{Z}}{d\tau }=\mathbf{A}\left(\tau \right)\mathbf{Z}, $$

(10)

where A(τ) is a matrix obtained from the Jacobian of the within-season system (1), evaluated along the trajectory $\overline{{\mathbf{x}}_{\mathbf{df}}}\left(y;\tau \right)$, where all rows and columns related to the susceptible compartments have been removed:

$$ \mathbf{A}\left(\tau \right)=\left[\begin{array}{cccccccc}\hfill {A}_{11}\hfill & \hfill {d}_{OC}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {B}_S\hfill \\ {}\hfill {d}_{CO}\hfill & \hfill {A}_{22}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill \left(1-\varphi \right){f}_B\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -{\mu}_F-h\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {F}_S\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \left(1-\varepsilon \right)h\hfill & \hfill {A}_{55}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \varepsilon\;h\hfill & \hfill \gamma \hfill & \hfill -{\mu}_F-\zeta \hfill & \hfill 0\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill {\pi}_B^{*}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -{\mu}_Z\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {\pi}_F^{*}\hfill & \hfill \kappa\;{\pi}_F^{*}\hfill & \hfill 0\hfill & \hfill -{\mu}_Z\hfill \\ {}\hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \end{array}\right], $$

where A ₁₁ = r(1 − ρ B _s ) − (d _CO  + ψ), A ₂₂ = r _O (1 − ρ B _s ) − (d _OC  + ψ), A ₅₅ = − μ _F  − a − γ. Equation (10) must be integrated from τ ₀ to τ ₁, with initial condition Z(τ ₀) = I (identity matrix). We thus have V = Z(τ ₁).

Matrix W can be obtained from the Jacobian matrix of the between-season system (4), by disregarding all rows and columns related to susceptible compartments. W is a sparse matrix of order 8 with the following non-null entries: W ₁₁ = σ _B , W ₁₂ = σ _O , W ₁₃ = ν, W ₆₄ = p _E σ _F , W ₆₅ = p _I σ _F , W ₆₆ = σ _F .

The stability of the disease-free trajectory of system (1) and (4) corresponds to the stability of the fixed point x _I  = 0 of the Poincaré map defined in (9), i.e. x _I (y + 1; τ ₁) = VW x _I (y; τ ₁). Therefore, the necessary and sufficient condition for the exponential instability of the DFT reads

$$ {\lambda}_{max}= \max \left|\lambda \left(\mathbf{V}\mathbf{W}\right)\right|>1, $$

(11)

where λ are the eigenvalues of matrix VW. When condition (11) is met, the parasite can invade the system. We therefore focus our analysis on λ _max , i.e. the maximum modulus of the eigenvalues of VW. For a given set of parameters, the corresponding value of λ _max can be computed numerically.

Sensitivity analyses

To understand how model parameters affect parasite invasibility and long-term PKD impact, we performed sensitivity analyses. Specifically, we investigate the effect of parameters on the value of λ _max and on the PKD-induced fish loss, here estimated as the percentage of fish at the end of the season with respect to the fish population size if the disease were absent. Computations are run by varying two focus parameters at a time while keeping the others at their reference value, as specified in Table 3. Temperature-dependent parameters are expressed via the functional forms of Fig. 2. The effect of these parameters is explored by varying their value at 15 °C, while keeping their minimum value at 0 °C (at 25 °C for d _OC ) constant. For parabolic functional forms, the null derivative at their minimum is also kept constant. The sinusoidal signal of water temperature is derived from the time series shown in the top panel of Fig. 3. Seasons are assumed to start on April 1 and last for 200 days. The effect of water temperature is also investigated by varying both the yearly mean value and the relative mid-amplitude of the sinusoidal signal, while using the relationships of Fig. 2.

Fig. 3

Model simulation. Evolution of the state variables of model (1) and (4) in a 12-years-long simulation. Seasons start on April 1 and last for 200 days. Black solid lines represent the time evolution of the overall bryozoan and fish populations; the dashed line indicates fish abundance if PKD is absent. Initial condition: B _S  = 0.495ρ ^-1, B _C  = 0.005ρ ^-1, F _S  = 0.99ξ ^-1, F _C  = 0.01ξ ^-1; other state variables are null. Temperature-dependent parameters are taken from Fig. 2; other parameters are set to their reference values (reported in Table 3). The water temperature time series used for the simulation is reported in the top subplot

Results

Model simulation

Figure 3 shows an example of model simulation where it is assumed that at the beginning of the first season 1 % of the bryozoans are covertly infected and 1 % of the fish belong to the carrier class. The remaining fractions of the host populations are susceptible. This model setting mimics the invasion of the parasite in a fully susceptible system. The total bryozoan population tends to reach the carrying capacity ρ ^-1 towards the end of the season. Overt stages of infection are absent in early spring but peaks are observed in summer. The irregular shape of the curves for B _C and B _O mirrors the variability of water temperature. Regarding fish, peaks of infection are observed during summer. After a few years, prevalence reaches about 80 % at the end of the summer, where few susceptible fish are present and most of the survived individuals belong to the carrier class. This is in agreement with the fact that, although young-of-the-year fish surviving the infection are not likely to develop clinical PKD in the following year, the parasite can remain viable in the host for several seasons after initial exposure [16, 43, 54]. In this model simulation, some 28 % of the fish population which is alive at τ ₀ dies during the season because of PKD. About ten years after the initial invasion, the disease becomes endemic and the seasonal state variables’ trajectories remain virtually unchanged. Overall, the model reproduces patterns of disease spread in bryozoan and fish populations which are in good qualitative agreement with evidence from the literature and field observations (see [19]).

Analysis of the within-season model

The analysis of the within-season model (Additional file 1) reveals two types of instability of the DFE. The first is related to the sign of the quantity

$$ \mathcal{T}={\psi}^2+\left[{d}_{CO}+{d}_{OC}-\left(r+{r}_O\right)\left(1-\rho {B}_S\right)\right]\psi -\left(1-\rho {B}_S\right)\left[{r}_O{d}_{CO}+r{d}_{OC}+r{r}_O\left(1-\rho {B}_S\right)\right], $$

which only depends on parameters regarding the bryozoan sub-model (note that the biomass of susceptible bryozoan B _S is assumed as a model parameter for this analysis). When $\mathcal{T}<0$, the DFE is unstable and the parasite can spread in the bryozoan population even in the absence of the fish host. Notably, if bryozoan cannot recover (ψ = 0), $\mathcal{T}$ is always negative, as long as the bryozoan biomass is lower than the carrying capacity. When $\mathcal{T}>0$, the parasite needs to cycle between the two hosts to possibly invade the system and the instability of the DFE occurs when the reproductive number

$$ \mathcal{R}=\frac{B_S{F}_S{d}_{CO}h{\pi}_B^{*}{\pi}_F^{*}\left[\varepsilon \kappa a+\kappa \gamma +\left(1-\varepsilon \right)\zeta +\left(1+\varepsilon \kappa -\varepsilon \right){\mu}_F\right]}{\mathcal{T}\left(h+{\mu}_F\right)\left(a+\gamma +{\mu}_F\right)\left(\zeta +{\mu}_F\right){\mu}_Z^2} $$

is greater than unity. Starting from the reference parameter set, $\mathcal{R}$ increases as parameters F _S , d _CO , h, π ^* _B , π ^* _F , increase, and decreases as parameters ρB _S , d _OC , a, ψ increase (see Additional file 1). Figure 4 illustrates how $\mathcal{R}$ evolves during a season, as temperature affects parameter values, and the population sizes of fish and bryozoan follow the DFT. The effect of temperature dominates and maximizes the reproductive number, and the related risk of an outbreak, during the warmest period. With the reference parameter set, $\mathcal{T}>0$ and $\mathcal{R}$ > 1 throughout the season: the parasite requires the fish host to spread into the system and the DFE is always unstable. Figure 4 also shows how reducing one of the transmission parameters (π ^* _B or π ^* _F ) or the fish population size F _S by a factor 10 and 20 can curb the reproductive number below unity during the coolest periods and throughout the season, respectively.

Fig. 4

Short-term risk of invasibility. Reproductive number $\mathcal{R}$ as a function of time. Water temperature is approximated as a sinusoidal signal (as for the computation of the DFT) and temperature-dependent parameters are taken as in Fig. 2. F _S (τ) and B _S (τ) follow the DFT. The red line refers to the reference parameter set reported in Table 3. The blue and the green lines show the effect of reducing one of the transmission parameters (π ^* _B or π ^* _F ) or the fish population size F _S by a factor 10 and 20, respectively

Sensitivity analyses of the full model

Parasite invasion

An exploration of the values of λ _max over wide ranges of the parameter space is reported in Fig. 5. Numerical results show high invasibility of the system (λ _max  > 1) for wide ranges of realistic parameters. The DFT becomes stable if one of the contamination rates (π ^* _B and π ^* _F , Fig. 5a) or if the characteristic sizes of bryozoan or fish populations (ρ ^-1 and ξ ^-1, Fig. 5b) are small. Recovery mechanisms (ψ and ζ, Fig. 5f) can promote the stability of the DFT: in particular, the DFT is stable if infection-clearing in bryozoans is fast enough, whereas if only ζ is increased, the persistence of the infection in bryozoans hinders the stability of the DFT. The DFT is predicted to be stable also if the rate of transition from overt to covert infection d _OC is considerably higher than the antithetic rate d _CO (Fig. 5c), although this circumstance occurs for unlikely values of these parameters. Finally, stability of the DFT is observed for extremely low values of the fish reproduction and recovery rates (η and γ, Fig. 5e).

Fig. 5

Sensitivity analysis: Parasite invasion. Values of λ _max as a function of model parameters π _B ^*​ vs. π _F ^*​ (a); ρ vs. ξ (b); d _OC vs. d _CO (c); a vs. h (d); γ vs. η (e); ψ vs. ζ (f); temperature (g). Pink dashed lines identify feasible parameter ranges. Black dots refer to the reference parameter set. For the sake of clarity, all rates are expressed as mean times (i.e. by their inverse). With regards to temperature-dependent parameters, their value at 15 °C is displayed. Black solid lines in g identify levels of mean temperature during the warm season

When the DFT is unstable, the maximum modulus λ _max of the eigenvalues of matrix VW can be interpreted as the rate at which model trajectories depart from the disease-free trajectory. λ _max thus represents a measure of how fast the parasite spreads into a susceptible population. When λ _max is slightly greater than unity, the outbreak develops slowly and stochastic effects (e.g. due to the demographic stochasticity of fish and bryozoan infected populations), which are not accounted for in the current model formulation, may lead to the extinction of the disease. When the DFT is stable, the closer λ _max is to unity, the slower perturbations to the DFT fade out and the system returns to a disease-free state. Therefore, λ _max can also be seen as an approximate estimate of the actual risk of invasion. In this perspective, Fig. 5 provides information on the role of the model parameters in promoting the establishment of PKD in fish populations. π ^* _B and π ^* _F have analogous effects in enhancing the risk for an outbreak. ρ ^-1, ξ ^-1 and the fish baseline fertility rate η are also positively correlated with λ _max . Concerning the effects of transmission and mortality rates, the highest values of λ _max are found when the disease development rate h is maximum and PKD-caused mortality a is minimum (Fig. 5d). This condition maximizes the number of non-fatal infection in fish, thereby causing an increase in the release of spores Z _F . As for the thermal regime, the mean water temperature during the warm season stands as the main factor controlling λ _max (Fig. 5g).

PKD-induced fish loss

Figure 6 shows the results of the sensitivity analysis of PKD impact on fish population size. The residual population size is larger (i.e. PKD impact is reduced) for low values of the contamination rates and the bryozoan carrying capacity (Fig. 6a and b), whereas no sensitivity to high values of these parameters is observed. The bryozoan baseline growth rate (Fig. 6b) has no effect on PKD-induced fish loss. Fish population is preserved when d _OC is high and d _CO is low, while the opposite case does not result in a severe population decay (Fig. 6c). Expectedly, there is a positive correlation between PKD-induced fish loss and PKD-induced mortality rate a (Fig. 6d). However, if a is extremely high, the immediate death of the host limits the transmission of the disease and thus its impact, a common feature of epidemiological models. PKD impact is enhanced by the rate of disease development h. Fish population decay is also enhanced when both ε and γ are low, as this condition maximizes the abundance of acutely infected fish (Fig. 6e); as previously pointed out, if γ is extremely low, this effect is balanced by the augmented relative importance of a, which weakens PKD impact. As the recovery time of bryozoans, ψ ^-1, decreases, PKD impact monotonically decreases (Fig. 6f). On the other hand, PKD impact reaches a peak for intermediate rates of fish recovery (ζ ≈ 0.005 d^-1), presumably due to concomitant high level of contamination and increased availability of susceptible fish.

Fig. 6

Sensitivity analysis: PKD-induced fish loss as a function of model parameters π _B ^*​ vs. π _F ^*​ (a); ρ vs. r (b); d _OC vs. d _CO (c); a vs. h (d); γ vs. ε (e); ψ vs. ζ (f); temperature (g). Simulations are run until convergence (100 seasons), with each season lasting for 200 days. Colors refer to the percentages of the population size at the end of the last season with respect to the same quantity calculated along the disease-free trajectory [computed as in Eqs. (7) and (8)]. For example, 50 % means that at the end of the last season the population size is half of the population that would have survived if the disease were absent. State variables at the beginning of the first season are set as in the model simulation of Fig. 3. Pink dashed lines identify feasible parameter ranges. Black dots refer to the reference parameter set. All rates are represented as mean times. With regards to temperature-dependent parameters, their value at 15 °C is displayed. Black solid lines in g identify levels of mean temperature during the warm season

As expected, higher water temperatures promote higher fish mortality (Fig. 6g). In particular, PKD impact is mostly related to the mean temperature during the warm season, whereas the amplitude of the temperature sinusoidal signal plays a lesser role. For instance, assuming a relative mid-amplitude of 0.5, an increase of 5 °C (from 10 to 15 °C) in the seasonal mean water temperature would produce an overall reduction of 26 % of the total fish population.

Discussion

According to the modelling framework developed in this work, the disease-free trajectory is found to be unstable over wide ranges of the ecological and epidemiological parameters. This property means that the introduction of the parasite in a fully susceptible community of salmonids and bryozoans would very likely lead to a PKD outbreak and to long-term parasite establishment. The ability of the parasite to cycle between covert, non-virulent, and overt, transmissive phases in the bryozoan host allows the parasite to produce large numbers of transmissive stages without compromising the susceptible host population. Similarly, long-term infections in trout are also responsible for continuous infections of naive bryozoan hosts and promote parasite persistence. This result holds for the deterministic model (1) and (4). Stochastic effects, currently not accounted for in the model formulation, could actually prevent the invasion. Our findings obviously depend on the current knowledge of the transmission modes of PKD (and on how it has been translated into mathematical terms), and the lack thereof. Indeed, the analyses of the within-season and the full models show how recovery mechanisms for both bryozoans and fish reduce the risk of outbreak and can make the disease-free system stable. However, knowledge and investigations on these critical processes are still scant. While the existence of a recovery dynamic for bryozoans is currently unknown (but see [37]), some experiments conducted on rainbow trout [55] actually revealed hints of a possible recovery from infection in salmonids. Further studies on this topic are needed to better elucidate the mechanisms underlying the persistence of T. bryosalmonae in fish and bryozoan communities. Overall, these results underline the reasons for the emerging status of PKD throughout Europe - the extremely successful invasion mechanism of the parasite may have facilitated its historical spread in European streams, and the current changes in climate and other environmental conditions are causing it to proliferate.

Our results highlight how both outbreak risk and the long-term persistence of the parasite are critically enhanced by warmer water temperature. The design of possible control strategies (e.g. the control of the population of one of the host or the interactions among them) should thus take into account that warmer periods pose a higher risk of a PKD outbreak. On the other hand, the effect of temperature offers the opportunity to design alternative intervention strategies aimed at controlling stream water temperature through e.g. tree shading or the selective release of cold water from upstream reservoirs. The analysis of the within-season model also shows that, for certain parameter combinations, the parasite can initially spread even in the absence of the fish host. This possible behavior is determined by the fact that the model assumes that infected bryozoan grow and that T. bryosalmonae can simultaneously proliferate inside them. This specific dynamics should be better scrutinized with additional experimental studies because the possibility of epidemics hosted solely by the bryozoan population has relevant implications for control strategies. It implies, for instance, that strategies focusing only on the fish population or on limiting the interactions between the two hosts might not be effective, under certain conditions, to prevent the invasion of the parasite.

Simulations and sensitivity analyses provided insights about the role and importance of the different parameters in parasite invasion and outbreak severity. In particular, the analysis of PKD-induced fish loss highlights the crucial role of the parameters that govern the transition between epidemiological classes in determining the impact of PKD outbreaks when the disease is endemic. While the carrying capacity of bryozoans and the contamination rates are the main factors controlling the stability of the disease-free trajectory, they appear less relevant in determining PKD impact in endemic settings. The underlying reason for this result is that these parameters solely control the rate at which susceptible fish are exposed to the parasite, whereas the infection develops over time scales ruled by the parameters governing the transition between epidemiological classes. A promising aspect is the fact that such parameters and their temperature dependence could be estimated in laboratory experiments by exposing hosts to the parasite and by monitoring the development of their infection status. Similar studies have already been carried out (see corresponding references in Table 3). More are needed, however, to better elucidate the role of temperature on PKD dynamics. It should be noted, in fact, that most of the available literature on PKD is based on rainbow trout (Oncorhynchus mykiss), which is known to be a dead-end host for T. bryosalmonae, as this species does not allow the completion of the parasite cycle by further infecting bryozoans [56]. Conversely, only few studies have focused on brown trout [43, 56, 57].

While it could be relevant to experimentally assess how the release of spores (by both infected fish and bryozoans) varies with temperature, the estimation of the actual value of the release rate may be less important. Indeed, the analysis of the model reveals that transmission dynamics are controlled by the product between contamination and exposure rates. The latter can hardly be estimated under field conditions as they depend on the probability of a successful encounter between viable spores and hosts. Therefore, efforts to precisely estimate contamination rates would be frustrated by the large uncertainties associated to exposure rates. The proposed rescaled model (3) features two parameters (one per host) that represent the product between contamination and exposure rates. Such parameters are key to disease dynamics, as discussed above. However, reasonable values can hardly be estimated a priori, and would need to be calibrated for each case study by contrasting model simulations with experimental or field data.

Knowledge and literature on bryozoans are rather scant compared to the vast and traditional literature on population dynamics and habitat distribution of salmonids. To improve our understanding of PKD as an emerging disease, such knowledge gap must be filled. In particular, knowledge of bryozoan habitat suitability needs be improved for an effective mapping on the risk of PKD invasion. This task could be achieved by using species distribution models [58] to relate the presence/absence of bryozoans to environmental variables (e.g. hydrological conditions, characteristics of the river-bed and banks, water quality, temperature). The few existing studies (e.g. [59]) focus on lakes, while large-scale studies of habitat suitability in streams and rivers are not available yet. Furthermore, many aspects of the relationship between infection status and bryozoan proliferation are still to be elucidated: for example, population genetic diversity of bryozoans may be linked to their ability to resist to infection; infective stages may propagate in partially infected bryozoan zooids even in the absence of further exposure to T. bryosalmonae.

Our results need be accompanied by an assessment of the limits of the model. It is acknowledged that young-of-the-year fish are highly likely to contract PKD when exposed to the parasite for the first time [60]. However, for the sake of simplicity, the age structure of the fish population has not been accounted for in the current model formulation. The fish population can indeed be split into age-specific classes (e.g. juveniles and adults), each with its own epidemiological compartments with different mortality, recovery and infection rates. To better elucidate this aspect, experimental investigations are currently being conducted to assess possible age-structure effects on T. bryosalmonae spore load in freshwater systems (Strepparava, personal communication). Our epidemiological model also relies on the simplifying hypothesis that all statoblasts hatch at the beginning of the following season. Hence, another aspect needing further investigation is the possibility of the existence of a statoblast bank (as suggested by Freeland et al. [61]). A statoblast bank could form on the bed of the water body if statoblasts were able to survive and maintain their infection status for more than one season. Dormant infected statoblasts could thus trigger new PKD outbreaks even if the disease were no longer present in fish and bryozoans.

Although the described transmission processes imply proximity between spores and hosts, several mechanisms allow long-distance spreading of PKD along river networks and lakes, among which bryozoan fragmentation, buoyancy of surfaces with attached colonies, fish migration, and hydrodynamic dispersal of T. bryosalmonae spores and bryozoan statoblasts may play a remarkable role [19]. To understand the role played by spatial processes in PKD transmission, the local model presented herein should be extended by considering the hydrological connectivity among different river reaches. In particular, modelling efforts at the river network scales (corroborated by field analyses) should try to understand whether PKD transmission is a spatially diffuse process or rather concentrated in transmission hot-spots. Such hot-spots could be represented by favorable habitats for both salmonids and bryozoans. The proximity of the hosts in these habitats could promote high transmission rates that may in turn be able to sustain the infection at the network scale through dispersal of fish and hydrodynamic transport of parasite spores. The possible identification of such hot-spots could offer alternative strategies to control the disease in the wild, like the confinement or the removal of one of the hosts in the isolated hot-spots.

Conclusions

We have developed a discrete-continuous hybrid epidemiological model that incorporates the current knowledge about T. bryosalmonae life-cycle and PKD transmission in wild salmonid populations. Our analysis showed that, for realistic parameter ranges, parasite invasion and long-term establishment are very likely, leading to a high risk of PKD outbreaks in susceptible environments. When the disease is endemic, the impact of PKD outbreaks is shown to be mostly controlled by the rates at which the disease develops in the fish population. These parameters can be estimated via experimental and field studies. Further studies are needed to gain additional insights on key disease dynamics, in particular concerning possible recovery in both hosts.