Introduction

Bivalve fisheries management, population restoration and aquaculture industries have recognized losses to predation as an important factor determining their productivity (Murray et al. 2007; Marino et al. 2009; Agüera et al. 2012; Wilcox and Jeffs 2019). Predators with a capacity to top-down control prey populations are capable of decimating populations of wild bivalves and aquaculture stocks (Menge 1982; Witman et al. 2003). Seastars of the genus Asterias are well-known keystone predators with the capacity to exert a control on their prey population, consisting mainly of bivalves (Menge 1982). This aspect has led them to be considered pests by shellfish producers and fishers (Agüera 2015) and to be catalogued among the most devastating invasive species (Byrne et al. 1997; Agüera and Byrne 2018). The impact of Asterias spp. in their prey is mediated by the capacity of seastar populations to rapidly fluctuate according to the presence of resources. Outbreaks, swarms or feeding fronts that eradicate prey populations are commonly described for Asterias rubens (Sloan 1980; Witman et al. 2003).

The impact of keystone predators on species of commercial value has to be considered in the management of shellfish resources including the development of mitigating actions to reduce or control their impact (Agüera 2015). Assessment and mitigation of the impact of A. rubens in northern Europe shellfish aquaculture and fishing have been the driver of several projects (Hancock 1958; Dolmer 1998; Gallagher et al. 2008; Agüera 2015). This interest has resulted in a large number of studies on A. rubens growth, feeding, thermal response, osmotic stress, reproductive dynamics and larval development (Barker and Nichols 1983; Nichols and Barker 1984; Wong and Barbeau 2005; Agüera et al. 2012). These results offer a description of the biology and physiology of this species, pointing out the particularities that make of A. rubens a species difficult to manage, such as its indeterminate growth, its highly plastic body continuously adapting to available prey, its resilience to adverse food and environmental conditions and its high reproductive output. Managing the impact of A. rubens on shellfish populations requires the capacity to use the existing knowledge on its ecology, physiology and population dynamics. However, the use of currently available descriptive data and statistical models is not enough to assess the potential impact of extremely plastic A. rubens populations under variable environments. There is a need for tools capable of assessing A. rubens individual and population growth, reproductive output, predation pressure exerted and other key traits under varying biotic and abiotic ecosystem components.

The metabolic processes under the varying environments are the underlying mechanisms behind species abundance, distribution and life history (Kearney and Porter 2009). Bioenergetic models assess these mechanisms integrating available knowledge to better understand and describe the processes behind A. rubens plasticity under varying biotic and abiotic environmental conditions. Dynamic energy budget (DEB) theory provides first-principles models describing the processes of matter or energy uptake and their uses in maintenance, development and reproduction. DEB theory is unique as it integrates the entire life cycle of a species as a function of its physical environment and food availability, providing mechanisms underlying the effects of environmental variability on the species population dynamics (Kooijman 2010). Due to these characteristics, DEB models can be used to predict the performance, abundance, resource requirements and potential distribution of individuals and populations under dynamic environments or under future climate scenarios (Kearney et al. 2010). DEB models have been used as an important tool to describe population dynamics of key species for aquaculture (van der Veer et al. 2006; Saraiva et al. 2012; Maar et al. 2015), fisheries (Pecquerie et al. 2009; Jusup et al. 2011), conservation (van der Meer and Kooijman 2014; Agüera et al. 2015c; Augustine et al. 2017; Marn et al. 2017) and invasion (Sarà et al. 2013; Agüera and Byrne 2018; Gourault et al. 2019). In the present context, DEB models are a powerful tool, which can assist in the management of this pest species.

In the present study, we developed and calibrated a DEB model for A. rubens. This DEB model integrated the available information on the biology and eco-physiology of A. rubens. Although there exists a vast literature on A. rubens, available data were insufficient. To improve the data available we performed a number of eco-physiological experiments exposing A. rubens to controlled prey availability treatments under a natural varying temperature. These experiments focused on the development of the functional response of A. rubens when feeding on its preferred prey, the blue mussel (Mytilus edulis) with the aim to provide a model capable of assessing the prey resources required by this predator and the importance of mussel beds to the biology of A. rubens. These experiments were repeated in different seasons and at different developmental stages allowing us to obtain a comprehensive dataset to estimate the parameters of a DEB model for the whole life cycle of A. rubens. Experimental data were combined with field data from a Danish estuary. An independent experimental dataset was used for the validation of estimated parameters and further exploration of the model was done by performing simulations against field-collected data on growth and reproductive condition.

Our aim was to develop this model as a tool to be used for management of the populations of A. rubens in Northern Europe, where this species is a known pest in shellfish production and also where alternative uses of this species have resulted in the development of a fishery in Denmark (Petersen et al. 2016). This DEB model can be used to estimate food demand, consumed food, growth, reproductive output and other quantities that could be readily integrated into models for the management of A. rubens populations in Limfjorden but also at other locations.

Materials and methods

A series of experiments were performed with the aim of recording the effect of controlled food availability and natural temperature variability on predation, growth, metabolic rate and reproductive output of A. rubens. These experiments provided a complete overview of A. rubens bioenergetics at different seasons and developmental stages. The use of a wide set of observations covering several biological traits and developmental stages aimed to improve the applicability of the yielded DEB parameters. Parameters obtained from more complete datasets are more robust when used under different scenarios and applications (Lika et al. 2011b). Independent datasets from other experiments and from field observations were used to validate model parameters.

Assessing A. rubens performance as a function of food and temperature

Feeding and growth under variable temperature and different food availability

Two growth experiments were performed in autumn 2016 and spring 2017. Experiments were planned with the aim of assessing growth at different food levels while covering the range of temperatures that these animals experience in the field. The experiments were carried out at the Danish Shellfish Centre (DTU Aqua, Nykøbing Mors, Denmark). The setup consisted of containers in an open circuit. Water was pumped directly from the adjacent estuary and passed through the containers in an upwelling direction (inflow at the bottom of the container, outflow from the top). Each container was divided into three cells. A total of 15 containers were used, i.e. 45 cells. Containers were sorted in 5 different food treatments 0, 2, 4, 8 and 16 mussels per cell (equivalent to a mussel density of 0, 1.6, 3.2, 6.4 and 12.8 g AFDW m−2 respectively), all 3 cells in the same container (pseudoreplicates) contained the same food treatment and each food treatment was replicated in 3 different containers (see Supplement 1, Figure S1 for a diagram of experimental set-up). During the experiments, food levels were kept constant by replacing eaten mussels daily. One seastar of an arm length (AL, distance measured from the centre of the disc to the edge of the longest arm) of 3.5 cm was added to each cell. Mussels were selected to provide a constant prey size of approximately 2 cm shell height (distance from the umbo to the opposite shell edge). Both seastars and mussels were hand-picked from long lines on a mussel farm in Sallingsund, Limfjorden (Denmark).

The autumn experiment was run for 63 days (September 27th to December 6th) and the spring experiment for 49 days (May 10th to June 28th). At the start of each experiment seastar set in the cells were measured for arm length and weighed for wet weight. A sample (n = 20) was taken from the same stock and within the same size range to determine initial organ indices and ash-free dry weights (AFDW). During the experiment eaten mussels were counted daily, their shell length measured and replaced by new ones. Every week a sample from the mussel stock (n = 25) was taken to assess their AFDW to shell length relationship. Every week all seastars in the cells were measured for arm length and weighed for wet weight. Water temperature was recorded daily. At the end of each experiment seastars in the cells were dissected to determine their organ indices and AFDW.

The experiments previously described were replicated during 2014. However, some differences apply that justify the replication to improve the detail in certain aspects of the data recorded (e.g. constant food and potential differences due to seasonality). The 2014 experiment used forty containers, i.e. 120 cells. Containers were sorted into three different size categories for both seastars and offered mussels: small (5 cm AL, fed with 2 cm mussels), medium (7 cm AL, fed with 3.5 cm mussels) and large (9 cm AL, fed with 5.5 cm mussels). Mussels were given in 4 different food treatments 0, 1, 5, 10 mussels per cage. Every 5 days eaten mussels were counted and all mussels replaced. Both seastars and mussels were collected on the same mussel farm in Sallingsund, Limfjorden (Denmark). Food levels were not constant but averaged over 5 days considering the number of mussels eaten during that time. In the scale of the 2016–2017 experiment the food treatments approximated to 0, 0.2, 1 and 2 mussel per cage. This experiment considered a larger size range of seastars and over a longer-term, as it lasted 130 days (from September 5th 2014 to January 13th 2015).

Data obtained from these three experiments included: A. rubens growth in arm length and weight at different temperatures and size, and feeding rates at different food levels, temperatures and size.

Oxygen consumption rates at different food levels and constant temperature

An experiment was performed to determine the oxygen consumption rate (OCR) of A. rubens at two different food levels. For this experiment, animals were kept in a temperature-controlled cabinet in 3 L containers at a constant temperature of 18 °C. Each container contained 10 specimens belonging to one of two size classes (1.3 and 2 cm AL). Two containers were used per size class, in total 40 seastars were used. One container of each size class was provided with an excess of mussels (10 mussels per seastar) while the seastars in the other container were starved. Every week close chamber respirometry was performed to measure OCR in 6 starved and on 6 fed animals combining individuals from different size classes. OCR was measured as consumed mg O2 day−1 per individual. After three weeks, OCR of starved individuals reached a stable minimum. A week later OCR was measured on 12 feeding individuals and on 12 starved individuals. OCR was measured between 100 and 70% oxygen saturation using a Presens Oxy-4 oxygen meter. The data from the last set of respiration measurements were used to assess the metabolic rate at size for feeding and starved animals.

Larvae growth at constant food and temperature

There are several sources of information on the larval development of A. rubens, however, most of these studies had a morphological focus and failed to properly describe the larvae rearing procedure or to follow the growth of the larva with detail enough for the data to be used in a DEB model parameterisation. To fill this gap, A. rubens larvae were reared from eggs fertilised in the laboratory. Adults were collected in Limfjorden, their sex determined through visual inspection of the gonads from small incisions. Gonads from three females were dissected out and placed in filtered seawater with 1 µM 1-methyladenine to induce maturation (Strathmann 1987). Sperm was obtained from dissected testes of three males by mechanical extraction and diluted in filtered seawater and used to fertilise the eggs. Fertilised eggs were rinsed to eliminate the excess of sperm after scoring fertilisation and separated in 10, 500 mL glass beakers. These beakers were kept at 15 °C and contained a density of 1 egg or larvae per mL. Samples were observed daily to confirm development stages. When larvae reached the early bipinnaria stage (start of exogenous feeding, i.e. DEB birth) healthy-looking larvae were transferred to five 1L glass beakers and kept at a density of 0.25 larvae per mL. From that moment, larvae were fed ad libitum with a mix of 3 species of microalgae (Tisochrysis lutea, Pavlova gyrans and Rhodomonas salina at a concentration of 7000 cells mL−1 of each). Water was changed every two days and new food was added then. Samples of at least 10 individuals were taken from each beaker at day 0, 2, 4, 6, 8, 10, 12, 14, 19, 26 and 32. The larvae were anaesthetized by exposing them to an isosmotic solution of 10% magnesium chloride (MgCl2) in seawater. Relaxed larvae were put in a welled slide and photographed in a stereomicroscope. The software ImageJ (Schneider et al. 2012) was used to process the pictures and measure the larvae length to the 1 µm. A total of 50 larvae were measured at each sampling. When larvae were already developing the primordia and getting ready for settlement the experiment was terminated.

To approximate egg quality, a sample of unfertilised eggs was taken and processed following Hoegh-Guldberg (1994) procedure to determine their AFDW. A total of 5 samples of approximately 1500 eggs each were taken in centrifuge tubes and let to settle, excess of seawater was then removed, eggs were then washed twice in an isosmotic ammonium formate solution and transferred to tared aluminium crucible. Samples were left to dry at 65 °C for three days, weighed for dry weight and burnt at 550 °C for 4 h to assess ash content.

The DEB model for Asterias rubens

Experimental results were used to parametrise a Dynamic Energy Budget (DEB). The DEB theory (Kooijman 2010) describes the fluxes of energy and mass of the organism throughout life. The DEB model divides the energy and mass between four state variables: reserves (E), structural volume (V), maturity (EH) and reproduction buffer (ER) and describes the fluxes between these state variables. A detailed description of the standard DEB model can be found in the following work (Kooijman 2010; Agüera and Byrne 2018) (Fig. 1).

Fig. 1
figure 1

Schematic representation of the DEB model for Asterias rubens. Left: arrows represent energy fluxes (J day−1) that drive the dynamics of the four state variables (boxes). Energy enters the organism as food (X), which is assimilated at a rate of ṗA into reserves (E). Mobilization rate (ṗC), regulates the energy leaving the reserves to cover somatic maintenance (ṗM), structural growth (ṗG), maturity maintenance (ṗJ), maturation (ṗR) (immature individuals) and reproduction (ṗR) (mature individuals). κ is the proportion of the mobilised energy diverted to ṗM and ṗG, while the rest is used for ṗJ and ṗR

Full size image

Some aspects of A. rubens life history and biology depart from the standard DEB model assumptions. Therefore, in this work the standard DEB model has been extended and modified to accommodate these deviations. In DEB models several parameters (Table 1) are related to structural volume or surface. This structure relates to physical measurable dimensions through conversion coefficients, such as the shape coefficient (δM). The shape coefficient is constant as DEB theory assumes isomorphism during its whole life (Kooijman 2010). A. rubens undergoes a larval stage during which the individual grows exponentially until it settles and metamorphoses into the final adult (Gemmill 1914). This larval stage implies a different animal shape and therefore the model used in this study accounted for two different shape coefficients, one for the post-metamorphic stages and another one for the larvae (Table 1). The exponential growth during the larval stage relates to metabolism acceleration resulting from a change in the size-specific ingestion rate, assimilation and energy mobilization (Kooijman 2014). The DEB model of A. rubens was also modified to incorporate an acceleration factor to account for the changes in the parameters related to the exponential growth (DEB-abj model).

Table 1 DEB parameters values for Asterias rubens
Full size table
Table 2 Zero-variate data used for the estimation of the DEB model parameters. MRE: mean absolute relative error
Full size table

The gonad development of A. rubens and its relationship with the pyloric caeca is similar to that observed and described for A. amurensis (Agüera and Byrne 2018). A. rubens gonad grows exponentially over the winter and spring (Jangoux and Vloebergh 1973; Oudejans and van der Sluis 1979b). For this DEB model, it was assumed that A. rubens energy reserves are stored mainly in the pyloric caeca linking this organ to the DEB state variable E. Moreover, the reproductive buffer was assumed to be also stored temporally in the pyloric caeca. The reproductive buffer stored in the pyloric caeca is mobilized to the gonad during the gonad growth season to support gonad exponential growth. The set of reproductive buffer handling rules developed for A. amurensis by Agüera and Byrne (2018) were applied here for A. rubens. These reproduction buffer handling rules needed the estimation of two extra parameters α and rg (Supplement 1, Figure S2). The link between the energy reserves with the pyloric caeca and reproduction buffer described by Aguera and Byrne (2018) can also be used to propose a metabolic strategy of A. rubens when facing starvation for different periods of time. Asterias rubens can withstand long periods of starvation. During starvation A. rubens individuals lose weight but also shrink in arm length The present model assumed a priority of maintenance over reproduction, and the use of the reproductive buffer stored in the pyloric caeca over the gonad tissue during prolonged starvation.

Parameter estimation: The covariation method

The covariation method was used to estimate A. rubens DEB parameters using the toolbox package DEBtool (available at http://www.bio.vu.nl/thb/deb/deblab/debtool) in a MATLAB® 2018a. The covariation method used the experimental and field observations at different life stages and under contrasting conditions of food and temperature. These observations were approximated using expressions derived from DEB theory and using the DEB parameters for the species (Jusup et al. 2017). Approximation of parameters values is done using a Nelder-Mead numerical optimisation to minimise the difference between observed and predicted values based on a weighed least-squares criterion (Lika et al. 2011a, b). The covariation method has been used for the DEB parameter estimations in hundreds of species proving to be a robust method. The covariation method required starting parameter values to initiate the estimation of the parameters for A. rubens. The parameters previously estimated for the closely related Asterias amurensis (Aguera and Byrne 2018) were used as starting parameter values. The covariation method was also used to determine the parameters describing the reproductive buffer handling rules simultaneously with the other parameters.

The dataset from the feeding experiments was divided. Data on growth, feeding rate and temperature from the shorter experiments performed in 2016 and 2017 were used for parameter estimation, these datasets offered more detailed data on feeding rate at the size and at food density and during two contrasting seasons without large variations in temperature within each of them. The data from the experiment performed in 2014 was used for model validation as it comprised a long-term simulation with a larger size range.

Gonad growth and pyloric cycle

Data on A. rubens gonad growth and pyloric index annual variability in the Limfjorden were not complete enough to model gonad growth and the pyloric cycle. Instead, data from the research conducted by Oudejans and van der Sluis (1979a) and Oudejans et al. (1979) were used to model gonad growth and pyloric cycle and to parametrise the rules for the handling of the reproductive buffer (Agüera and Byrne 2018). These works reported monthly values of gonad and pyloric indices for a population from the Wadden Sea (the Netherlands) during a year. Literature data were complemented with temperature measurements taken during the same time at a nearby location (Texelstroom, http://www.waterbase.nl). These data were used in the covariation method to fit the equations described by Agüera and Byrne (2018) for A. amurensis. A function script was added to the already available DEBtools toolbox to be used for the determination of the parameters α (which determines the initial size of the gonad for the exponential growth) and rg (the exponential intrinsic growth rate of the gonad) determining the transfer of energy from the reproductive buffer stored in the pyloric caeca to the gonads during the gonad growing season (Agüera and Byrne 2018, Supplement 1 Figure S2). The length at puberty (Lp) was determined as the smaller individual with gonads in May 2017. May is within the spawning season and gonads are present (spent or spawning) in all individuals mature enough to produce them.

Parameter validation: experimental and field data

The performance of the estimated parameters was explored on an independent set of data. Model simulations were used to assess the model performance by comparing these simulations with the results from the experiment performed during the year 2014–2015. Further validation consisted of the application of the model parameters to simulate wild populations and comparing the model with empirical data determined gonad and pyloric cycles.

Field sampling for validation of somatic and gonad growth of a field population

Between 2014 and 2015 four sites in Løgstør, Limfjorden, Denmark were sampled every 1 or 2 months using a mussel dredge. All individuals or subsamples depending on the numbers were processed: length distribution, wet weight and gonad and pyloric indices were recorded for these field populations. Organ indices were defined as the ratio of the organ fresh weight by the fresh weight of the individual. Field data on gonad and pyloric index from Limfjorden comprised 5 sampling events between August 2014 and March 2015. The parameters estimated in this study were used to predict the changes in growth, gonad and pyloric indices for these populations between sampling points. There were no data of the available food (mussels or other bivalves), therefore the food level was determined from the animal condition in August 2014, assuming that at this time of the year the pyloric caeca contains only reserves which were in equilibrium with available food at that moment. Seawater temperature was obtained from a logging station in Løgstør, nearby the sampling locations (Danish nationwide monitoring program, NOVANA).

Model applications

The resulting DEB model can be used to estimate food demand (maximum amount of prey that could be consumed at given environmental conditions if food availability was not restricting it), consumed food (prey consumed at given environmental conditions and food availability), growth, reproductive output and reproduction status and other quantities that could be readily integrated into models.. The parameters obtained in this study were used to run a simulation using data from ongoing monitoring reports on A. rubens population biomass and length distribution in Limfjorden. These simulations were used to approximate the impact of this population on the local mussel stock and A. rubens population reproductive output.

Results

Experiments assessing A. rubens performance as a function of food and temperature

The experiment during autumn 2016 exposed A. rubens to a decreasing temperature from 16 to 4 °C with the mean temperature during the experiment being of 8.4 ± 3 °C (all data reported as mean ± standard deviation) (Supplement 1, Figure S3). During this time fed seastars grew from the initial AL of 3.4 ± 0.16 cm and 7 ± 0.72 g wet weight to a final 3.8 ± 0.38 cm AL and 14.4 ± 3.53 g (Fig. 2a, b). Non-fed seastars shrank from 3.2 ± 0.12 cm AL to 2.9 ± 0.25 cm by the end of the experiment and losing about 2.3 ± 0.5 g in wet weight (Fig. 2a, b). Mussels consumed per seastar were stable during the first 40 days and dropped significantly between days 40 and 50 (Fig. 2c). This drop was not captured by the model, suggesting that it was not related to temperature or food. The experiment during spring 2017 showed a different pattern. The temperature increased from 11 °C at the beginning of the experiment to a maximum of 18 °C (Supplement 1, Figure S3). Seastars grew from the initial 3.5 ± 0.13 cm and 8 ± 0.85 g wet weight to 5 ± 0.51 cm AL and 32 ± 8.04 g wet weight by day 49. Non-fed seastars shrank from 3.5 ± 0.11 to 3.2 ± 0.15 cm AL and lost 2.9 ± 0.25 g of wet weight (Fig. 2d, e). Feeding rate increased consistently during the experiment and reached a fivefold increase by the end of the experiment (Fig. 2f). In both experiments growth trajectories at different treatments were not clearly different, the same applies to feeding rate were observed variability was very large (Fig. 2c, f).

Fig. 2
figure 2

Model output and experimental observations. a–c Autumn experiment (2016), d–f Spring experiment (2017). a, d Wet weight at time during the experiment, b, e arm length at time during experiments and c, f daily feeding rate during experiments. Dots and bars are observations mean and standard deviation (n = 9), lines and shaded areas are predictions means and 95% confidence intervals. DEB model predictions were simulated from the experimental individuals starting values (weight and length, reserves and reproductive buffer at experimental day 1)

Full size image

There is an important impact of assimilation in the respiration rate of A. rubens at a temperature of 18 °C. At this temperature respiration rate of A. rubens with an arm length from 1.6 to 2.85 cm ranged between 5 and 12.3 mg O2 day−1 for A. rubens fed ad libitum (f = 1). Starved A. rubens (f = 0) within the same size range exhibited a respiration rate ranging from 1.25 to 4.5 mg O2 day−1 (Fig. 2).

The fertilized eggs of A. rubens hatched 24 h after fertilization and reached the early bipinnaria stage after 3 days at 15 °C. Larvae were fed and measured from day 4 (352 ± 32 μm length, n = 50). Larvae reached settlement competency (late brachiolaria stage, 2.92 ± 0.25 mm length, n = 50) 32 days after fertilization. Mean ash-free dry weight of A. rubens eggs was 1.002 ± 245 μg (n = 5).

DEB model

The covariation method successfully converged on the combination of experimental and literature data over the full life cycle of A. rubens. The parameters fitted the data with a mean relative error (MRE) of 0.180 and symmetric mean squared error (SMSE) of 0.182. These values indicated a small deviation (less than 20%) between observations and model fitted values. These error values are within the range of DEB model parameters fit obtained for thousands of species and reported at https://www.bio.vu.nl/thb/deb/deblab/add_my_pet/

Compared against the DEB parameters for the generalised animal, A. rubens parameters showed a larger maintenance rate [(ṗM), 21.7 vs 18 J day−1 cm−3] and energy conductance (\(\dot{\nu }\), 0.055 vs 0.02 cm day−1). Compared with the sibling species A. amurensis, whose parameters were used as starting values, A. rubens shown a larger maintenance rate [(ṗM), 21.7 vs 15.3 J day−1 cm−3], energy conductance (\(\dot{\nu }\), 0.055 vs 0.033 cm day−1) and the fraction of energy allocated to somatic maintenance and growth (κ, 0.8 vs 0.78) (Table 1).

When evaluating the model output against the datasets used for the parameter estimation procedure. The estimated DEB parameters provided a good fit to the time series on growth and feeding for the experiments performed in autumn 2016 and spring 2017. The model predictions of growth (MRE = 0.092 and 0.076 for wet weight and AL in spring and MRE 0.102 and 0.153 for wet weight and AL in autumn, Fig. 2) were better than predictions of feeding rate (MRE 0.232 in Spring and 0.298 in autumn) which is also a consequence of the observed variation in feeding rate in both experiments. In general, the fit is better for the spring dataset, a consequence of the inability of the model to describe the drop of feeding rate occurring after day 40 during the autumn 2016 experiment (Fig. 2c). The inability to predict this drop-in feeding rate, that was not related to food or temperature, resulted in overestimation in the predictions of AL and wet weight after day 40 (Fig. 2a, b). The model underestimated the respiration rate (Fig. 3) for seastars fed ad libitum (MRE = 0.333). The model fit improved for the non-fed seastars (MRE = 0.201).

Fig. 3
figure 3

Model output and experimental observations of respiration rate at two food levels (‘ad libitum’ and starved). Dots are observations of individual respiration rates. Lines are the model predictions at the given food levels and at a temperature of 15 °C

Full size image

Growth during the larval stage was properly described by the model (Fig. 4, MRE = 0.09) in agreement with the results for the zero-variate data on age and size at birth and metamorphosis (Table 2).

Fig. 4
figure 4

Model output and experimental observations of larvae length. Dots and bars are observations ± standard deviation (n = 50). Line is the model prediction at experimental temperature of 15 °C and ‘ad libitum’ food level

Full size image

The model provided a good prediction of the gonad cycle observed by Oudejans and van der Sluis (1979a, b) with a MRE of 0.179 when simulating the gonad index. However, when predicting the changes in the pyloric index for the same population the fit was worse (MRE = 0.333) failing to fully describe sudden changes in the pyloric index during the year (Fig. 5). Both gonad and pyloric indices observations exhibited large variability (Fig. 5) probably due to differences in individuals’ size, access to food, reproductive state, etc. When simulating for a population with the size distribution observed by Oudejans and van der Sluis (1979a, b) simulations also exhibited a large variability encompasing the error of the field observations (Fig. 5, shaded areas).

Fig. 5
figure 5

Data and model output for gonad and pyloric indices variation during a year. Dots and bars are data and standard deviation from Oudejans and van der Sluis (1979a, b). Lines and shaded areas are mean model predictions and standard deviation for a population of n = 1000, with the same length distribution and subject to the same food level and environment temperature than the population sampled by Oudejans and van der Sluis (1979a, b)

Full size image

Validation results

The experiment performed during 2014 and 2015 exposed A. rubens to a similar range of temperature that the autumn 2016 experiments, decreasing from 17 to 4 °C, with a temperature close to the lowest maintained for 45 days (Supplement 1, Figure S4). The results showed larger differences between feeding treatments with growth and feeding for the treatment with 10 mussels growing faster (Fig. 6). Treatments with 5 and 1 mussel per cage showed a similar output, but with a larger variability related to the combination of food treatment and initial seastar size (Fig. 6). Smaller seastars grew faster increasing their mean weight by 2.5-fold during the experiment and growing 2 cm in AL in the 10 mussel treatment, while the large seastars duplicated their average weight and increased 1 cm their AL in the same food treatment. All starved seastars shrank and lost weight, however, there were differences between size with large seastars shrinking 5% in AL and losing 20% of their weight while the small seastars shrank 8% in AL and lost 15% of their weight (Fig. 6a, b). Observed patterns were therefore comparable to those observed in the dataset used for model parameter estimation. This experiment provided a combination of seastars size, food availability and environmental conditions adequate to validate the obtained parameters.

Fig. 6
figure 6

Model exploration and validation. DEB model predictions for independent experimental data (2014–2015 experiment). Dots and bars are observations mean and standard deviation (n = 26), lines and shaded areas are predictions means and 95% confidence interval. DEB model predictions are done from the experimental individuals starting values (weight and length, reserves and reproductive buffer at experimental day 1). Food treatments are # mussels per cage, however in this experiment mussels were replaced every 5 days, therefore the food scale in this experiment is not comparable to Fig. 2 treatments

Full size image

The parameters estimated for A. rubens provided a good fit when used on this independent experiment dataset. Model simulations were able to describe growth trajectories (MRE = 0.245 for wet weight and MRE = 0.109 for AL) and feeding (MRE = 0.379) for 125 days (Fig. 6).

Field samples from the Limfjorden showed a fast increase in the gonad index between January and March while the changes pyloric index showed an increase from August reaching a maximum in November and changing little until March. The application of the model to a population with the observed size distribution under constant food level resulted in simulations fitting the processes of gonad growth and the reserve dynamics in the field (Fig. 7).

Fig. 7
figure 7

Gonad and pyloric caeca index from field locations in Limfjorden. Dot and bars are observations and standard deviation at three sites in Løgstør, Limfjorden, Denmark. Line and shaded area are DEB model predictions for a population with the same size distribution (n = 400) observed in the field and under the same temperature and food conditions

Full size image

Model applications

The standing stock of A. rubens at March 2018 was of 41,000 t (Nielsen et al. 2018). This population had an average arm length of 7 cm and seastars were within their reproductive season. The application of the DEB model parameters obtained in this study predicted that this population released about 595 trillion eggs out to Limfjorden waters during its 2018 reproductive season, assuming a sex ratio 1:1 and that all the population was reproductively mature. Moreover, at that time of the year this standing stock have a daily food demand of 67 t of mussels (considering wet weight including the shell).

Discussion

Asterias rubens is a voracious consumer of mussels and other bivalves. As such it is considered a pest by mussel farmers and bivalve fisheries all over its distribution range (Hancock 1958; Barbeau et al. 1996; Dolmer 1998; Agüera 2015). This situation prompted a vast series of works that aimed to assess several biological and physiological characteristics of this species, with a focus on its capacity to impact the natural and cultured population of bivalves (Barbeau and Scheibling 1994; Agüera 2015). However, the management of this pest requires tools capable to assess the performance and impact of this keystone species under varying environmental conditions. In this study, we developed and calibrated a Dynamic Energy Budget model, a full life cycle eco-physiological model describing the functional response of A. rubens and the energy allocation to growth and reproduction.

The modelling of Asterias rubens functional response resulted in a very low half saturation parameter XK, which indicates that they fed at their maximum even at a density of 1.6 g AFDW of mussel flesh per m−2. As shown during the validation, mussel availability lower than that (treatments 1 and 5 of 2014–2015 experiment) already limited feeding and therefore hinder growth. However, these feeding rates were already large with a seastar of just 3.5 cm AL consuming a mean of 1 mussel of 25 mm per day (at 16 °C) and a seastar of 9 cm AL consuming almost 2 mussels of 55 mm per day (at 16 °C). A. rubens growth and feeding changed with the season, mainly driven by changes in environmental temperature. As such our experimental animals double their weight in 20 days during spring, while it took them almost 60 during winter. Similarly feeding rates were higher during the spring experiment, resulting from both larger individuals and increased proportional feeding rate. Both growth and feeding rate results are comparable to previous studies on other locations along with the distribution that shows the strong dependency of A. rubens growth and performance related to temperature (Agüera et al. 2012; Agüera 2015). Other factors not considered in these experiments are also known to affect predation rate by A. rubens such as salinity (Agüera et al. 2015b) or other stressors (Agüera et al. 2015a). The decrease of feeding rate during the experiment in autumn 2016 (Fig. 2c) was not explained by the temperature sensitivity captured by the model. However, salinity changed abruptly from 28 to 24 at day 40 (Supplement 1, Figure S4), and fluctuated between 24 and 27 during the rest of the experiment. Sudden changes in the salinity of 4 units are known to decrease feeding in A. rubens temporally (Agüera et al. 2015b). Further development of the model should aim to incorporate the impact of salinity in the bioenergetics of A. rubens.

As is typical for seastars, A. rubens exhibit a low metabolic rate (Propp et al. 1983; Saranchova and Lukanin 1989) with a considerable contribution of the specific dynamic action related to food assimilation (Fig. 2). This low metabolic rates also relate to the resilience to the prolonged starvation characteristic of echinoderms (Monaco et al. 2014; St-Pierre and Gagnon 2015) as they relate to low maintenance costs.

One key aspect of asteroids is their capacity to undertake long periods of starvation, in these periods the animals may shrink consuming both reserves and structure to maintain themselves (Monaco et al. 2014). The present model aimed only to describe short-time starvation (Figs. 3, 6), in this case a simple assumption of re-mobilization of structure (V) when reserves were too small to maintain the animal. This approach fitted well the observed changes in weight during short starvation. A. rubens survived more than 125 days without any prey in our experiments. The described mechanism does not fit the observed changes in arm length which suggests that under starvation arm-length may not be an accurate descriptor of structural length and volume. Indeed previous research point out to changes in the general shape of A. rubens under starvation with changes not only in arm-length but also in arm thickness (Briggs 1983) which implied a change in the shape coefficient (δM) due to fail to preserve isomorphism, a fundamental assumption of DEB theory (Kooijman 2010). Studies on other asteroids have suggested that they can reduce their metabolic demand during starvation reducing the need to remobilize structure (Monaco et al. 2014). The parameter (ṗM), related to maintenance of structure, assumes an average activity levels, i.e. maintaining of the structure under normal use (food search, intra and inter-specific interactions, etc.) (Kooijman 2010), it is possible that after a certain starvation period A. rubens and other asteroids reduce their activity levels reducing their need of energy entering a resting/dormancy stage. This resting stage may relate to the ‘waiting stage’ described by Nauen (1978). Newly settled juveniles may have difficulties to access food due to competition with larger individuals (Nauen 1978) or because present prey is outside of their predation size range (O’Neill et al. 1983). Under these difficulties, juveniles may enter in a prolonged resting state or feed suboptimally without growing until they have better access to prey. We did not have the appropriate data to aim for the parametrization of a (ṗM) under starvation without risking a fit for the wrong reason on the rest of the parameters. Further research on the capabilities of asteroids to undertake long starvations periods is therefore needed.

Model applications

The main aim of this work was to provide a tool that could be used to manage A. rubens populations. The presented DEB model could be used to estimate predation rate, growth, short-term resilience, temperature sensitivity of metabolic processes, reproductive output, larval development, among many other traits related to the bioenergetics of A. rubens. The simulations performed to validate the model already relate to some of these applications and may serve to exemplify them. Furthermore, we performed a simple application of this DEB model to estimate food demand and reproductive output that could be readily used for the management of A. rubens populations in Limfjorden. These values provide an insight into the large impact of A. rubens on the bivalve stock within Limfjorden and why it is necessary to control its population to alleviate the competition of this species with fishing and aquaculture activities. Simulations also pointed towards a need to reduce A. rubens populations before the reproductive season to reduce the number of propagules in the water and potential recruitment. Seastar fishery in the Limfjorden is currently occurring in the period between November and May before reproduction and when biochemical composition is of more interest for its uses as a resource (van der Heide et al. 2018). Seastar fishery landings in the Limfjorden have fluctuated between 190 to 1700 t yea−1 for the past 7 years. This is lower than the intended 10,000 t year−1 anticipated as supply for a newly opened seastar meal factory in Skive, Denmark. In 2018, 2.45% of the standing stock was fished, representing a potential reduction of 15.5 trillion eggs and saving 1.39 t of mussel flesh daily according to the model. Although, seastars in the Limfjorden also the prey of other important commercial species of bivalve such as cockles (Cerastorderma edule), the main concern for the fishermen has been the protection of mussel beds and especially of relayed beds in the on-bottom culture part of the mussel production. With a strategy to reduce the fertile population, it is the strategy of the fishermen to reduce predation pressure by the seastar and at the same time create an income by fishing seastars.