A stratified compartmental model for the transmission of Sparicotyle chrysophrii (Platyhelminthes: Monogenea) in gilthead seabream (Sparus aurata) fish farms

The rapid development of intensive fish farming has been associated with the spreading of infectious diseases, pathogens and parasites. One such parasite is Sparicotyle chrysophrii (Platyhelminthes: Monogenea), which commonly infects cultured gilthead seabream (Sparus aurata)—a vital species in Mediterranean aquaculture. The parasite attaches to fish gills and can cause epizootics in sea cages with relevant consequences for fish health and associated economic losses for fish farmers. In this study, a novel stratified compartmental epidemiological model of S. chrysophrii transmission was developed and analysed. The model accounts for the temporal progression of the number of juvenile and adult parasites attached to each fish, as well as the abundance of eggs and oncomiracidia. We applied the model to data collected in a seabream farm, where the fish population and the number of adult parasites attached to fish gills were closely monitored in six different cages for 10 months. The model successfully replicated the temporal dynamics of the distribution of the parasite abundance within fish hosts and simulated the effects of environmental factors, such as water temperature, on the transmission dynamics. The findings highlight the potential of modelling tools for farming management, aiding in the prevention and control of S. chrysophrii infections in Mediterranean aquaculture.


Introduction
Aquaculture has experienced steady growth in global fish production in recent decades and is expected to expand further in the coming years, also in view of the world's growing population [1]. In Europe, aquaculture contributed to 10% of the total fish production in 2020 [1]. However, Mediterranean aquaculture is also known to bear negative impacts on the environment and the marine habitat [2], impairing water and bottom sediment quality [3][4][5]. In addition, the expansion of fish farming in cages has been linked to the spread of infectious diseases, pathogens and parasites [6].
Gilthead seabream is the most produced species in European marine finfish aquaculture (36% of total production, 88 000 tons yr −1 ), with an estimated economic value of €473 million [7]. However, one of the challenges facing seabream farmers in the region is the emergence of infections of Sparicotyle chrysophrii, a flatworm parasite of the phylum Platyhelminthes. This ectoparasitic monogean flatworm attaches to fish gills and can cause lethal epizootics in sea cages [8] due to its blood-feeding activity. Severe infection can lead to hypoxia and anaemia, which in turn may induce a state of lethargy. The parasite also causes haemorrhages and necrosis [9][10][11]. Moreover, infected fish are more susceptible to secondary bacterial, viral or parasitic infections due to their compromised immune system [12]. The outbreak of S. chrysophrii can thus have important economic consequences for seabream farming.
S. chrysophrii is a hermaphrodite parasite. Once it releases eggs, these can either directly remain attached to the host or be released into the environment. In an intensive fish-farm setting, free-floating eggs can attach to cage nets, thus spatially constraining the transmission cycle and likely increasing infection rates. After hatching, the oncomiracidia (free-swimming ciliated larvae) can survive in the water column up to a couple of days [13] while seeking for a fish host. Previous studies showed how infections depend also on environmental factors, in particular on water temperature [8,14]. However, the literature reports contrasting effects of this variable on the parasite life cycle: warm temperature promotes fast parasite development, accelerating the spreading of infections [14], although disease outbreaks have been reported also in winter, when gilthead seabream are more immunodepressed, and, accordingly, more susceptible to the infection [8]. Furthermore, it has been reported that the parasite life cycle could be influenced by fish size and age [8].
Fish farmers routinely perform parasite surveillance and treatment when necessary. Surveillance is typically conducted through gill inspection focused on S. chrysophrii detection, while treatment commonly relies on anti-parasitic baths (usually with formalin, where permitted, or other suitable compounds) [15,16]. Other intervention strategies to limit parasite transmission are based on the substitution or cleaning of cage nets to avoid the trapping of parasite eggs. Alternative strategies against S. chrysophrii consist in functional feeding [17][18][19] aimed at preventing and/or treating against infection through dietary additives. Regardless of the type of action taken to treat the fish, to be effective, control efforts have to take into account the parasite life cycle, environmental conditions and infection dynamics [13,14,20].
In this context, epidemiological modelling potentially represents a valid tool to help fish farmers understand parasite transmission and, accordingly, design control measures. To our knowledge, no modelling studies exist on S. chrysophrii transmission. For this purpose, in the following, we propose and analyse a novel epidemiological model that investigates S. chrysophrii transmission within gilthead seabream farms. We propose a stratified compartmental model, taking into account parasite life cycle, environmental variables and aquaculture practices. This kind of approach has already been adopted in other studies on parasite infection dynamics [21][22][23][24] and can be adapted depending on the object of study and the epidemiological dynamics involved. We then applied the model to data collected in a seabream farm managed by Cromaris (Bisage, Croatia). pathogenicity and fertility of adult parasites, we differentiate the parasites in two classes: juveniles (attached larvae) and adults. As a result, the fish population is actually subject to a double stratification based on the abundance of juvenile and adult parasites hosted by each individual. Specifically, we define X j,a as the number of fish with respectively j juvenile and a adult parasites attached, with j = 0, 1, 2, …, J and a = 0, 1, 2, …, A (where J and A are the maximum number of juvenile and adult parasites, respectively). Approximating X j,a as a continuous function of time, X j,a (t), its dynamics for 0 < j < J and 0 < a < A can be written as follows: ð2:1Þ Equation (2.1) accounts for four input and six output fluxes, which are represented in figure 1 and detailed in the following. A parasite larva that attaches to a fish host, with j − 1 juvenile larvae attached, causes an instantaneous increase of X j,a and a related decrease of X j−1,a . The rate of such process is F = βM/V, where β is the exposure rate and M is the number of oncomiracidia within a control volume V. A juvenile parasite becomes adult at a rate r. Thus, in a fish hosting j + 1 juveniles, the rate at which one among the j + 1 juveniles become adult is r( j + 1), which is thus the rate of the transition from X j+1,a−1 to X j,a . Juvenile parasites die at rate μ j , whilst adults die at a rate μ a . Similarly to the rationale introduced for parasite maturation, the rate at which one of the j juvenile (or a adult) parasite dies is μ j j (or μ a a). The latter rates accounts for the transitions X j+1,a → X j,a and X j,a+1 → X j,a reported as input terms in equation (2.1). Output fluxes from the state variable X j,a are due to fish mortality, at a rate μ F , and fish excess mortality due to the parasite. The latter is assumed to be linearly proportional to the adult parasite abundance, aμ F 0, while extra-mortality possibly linked to the presence of juvenile stages is considered to be negligible. The remaining four negative fluxes reported in equation (2.1) (with overall rate: F + jr + jμ j + aμ a ) represent the output fluxes corresponding to the four input fluxes described earlier. The fluxes involving state transitions are described in table 1. If we imagine that the state variables X j,a are arranged in a matrix, where rows and columns represent juvenile and adult parasite abundance, respectively, special care has to be devoted to the equations representing the dynamics of the state variables at the corners of the matrix, i.e. having the maximum or minimum number of juveniles and adults ( j = 0, J and a = 0, A) (equations (2.2)-(2.5)): royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 221377 and and at the edges (equations (2.6)-(2.9)): which represent the boundary conditions for equation (2.1). The description of the parasite life cycle is completed by modelling two additional state variables: the number of parasite eggs in the system, E, and the number of oncomiracidia, M. The equations regulating their dynamics are as follows: Each adult parasite can release eggs with a shedding rate f, which can survive and hatch with probability ρ, for a total amount of rf P J j¼0 P A a¼1 aX j,a eggs released per unit time by the fish population. Eggs hatch into oncomiracidia at rate ξ, and eventually find a host, thus being transferred to the fish compartment (second term of the right-hand side of equation (2.11)), or are removed from the cage. The latter process occurs at a rate μ M , which accounts for both oncomiracidia mortality and possible exit from the cage.
In the absence of evidence of density-dependent effects, we assume that parasite growth, infection and death rates are independent of the parasite abundance (i.e. r, μ a , μ j , μ M , f and ξ do not depend on a and j).

Data
We applied the model to data collected in six cages (in the following referred to as cages 1-6) from a seabream farm managed by Cromaris d.d. The cages are located in Bisage, Croatia (44°01 0 30 00 N, 15°13 0 10 00 E), and have a unit volume of about 225 m 3 each. The experiment started in February 2021 with about 10 4 fish per cage with an average weight of 8.5 g. A repeated parasitological survey was run between February and November 2021: each month, 30 fish were collected from each cage (1800 total fish sampled) and, for each fish, all eight arch gills were examined to count the number of attached adult parasites. The fish stock size was re-evaluated each month during the experiment. At the end, the number of fish in the cages was approximately 9000, with an average weight of about 287 g. At each sampling, water temperature, T (°C), was also measured. Table 1. Compartment transitions and rates involving as a target the element X j,a of the fish matrix, X. Notice that fish death rates, μ F and μ F 0, are not included, since they do not imply a transition between compartments, causing rather a net loss. event event rate transition parasite transmission F X j−1,a → X j,a death of an adult parasite (a + 1)μ a X j,a+1 → X j,a death of a juvenile parasite ( j + 1)μ j X j+1,a → X j,a juvenile parasite maturation ( j + 1)r X j+1,a−1 → X j,a royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 221377

Model setup and calibration
To run the model, we set the maximum number of adult and juvenile parasites, respectively, as A = 30 and J = 15 individuals. These values do not affect the model output as long as the number of fish that reach the maximum intensity of infection remains close to zero during the simulation. We checked a posteriori that such condition was verified in our simulations. The initial population was set to H = 10 4 fish like in the experimental cages. We assumed that the fish mortality parameters, μ F and μ F 0, the exposure and juvenile maturation rates, β and r, and the hatching rate, ξ, possibly vary with temperature following a widely used power-law formulation. Accordingly, we set their reference values at 20 C (i.e. μ F20 , μ F 0 20 , β 20 , r 20 and ξ 20 , respectively) so that where the θ parameters express the sensitivity of each rate to water temperature. As initial conditions, we assume that M 0 oncomiracidia are present in the cage while all other state variables are null, except for the uninfected fish abundance. We thus set, as of February 2021 (t = 0): Parameter estimation relies on a Bayesian framework and is based on Markov chain Monte Carlo (MCMC) sampling. We used the DREAM algorithm (Differential Evolution Adaptive Metropolis) [25], which is an implementation of MCMC that runs multiple chains simultaneously, in order to efficiently explore the parameter space. In particular, we adopted the DREAM ZS version [26]. We assign a range of values where the algorithm is allowed to explore parameters up to convergence (a total of Oð3 Â 10 5 Þ iterations were run). We also assigned prior Gaussian marginal distributions to those parameters for which literature information is available (table 2).
To compute the likelihood of the observation we proceeded as follows. From the model simulation, we estimated the time-dependent probability that a randomly sampled fish hosts a adult parasites  Terming n a (m i ) the number of fish with a adult parasites (a ¼ 0, . . ., A) in the monthly sample m i of N = 30 fish, the log likelihood of a sample can be calculated using a multinomial probability distribution: n a ðm i Þ log p a ðtÞ: ð3:2Þ Finally, the log likelihood of the samples taken at one cage through the growing season is expressed as follows: Lðm i Þ: ð3:3Þ To better identify the parameters related to fish mortality, we included also data about the number of surviving fish in the likelihood computations. In particular, we assumed that the abundance of fish in a certain month follows a Gaussian distribution with an average equal to the abundance predicted by the model and a standard deviation equal to the one observed among the six cages. The six cages have slightly different epidemic trajectories. These differences could be attributed to slightly different environmental conditions such as net fouling and the related dissolved oxygen concentration. Moreover, different fish biomass growth trajectories, as observed in the data, could affect parasite susceptibility and fish mortality. We therefore allow some specific parameters to possibly assume different values across the different cages. Specifically, they are the mortality rate μ F,20 , the exposure rate β 20 , and the initial concentration of oncomiracidia M 0 . The remaining parameters are instead assumed to the the same for all cages.

Results
An example of observed distribution of parasite abundance among the sampled fish in a single cage is shown in figure 2. Analogous figures for the remaining five cages are shown in appendix A. At the beginning of the sampling period, all fish have zero parasites. Then, parasite abundance increases, and in summer, the mode of the parasite abundance distribution in the sampled fish reaches one. Figure 3 provides an overview of fish and parasite population statistics. At the end of the sampling period, about 95% of fish survive in all cages, with an observed increase in mortality between July and August. The prevalence, i.e. the fraction of fish with at least one parasite, varies markedly among cages. It is highest in cage 1, in which the share of infected fish reaches more than 60% in April, while the other cages reach high values more gradually, in late spring or early summer. In cages 5 and 6, infected fish prevalence does not exceed 50%. The mean parasite abundance (i.e. the mean number of parasites in the sampled fish [27]) shows a pattern similar to the prevalence, with mean abundance at the end of the period decreasing going from cage 1 to 6. Figure 3 reports the mean intensity, defined as the mean number of parasites in sampled infected fish [27], where the latter are fish hosts with at least one parasite. In cages 5 and 6, infected fish typically have a single parasite, while the mean intensity reached almost five in cage 1. The developed model was able to reproduce the distribution of parasite abundance (figures 2 and 5-9) as well as the main patterns of the summary statistics presented in figure 3. The marginal posterior distributions of the estimated parameters (figure 4) enable us to gain further insights into the processes controlling parasite dynamics and to highlight some model limitations. For some parameters, namely, the pathogen juvenile and adult mortality μ j and μ a , the egg hatching rate ξ 20 , and the oncomiracidia loss rate μ M , the posterior distribution shows a good overlap with the prior. This result implies that the dataset does not contain sufficient information to further characterize these process rates, and, in turn, that without prior information, these parameters could not be properly identified. Baseline fish mortality μ F shows some variability across cages, which leads to the slightly different observed survival trajectories (figure 3), but with a consistent overlap. Results also show that the differences in parasite abundance distribution among the different cages are explained by differences in the exposure rate β 20 and initial conditions M 0 . We also note that the royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 221377 posterior sample shows a negative correlation between the exposure rate β 20 and the initial condition M 0 (average correlation coefficient among cages equal to 0.51). This result implies that it is difficult to estimate both parameters, which control the concentration of oncomiracidia in the water, based solely on observations of parasite abundance in fish. Correlation coefficients for the other pairs of parameters are all less than 0.3.   Regarding the effect of temperature on crucial parasite and fish process rates, results indicate that higher temperature enhances fish mortality (median u m F ¼ 1:10) and egg hatching rate (median θ ξ = 1.06). By contrast, for the exposure and parasite maturation rates, the posterior distribution does not allow one to identify a clear temperature effect (i.e. temperature sensitivity parameter close to unity: median θ β = 0.98, median θ r = 0.99). As an example, to properly interpret temperature sensitivity values, θ = 1.1 implies an increase of 10% of the relevant parameters for every increase of 1 C of water temperature.
Finally, to understand the modelled contribution of parasite infection to overall fish mortality, we run model simulations using the estimated parameters but setting μ F 0 20 to zero (no parasite-related mortality). Results indicate that the percentage of deaths related to the parasite may range from 6% in cage 1 to 0.1% in cage 6.

Discussion
In this work, we have proposed a novel mathematical formulation to study the transmission dynamics of S. chrysophrii in S. aurata aquaculture farms, accounting for the whole life cycle of the parasite. The novelty of the proposed model relies on the double stratification of parasite abundance within the fish host population based on the number of juvenile and adult parasites attached to the gills of each fish. By using this model, we were able to account for the higher pathogenicity of adult parasites, as well as for the link between fish mortality and (adult) parasite abundance. We also included in the model the role of water temperature, which has been found to affect the parasite transmission [8,14]. We applied the model to data obtained from a controlled parasitological survey conducted in S. aurata farms located in Bisage (Croatia) and managed by Cromaris d.d. We estimated model parameters using a Bayesian approach that relies on MCMC. Results show that our model reproduces well the available data and the temporal patterns observed in infection dynamics. The results also emphasize the importance of combining modelling studies with laboratory and field experiments to enhance our understanding. In fact, parameter estimation did not provide additional information on some process rates that regulate the parasite's life cycle, as the posterior distribution overlapped with the prior distribution. Therefore, for the application of the model, it is important to rely on information obtained independently through experimental studies. This will help to improve the accuracy and reliability of the model, particularly in cases where parameter estimation cannot provide sufficient information.
We tested the possible temperature effects on critical parameters that regulate the transmission cycle of S. chrysophrii. Our estimation procedure revealed the egg hatching rate ξ as temperature dependent, which is consistent with experimental findings [14]. There was no clear evidence to support the temperature dependence of the parasite maturation rate r and the exposure rate β. However, it should  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 221377 be noted that this could be attributed to limitations in the specific modelling and estimation setup, and further experimental and modelling studies may provide clearer evidence. Compartmental epidemiological models, like the one formulated in this study, have long been shown to be flexible, yet parsimonious, tools to model a wide array of parasite dynamics and life cycles [21,28,29]. Macroparasitic infections pose specific challenges because the number of parasites within each host matters and should thus be properly accounted for in the model formulation. To this end, stratified compartmental models have been proposed [21][22][23]. So far, this class of models was developed and tested on human parasitic diseases like schistosomiasis [23,24]. The main challenge in these former applications was that the number of parasites per host could not be directly observed, and therefore the full potential of these tools could not be properly tested. In this regard, the application presented herein offers the ideal setting to test this class of models. Indeed, in farms, the fish population is confined and monitored, and the monthly samples of 30 fish per cage, whose eight arc gills were all examined to count the number of adult ectoparasites, allowed us to reconstruct the distribution of parasite abundance within hosts, which is directly related to the state variables of stratified compartmental models.
One crucial benefit of compartmental models, versus, say, purely statistical models, is that the mechanistic representations of the key processes controlling the parasite life cycle allow one to test the effect of alternative intervention strategies to control the transmission. For instance, the proposed model could be readily adapted to represent treatments like anti-parasitic bath that periodically (or triggered by the surveyed parasite abundance) abate the parasite population in the host. Replacement or cleaning of the cage nets could be represented in the model as a sudden decrease of egg abundance. Net clogging could be represented as a reduction of the flux of eggs that leave the cage. Different intervention strategies, or combinations thereof, could be tested, possibly in an optimal control mathematical framework [30,31], to select the most efficient, accounting also for the ensuing costs.
Infection with S. chrysophrii has been reported to reduce the appetite of seabream and affect its metabolism, which in turn affects its growth trajectory. This variable is crucial in fish farming, and a promising avenue for future development is to couple the proposed epidemiological model with a metabolic growth model (e.g. [32][33][34]). Average fish size or size distribution is commonly monitored in fish farms, and when combined with a well-designed S. chrysophrii surveillance system like the one presented here, these data can be used to estimate the parameters of a coupled model that also includes growth dynamics. Such a tool could potentially account for the differential susceptibility of fish of different sizes and the effect of fish biomass density on S. chrysophrii transmission, features that cannot be addressed in the current model formulation. We believe that such a tool could more comprehensively account for the full spectrum of the impacts of this parasitic infection and provide a valuable tool for precision fish farming [34].
All authors gave final approval for publication and agreed to be held accountable for the work performed therein.      royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 221377