Growth and Age Validation of the Thornback Ray (Raja clavata Linnaeus, 1758) in the South Adriatic Sea (Central Mediterranean)

Raja clavata is the most widespread and landed skate species in the Mediterranean Basin. Despite its diffusion and economic importance, several aspects of its life history, such as age and growth, are poorly understood. This study evaluated the species’ growth in the South Adriatic Sea (Central Mediterranean Sea) and for the first time attempted an age validation through a tagging experiment. Thin sectioning of vertebral centra proved to be a more accurate preparation method in terms of age estimation precision and reproducibility than whole vertebral centrum staining (cobalt nitrate and ammonium sulfide technique). Marginal analysis showed a clear seasonal pattern, confirming the hypothesis of a single annulus deposition per year. A total of 291 vertebral centra were sampled and used for age estimation purposes. The oldest female was estimated to be 12 years old [total length (TL) = 89 cm], while the oldest male was aged 8 years (TL = 79.9 cm). Females were also found to be characterized by a slightly wider longevity range (ωL = 11.5, ωU = 16.8 years) than males (ωL = 7.8, ωU = 11.2 years). The von Bertalanffy growth curve fit the age and length data more accurately than the Gompertz and logistic models. Eighty-three thornback rays were tagged and released, of which two were recaptured. In both recaptured specimens, oxytetracycline marks were clearly visible. The band deposition after oxytetracycline injection and growth during the freedom period (about 1 year) were consistent with the age estimation method and criteria used and with the obtained growth results. Thus, the analysis of the vertebral centra extracted from the two recaptured specimens confirmed the hypothesis of the deposition of a single annulus per year and in general the age estimation criteria used in this study.

Raja clavata is the most widespread and landed skate species in the Mediterranean Basin. Despite its diffusion and economic importance, several aspects of its life history, such as age and growth, are poorly understood. This study evaluated the species' growth in the South Adriatic Sea (Central Mediterranean Sea) and for the first time attempted an age validation through a tagging experiment. Thin sectioning of vertebral centra proved to be a more accurate preparation method in terms of age estimation precision and reproducibility than whole vertebral centrum staining (cobalt nitrate and ammonium sulfide technique). Marginal analysis showed a clear seasonal pattern, confirming the hypothesis of a single annulus deposition per year. A total of 291 vertebral centra were sampled and used for age estimation purposes. The oldest female was estimated to be 12 years old [total length (TL) = 89 cm], while the oldest male was aged 8 years (TL = 79.9 cm). Females were also found to be characterized by a slightly wider longevity range (ω L = 11.5, ω U = 16.8 years) than males (ω L = 7.8, ω U = 11.2 years). The von Bertalanffy growth curve fit the age and length data more accurately than the Gompertz and logistic models. Eighty-three thornback rays were tagged and released, of which two were recaptured. In both recaptured specimens, oxytetracycline marks were clearly visible. The band deposition after oxytetracycline injection and growth during the freedom period (about 1 year) were consistent with the age estimation method and criteria used and with the obtained growth results. Thus, the analysis of the vertebral centra extracted from the two recaptured specimens confirmed the hypothesis of the deposition of a single annulus per year and in general the age estimation criteria used in this study.

INTRODUCTION
Elasmobranchs represent a considerable proportion of bycatch landings in the Mediterranean Basin (Food and Agriculture Organization [FAO], 2016) accounting for 2% of small-scale fishery and 0.8% of bottom trawl total landings in the South Adriatic Sea (Maiorano et al., 2019). The thornback ray (Raja clavata Linnaeus, 1758) is the most important species in terms of landing (Food and Agriculture Organization [FAO], 2018; International Council for the Exploration of the Sea [ICES], 2018), accounting for up to 38% of the total elasmobranch landings in Italian fisheries (Food and Agriculture Organization [FAO], 2016). In the South Adriatic Sea, the thornback ray is caught mostly by bottom longlines targeting the European hake and bottom trawls with mixed target fishes (Merluccius merluccius and Mullus barbatus) and crustacean decapods (Parapaeneus longirostris and Nehrops norvegicus) (Maiorano et al., 2019). Living in a wide bathymetric range (Krstulović-Šifner et al., 2008;Marongiu et al., 2017;Follesa et al., 2019b), the thornback ray is accessible to different fishing métiers (International Council for the Exploration of the Sea [ICES], 2018). Moreover, like most elasmobranchs, R. clavata exhibits a k-selected strategy, which is characterized by slow growth rate, long life spans, late-age sexual maturity, and low fecundity (Gallagher et al., 2005;Whittamore and McCarthy, 2005;Kadri et al., 2014), which makes the species particularly vulnerable to fishing exploitation. Indeed, several studies have shown that k-selected species with low intrinsic population regeneration potential are the most impacted by harvesting (Stevens et al., 2000;Ferretti et al., 2010;Dulvy et al., 2014).
The thornback ray has been listed as near threatened by the International Union for Conservation of Nature (Ellis et al., 2016); therefore, sustainable management of this resource seems to be necessary. However, stock status assessment, which is essential for developing management plans for cartilaginous species, and particularly for R. clavata, requires basic knowledge of their life history traits, such as growth rate, maturity, and fecundity, as fundamental data input. Most stock assessment models, especially the analytical ones, require robust data on the demographic structure of the exploited stocks. One of the first steps for running such models is the conversion of catches by length into catches by age. This routine is performed by means of an age slicing procedure that uses parameters from the von Bertalanffy growth function (VBGF) or age-length keys. For this reason, inappropriate or uncertain growth parameters or age-length keys for converting a size distribution into an age structure can lead to unreliable scientific advice (Scientific, Technical and Economic Committee for Fisheries [STECF], 2017). The evaluation of the stock status is in general a difficult work, in particular this is true for the elasmobranchs fish due to the lack and/or poor quality of the biological information for these species (Musick and Bonfil, 2005;Ellis et al., 2008;Cortés et al., 2015). An age overestimation could lead to an erroneous stock assessment, constructing a misleading scenario of a population composed of older individuals than in reality, which is thus seemingly affected by lower fishing mortality than in the real world. In the opposite scenario, a population would be considered younger, with an overestimation of its fishing mortality (Campana, 2001). Moreover, age errors also impact on the estimation of natural mortality and maturityat-age data and consequently the recruitment strength and spawning stock biomass (Maunder and Piner, 2015;Olsen et al., 2011). Ultimately, growth and age data affect mostly the short-term predictions of the stock status and the related management measures (Punt et al., 2008;Eero et al., 2015;Hüssy et al., 2016).
Validated age and growth parameters of R. clavata are still not available in many areas, especially in the Mediterranean Basin (Kadri et al., 2014). Some studies have performed age analyses typically using hard structures such as vertebral centra and dermal denticles (Serra-Pereira et al., 2008;Campana, 2014;Kadri et al., 2014;Bellodi et al., 2019). In both cases, age estimations are based on the annual patterns of the growth area in the form of a succession of opaque and transparent zones Campana, 2014;Carbonara and Follesa, 2019). One of the main problems of this kind of analysis is the presence of false rings and/or multiple bands on the hard structure, which are particularly pronounced in older specimens (Campana, 2014;Bellodi et al., 2019), reducing both the precision and the accuracy of the age data (International Council for the Exploration of the Sea [ICES], 2011), thus negatively affecting the stock exploitation analysis (Scientific, Technical and Economic Committee for Fisheries [STECF], 2017). In these contexts, it is highly recommended that age validation analysis be performed (International Council for the Exploration of the Sea [ICES], 2011[ICES], , 2013. In terms of effectiveness, direct age validation studies, such as mark-recapture experiments with chemically marked otoliths, present a high level of reliability (Campana, 2001; International Council for the Exploration of the Sea [ICES], 2011, 2020). However, mark-recapture surveys are generally expensive and time-consuming, and the success of tagging programs largely depends on tagging-related mortality, tag loss, and tag detection and reporting rates (Björnsson et al., 2011). Moreover, the injected chemical may not always leave a clear mark on the hard structure. This limitation is amplified in elasmobranchs because the chemical marker binds to calcium in bone structures, which is typically low in sharks and rays (Campana, 2014). However, structures such as caudal thorns, spines, neural arches, and vertebral centra contain systematic deposits of calcium phosphate and have often been proven reliable for age estimations in the context of mark-recapture experiments (Cailliet et al., 2006;McFarlane and King, 2006). Thus, when chemically marked individuals are recaptured, the number of opaque/transparent areas formed within the hard structure between release and recapture can be directly related to the time of release (Cailliet and Goldman, 2004;McQueen et al., 2018). Semi-direct validation studies such as marginal analyses represent the most commonly applied methodology for elasmobranch species (Cailliet and Goldman, 2004) and especially R. clavata (Gallagher et al., 2005;Kadri et al., 2014). Although this approach does not validate the absolute age, it is a useful tool for understanding the yearly deposition pattern of bands laid down on the hard structure (Campana, 2001).
This study combined a mark-recapture experiment and marginal analysis as a validation study to provide reliable age and growth data for the thornback ray in the South Adriatic Sea (Central Mediterranean). The research effort was focused on a thorough investigation of the complex life history dynamics of this valuable but vulnerable bycatch species, aiming to present solid data that can contribute to more effective management plans.

MATERIALS AND METHODS
In the period 2014-2018, thornback ray specimens were collected through commercial landings and discard monitoring (Data Collection Framework; EU Reg. 1004/2017) in fishing ports along the Italian South Adriatic coasts [Geographical Sub-Area (GSA) 18; GSA sensu Food and Agriculture Organization of the UN-General Fisheries Commission for the Mediterranean; Figure 1]. Additional samples were obtained from the Mediterranean International Trawl Survey (MEDITS; 2014-2018) (Spedicato et al., 2019) in the South Adriatic Sea, including Albanian and Montenegrin waters (GSA 18). The following data were collected from the collected specimens: total length (TL; in centimeters) to the nearest 0.1 cm, total weight (TW) to the nearest 0.1 g, sex, and maturity stages (AAVV, 2017;Follesa et al., 2019a).
A total of 8 specimens, ranging from 13.6 to 28.3 cm TL and 11.4-110.6 g TW, were sampled from the monitoring of discard, 274 specimens, ranging from 26.7 to 89 cm TL and 90-4180.9 g TW, from the monitoring of landings and 9 specimens, ranging from 32 to 76.6 cm TL and 745.3 and 2940.8 g TW, from the MEDITS survey, respectively.

Age Analysis
From each specimen, about 10 vertebrae were extracted from the thoracic area (immediately posterior to the scapular origin of the vertebral column), where the vertebral centra are bigger and thus easier to analyze (Campana, 2014;Bellodi et al., 2017Bellodi et al., , 2019. Before age analysis, the vertebrae were stored at −20 • C. Subsequently, the vertebral centra were defrosted, soft tissue (muscle, skin, and neural material) was removed by immersion in a 5% sodium hypochlorite solution (Goldman, 2006;Bellodi et al., 2017;Porcu et al., 2020), and the neural and hematic arches were removed with a scalpel. The centra were then measured considering the vertebral length (VL) and vertebral diameter (VD) to the nearest 0.1 mm (Goldman, 2006; Figure 2A) and embedded in epoxy resin (Prochima E-30).
Thin sections were cut with a cutting machine (IsoMet 1000) with a circular diamond blade (Bellodi et al., 2019;Porcu et al., 2020). The sections had an initial width of 0.5 mm. If necessary, they were thinned to increase their readability by grinding with sandpaper (grit 800). The final width ranged between 0.3 and 0.5 mm. The sections were then immersed in a clarification medium (sea water) and examined under a microscope with reflected light against a black background. Under these conditions, the opaque and transparent bands appeared white and dark, respectively ( Figure 2B). The length of the corpus calcareum (CCL) and the distance from the vertebra center to each transparent band (B 1 , B 2 ,. . .., Bn) were measured ( Figure 2B). The birthmark (BM) in each vertebral centrum was identified as the first clear mark corresponding to a change of angle in the CCL shape (Sulikowski et al., 2003; Figure 2B). Additionally, a subsample of 120 whole vertebral centra were stained using the cobalt nitrate and ammonium sulfide technique (Hoenig and Brown, 1988;Goldman, 2006) and observed with reflected light under a stereomicroscope (Leica S9D). Using this technique, the opaque winter bands appeared as darker circles ( Figure 2C).
The age and growth estimations were performed by band count, with one transparent band followed by an opaque one assumed to be an annulus (Sulikowski et al., 2003). An age was assigned to each vertebral section according to the age estimation scheme proposed by Carbonara and Follesa (2019) with a yearly resolution, accepting the birthday as January 1. Each section's age was independently estimated by two experienced readers blinded to the specimen's sex and size. Only concordant readings were included in the age analysis. The transparent band counts started from the first band pair that was easily recognizable after the birthmark (BM) ( Figure 2B). The transparent bands in the vertebral sections were counted using a microscope (Nikon Eclipse E400). The aging accuracy and the inter-reader agreement of the two preparation methods (thin sectioning and whole vertebral centrum staining) were evaluated by the coefficient of variation (CV) and percentage of agreement (PA) (Eltink, 2000; International Council for the Exploration of the Sea [ICES], 2011) as follows: where R is the number of times each fish's age is estimated, X ij is the ith age estimation of the jth fish, X j is the mean age calculated for the jth fish, and n diff is the difference in age estimation between the two readers. Moreover, the index of average percent error (IAPE) (Beamish and Fournier, 1981) was calculated: where N is the number of samples whose age was estimated, R is the number of readings, X ij is the ith age estimation of the jth fish, and X j is the average age calculated for the jth fish. The lower (ω L ) and upper (ω U ) species longevity ranges were calculated following Barnett et al. (2013). The lower range was obtained as while the upper range was calculated as where Tmax is the maximum observed age, E is the IAPE value calculated using the empirical data of the greatest 20% of the age groups observed, and c is an arbitrary constant (c = 1.4) to account for the possibility that the absolute maximum age was not observed.

Marginal Analysis
The pattern of annulus deposition on the vertebrae was analyzed by a semi-direct qualitative marginal analysis approach (Campana, 2001(Campana, , 2014. The marginal analysis considered the monthly evolution of the edge types (transparent or opaque) of the hard structures. The two edge types were defined when over 75% of the corpus calcareum margin appeared transparent or opaque. Vertebral sections in which about 50% of the edge was opaque or transparent were not included in the analysis .

Back-Calculation
On the basis of the obtained relationships between the TL and morphological measurements of vertebral centra (VL, VD, and CCL), the fish lengths at which different transparent bands were deposited were calculated using the back-calculation method. In particular, the linear-modified Dahl-Lea method (Francis, 1990) was applied, previously used in other Rajidae species (using thin vertebral sections) with reasonable biological accuracy for growth modeling (Sulikowski et al., 2007). The following formula was used: where a and b are the linear fit parameter estimates, Li is the length at ring i, Lc is the length at capture, CCc is the corpus calcareum at capture, and CCi is the corpus calcareum at ring i.

Tagging Experiment
During seven fishing trips (bottom longlines targeting European hakes) between 2014 and 2015, 132 specimens of R. clavata were captured, of which 83 (42 males and 41 females) were tagged externally with spaghetti tags and marked chemically with an intramuscular injection of 25 mg of oxytetracycline per kilogram of body weight (Tanaka et al., 1990;Gelsleichter et al., 1998; Figure 3). The condition of caught thornback rays was classified using the scale presented in Table 1 (Benoît et al., 2010;Dapp et al., 2016). To maximize the chance of survival and recapture, only thornback rays classified as condition 1 (Righton et al., 2006) were tagged. The tagging procedures were performed immediately after capture and lasted 4 min on average. The specimens were measured (TL and TW), tagged with the spaghetti tag, injected with oxytetracycline, and released back into the sea. In recaptured specimens, thin vertebral centrum

Statistical Analysis
To analyze the growth pattern, four different growth models were tested: the three-parameter von Bertalanffy growth function (3VBGF; von Bertalanffy, 1938) L = L ∞ (1 − e −K(t−t 0 ) ), the twoparameter von Bertalanffy growth function (2VBGF; Fabens, 1965 , the Gompertz sigmoid curve (Winsor, 1932) TL = L ∞ e e −K(t−I) , and the logistic curve (Richards, 1959) TL = L ∞ /1 + e −K(t−I) , where L ∞ is the species' theoretical maximum length, K is the growth coefficient, t is the observed age, t 0 is the hypothetical age of an individual with TL = 0, β is ( L ∞ − L 0 ) L −1 ∞ , L 0 is the individual's length at birth, and I is the age at the curve inflection point. Each growth model's fitness for the collected age data was evaluated by the Akaike information criterion (AIC; Akaike, 1974;Haddor, 2001).
The linear relationships between the TL and the measurements of the vertebral centra were tested using analysis of variance for regression. The obtained growth curves were compared using Chen's test (Chen et al., 1992).

Age Analysis
A total of 291 R. clavata (168 females and 123 males) were included in the age analysis. The TL ranged between 14.1 and FIGURE 4 | Length frequency distribution of the R. clavata specimens by sex, used for the age analysis.
FIGURE 5 | Relationship between age and total length by sex. The two recaptured specimens (cod. 1 and cod. 2) are also indicated.
89.0 cm in females and between 13.6 and 79.9 cm in males (Figure 4). The age ranged between 0 and 12 years in females and between 0 and 8 years in males (Figure 5). The IAPE values calculated from the age data were 3.42% for females and 3.87% for males. Figure 6 shows a comparison between the reading results of the two preparation methods (thin sectioning and silver nitrate staining) in terms of PA and CV (accuracy). The overall PA values were 90.7% for thin sectioning and 71.8% for silver nitrate staining. The CV values were 4.5 and 14.1%, respectively.

Marginal Analysis
The monthly margin evolution showed a prevalence of opaque edges (>50%) between June and October and a prevalence of transparent edges from December to May (Figure 7). This pattern of annulus deposition revealed that one transparent band is followed by one opaque band yearly.
The growth model parameters for females, males, and the sexes combined with their respective AIC values are displayed in Table 2. The analysis showed that 3VBGF was the model that best described the age-at-length data. The statistical comparison between females and males in the 3VBGF model showed no significant differences (Chen's test, F obs > F crit ).
To estimate the longevity range, the IAPE values of the greatest 20% of the age groups were used for both females and males (3.42 and 3.87%, respectively). The longevity in females ranged between 11.5 (ω L ) and 16.8 (ω U ) years with a T max of 12 years, whereas the longevity in males ranged between 7.8 (ω L ) and 11.2 (ω U ) years with a T max of 8 years.

Back-Calculation
The centrum morphometric descriptors (VL, VD, and CCL) and fish TL were significantly linearly correlated in both sexes (linear regression p < 0.05; Figure 8). A comparison between the sexes showed no significant differences (ANCOVA p > 0.05).
For the back-calculation analysis, the linear relationship between the TL and CCL of the combined sexes was used, as the difference between females and males was statistically insignificant. The results of the back-calculation analysis are reported in Table 3. Given that the transparent bands were laid down during the winter (see section "Marginal Analysis"), the first transparent band was considered age 1, the second age 2, and so on. In this way, the age/length data were achieved, and the growth parameters obtained by back-calculation were as follows: L 8 = 105.38 cm; K = 0.14 years −1 ; t 0 = 1.65 years. The comparison between 3VBGF from back-calculation and the readers' age estimates of the sexes combined showed no statistically significant differences (Chen's test, F obs > F crit ).

Tagging Experiment
A total of 132 R. clavata specimens were collected at depths ranging between 106 and 385 m during seven longline fishing trips (Figure 9). Most (63%) specimens were classified as condition 1, while conditions 2 and 3 accounted for 26 and 11%, respectively. In total, 41 females and 42 males were marked with spaghetti tags and oxytetracycline injections (cf. Figure 3). L∞, maximum asymptotic length; K, growth coefficient; t 0 , theoretical age at which the total length equals zero; IP, inflection point; AIC, Akaike information criterion.
Frontiers in Marine Science | www.frontiersin.org     The TL of the tagged specimens ranged between 43.4 and 86.7 cm (Figure 10). Two specimens were recaptured after 11 (cod. 1 captured November 2014, recaptured October 2015) and 12 (cod. 2 captured November 2015, recaptured December 2016) months (Figure 9), having spent 328 and 417 days of freedom, respectively. The distances between the release and recapture points were 2.26 km for cod. 1 and 1.52 km for cod. 2. Their differences in TL were 10.4 cm (from 50.5 to 60.9 cm) and 7.4 cm (from 62.3 to 69.7 cm), respectively, with average daily increments of 0.32 and 0.18 mm, respectively. In both specimens, the oxytetracycline mark was clearly visible, revealing a completed annulus (one transparent and one opaque ring) (Figure 11).

DISCUSSION
Thornback ray age analysis may be performed using several hard structures (e.g., vertebral centra) and preparation methods (Campana, 2014;Bellodi et al., 2019). The accuracy of age reading linked to each hard structure and preparation method represents the first step to establishing a protocol for age estimations Vitale et al., 2019), deepening in this perspective not only the influence of each processing technique on the data quality but also its cost and time requirements (International Council for the Exploration of the Sea [ICES], 2011). The hard structures mainly used in R. clavata age analyses are the dermal thorns and vertebral centra (Serra-Pereira et al., 2008;Campana, 2014;Kadri et al., 2014;Bellodi et al., 2019). The dermal thorns are usually read after staining. However, staining techniques often have several limitations; for example, the bands remain clear for limited periods, and staining techniques are often time-consuming and expensive (Serra-Pereira et al., 2008). Dermal thorns have additional disadvantages related to the readability of the different types of thorns (Gallagher et al., 2002(Gallagher et al., , 2006Gravendeel et al., 2002). Caudal thorns present a high shape variability, which affects their readability (Gallagher et al., 2002;Serra-Pereira et al., 2008). Moreover, they may be lost and replaced in the course of the fish's life history (Meunier and Panfili, 2002), representing an important obstacle to age analysis. Nevertheless, caudal thorns can sometimes be a useful (or even the only available) hard structure. The tail can be sampled without purchasing the entire specimen, and although they may be less readable than small thorns, thorns with large bases can hardly be lost (Serra-Pereira et al., 2008).
In comparison to thorns, vertebral centra have the obvious advantage that they cannot be lost. However, the readability of this structure linked to the preparation method is an important aspect of delineating the standardization steps for age analysis of this species Vitale et al., 2019). Campana (2014) evaluated the precision of the two most common preparation methods (whole vertebra staining and thin sectioning) for age estimations of batoids. The whole vertebral centrum cobalt nitrate and ammonium sulfide staining technique (Hoenig and Brown, 1988;Goldman, 2006) involves several passages through different chemical solutions to obtain the dark bands corresponding to a more intense deposition of calcium (calcium phosphate). On the other hand, in thin vertebral sectioning, seasonal differences are highlighted as bands with different transparency due to the seasonal fluctuations of calcium deposition. Our results showed that thin sectioning of unstained centra achieved higher accuracy (PA and CV) in all age groups. These results are in line with findings regarding other Rajidae species, such as R. brachyura (Porcu et al., 2015) and Dipturus oxyrinchus ). However, the success of each age analysis preparation method for cartilaginous species greatly depends on the overall amount of calcium in their cartilaginous structures, which in turn is linked to the general environmental conditions and diet (Flik et al., 1995;Allen et al., 2011). The best preparation methods depend on the species and specimens and the area (Clement, 1992;Wischniowski, 2008;Başusta et al., 2017;Bellodi et al., 2019). Therefore, it is difficult to define a general rule for age analysis preparation methods (Campana, 2014).
The growth model in fish is typically curvilinear, following a logistic or sigmoidal pattern (Gompertz, 1825;von Bertalanffy, 1938;Ricker, 1975). The strategy in both logistic and sigmoidal models is the same, with a growth rate increase in the juvenile (concave) and late-juvenile (sigmoidal) life stages, and a reduction later, often after sexual maturation (Sebens, 1987;Higgins et al., 2015;Manabe et al., 2018). After sexual maturation, a balance is often achieved between energy allocated for body growth and energy used for reproduction and other vital functions (Kozlowski, 1996). This trade-off may explain the two opposite life strategies of the main fish species: r-selection, whereby organisms exhibit sexual maturity at younger ages, small body sizes, and short life spans, while k-selected organisms present sexual maturation at older ages, large body size, and long life spans (MacArthur and Wilson, 1967).
In fish growth pattern modeling, generally (Katsanevakis, 2006) and in particular for Batoids species the comparison of several model is suggested (Cailliet and Goldman, 2004;Fisher et al., 2013;Mejía-Falla et al., 2014;Porcu et al., 2020). Thus, approaching the growth modeling without a priori "true model" seems more reliable (Burnham and Anderson, 2002;Katsanevakis, 2006) indeed, several species follow different growth trajectories (Katsanevakis and Maravelias, 2008). The estimated growth parameters are obviously model-dependent; therefore, ignoring the model selection uncertainty may introduce substantial bias (Katsanevakis and Maravelias, 2008). This has serious implications for stock status assessments (Carbonara et al., 2019b). The AIC allows the selection of the "best model" (or models) from a set of candidate models (Katsanevakis, 2006), minimizing the risks associated with an a priori selected growth model (Katsanevakis and Maravelias, 2008).
The 3VBGF is the most widely used model for fish, including ray species Campana, 2014). However, it requires a balanced sample, including both small and large specimens. It is not accurate when the sample size is small or lacking small or very large specimens (Cailliet and Goldman, 2004;Higgins et al., 2015). The 2VBGF can improve data fitting mainly when small specimens are absent from the sample (Gutteridge et al., 2013;Porcu et al., 2015). The Gompertz and logistic curves are sigmoid growth functions often used to describe growth in the larval/juvenile stage or in fish that present a low initial growth rate (Campana and Jones, 1992;Quist et al., 2012). These two models differ from each other in the inflection areas: in the logistic curve, the areas above and below the inflection point are symmetric, whereas in the Gompertz curve, they are asymmetric (Quist et al., 2012). The 3VBGF was the best fit for our dataset for both sexes, making the most biologically reasonable estimates for L 8 , higher of the L max (Gallagher et al., 2005;Serra-Pereira et al., 2008;Kadri et al., 2014).
In this study, the difference in growth rates between the sexes was statistically not-significant, even though females appeared to reach larger sizes and older ages. This is in line with findings from the Portuguese coasts (Serra-Pereira et al., 2008) and Welsh coastal waters (Gallagher et al., 2005;Whittamore and McCarthy, 2005). Conversely, significant differences were found between the sexes in Tunisian waters (Kadri et al., 2014). In our study, the similar growth rates of the two sexes were also observable in the relationship between the morphological vertebral centrum measurements and the TL. Although a comparison of growth curves from back-calculation and direct age analysis is not considered a validation method (Campana, 2001), the absence of statistically significant differences suggests that the age criteria and age scheme were consistent and used coherently .
A comparison of our estimated growth rate with those reported in the literature reveals similar patterns of a fast growth rate in the first years and a slowdown from about the fourth year, which coincides with reaching the size at first maturity (Cannizzaro et al., 1995;Gallagher et al., 2005;Kadri et al., 2014). In particular, growth curve comparisons between different studies reveal a difference in TL in the order of 7-9 cm in the central part of the curve for each age group (e.g., from 2 to 6 years of age) ( Table 4). This is reasonable, considering the wide geographic area of origin. However, several other factors may be responsible for higher variability of age estimation, apart from the geographic area  and genetic origin (Matić-Skoko et al., 2018): differences in sampling methodologies (Coggins et al., 2013); age schemes/criteria (International Council for the Exploration of the Sea [ICES], 2008;Hüssy et al., 2016;Carbonara et al., 2018); preparation methods (Smith et al., 2016); and reader experience (Kimura and Lyons, 1991;Carbonara et al., 2019b).These combined factors can thus lead to considerably greater differences than those observed in R. clavata. Then, it could be needful, taking in account the results obtained by the present work, considering the use of thin sectioned vertebra as suggested preparation method.
To our knowledge, there are no longevity data available for the thornback ray. The typical method for longevity estimations based on Taylor's (1960) formula employs the K value from the VBGF, which can easily lead to an overestimation (e.g., Ebert et al., 2009;Barnett et al., 2013;Mejía-Falla et al., 2014). To avoid this risk, in this study, longevity was estimated according to the method proposed by Barnett et al. (2013), which provides a lower and upper range calculated on the basis of the maximum observed size and the reading precision (IAPE) of the oldest age group. Females showed higher longevity than males, probably due to the difference in maximum observed age and TL between females and males.
Marginal analysis is considered a reliable method for investigating the deposition pattern of opaque and transparent bands (Campana, 2001). In our study, we found a prevalence of opaque bands between June and October, which is comparable to that in the Irish Sea, where the opaque band deposition is prevalent between July and November (Gallagher et al., 2005). This element not only allowed us to assign an annual frequency to the succession of bands but also agreed with the seasonal deposition patterns linked to the specific environmental conditions .
The tagging experiment was performed on the east side of the South Adriatic Sea in a bathymetric range of 100-500 m, according to the species' distribution in the Mediterranean (Follesa et al., 2019b). The condition of the specimens at the moment of capture was good, which seems an important result in itself, not only for the tagging experiment but also from a management perspective. The thornback ray is considered a bycatch species in the investigated area (Carbonara et al., 2019a;Maiorano et al., 2019). Thus, the release of the caught specimens may be considered a reasonable management measure in terms of survival chances (Bell and Lyle, 2016;Ellis et al., 2017Ellis et al., , 2018. The survival rate seems to be closely associated with the condition at capture for all ray species (Braccini et al., 2012), and especially for thornback rays captured by longlines (Ellis et al., 2018).
This study is the first to conduct a thornback ray tagand-release experiment in the Mediterranean Sea. Analogous experiments have been conducted mostly in NE Atlantic Ocean (North Sea, English Channel, and Irish Sea), allowing assessments of various aspects of the life history of R. clavata in terms of migration, stock unit, and size and sex distribution (Hunter et al., 2006;Ellis et al., 2008Ellis et al., , 2011Bird et al., 2020). The recapture rates in those studies were considerably higher than that in our study (2.4%), ranging between 6.5% (Ellis et al., 2008) and 51% (Hunter et al., 2006). However, the numbers of tagged fish in those studies were significantly larger: 197 in Hunter et al. (2006) and 4,152 in Ellis et al. (2008). Moreover, the tag-release area in our study represents a fishing ground for several fleets from several ports and countries (Lembo and Spedicato, 2011); therefore, it was difficult to intercept any recaptures.
Previous studies demonstrated that the thornback ray is capable of covering distances of several hundred kilometers (Hunter et al., 2006;Bird et al., 2020). The proximity of the recapture site to the release site in this study may not necessarily reflect restricted spatial movements for this species. Skates often exhibit a philopatric behavior (repeated return migrations) (Flowers et al., 2016), which is difficult to observe from markrecapture data (Hunter et al., 2005). However, in our study, the sampled specimens were released and recaptured in about the same time of the year.
In both recaptured specimens, the mark of oxytetracycline was clearly visible, showing a completed annulus with a transparent and an opaque band. The oxytetracycline mark in the vertebra is in line with the deposition pattern described here, as well as with the bands (opaque and transparent) subsequently laid down up to the margin. Even though we only had two recaptures, the deposition of one annulus per year and the band count can be considered valid (Campana, 2001). Moreover, the age data from the recaptured specimens are consistent with the growth curve obtained by the age analysis.
Although a larger number of OTC-injected recaptured thornback rays with greater time at liberty would be desirable in this type of study (Cailliet et al., 2006), the tag-recapture with chemical mark, that is considered one of the most power validation method (Campana, 2001;Cailliet and Goldman, 2004), allows to obtain reliable results also in studies with a low number of recapture. Indeed, the low number of specimens recaptured is not an atypical condition for age and growth studies on elasmobranchs that apply this age validation technique (Wintner and Dudley, 2000;Natanson et al., 2002;  L∞, maximum asymptotic length; K, growth coefficient; t 0 , theoretical age at which the total length equals zero. Frontiers in Marine Science | www.frontiersin.org Skomal and Natanson, 2003). Moreover the results obtained from the recapture are in line with the marginal analysis and vertebra thin section readings. Further studies with larger numbers of marked and recaptured thornback rays are needed to obtain complete validated growth data in light of the encouraging results of this study. The availability of the growth parameters validated for the thornback ray will allow to produce diagnosis of the state of the stocks and the consequent management measures for most of the Mediterranean areas where they have not yet been implemented, especially for the study area: the South Adriatic Sea.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

ETHICS STATEMENT
The animal study was reviewed and approved by the Committee on the Ethics of Animal Experiments of COISPA (Italian Ministry of Health 15/2015-UT).

AUTHOR CONTRIBUTIONS
PC contributed to experimental design, experiment development, data analysis, results processing, and wrote the manuscript. AB participated in data analysis, results processing, and manuscript revisions. MP took part in data analysis and experiment development. AM, CP, and MD contributed to data analysis and manuscript revisions. WZ participated in results processing and manuscript revisions. RC took part in manuscript revisions. LS contributed to experimental design and revision. MF contributed to data analysis, results processing, and manuscript revisions. All authors contributed to the article and approved the submitted version.

FUNDING
The tagging experiment was funded by RitMARE project (SP2_WP1_AZ4_UO08_D03). Moreover the samples from commercial fishery (landing and discard) are collected under the Data Collection Framework supported by the Italian Ministry of Agriculture, Food and Forestry Policy (MiPAAF) and by the European Commission (EU Reg. 199/2008).