Measurement of the Growth of the Main Commercial Rays (Raja clavata, Raja brachyura, Torpedo marmorata, Dipturus oxyrinchus) in European Waters Using Intercalibration Methods

Simple Summary Poor quality of biological information, such as age and growth parameters, could be a source of variability with a significant impact on stock assessment results. Concerning the ageing process, variability is frequently linked to differences in the interpretation of calcified structures. The evaluation of precision and accuracy therefore represents a keystone in the ageing procedure. Abstract The intercalibration of age readings represents a crucial step in the ageing procedure; the use of different sampling methods, structures, preparation techniques, and ageing criteria can significantly affect age and growth data. This study evaluated the precision and accuracy of ageing for the most important North Atlantic (NA) and Mediterranean (M) ray species, Raja clavata, Raja brachyura, Torpedo marmorata, and Dipturus oxyrinchus, through exchange exercises carried out by readers from different laboratories. In addition, growth parameters were estimated from the obtained data. A total of 663 individual batoids were analysed. R. clavata and R. brachyura samples were obtained from both the NA and the M, while vertebral centra of T. marmorata and D. oxyrinchus were only available for the M. High reading variability was observed for all four evaluated species in terms of CV, APE, and PA. D. oxyrinchus and T. marmorata showed relatively slow growth and the von Bertalanffy model with fixed t0 and Gompertz’s model were, respectively, the most precise models for each of these species. In R. brachyura, females had a faster growth rate compared to combined sexes. The vbt0p proved the most precise model for describing growth in this species, and no statistical differences were found between the NO and the M. For R. clavata, the best-fitting model was the vbt0p for females and males in the NO and for females from the M, while the best-fitting model for males from the M and sexes combined for both areas was log.p. Distinct growth patterns were observed between the two study areas.


Introduction
Elasmobranchs are globally indicated as one of the most highly threatened vertebrate groups [1], due to their life history traits, which are typical of k-selected life strategy species [2].These traits make these fish highly susceptible to anthropogenic impacts and in particular to fishing pressure [1,3].Several shark and skate species are commercially valuable for their fins, meat, liver oil, gill rakers, or leather and are an important food resource [4].Sharks and rays were once considered a lower-value bycatch of more profitable fisheries stocks, such as tuna, cod, and shrimps [5].The rising demand for marine products, coupled with the decline of valuable target stocks has, however, resulted in rising catches and retention of these taxa [6].Sharks and skates are found as bycatch species in fisheries worldwide, and the Mediterranean and Atlantic Ocean are areas with a high level of bycatch for both sharks and rays at a global scale [7].Fishing pressure thus puts cartilaginous fish at a higher risk, which emphasises the need for effective conservation and management measures.Shark and ray fisheries have only recently been subject to management strategies, and attention to their need for conservation has grown [8].The International Union for the Conservation of Nature's (IUCN) Red List of Threatened Species estimates that about a quarter of elasmobranch species are threatened with extinction (i.e., assessed or estimated to be Vulnerable, Endangered, or Critically Endangered), and overfishing is the principal threat behind elasmobranch population declines [1].Despite the advances in shark and ray fisheries management, there are concerns that elasmobranch fisheries are following the pattern of unregulated fisheries, resulting in wild population decline and collapse of stocks [9].Indeed, shark and ray landings increased from 1950 (the first year of data collection) to the peak year at the beginning of the 2000 s, then subsequently declined in the following years [10].
Conservation and fisheries management measures require a solid assessment of population status, which must be based on reliable information regarding species life history traits.Age and growth data are essential to obtain mortality data and productivity estimations, which are crucial in stock assessments [11].Describing growth parameters has been historically more difficult in cartilaginous fish than bony fish [12].Ageing in bony fish involves the use of calcified structures, such as otoliths, which are absent in elasmobranch species.Cartilaginous structures in elasmobranchs, such as vertebrae, have low calcification levels and often require additional enhancements of growth band visibility with the use of staining techniques.These characteristics make ageing of elasmobranchs a rather complex process [13].
Inaccurate population assessments resulting in stock collapse have been caused in some cases by poor-quality ageing data [14].A growing body of research has therefore been conducted to increase the accuracy of age data, particularly within the European Union Data Collection Framework [15].Variations in age data between institutes may occur due to sampling methods (e.g., commercial fishing or scientific surveys) [16], the use of different structures (vertebral centra, spines, scales, etc.), preparation techniques [17], and the ageing criteria used [11,18].Variable levels of fishing pressure [19,20] and spatial variations related to environmental conditions or genetic factors [21,22] could be the basis of variation in growth patterns in fish stocks, including contiguous populations.The fitting of growth models to age data could also be a source of variability with significant impact on stock assessment results.Finally, the level of reader experience can be a very important additional source of variability [23].All these factors can compromise both precision and accuracy of age data and consequently the analysis of the level of population exploitation.Unreliable scientific advice can result from using incorrect growth parameters or age-at-length keys to translate size distribution into age structure.If the age of a population is overestimated, the stock assessment will incorrectly predict that fishing mortality will be lower because the population will be composed of older specimens.Conversely, if the age is underestimated, fishing mortality will be overestimated, and the population will appear to be composed of younger specimens [24].Additionally, age and growth have an impact on how data on natural mortality and maturity at age are estimated.As a result, these measures also impact how recruitment strength and biomass of spawning stock are estimated.Ultimately, the most significant impact is related to short-term stock status forecasts and the corresponding management actions [18].
In this context, an age calibration exercise was initiated on the most important North Atlantic and Mediterranean ray species from commercial (Raja clavata; Raja brachyura) [10,25] and conservation (Torpedo marmorata; Dipturus oxyrinchus) [26,27] perspectives.The exercise aimed to assess the precision of age and growth data for these species.The present study, resulting from the exercise, aims to provide reliable age and growth data for R. clavata, R. brachyura, T. marmorata and D. oxyrinchus in European waters (North Atlantic Ocean and Mediterranean Sea).These data will then contribute to more effective management plans.As some recent papers stated that alternative growth models could offer a better fit for cartilaginous fish age-at-length data, in particular for batoids (e.g., [28,29]), we applied a number of alternative growth models to age-at-length data in addition to the frequently used von Bertalanffy function.

Sampling
Individuals were sampled in two areas in the North Atlantic Ocean (North Sea, ICES area 27.4; Eastern Channel, ICES area 27.7.d)and in two areas in the Mediterranean Sea (Ligurian and North Tyrrhenian Seas, GSA 9; western part of Sardinia, GSA 11.1) between 2010 and 2020 from commercial sampling and during scientific surveys.All individuals were taken to the laboratory for accurate measurements.Each individual was measured to the nearest mm for total length (TL) and to the nearest g for total weight (WT).Finally, the sex of each specimen was recorded.

Ageing Procedures
Vertebrae were excised from the spine during dissection and subsequently used to estimate age data.The preparation methods used varied between institutions, where vertebrae sampled from the North Atlantic Ocean were stained using varying concentrations of crystal violet and read as whole structures, while vertebral samples from the Mediterranean Sea were sagittally sectioned and left unstained (Figure 1).A number of vertebrae were collected for each individual.One vertebra per individual was photographed using a binocular microscope under transmitted light.
Alternating translucent and opaque bands were visible in the vertebra of these elasmobranchs.It was assumed that each annual growth ring consisted of one opaque and one translucent band, as is standard in temperate fish sclerochronological analyses.The age was therefore expressed in consistent age groups, e.g., a fish in age group 0 lived between 1 day and 364 days (i.e., between hatching and before 1 year), as recommended by international expert groups [30][31][32].

Ageing Data Precision
To limit interpretation error and reading bias, each individual was analysed by eleven European readers from Italy, Greece, Belgium, the Netherlands, and France during the European exchange in 2022 to evaluate precision.Precision is defined as the reproducibility of repeated measurements on a given scale, whether or not measurements are accurate [30].Precision was measured from the average percent error (APE), the percentage agreement (PA), and the coefficient of variation (CV).The formula presented by Beamish and Fournier [33] was used to calculate APE: where Xij is the ith age determination of the jth fish, Xj is the average age calculated for the jth fish, and R is the number of times each fish was aged.CV and PA within one year (+/−1 yr) were proposed by [24]: where R is the number of times each fish is aged, Xij the i(th) age determination of the j(th) fish, Xj is the mean age calculated for the j(th) fish, and ndiff is the difference in age determination between the readings of two readers.

Growth Model Estimation
Non-linear growth models were fitted to length-at-age data.Mean body growth patterns of the commercial ray species sampled were described using four different growth models including the following: the unconstrained von Bertalanffy model [34] (vbp):

Ageing Data Precision
To limit interpretation error and reading bias, each individual was analysed by eleven European readers from Italy, Greece, Belgium, the Netherlands, and France during the European exchange in 2022 to evaluate precision.Precision is defined as the reproducibility of repeated measurements on a given scale, whether or not measurements are accurate [30].Precision was measured from the average percent error (APE), the percentage agreement (PA), and the coefficient of variation (CV).The formula presented by Beamish and Fournier [33] was used to calculate APE: where Xij is the ith age determination of the jth fish, Xj is the average age calculated for the jth fish, and R is the number of times each fish was aged.CV and PA within one year (+/−1 yr) were proposed by [24]: where R is the number of times each fish is aged, Xij the i(th) age determination of the j(th) fish, X j is the mean age calculated for the j(th) fish, and n diff is the difference in age determination between the readings of two readers.

Growth Model Estimation
Non-linear growth models were fitted to length-at-age data.Mean body growth patterns of the commercial ray species sampled were described using four different growth models including the following: the unconstrained von Bertalanffy model [34] (vbp): the von Bertalanffy model with forced t 0 = 0 (vbt0p): the Gompertz model [35] (vbL1p): the logistic model [36] (log.p): where TL 1 , TL t , andTL ∞ are, respectively, the length at age 1, at age t and the asymptotic length, and K is the rate at which the asymptote is reached, also called the growth coefficient.

Data Analysis
For each individual, the total length and the age group were estimated according to the sex and/or geographical sampling area.With these all-individual data, all growth models were tested and the best growth model was identified as the one that minimised the small sample bias-corrected form of the Akaike Information Criterion (AICc) [37,38].The AICc balances the trade off between the quality of fit and the number of parameters used [39] while accounting for small sample bias, and is defined as follows: where n is the sample size, k is the total number of parameters of the model, and TL is its likelihood.Fish growth was estimated using the growth performance index (ϕ) [40]: Growth performance index was more appropriate for growth comparison versus comparison of TL ∞ and K individually, as these two parameters are highly correlated [41].
The lifespan (t max ) was estimated from the empirical relationship with growth rate k [42] as follows: [43] was used to look for potential differences in growth between areas.

Sample Composition
A total of 663 individual batoids were analysed in the present study.Table 1 reports the specific data of the sample composition.R. clavata and R. brachyura samples were obtained both from the Atlantic Ocean and the Mediterranean Sea, while vertebral centra of T. marmorata and D. oxyrinchus were only available for the Mediterranean basin.R. clavata was the most sampled species, with 224 females (131-955 mm TL) and 204 males (209-900 mm TL) (Table 1), followed by R. brachyura for which a total of 115 samples were collected (54 females 175-990 mm TL; 61 males 70-955 mm TL).A total of 60 specimens of T. marmorata and 61 specimens of D. oxyrinchus were analysed from the Mediterranean basin.TLs of these species ranged between 127 and 557 mm and between 220 and 1120 mm, respectively, for sexes combined (Table 1).

Ageing Precision
The ageing precision evaluation returned relatively high reading variability for all four evaluated species.CV and APE values ranged between 47 and 49% and between 33 and 37%, respectively, for R. brachyura and T. marmorata.For D. oxyrinchus and R. clavata, slightly more precise results were observed (CV 30-34%; APE 21-26%) (Table 2).Nonetheless, the PA was similar for all evaluated species (44-52%) (Table 2).Figure 2 shows the CV and PA values, specifically for each modal age of the four examined species against the standard deviation.In general, PA is observed to be higher in younger age classes, showing values close to 75% (Figure 2D, D. oxyrinchus) or higher (Figure 2C, T. marmorata), while in the two species belonging to the genus Raja this value seems more stable around the 50% in all modal ages, with an inflection observed for specimens older than 5-6 years (Figure 2A,B).CV values also show this trend where lower reader variation is observed for younger modal ages and this variation increases with modal age (Figure 2).Logically, the reading standard deviation seemed to follow the opposite tendency, with higher values for older modal ages (Figure 2).

Growth Parameters
Due to the relatively small sample number of D. oxyrinchus and T. marmorata, it was only possible to estimate growth for the sexes combined.The AICc indicated that the mos precise models were the von Bertalanffy model with fixed t0 and Gompertz's model in fitting the observed data for D. oxyrinchus and T. marmorata, respectively (Table 3).Both species appeared to grow relatively slowly (D. oxyrinchus k = 0.101; T. marmorata k = 0.175 and to be capable of a long lifespan, with estimations of up to 30 years for D. oxyrinchus and up to 17 for T. marmorata (Table 3).
The age-at-length data were also insufficient to model the growth of male R brachyura, thus the species growth pattern was estimated only for females and combined sexes both in the Mediterranean Sea and the Atlantic Ocean (Table 3).In both study areas the vbt0p was the most precise model in describing the species growth for sexes combined according to the AICc, while for females this was the vbp model (Table 3).When considering combined sexes, R. brachyura appears to be a relatively slow growing species however females showed a faster growth rate (k = 0.397 in Atlantic Ocean, k = 0.429 in Mediterranean Sea).The observed ages ranged between 0 and 8 years in the Atlantic Ocean and between 0 and 10 years in the Mediterranean Sea.This resulted in a higher estimated lifespan for this species in the Mediterranean Sea (18 years) compared to in Atlantic waters (12 years) (Table 3).The Chen test comparison of the obtained vbp growth curve did not show a statistical difference between the areas (Fobs < Fcrit), thus this

Growth Parameters
Due to the relatively small sample number of D. oxyrinchus and T. marmorata, it was only possible to estimate growth for the sexes combined.The AICc indicated that the most precise models were the von Bertalanffy model with fixed t 0 and Gompertz's model in fitting the observed data for D. oxyrinchus and T. marmorata, respectively (Table 3).Both species appeared to grow relatively slowly (D. oxyrinchus k = 0.101; T. marmorata k = 0.175) and to be capable of a long lifespan, with estimations of up to 30 years for D. oxyrinchus and up to 17 for T. marmorata (Table 3).
The age-at-length data were also insufficient to model the growth of male R. brachyura, thus the species growth pattern was estimated only for females and combined sexes both in the Mediterranean Sea and the Atlantic Ocean (Table 3).In both study areas, the vbt0p was the most precise model in describing the species growth for sexes combined according to the AICc, while for females this was the vbp model (Table 3).When considering combined sexes, R. brachyura appears to be a relatively slow growing species, however females showed a faster growth rate (k = 0.397 in Atlantic Ocean, k = 0.429 in Mediterranean Sea).The observed ages ranged between 0 and 8 years in the Atlantic Ocean and between 0 and 10 years in the Mediterranean Sea.This resulted in a higher estimated lifespan for this species in the Mediterranean Sea (18 years) compared to in Atlantic waters (12 years) (Table 3).The Chen test comparison of the obtained vbp growth curve did not show a statistical difference between the areas (Fobs < Fcrit), thus this species follows a similar growth pattern in both the Atlantic and Mediterranean regions (Figure 3; Table 3).Table 3. Summary of the growth modelling results obtained for each species for females (F), males (M), and sexes combined (F + M) in the different geographical areas investigated.N is the number of analysed specimens, with information on size composition (TL max and TL min, in mm) and observed age range (age max and age min, in years).Growth model indicates the best-fitting model to the observed age-at-length data according to the AICc (vbp = unconstrained von Bertalanffy growth model, vbt0p = von Bertalanffy model with forced t 0 = 0, vbL1p = Gompertz' growth model, log.p = logistic growth model); TL ∞ is the asymptotic length (mm); k is the growth coefficient; TL1 and t 0 are, respectively, the length at age 1 and the theoretical length at time 0. The estimation of the lifespan and the growth performance index (φ) are also reported.Finally, the vbt0p was the best-fitting model to the observed age-at-length data fo clavata females and males in the Atlantic Ocean and for females from Mediterranean S The logistic model (log.p)returned the best-fitting results (AICc) for both males fr Mediterranean Sea and both sexes combined for both the investigated areas (Table 3) clavata seems to be a relatively slow-growing species, with males that appear to be capa of growing faster yet show a shorter lifespan than females.In contrast to R. brachyura, t distinct growth patterns were observed for R. clavata between the Atlantic Ocean and Mediterranean Sea (Chen test, Fobs > Fcrit).The Atlantic Ocean population produ higher growth rates compared to the Mediterranean population (Figure 3).

Species
The age-at-length data for all species are plotted in Figure 4 with the growth mo that provided the best-fitting results for the observations.Finally, the vbt0p was the best-fitting model to the observed age-at-length data for R. clavata females and males in the Atlantic Ocean and for females from Mediterranean Sea.The logistic model (log.p)returned the best-fitting results (AICc) for both males from Mediterranean Sea and both sexes combined for both the investigated areas (Table 3).R. clavata seems to be a relatively slow-growing species, with males that appear to be capable of growing faster yet show a shorter lifespan than females.In contrast to R. brachyura, two distinct growth patterns were observed for R. clavata between the Atlantic Ocean and the Mediterranean Sea (Chen test, Fobs > Fcrit).The Atlantic Ocean population produced higher growth rates compared to the Mediterranean population (Figure 3).
The age-at-length data for all species are plotted in Figure 4 with the growth model that provided the best-fitting results for the observations.

Discussion
This paper presents the first attempt at an intercalibration of age readers for elasmobranch species at a European level.For the first time, eleven international readers from five European countries took part in an age reading exercise that involved over 600 calcified structures extracted from four different batoid species.Vertebral preparation methods varied between institutions, although only vertebral sections were available for T. marmorata and D. oxyrinchus.Both vertebral sections and whole vertebrae were obtained for R. brachyura and R. clavata, which were collected in the Mediterranean Sea and the North-eastern Atlantic Ocean, respectively.Vertebrae were therefore analysed with the same preparation method for each geographical sampling area.
In consideration of the large number of scientists involved, the results obtained from the analysis of the ageing precision and reproducibility, although appearing far from the thresholds usually considered acceptable for elasmobranch ageing studies (sensu [12]), can be considered encouraging.Indexes such as the CV, the APE, and particularly the PA can easily be negatively affected by a high number of readers.Additionally, readers from different countries, while having a good experience level in interpreting hard structures, were also asked to read structures prepared following different protocols from those to which they were accustomed.This difference potentially played a role in the age reading variability observed in this study.
It is well known that reader experience is the most important factor affecting ageing precision.This has been confirmed in other age calibration studies, which compared other potential sources of bias such as the identification of first annulus or the interpretation of possible false rings [23].The present study also seems to confirm this assumption as, despite the application of different structure preparation methods, the best outcomes in terms of reading precision were obtained for R. clavata.This species is analysed in all the laboratories involved and ageing was consequently more familiar to the readers.The intercalibration results obtained should therefore be considered encouraging, and future reading exercises and workshops must be endorsed.
Although the growth modelling of the analysed species was not the principal purpose of the age reading exchange, the comparison of the obtained growth parameters with those reported in the extant literature (Table 4) revealed, in most cases, no major differences.The growth patterns of the two species investigated in both the Atlantic and the Mediterranean Sea, R. clavata and R. brachyura, appeared similar to previous observations in the two areas (Table 4).The main differences may be ascribed to the growth model selected, as only the common von Bertalanffy function was considered in many studies, or to the age estimation method (e.g., tagging [44]) and the hard structure used (e.g., caudal thorns [45,46]).The logistic, the vbp, and the vbt0p models provided the best fit to the length-at-age data following model selection with AIC.This is in accordance with the study of Thys et al. [29] where the best-fitting models were the logistic and vbp model (vbt0p not tested).Growth patterns for T. marmorata appeared to be in line with the literature [47,48] indicating this batoid as a slow-growing and long-lived species.Conversely, the calculated growth parameters for D. oxyrinchus, while comparable to those in Sardinian [49] and Tunisian waters [50], appeared rather different from those estimated by Yigin and Ismen [51].The data from the present study indicated a much faster growth rate and an asymptotic length of almost half the size previously reported.
It is also interesting to note that the different growth rates determined in the analysed species seem unlikely to be related to their trophic level.In fact, the two species that showed the faster growth rates, R. brachyura and R. clavata, are recognized as a specialist bony fish predator and as a generalist feeder, respectively [52].Similarly for the two slowgrowing species, T. marmorata specializes in hunting fish [48] and D. oxyrinchus is more generalist [51].Nonetheless, it could be interesting to investigate further how interactions with the environment could affect the growth of these species, and future studies on this aspect should be endorsed and welcomed.
The present study also observed different growth patterns for R. clavata and R. brachyura caught in the Atlantic Ocean and the Mediterranean Sea.R. brachyura produced similar growth patterns for both investigated areas, while R. clavata appeared to be capable of growing faster and larger in the Atlantic Ocean compared to the Mediterranean Sea.Although different preparation methods were used between the different areas for both species, i.e. whole vertebrae from the Atlantic and sectioned vertebrae from the Mediterranean, differences in growth patterns between areas were only observed for one species, R. clavata.This therefore suggests that these divergences in growth are not caused by the preparation method but may be linked to other factors such as different environmental conditions (e.g., water temperature, prey availability, nutrient levels, pollution, etc.) [21,53], fishing pressures [54], or strong regional genetic differentiation between Atlantic and Mediterranean populations of R. clavata [55].This is the first study investigating growth of these skates and including samples collected from the two different areas, and represents a first step towards a better understanding of the factors that may influence growth patterns of these batoids.

Conclusions
This study was the result of an international exchange with eleven readers representing five European countries.Age and growth parameters were successfully evaluated for four batoid species, R. clavata, R. brachyura, T. marmorata, and D. oxyrinchus, sampled in European waters, namely the Atlantic Ocean and the Mediterranean Sea.Although the precision of the age readings in this study is relatively low, the results are still encouraging considering the large number of age readers.Precision was higher at lower ages and decreased with age for all species.Alternative growth models were used to describe the age-at-length data where different models performed better than others, depending on the species, sex, and location.
The outcomes of this research, while preliminary, emphasise the need for intercalibration events involving large numbers of different laboratories and scientists from different countries.In this way, it is possible to increase ageing data quality for these ecologically important species, while providing solid inputs for their stock evaluation and management.

Figure 1 .
Figure 1.Summary of the preparation methods used for the analysed species in the Atlantic Ocean and Mediterranean Sea.Dots and lines represent the structure interpretation by different readers.

Figure 1 .
Figure 1.Summary of the preparation methods used for the analysed species in the Atlantic Ocean and Mediterranean Sea.Dots and lines represent the structure interpretation by different readers.

Figure 3 .
Figure 3.Comparison of the estimated growth curves and the observed age-at-length data obtai for R. clavata (left) and R. brachyura (right) in the Atlantic Ocean and Mediterranean Sea.

Figure 3 .
Figure 3.Comparison of the estimated growth curves and the observed age-at-length data obtained for R. clavata (left) and R. brachyura (right) in the Atlantic Ocean and Mediterranean Sea.

Figure 3 .
Figure 3.Comparison of the estimated growth curves and the observed age-at-length data obtained for R. clavata (left) and R. brachyura (right) in the Atlantic Ocean and Mediterranean Sea.

Figure 4 .
Figure 4. Age-at-length data for each species with the estimated growth curves from each best-fitting model: von Bertalanffy with constrained t 0 (vbt0p) curves for R. brachyura and D. oxyrinchus, logistic curve (log.p) for R. clavata, and Gompertz's curve (vbT1p) for T. marmorata.

Table 1 .
Sample composition per species with number of specimens (N), mean TL with standard error (SE), size range (in mm), mean age (in years) with SE, and age range for females (F), males (M), and sex combined (F + M).

Table 2 .
Summary of the ageing precision results, with vertebral centra preparation method, for each elasmobranch investigated species.

Table 4 .
Biogeographic comparison of the biological parameters of the investigated elasmobranch species with sampling details (number of samples, observed maximum TL and Age), and details of the growth model (vbp = unconstrained von Bertalanffy growth model, vbt0p = von Bertalanffy model with forced t 0 = 0, vbL1p = Gompertz' growth model, log.p = logistic growth model), with the type of data and the parameters of the growth model and growth performance index (φ).