A METHOD TO IDENTIFY BIMODAL WEIGHT–LENGTH RELATIONS: POSSIBLE ONTOGENETIC DIET AND/OR METABOLISM SHIFT EFFECTS IN ANGUILLA ANGUILLA (ACTINOPTERYGII: ANGUILLIFORMES: ANGUILLIDAE)

2018. A method to identify bimodal weight–length relations: Possible ontogenetic diet and/or metabolism shift effects in Anguilla anguilla (Actinopterygii: Anguilliformes: Anguillidae). Acta Ichthyol. Piscat. 48 (2): 163–171. Background. The power function W = a · L b is commonly used to describe the weight–length ( W – L ) relation (WLR) of fish. Smaller/younger specimens may present different WLR from larger/older ones, introducing errors in the derivation of WLR of the total population. This difference appears through a breakpoint in the log–log plot of W – L data and can be justified due to biological factors or due to errors in the sampling procedure. The aim of the study is to propose a bilinear model (LinBiExp) that identifies the specific coordinates of the breakpoint in the log-transformed W – L measurements. Materials and methods. The analysis was performed using 2627 W – L measurements of European eel, Anguilla anguilla (Linnaeus, 1758), from the Comacchio Lagoon (Italy). The bilinearity produced by LinBiExp model was verified through comparison of slopes and intercepts (ANOVA) of the two linear segments and through the 95% intervals of the highest posterior density (HPD) distribution of breakpoint coordinates estimated by bootstrap regression of LinBiExp. Additionally, gut content analysis was performed in order to detect any diet shift in order to justify the existence of the breakpoint. Results. LinBiExp identified the breakpoint coordinates ( L t , W t ) = (28.9 cm, 35.9 g). The ANOVA showed that there was a statistically significant difference between the slopes and between the intercepts of the two linear segments at 99.9% confidence level. The 95% HPD intervals of L t and W t were 28.4–29.4 cm and 34.5–38.0 g, respectively, based on 10 000 bootstrap estimates. The gut content analysis showed inclusion of other fish preys in the diet of eels when their weight and length exceeded the coordinates of the breakpoint in W – L data. Conclusion. The estimated breakpoint for the specific dataset was justified by the possible interrelation of ontogenetic diet shift with other metabolic processes (e.g., beginning of sexual maturation). The study showed that the LinBiExp function can be a valuable tool for detecting the absolute coordinates of a breakpoint in log-transformed W – L data, while the presented methodology can increase the robustness of weight–length analysis of fishes using the typical power function.


INTRODUCTION
The knowledge of fish weight-length (W-L) relation (WLR) is one of the most important elements for the description of body proportions of fish species. The typical power function (W = a · L b ) is the most widely used function for the description of WLRs, facilitating weight predictions from length measurements (Froese 2006). Furthermore, an important contribution of the typical power function is that its coefficients (a, b) for various fish species are stored in FishBase (Froese and Pauly 2018). These coefficients have also been used to analyse growth allometry and growth condition factors, which show the degree of well-being of a species in its habitat (Knights 1982, Bolger and Connolly 1989, Lima-Junior et al. 2002, Gomiero and de Souza Braga 2005, Acarli et al. 2014. The most extensive review about different types and methods for obtaining WLRs was provided by Froese (2006) while additional guidelines were proposed in a more recent editorial note by Froese et al. (2011). Froese (2006) noted that the derivation of WLRs for fish species may seem a simple technical procedure but there are many issues associated with their derivation, which need special attention. One of the most important issues was that small specimens from very young age classes might present different WLR from the older ones introducing errors in the final WLR of the total population. For this reason, a recommendation was the exclusion of the young specimens for improving fitting performance. The different WLR of the young specimens can be explained by justifications such as (Stergiou and Fourtouni 1991, Froese 2006, Simonović et al. 2011, Milardi et al. 2014): • The young specimens present different WLR from the larger/older ones due to ontogenetic shifts in diet while this difference can also be regulated by the availability of the different food types required at different growth stages. • Usually the amount of data from very young age classes is small and not representative due to the efficacy reduction of fishing techniques for capturing smaller specimens (gear/net selectivity. • Shifts in growth performance due to metabolic changes associated to sexual maturation. Taking into account the above, the safer procedures to obtain the WLR are either to use the typical power function after the removal of young specimens or to use a more complex model to describe the complete data. The method for removing the data of younger specimens is either empirical (e.g., from visual inspection) or it is based on statistical outliers identification procedures. If indeed, the young specimens have a different WLR, outliers' identification methods do not usually remove all the observations of young specimens, which present different WLR response, but may indicate as outliers, specimens of the larger classes. Thus, a robust method to identify a specific threshold of length or weight in order to remove the "problematic" young specimens is still missing. Such a method would allow the users to derive a more robust WLR for the older/larger specimens using the typical power function.
The aim of this study was: • To propose a bilinear model that can describe the complete dataset of weight-length measurements identifying the specific coordinates of any existing breakpoint in the log-log plot of weight-length data. • To use the specific breakpoint for dividing the data in parts, which would improve the fitting efficacy of the typical power function. The proposed methodology is applied using as an example a large dataset of weight-length data of the specimens of the European eel, Anguilla anguilla (Linnaeus, 1758), captured in the Comacchio Lagoon (Italy). Additional gut content data from the specimens were used to evaluate the hypothesis that the observed breakpoint in the log-log plot of weight-length data is related to an ontogenetic diet shift while additional justifications were given in order to justify the interrelation of diet shift with other shifts in metabolic processes.
Silver eels were sampled at fishing screens called "lavorieri" during the period of seaward migration (November-December, while silver eel migration has not been observed in the period January-October). The screens permit the entry or escape of elvers in the lagoon but entrap all silver eels when they begin their migration, thus, the eel catch in the screens represents ~100% of the migrating silver eel population of the Comacchio Lagoon. The total yield of silver eels of 2011 (for the fishing area of 8470 ha) was 3811.5 kg with a total silver eel abundance equal to 0.45 kg · ha -1 (Aschonitis et al. 2017a). The value of silver eel abundance of 2011 was among the 1% lower annual records starting from 1781. The last thirty years the silver eel abundance in Comacchio showed a continuous decline, and after 2010, has been reduced more than 96% compared to the mean value (~14.5 kg · ha -1 ) of 1781-1980(Aschonitis et al. 2017b. Yellow eels were caught by a set of 20 trap nets evenly distributed over the entire lagoon area. The fishing gear was a modified fyke net, locally called "cogollo", which is typically used for eel fishing in the shallow lagoons of the northern Adriatic. It consists of a leader, 50 × 1.5 m, that directs the fish toward two conical trap nets positioned at its distal ends. The structure consists of 8 × 8 mm mesh, large enough to prevent blocking by macroalgae, periphyton and detritus, but small enough to prevent loss of small age class specimens. The nets were monitored every two days in September and October, a period when metamorphosis to silver eels is considered complete, preventing overlap between yellow and silver eel population groups (van Ginneken et al. 2007).
Yellow and silver eels were counted and anaesthetized with ice in order to measure length and weight. Subsamples of 573 yellow and 366 silver eels (939 in total) were randomly selected for age and sex determination (Castaldelli et al. 2014, Aschonitis et al. 2017a) while all the rest eels were released. The age of specimens was determined by double reading after grinding and polishing the otoliths (Anonymous 2009). Sex was determined by macroscopic examination of the gonads and when specimens were smaller than 35 cm, then microscopic examination of the gonads was performed (Colombo andGrandi 1996, Tesch 2007). The females represented 97.5% of the collected specimens, the males represented 1.4%, and the share of sexually undifferentiated specimens was 1.1% (Aschonitis et al. 2017a). The extremely high percentage of females is attributed to the extremely low abundance of the population in the lagoon (Aschonitis et al. 2017a(Aschonitis et al. , 2017b (it is already known the strong negative correlation of European eel feminization rate with the abundance of the population) (Roncarati et al. 1997, Krueger and Oliveira 1999, Han and Tzeng 2006, Aschonitis et al. 2017a, 2017b. Male and female eels present different WLR and so the error, which might be introduced by mixing males and females in WLR analysis, is negligible in this dataset. Thus, the analysis that will follow, mainly describes the characteristics of a mixed yellow and silver female European eel population. The sub-dataset of mixed yellow and silver female individuals (925 in total) of known age and length have already been used by Castaldelli et al. (2014) for the derivation of the length-age relation using von Bertalanffy function: where L is length [cm], X is age [years], with L ∞ = 155.94, L 0 = 7.5, k = 0.087 and R 2 = 0.97. The total observed age classes of specimens subjected to age analysis using otoliths were 11 (from 0+ to 10+ years, the + accounts for ~0.5). The small number of sexually undifferentiated European eels was restricted to the 0+ and 1+ age classes (Aschonitis et al. 2017a).
Since ontogenetic diet shifts have been documented as one of the reasons for the existence of such breakpoints Fourtouni 1991, Froese 2006), additional observations about the gut content of eels were also obtained. Observations of the gut content were performed based on 62 specimens selected from the initial dataset, with length ranging between 11 and 83 cm. Identification of prey type was performed using a stereo-microscope. The initial number of analysed specimens was much larger, but their majority was excluded because their gut content did not allow robust identification of prey due to almost complete digestion. Prey was identified to the species level and gut content of the selected specimens was expressed as a percentage of biomass per food/prey type. Weight-length relations. The analysis of WLR was performed using the typical power function and the linearized biexponential function (LinBiExp) proposed by Buchwald (2007).
The typical power function and its respective logtransformed form are the following: where W is the weight [g], L is the length [cm] and a and b are regression coefficients. The common procedure for analysing WLRs using Eq. 2 requires the initial use of simple regression for Eq. 2b using the log-transformed W and L variables for identification and removal of outliers (Froese 2006). Then, the fitting procedure is repeated without outliers using either non-linear regression with Eq. 2a or linear regression with Eq. 2b. The LinBiExp provides smooth and fully parametrizable transitions between two linear segments maintaining a clear connection between them. The LinBiExp is fitted on the log-transformed W, L variables and is given by the following function (Buchwald 2007): where W′ and L′ are the log-transformed weight [log(g)] and length [log(cm)], a 1 and a 2 are coefficients that regulate the slopes of the two linear segments, c is a parameter for adjusting the smoothness/abruptness of the transition and the form of angle between the two linear segments, d is a constant for shifting the curve along the vertical axis (log-W axis), and L t ′ is a constant that defines the break point between the two linear segments at horizontal axis (log-L axis). The transition between the two linear segments does not require a sharp break-point, it can take place along a smooth, continuously differentiable, curved portion of adjustable width (with the deviation from linearity having an exponential character). Nevertheless, a model should be considered bilinear only if it shows linearity at both ends of its considered range (Buchwald 2007). Positive values of c coefficient indicate that the angle above the two linear segments is < 180 o while negative c values indicate that the respective angle is > 180 o . The larger the absolute value of c is, the smoother is the transition between the two linear segments. The c coefficient is generally sensitive and a large variation of its values only in a positive or a negative range just affects the steepness in the transition between the two linear segments. For this reason, a high statistical significance for the case of c coefficient is not required. The angle is regulated by the a 1 ÷ a 2 ratio and when it is equal or close to 1 indicates that bilinearity hardly exists. Thus, the smaller the ratio is, the higher is the degree of bilinearity. It has also to be mentioned that the only difference between the two exponent factors inside Eq. 3 is the value of a 1 and a 2 coefficient and thus exchanging their values does not change the curve (e.g., for a 1 = 1 and a 2 = 4 or a 1 = 4 and a 2 = 1, the curve is the same).
The coordinates of the breakpoint in the log(W)log(L) plot are provided by the following functions: for log-L axis: x΄ = L ΄ t (4a) and for log-W axis: The respective coordinates of the breakpoint in regular weight-length curves after removing the logarithmic transformation are provided by: Steps of analysis. The first step of the analysis included simple non-linear regression (Simple-NLR) of LinBiExp (Eq. 3) using the full dataset for detecting outliers. The fitting analysis together with outliers detection was performed using StatGraphics Centurion XV software (Statpoint Technologies, Warrenton, VA, USA). Those observations with studentized residuals greater than 2 in absolute value were considered outliers (standard procedure of the aforementioned software). The new dataset without outliers was used again to perform Simple-NLR of Eq. 3 for the final calculation of its coefficients and the respective coordinates of the breakpoint.
The second step of the analysis was performed for assessing the robustness of bilinear response. This step included two procedures that complement each other.
In the first procedure, the analysis aimed to investigate the robustness of a 1 ÷ a 2 ratio and the robustness of breakpoint coordinates (L t , W t ) provided by LinBiExp (Eq. 3) based on their observed variation when Eq. 3 is subjected to bootstrap non-linear regression (Boot-NLR). Boot-NLR is based on the generation of a large number of new datasets by randomly sampling data with replacement (Efron and Tibshirani 1994) and it is considered among the most robust methods for assessing the variability of regression coefficients. The Boot-NLR was performed by applying the "nls.lm" function of the {minpack.lm} package (Elzhov et al. 2016) in R software. The "nls.lm" function uses the Levenberg-Marquardt non-linear leastsquares algorithm. The Boot-NLR procedure was applied for 10 000 iterations that led to a respective number of (a 1 , a 2 , c, d, L t ′) solutions. The 10 000 bootstrap sets of coefficients were then used to assess the respective values of a 1 ÷ a 2 ratio and the coordinates of the breakpoint (L t , W t ) based on Eq. 4 and Eq. 5 for each bootstrap set of variables. The range of bootstrap estimations of the regression coefficients but also of a 1 ÷ a 2 , L t , and W t , was defined by the 95% confidence interval, which was estimated based on the probability distribution of their 10 000 estimations. This method was applied in order to estimate the values of the highest posterior density (HPD) distribution that indicates the 2.5% and 97.5% thresholds (HPD thresholds), which contain the central 95% of the HPDD distribution. The probability interval was computed using the "p.interval" of {LaplacesDemon} package (Bernardo 2005) in R software.
In the second procedure, the analysis aimed to investigate if the slopes (b) and intercepts log(a) of Eq. 2b are statistically different between the two linear segments defined by the breakpoint of Eq. 3. For this reason, the clean dataset of observations, after the removal of outliers from step 1, was divided in two subsets using the breakpoint coordinates that were also obtained from step 1. Thus, any observed pair of (W, L) with W < W t and L < L t was included in subset 1 (smaller specimens) and the rest were included in subset 2 (larger specimens). Comparison of the linear regressions based on the log-transformed variables was performed for both subsets using Eq. 2b, while ANOVA was used to compare the statistical difference of their slopes (b) and intercepts log(a).
The third step included the examination of possible hypotheses related to the existence of breakpoint in the log-WLRs. The very small size of nets used for capturing eels (8 mm) was chosen in order to minimize the possible errors of nets selectivity for young specimens. Thus, any possible occurrence of a breakpoint in the log-WLR of this dataset could be attributed to other reasons (e.g., environmental or biological factors). Three possible hypotheses of biological background for the existence of the breakpoint in the specific dataset were examined: • A shift of metabolic activity to achieve higher elongation rates at the young age classes.
• A shift of metabolic activity to achieve sexual maturation. • An ontogenetic diet shift.

RESULTS
Step 1: Simple-NLR of LinBiExp for removing outliers and determination of the breakpoint. The minimum and maximum observed weights in the dataset of 2627 specimens were 8.4 and 2371.5 g, respectively, and the minimum and maximum observed lengths were 11.0 and 105.0 cm, respectively. Simple-NLR of Eq. 3 was performed based on the full dataset (n = 2627) of log-transformed variables for detecting outliers (Fig. 1A). This procedure identified 143 outliers (Fig. 1A), which were removed from the dataset. Simple NLR of Eq. 3 was applied again based on the new clean dataset (n = 2484) and the results are given in Fig. 1B. The log-transformed coordinates of the breakpoint in the log(L) vs. log(W) plot (Fig. 1b) were equal to (L t ′, W t ′) = (1.461, 1.555) (Eqs. 4a, 4b), while their values without the log-transformation were equal to (L t , W t ) = (28.9 cm, 35.9 g) (Eqs. 5a, 5b).
Step 2: Robustness of bilinear response. The general statistics together with the 2.5% and 97.5% HPD thresholds for a 1 , a 2 , c, d, L t ′ but also for a 1 ÷ a 2 , L t , and W t using the 10 000 bootstrap estimates of Boot-NLR of Eq. 3 are given in Table 1. For a 1 ÷ a 2 , L t , and W t are also given the HPDD graphs (Fig. 2). Taking into account the results of Table 1, it is observed that almost all the parameters present a restricted range of variation according to HPD thresholds, which suggests a robust performance of the function. The only parameter that presents large variance is the c coefficient, which is very sensitive because it regulates the smoothness/abruptness of the transition between the two linear segments. This is not a serious problem since its bootstrap estimates fluctuate only in a positive range, while the parameters of a 1 ÷ a 2 , L t , and especially W t (which includes c according to Eq. 4b) show small variance, which is also verified by Fig. 2.
The clean dataset (n = 2484) was further divided into two subsets based on the breakpoint of Eq. 3 obtained from step 1. Eq. 2b was then fitted on the two subsets (Fig. 3) and ANOVA was used to compare their slopes (b) and intercepts log(a). The ANOVA showed that there was a statistically significant difference between the slopes (b) but also between the intercepts log(a) of the regression lines of the two subsets at 99.9% confidence level, respectively, verifying the robustness of the bilinear response. The individual regressions of the two lines showed R 2 = 0.58 (P < 0.001) and R 2 = 0.98 (P < 0.001) for subset 1 and 2, respectively.
Step 3: Possible hypotheses that justify the breakpoint of log-WLR. Regarding the first hypothesis about the regulation of metabolic activity to achieve higher elongation rates at the young age classes, it was found that when L is equal to L t ≈ 28.9 cm, eels have already reached approximately the 28% of their maximum observed length in comparison to the weight, which is still less than 2% of their maximum observed value. This performance may be related to a regulation of metabolic activity to achieve higher elongation rates at the young age classes, which presents a shift after reaching L t . The elongation rates of subset 1 and 2 can also be considered statistically different since the slopes b of Eq. 2b (Fig. 3) were also statistically different.
The fact that the elongation rate of subset 2 was smaller than the one of subset 1 may also be associated to a metabolic shift for achieving sexual maturation. This hypothesis can be justified based on the length-age relation (Eq. 1). According to Eq. 1, the range between the 2.5% and 97.5% HPD values of L t (28.4-29.4 cm, Table 1) corresponds to eels age of 1.7-1.9 years. The fact that the observed sexually undifferentiated eels of this population were all < 2 years old (Aschonitis et al. 2017a) Table 1 General statistics together with 2.5% and 97.5% HPD intervals for a 1 , a 2 , c, d, L t ′, a 1 ÷ a 2 , L t and W t using the 10 000 bootstrap estimates obtained from Boot-NLR of Eq. 3   Fig. 2. HPDD graphs of the a 1 ÷ a 2 , L t , and W t parameters (the central black portion of the graph describes the 95% central part of the distributions) for the European eel, Anguilla anguilla, of the Comacchio Lagoon that the high metabolic performance for achieving high elongation rates is reduced after L t on behalf of metabolic processes associated to sexual maturation. The third hypothesis of ontogenetic diet shift was investigated based on the gut content data. The gut content analysis is presented using five length classes (Fig. 4). The threshold value between the second and the third length class is equal to L t = 28.9 cm. The gut observations for the specimens with L < L t , showed that eels diet primarily consisted of small amphipods (Gammarus spp.) and secondly by small shrimps (Palaemon spp.). On the other hand, the analysis of specimens with L > L t , showed a diet shift by showing higher shrimps content in comparison to amphipods and inclusion of small fish species like anchovies, Engraulis encrasicolus (Linnaeus, 1758). The percentage of fish species in their diet was gradually increased with the length increase, including also other species like grass goby, Zosterisessor ophiocephalus (Pallas, 1814), or larger shrimps and crabs (Carcinus aestuarii). The existence of brackish/saltwater fish species like anchovy and grass goby in their diet is due to the saltwater conditions of the Comacchio Lagoon, which is connected to the Adriatic Sea.

DISCUSSION
The results of this study could not exclude any of the initial three hypotheses as causes of the breakpoint existence. On the other hand, it is likely that an interrelation and synergy among them exists as a result of both biological (metabolic processes) and environmental (prey type and availability) factors. The combined explanation about the specific breakpoint occurrence could be that the specimens are trying to optimize the exploitation of energy, which is derived by various prey types of different size, in this specific environment. At the initial stages of growth, individuals invest energy on somatic growth in order to perform faster the following two strategies: • Their protective mechanisms (swimming performance).
A justification for this behaviour is that it probably increases the defensive ability of the species from predators since the longer elvers swim faster and since they are still thin (small trunk section), they can easily hide in small holes of the bottom or inside the stones. • Their ability to capture larger preys of higher energetic content that will allow them to cover the higher energy demand of metabolic processes associated with sexual maturation. When their size reaches a specific threshold (breakpoint), which allows them to consume enough energy for boosting sexual maturation, then their metabolic processes change towards this latter purpose.
Fish consumption (especially anchovy), which started to appear in the gut content analysis when L became larger than L t (Fig. 4), probably plays a key role in European eel maturation. Anchovy belongs to the oily fish group, providing high amounts of oil and fatty acids (Üstün et al. 1996) that have extremely high nutritional value for eels, which also belong to the same fish group (Pike and Jackson 2010). The fact that young European eels invest more energy to increase their length in comparison to their transect diameter probably leads to a faster elongation of their stomach, which also has an elongated shape (see fig. 1 in Peters 1982). This may also help them to consume and digest more efficiently anchovies, which have an elongated body shape in spite of their small size (Whitehead et al. 1985). The process of elongation is, of course, a genetic attribute of eel species associated with their particular snakelike morphology but it could be hypothesized that this rate might be enhanced at the initial growth stages in order for a specific population to adapt itself to better exploiting the fish preys present in the environment. If this hypothesis is true, then it would indicate high plasticity of elongation rates associated to specific food types. European eels have been shown to be capable to adjust their morphology for similar reasons. For example, in a study by De Meyer et al. (2016), it was found that a controlled population of eels presented phenotypic difference in their head when they were fed either hard or soft diets. The authors found that hard feeders developed a broader head and a larger abductor mandibulae region that made them capable of stronger bites. On the other hand, soft feeders developed a sharper and narrower head, which could reduce hydrodynamic drag, allowing more rapid strikes towards their prey.
Furthermore, it should be emphasized that this is the first time that a breakpoint in the log-WLR of European eel is documented. This can be attributed either to the fact that other eel datasets were not adequately expanded in the small age classes for identifying the existence of a breakpoint or to the fact that the existence of the breakpoint in this dataset is a rare case based on the local preys' preferences of eels (Fig. 4). It has to be noted that dietary shifts of European eels with the gradual inclusion of fish preys may not be observed in other areas due to differences in the availability of different prey types and the high competition for fish preys with other predators. For example, Dörner et al. (2009) analysed the dietary habits of European eel in two lakes-one in Germany (Großer Vätersee) and the other one in Denmark (Vallum Sø). Despite the fact that Vallum Sø had significantly higher amount of small prey fishes and macrozoobenthic preys compared to Großer Vätersee, the eels of Vallum Sø were mainly fed with macrozoobenthic preys (e.g., chironomid larvae), while the eels of Großer Vätersee used fishes as the main food component. The above-cited authors concluded that the density of macrozoobenthos generally controlled the degree of piscivory. Another element that should also be considered is the possible seasonal variability in eel diet. The samplings of this study were performed during September-December, and thus seasonal effects could not be identified. On the other hand, Bouchereau et al. (2009) found seasonal variations in the diet of eels sampled in Languedocian Mauguio Lagoon (Gulf of Lion, France) but they did not consider the diet differences among different eel size classes.
Taking into account the aforementioned findings and comments, it is observed that there is a large knowledge gap regarding: • The feeding preferences of European eel which can vary between different size classes, different seasons, or different habitats. • The effects of dietary habits on its growth and morphometric characteristics in natural environments. Future research should focus on the aforementioned issues giving also special attention to the sampling procedures for improving the capture of smaller specimens in order to verify the breakpoint existence in the log-WLR of other eel populations. If other breakpoints can be found in other extensive eel datasets, excluding the possibility that they are related to the selectivity of fishing technique, then the breakpoints could be used as a biological and environmental indicator. Furthermore, the aforementioned analysis should be expanded to other fish species, since the breakpoint point was also observed in log-transformed weight-length data of other species (Simonović et al. 2000(Simonović et al. , 2011. Additionally, it would be of special interest the investigation of breakpoint existence in log-transformed WLRs of fish species, which are considered to have isometric or slight allometric growth. Finally, regardless of the reasons of the breakpoint existence, its observation and identification through LinBiExp (Eq. 3) are very important since it can be used to improve the WLR when it is described by Eq. 2. The coefficients of Eq. 2 are among the most important data of FishBase (Froese and Pauly 2018) and it is unknown if any of the thousand coefficients of many fish species were biased due to the existence of a breakpoint in the initial data.
of the Fisheries Bureau of the Emilia-Romagna Region for the support to the research in the context of long-term collaboration.