Multi-Model Approach on Growth Estimation and Association With Life History Trait for Elasmobranchs

Age and growth information is essential for stock assessment of fish, and growth model selection may influence the accuracy of stock assessment and subsequent fishery management decision making. Previous descriptions of the age and growth of elasmobranchs relied mainly on the von Bertalanffy growth model (VBGM). However, it has been noted that sharks, skates and rays exhibit significant variety in size, shape, and life history traits. Given this variation, the VBGM may not necessarily provide the best fit for all elasmobranchs. This study attempts to improve the growth estimates by using multi-model approach to test four growth models—the VBGM, the two-parameter VBGM, the Robertson (Logistic) and the Gompertz models to fit observed or simulated length-at-age data for 38 species (44 cases) of elasmobranchs. The best-fit growth model was selected based on the bias corrected Akaike’s Information Criterion (AICc), the AICc difference, the AICc weight, the Bayesian Information Criterion (BIC), and the Leave-one-out cross-validation (LOOCV). The VBGM and two-parameter VBGM provide the best fit for species with slow growth and extended longevity (L∞ > 100 cm TL, 0.02 < k < 0.25 yr–1), such as pelagic sharks. For fast-growing small sharks (L∞ < 100 cm TL, kr or kg > 0.2 yr–1) in deep waters and for small-sized demersal skates/rays, the Robertson and the Gompertz models provide the best fit. The best-fit growth models for small sharks in shallow waters are the two-parameter VBGM and the Robertson model. Although it was found that the best-fit growth models for elasmobranchs were associated with their life history trait, exceptions were also noted. Therefore, a multi-model approach incorporating with the best-fit model selected for each group in this study was recommended in growth estimation for elasmobranchs.


INTRODUCTION
Elasmobranch life histories and assessments are different than those of teleosts for a number of reasons. Most elasmobranchs are characterized as slow growth, late maturation, extended longevity, few offspring, and low mortality (Hoenig and Gruber, 1990;Winemiller and Rose, 1992;King and McFarlane, 2003). Elasmobranchs also take a long time to recover when subjected to high fishing pressure (Stevens et al., 2000). The life history traits, particularly the reproductive traits of elasmobranchs, are more complex than those of teleosts, which are mostly oviparous.
Similar to many marine mammals, elasmobranchs are among the ocean's apex predators, and their life history characteristics make them vulnerable to overexploitation. A collapse of the elasmobranch population could result in imbalances in marine ecosystems (Stevens et al., 2000). Age, growth, and reproduction parameters are crucial for accurate stock assessment and evaluation of their population dynamics (Cailliet et al., 1986;Cailliet and Goldman, 2004). Information on age and growth can be used in natural mortality, longevity, and yield-per-recruit estimates (Ismen, 2003). However, age and growth study of elasmobranchs is more challenging compared with that of teleost species for several reasons. Some deep-water species are difficult to age such as spiny dogfish (Squalus spp.) (Campana et al., 2006), and the periodicity of band pair deposition may vary in different life stages such as shortfin mako shark Isurus oxyrinchus (Kinney et al., 2016) which may lead to the uncertainty in age estimation. In addition, ageing based on hard parts of elasmobranchs is typically underestimated in larger and older individuals based on bomb radiocarbon analysis (Passerotti et al., 2010;Harry, 2018;Natanson et al., 2018). These studies showed that vertebrae may not be useful for age estimation in some species once asymptotic length is reached. The underestimated age may lead to an underestimation of longevity and biased estimation of mortality which may affect the result of stock assessment and management measures.
In the 1950s, Beverton and Holt (1957) first applied the von Bertalanffy growth model (VBGM) to fish population dynamics. However, VBGM may not necessarily provide the best fit for all elasmobranchs (Cailliet and Goldman, 2004). Therefore, selecting the most appropriate growth model is important in stock assessment and fishery management of elasmobranchs (Gelsleichter et al., 1998).
Four growth models are commonly used in description of age and growth of elasmobranchs: the VBGM, the two-parameter VBGM, the Robertson (Logistic) model, and the Gompertz model . Derived from the allometric relationship between metabolic rate and body weight, the VBGM has been widely used to describe the growth of fish (Haddon, 2001). The ideas underpinning this model are that energy transformation during growth can be expressed as the difference between anabolism and catabolism and that the growth rate decreases exponentially with age (Pütter, 1920). While the Robertson and the Gompertz models are S-shaped curves with inflection points occurring at an intermediate age when the growth rate starts to decrease (Wang and Zuidhof, 2004). The inflection point of the Robertson model occurs at 50% of L ∞ , but it occurs at approximately 37% of L ∞ for the Gompertz model (Winsor, 1932). Under these two models, growth rates increase with age to a maximum at the inflection points and then decrease thereafter (Ricker, 1975(Ricker, , 1979. Discrepancies in life history among elasmobranch species are likely to affect the result of selecting the best-fit growth model. In addition, energy allocation for elasmobranchs varies with habitat and reproductive strategies; this may result in differences in growth. Numerous examples exist that used VBGM to estimate the age and growth of elasmobranchs for example, the pelagic thresher Alopias pelagicus (Liu et al., 1999), whiskery shark Furgaleus macki (Simpfendorfer et al., 2000), and undulate ray Raja undulate (Coelho and Erzini, 2002). Although the VBGM is the most commonly-used growth model, others have also been employed. The two-parameter VBGM has been successfully applied in describing growth of the bull shark Carcharhinus leucas (Neer et al., 2005), the shortfin mako shark Isurus oxyrinchus (Natanson et al., 2006;Chang and Liu, 2009;Semba et al., 2009), and the oceanic whitetip shark, Carcharhinus longimanus (Joung et al., 2016) for which size-at-birth data are available. An S-sharp model (Gompertz model) has been used to estimate the growth of the Pacific electric ray Torpedo californica (Neer and Cailliet, 2001), pelagic stingray Dasyatis violacea (Mollet et al., 2002), and the Kwangtung skate Dipturus kwangtungensis (Joung et al., 2015a). While, the Robertson (Logistic) model has been used to fit the length-at-age data for the spinner shark Carcharhinus brevipinna (Carlson and Baremore, 2005), the dusky shark Carcharhinus obscures (Natanson et al., 2014;Joung et al., 2015b), and Mediterranean skate Raja polystigma (Porcu et al., 2020). In recent years, multi-models have been used in age and growth studies of sharks. The use of multimodels in fitting length-at-age data is considered preferable to using single growth model (Katsanevakis and Maravelias, 2008;Smart et al., 2016). Romney and Campana (2009) examined four skate species and concluded that the VBGM best fit the winter and thorny skate Amblyraja radiata, while the Robertson model best fit the little skate, Leucoraja erinacea, and the Gompertz model best fit the smooth skate. Ebert et al. (2007) concluded that the VBGM provided the best fit for the Aleutian skate Bathyraja aleutica but that the Bering skate B. interrupta was more accurately described by the Robertson model. Katsanevakis (2006) also concluded that various growth models best described the growth of different chondrichthyan fish. Chen (pers. comm., 2004) applied several growth models to a variety of teleost species and concluded that the Richards and Robertson models best fit fish with slender and long lateral profile, while the VBGM and Gompertz model best fit other species. However, despite of the examination on the influence of taxa, body size, and sex segregation on candidate model fits for elasmobranchs (Smart et al., 2016), no other study has attempted to examine the correlation of growth model fits and life traits for sharks and rays.
Given the influence of growth model selection on the results of stock assessment, in particular, age-structured models, the objectives of this study were twofold: first, to fit the length-at-age data using multi-model approach, selecting the best-fit model for each species; and second, to group species on the basis of the best-fit model, examine the life history traits for each of these groups, and discuss the possible correlation involved. The hypothesis of this study is that age and growth of life-history closely related species can be best described by a similar growth model. We hope that our findings can provide an important future reference for the selection of the most appropriate growth model for elasmobranchs.

Source of Data
This study collected and analyzed the length-at-age data of 38 species, including the observations of vertebral band counts of 8 species and the age-length key data of 30 species from the literature (Supplementary Table 1). Of which, sex-specific data were available for 6 species. These species fell into 6 orders and 12 families (Supplementary Table 2). Two species were from Hemiscylliidae and Rhincodontidae (Orectolobiformes), 3 species were from Odontaspididae, Lamnidae, and Alopiidae (Lamniformes), 19 species were from Triakidae, Carcharhinidae, and Sphyrnidae (Carcharhiniformes), 2 species were from Etmopteridae (Squaliformes), 1species was from Rhinobatidae, 10 species were from Rajiformes (Rajidae), and 1 species was from Dasyatidae (Myliobatiformes). Life history parameters and ecological information, including habitat information, reproductive strategy, fecundity, reproductive cycle (R c ), and size at maturity (L mat ), were collected through literature searches in FishBase 1 as well as from published scientific articles and gray literature.

Data Process
In addition to observed length-at-age data, age-length key data adopted from the literature comprised the following data sets: (1) detailed age-specific length distribution data that can be directly fit by the growth models, (2) age-specific mean length with standard deviation, and (3) age-specific length interval. For data sets 2 and 3, a simulation process was used to generate (mimic) individual observations. For data set 2, a normal random number generator in Excel was used to generate simulated observations based on the sample size, mean length, and standard deviation for each age. The simulated data set for each age was adopted when its mean length and standard deviation approximate to the observed values within a limit of 0.1 for the observed mean length and standard deviation. For data set 3, the length distribution of each age was assumed to be a uniform distribution, and simulated observations were generated from a uniform random number generator in Excel based on the sample size and the maximum and minimum length of each age. The simulated data set for each age was adopted when its mean length (the average of the maximum and the minimum length) approximated to the observed mean within a limit of 0.5 (see examples in Supplementary Data Sheet).
The literature reveals an inconsistency in the way that body length is measured. Total length (TL) of 26 species, fork length (FL) of 4 species, precaudal length (PCL) of 7 species, and disc width (DW) of 1 species have all been used. Size-at-birth data were available for 23 of the 38 species (Supplementary Table 1). As most length measurements were in TL, our analysis converted other length measurements to TL using linear relationships between TL and other length measurements.

Growth Models
Three commonly used growth models, the VBGM (Beverton, 1954), the Robertson (Logistic) (Robertson, 1923) and the Gompertz (Gompertz, 1825), were fit to the length-at-age data for all species. For those species where size-at-birth data were available, an additional model, the two-parameter VBGM (Fabens, 1965), was also used. The nonlinear least squares function (nls) in R program (R Core Team, 2013) was used to estimate the parameters of each model. The four growth models are described as follows: (1) VBGM (Beverton, 1954) (2) Two-parameter VBGM (Fabens, 1965) (3) Robertson (Logistic) model (Robertson, 1923) where L t is the length at age t, L ∞ is the asymptotic length, k is the growth coefficient, t is the age (year from birth), L 0 is the size at birth, t 0 is the age at length 0, c r and k r are the parameters of the Robertson model, and c g and k g are the parameters of the Gompertz model.

Model Selection
The goodness of fit of the four growth models was compared based on the bias corrected Akaike's Information Criterion (AIC c ), the AIC c difference ( AIC c ), and the AIC c weight (w i ) (Burnham and Anderson, 2002), the Bayesian Information Criterion (BIC; Schwarz, 1978), and the Leave-one-out cross-validation (LOOCV; Breiman and Spector, 1992). The model with the smallest values of the aforementioned criteria was selected as the best growth model. AIC c was expressed as: Hurvich and Tsai, 1989), where n is the total sample size, MSE is the mean square of residuals, and K is the number of parameters estimated in the growth model. The AIC c difference ( AIC c ) of each model was calculated as the difference between AIC c,i and the lowest observed AIC c value (AIC cmin ). The corrected Akaike weight (W i ) is expressed as a percentage, which is useful when there are only minor differences in AIC c values among the growth models (Burnham and Anderson, 2002). The BIC can be expressed as: BIC = n × ln(MSE) + K × ln(n) (Schwarz, 1978). Minimizing a CV statistic is a common way to do model selection such as LOOCV. The LOOCV statistic was expressed as: where e [i] = y i −ŷ [i] . for i = 1,. . .,n andŷ [i] is the predicted value obtained when the model is estimated after deleting the i th observation. The best fit model can be found out by cross validating (CV) and choosing the one that has lowest mean squares error (MSE).

RESULTS
The goodness of fit of the 4 growth models for 44 cases (38 species and 6 with both sexes data) are listed in Supplementary  Table 3. The best-fit model for each species was selected based on the minimum values of the AICc, BIC, and LOOCV. The best-fit models selected based on the three criteria were consistent for 32 of 38 species examined in this study. Only slight difference among the criteria was found for the other 6 species. The BIC values of the VBGM (the best-fit model) were slightly larger than those of two-parameter VBGM for the silky and sandbar shark. The values of LOOCV of twoparameter VBGM (the best-fit model) were slightly larger than those of the VBGM for the scalloped hammerhead and the whitespotted bamboo shark. While, the value of LOOCV of the Gompertz model (the best-fit model) was larger than that of the two-parameter VBGM model for the gummy shark Mustelus antarcticus (both sexes). The BIC value of the Gompertz model (the best-fit model) was slightly larger than that of the two-parameter VBGM for cukoo ray Leucoraja naevus (Supplementary Table 3).

VBGM as the Best-Fit Growth Model
The VBGM provided the best fit for 9 species/cases including 7 shark species: the pelagic thresher, female shortfin mako, blue shark, night shark Carcharhinus signatus, and tiger shark; and 2 skates: female Kwangtung skate Dipturus kwangtungensis, and blue skate D. batis ( Table 1). All are large-size sharks or skates except the female Kwangtung skate. The blue skate had the highest L ∞ (L ∞ = 475.9 cm TL), while the female Kwangtung skate had the lowest (L ∞ = 57.0 cm TL). The blue skate had the slowest growth rate (k = 0.024 yr −1 ), while the female Kwangtung skate had the fastest growth rate (k = 0.294 yr −1 ).

Two-Parameter VBGM as the Best-Fit Growth Model
The two-parameter VBGM provided the best fit for 13 of the 23 species with size at birth data, of which 11 were sharks and 2 were skates/rays, comprising 39% and 12% of the 26 species of sharks and 12 species of skates/rays, respectively analyzed in this study (Figure 1).
Only the smooth lantern shark, and the Atlantic sharpnose shark had an L ∞ < 100 cm TL. The remaining species had an L ∞ > 100 cm TL. The exceptionally large-sized whale shark fell into this group, with L ∞ = 1580 cm TL and k = 0.020 yr −1 ( Table 2). Apart from that, the sand tiger shark had the highest L ∞ (L ∞ = 299.5 cm TL), while the smooth lantern shark had the lowest (L ∞ = 53.1 cm TL). The Atlantic sharpnose shark had the fastest growth rate (k = 0.582 yr −1 ), while the spinner shark had the slowest (k = 0.081 yr −1 ).
little skate had the fastest growth rate (k r = 0.667 yr −1 ), while the dusky shark had the slowest (k r = 0.131 yr −1 ).

Gompertz Model as the Best-Fit Growth Model
The Gompertz model provided the best fit for 6 species (7 cases) ( Table 4), including the yellownose skate, male smooth lantern shark, winter skate, gummy shark (both sexes), male roundel skate, and cuckoo ray Leucoraja naevus. The cuckoo ray, male smooth lantern shark, and male roundel skate fell into the smallsize category (L ∞ < 100 cm). While the yellownose skate and winter skate fell into the large-size category (L ∞ > 100 cm). The male roundel skate had the fastest growth rate (k g = 0.448 yr −1 ), while the female gummy shark had the slowest growth rate (k g = 0.126 yr −1 ). In summary, sharks were best fit by the two-parameter VBGM (39%) and VBGM (25%). Large sharks were best fit by the VBGM (35%) and the two-parameter VBGM (35%), small sharks were best fit by the two-parameter VBGM (50%) and Robertson model (38%), and skates/rays were best fit by the Robertson model (50%) and the Gompertz model (25%) (Figure 2). For 6 species with sex-specific data, same best-fit models were found in both sexes for the gummy shark (Gompertz model) and sharpspine skate (Robertson model). While, different best-fit models between sexes were found for the spinner shark (twoparameter VBGM and Robertson model), smooth lantern shark (two-parameter VBGM and Robertson model), roundel skate L ∞ , asymptotic total length (cm); k r , growth coefficient of Robertson model (yr −1 ); c r , Roberson parameter.  L ∞ , asymptotic total length (cm); k g , growth coefficient of Gompertz model (yr −1 ); c g , Gompertz parameter.
Frontiers in Marine Science | www.frontiersin.org (Robertson and Gompertz model), and Kwangtung skate (VBGM and Robertson model) (Tables 1-4). Cailliet and Goldman (2004) intensively reviewed 115 publications on the age and growth of 91 species of chondrichthyans, and Cailliet et al. (2006) presented updated information on 28 new studies. However, most of these studies did not provide either length-at-age or age-length key data. Thus, only 38 species with either observed length-at-age or age-length key data were used in this study.

Uncertainties
As mentioned above, observed length-at-age data were available for only 10 of 38 species. For the remaining 30 species, figures were generated (simulated) from age-length key data. Because such simulations may not be representative of real observations, there may be inaccuracies in the growth parameter estimates. Some species were represented by a small sample size -the common stingray, sand tiger shark, cuckoo ray, etc. Because of this, and due to a lack of small or large specimens, the size-atage data may not cover the whole life history of the fish. As a result, estimated parameters might not accurately describe the growth over the entire life history. Cailliet and Goldman (2004) stated that growth parameter estimates are greatly influenced by a lack of very young or old individuals. The existence of length-at-birth information may therefore have a significant effect on the choice of growth model. The size at birth and the growth in neonate stage of fish cannot be well described by the VBGM, the Robertson, and the Gompertz models (Cailliet and Goldman, 2004). On the other hand, the two-parameter VBGM can better describe the growth in the early life stage (Cailliet and Goldman, 2004). In this study, the simulated observation data set was adopted only when its mean length and standard deviation approach to observed values within limit values. Several simulations were made for each species and although minor variations in growth parameter estimates were noted, this had no effect on the selection of best-fit growth model. As for the model selection, Brewer et al. (2016) concluded that the AIC performs as well as the BIC when the heterogeneity is small, but the BIC performs better if the unobserved heterogeneity between data sets is large. In the present study, as no large heterogeneity between data was found, the same best-fit growth model was selected for 32 of the 38 species using the three criteria. In addition, only slightly inconsistence in the minimum values of the three criteria was found for the other 6 species. Therefore, the best-fit growth model selection in this study was believed to be robust. Pardo et al. (2013) recommended three-parameter VBGF and discourage use of the two-parameter VBGF because it results in substantially biased growth estimates even with slight variations in the value of fixed L 0. That study was based on simulated data (large sample size for each age class) yet not the observed length at age data. However, in most age and growth studies including this study, the sample sizes of very small and very large fish are less than that of medium size fish. So, the conclusion made by Pardo et al. (2013) may not be necessarily true for all elasmobranchs. In addition, compared with teleost fish, the size at birth can be easier estimated for elasmobranchs particularly the viviparous or ovoviviparous species. We believe the two-parameter VBGM can better describe the growth in the early life stage of some large pelagic elasmobranch species as VBGM, Robertson, and Gompertz model usually overestimate or underestimate the L 0 .

The Basic Theory of Growth Equation
In this study, growth of the bull shark was best described by the two-parameter VBGM. According to Schmid and Murru (1994), the juvenile bull shark allocates most of its energy to catabolism and waste and little to growth. Conversely, the chain dogfish Scyliorhinus rotifer, similar to the species best described by the Robertson model in this study, allocates most of its energy to growth and reproduction (Duffy, 1999). This suggests that small-size species allocate the most energy to growth and reproduction, while the converse is true for large-size species. The fact that the VBGM or two-parameter best fit large sharks, while the Robertson and Gompertz models best fit small sharks may therefore be related to their energy allocation. Cailliet and Goldman (2004) suggested that the Gompertz model better describes changes in body weight over time than changes in length. However, this hypothesis cannot be tested because age-at-weight data are not available in the literature. Most sharks are torpedo-shaped and large, while most skates and rays are flat and small. The best-fit growth model might be related to the ratio of size-at-maturity and maximum observed size. Species for which the VBGM and two-parameter VBGM provide the best fit are mostly sharks that tend to be late-maturing (Table 5), e.g., the pelagic thresher shark, bull shark and blacktip shark. Species best fit by the Robertson and Gompertz models tend, in contrast, to be early-maturing, such as the common stingray, sharpspine skate, and cuckoo ray ( Table 5).

Other Growth Models
The four-parameter Richards growth model is a general form of the VBGM, Robertson, and Gompertz models and is considered superior to the three-parameter growth models (Quinn and Deriso, 1999). However, in this study, the lack of large specimens and the relatively small sample size for certain species may cause the inconvergence of iterations in parameter estimation by nonlinear procedures. Araya and Cubillos (2006) used a two-phase growth model (TPGM) to fit for the porbeagle shark Lamna nasus and leopard shark Triakis semisfaciata and concluded that the model provided a better estimate of elasmobranch growth than the VBGM. This finding was later supported by Braccini et al. (2007) and Mejía-Falla et al. (2014) in their study of the piked spurdog Squalus megalops and the round stingray Urotrygon rogersi. The TPGF is a four or five-parameter growth model; the additional parameter is the age at which transition between two phases occurs. Because more detailed age-length data are required for this model, it was not applied in this study. Rogers-Bennett and Rogers (2016) also proposed a twostep growth model to account for the variation in growth rate during different life stages (juvenile vs adult) for marine V, VBGM; V 2 , two-parameter VBGM; R, Robertson model; G, Gompertz model; L mat , size at maturity; L mat /L ∞ , ratio of size at maturity and asymptotic length; f, fecundity; R c , reproductive cycle; f/R c , annual fecundity. *: Reproductive cycle is assumed to be 1 year invertebrates. However, because the size at maturity information is not necessarily available for all species used in this study, this approach has not been used.

Estimation of Parameters
In this study, L ∞ estimates derived by the VBGM and twoparameter VBGM were larger than those derived by the Robertson and Gompertz models (Figure 3). A similar finding has been documented by Katsanevakis and Maravelias (2008). Therefore, it seems that L ∞ is closely related to growth model selection.
In this study, the VBGM provided the best fit for 9 species/cases, with estimated k values of 0.024 -0.294 year −1 . All but one of these are moderate or slow-growing species (k < 0.2 year −1 ). Recent ageing studies of the blue shark supported that VBGM is best-fit model for the blue shark with k < 0.2 year −1 (k = 0.130 year −1 , 0.128-0.164 year −1 in the South Atlantic and the South Pacific (Joung et al., 2017. The two-parameter VBGM was the best fit for 13 species/cases, with estimated k values of 0.020 -0.583 year −1 . Most of these are also moderate-or slow-growing species (k < 0.2 year −1 ), the exceptions being the Atlantic sharpnose shark Rhizoprionodon terraenovae, finetooth shark Carcharhinus isodon, blacktip shark C. limbatus, and common guitarfish Rhinobatos (k > 0.2 year −1 ). The Robertson model provided the best fit for 14 species (15 cases), with estimated k r values of 0.131-0.667 year −1 . Most of these are fastgrowing species (k r > 0.2 year −1 ), the exception being the dusky shark (k r < 0.2 year −1 ). The Gompertz model was the best fit for 6 species (7 cases), with estimated k g values of 0.126-0.448 year −1 . These included moderate-growing species, the yellownose skate, winter skate, and cuckoo ray (0.15 < k g < 0.25 year −1 ). The VBGM and two-parameter VBGM provided the best fit for slowto moderate-growing species, while the Gompertz model was the best fit for moderate-growing species, and the Robertson model was the best fit for fast-growing species (Figure 4).

Contingency of Fitting Models
For the 9 species best fit by the VBGM, the second best choice was the Gompertz model (78%) that without size at birth information followed by the two-parameter VBGM (22%) that with size at birth data. For the species best fit by the two-parameter VBGM (13 species), the second-best choice was the VBGM (92%) and the Gompertz model (8%), while for the species best fit by the Robertson model (15 species), the second-best choice was the Gompertz model (100%). For the species best fit by the Gompertz model (7 species), the second-best choice was the two-parameter VBGM, and the Robertson model (43% each). In their study of elasmobranchs, Katsanevakis and Maravelias (2008) proposed four growth models in order of fit, as follows: Logistic-Gompertz-VBGM-Power (where Gompertz is the best choice, and Logistic and VBGM are the second-best choices). They concluded that the VBGM provided the best description of growth among elasmobranchs and bony fish. Our study arrived at a similar order of growth models, namely Robertson (Logistic)-Gompertz-VBGM-two-parameter VBGM, although it should be noted that the previous study separated species into sharks, skates, and rays. Mollet et al. (2002) found that the best-fit for the growth of the pelagic stingray was the Gompertz model. In this study, the growth of skates/rays is best described by the S-shaped Gompertz or Robertson models.

Comparison With Literature Results
Of the 38 species analyzed in this study, 12 have been previously fit with more than one growth model in the literature. Of these, our study found the same best-fit growth model for 10 species. The remaining 26 species have previously been described using VBGM alone, but our study found that 21 of these species are better fit by an alternative model. Thorson and Simpfendorfer (2009) have suggested using AIC c , AIC c weight, and multi-model inference to obtain the most appropriate model to describe fish growth. In this study, the best-fit growth model for each stock was selected based on similar criteria, AIC c , AIC c , the AIC c weight, the BIC, and the LOOCVadditionally suggesting that the derived results are reasonable.

The Best-Fit Growth Models Associate With Life History Traits
Based on their life history traits, three groups of sharks have been identified using cluster analysis (Liu et al., 2015) as follows. Group 1: large size, extended life span, slow growth, e.g., the bull shark, scalloped hammerhead shark, and oceanic whitetip shark. These are similar to the species best described in this study by the two-parameter VBGM. Group 2: small size, short life span, rapid growth, e.g., the smooth dogfish and blacknose shark. These are similar to the species best described in this study by the Robertson model. Group 3: late-maturing, moderate life span, e.g., the pelagic thresher shark, blue shark and night shark.  These are similar to the species best described in this study by the VBGM. This study found that the Robertson and Gompertz models provided the best fit for skates and rays. Those best described by the Robertson model are characterized by small size and rapid growth, e.g., the thorny skate, common stingray and little skate. Those best described by Gompertz have the characters of small or large size and moderate growth, e.g., the winter skate, yellownose skate, and cuckoo ray.
As mentioned above, the Orectolobiformes, Lamniformes and large-sized species of Carcharihidae are best fit by the VBGM or two-parameter VBGM. Rajiformes and Myliobatiformes are best fit by the Robertson or Gompertz model, while the Robertson model best describes the growth of small-size species of Carcharihidae.
Most species for which the VBGM or two-parameter VBGM provided the best fit are viviparous (70% and 54% respectively), while all species best described by the Gompertz models are oviparous. Species best described by the VBGM or twoparameter VBGM have lower annual fecundity and mature later (higher L mat /L ∞ ) than those best described by the Robertson or Gompertz models (Table 5).
Although VBGM has been widely used in fitting age and length data, where alternative models have not been tried and evaluated, the derived age structure may be biased and inaccurate (Roff, 1980). This will cause further errors in the estimates of mortality, yield per recruit, and stock assessment. If the Robertson or Gompertz models better describe the growth of certain species, variations in different life stages can be considered, and stock assessment will be improved (Carlson and Baremore, 2005). Katsanevakis and Maravelias (2008) suggested that multimodel inference should be used to estimate growth for aquatic taxa. The authors concluded that the weighted means of the growth parameters derived from different growth models based on AIC c weights are better estimators than those from a single growth model. However, given the definitions of growth coefficient vary among different growth models, this approach may not be appropriate and more discussion is needed. Smart et al. (2016) examined 74 elasmobranch growth studies to test the hypothesis that the best growth model is related the shape, reproductive mode, or life history trait using general linear models. The authors rejected the hypothesis as VBGM best fit all taxa and reproductive modes. It is concluded that the best growth model is not associated with the shape or reproductive mode of sharks and rays and recommended multi-model approach should be applied to growth estimate for elasmobranchs. However, we found the best-fit growth models for elasmobranchs were associated with their life history traits but with some exceptions in this study. These findings meet the basic theory of growth equation. Given this, it is suggested that the best-fit model selected for each group in this study should be included as the candidate when using a multi-model approach in growth estimation for elasmobranchs in the future. Additional growth model such as Schnute and Richards (1990) model should also be considered.

CONCLUSION
The best-fit growth model for elasmobranchs depends on their size and life history characteristics. VBGM provides the best fit for large pelagic sharks that are late-maturing and of moderate longevity. These include the pelagic thresher and blue sharks. The two-parameter VBGM best fits large pelagic sharks that are slow-growing and have extended longevity, such as the bull and oceanic whitetip sharks. The Robertson model is the best-fit for fast-growing small sharks that inhabit deep water. For small sharks in shallow waters, the two-parameter VBGM and the Robertson model provide the best description. The Robertson model is also the best-fit for medium and small-size demersal skates and rays, which are fast-growing and of short longevity, such as the smooth dogfish and thorny skate. The Gompertz model best fits large or small, median-growing skates and rays, such as the yellownose and winter skates. For the whale shark, with its huge size, slow growth, extended longevity, late maturity, and prolonged reproductive cycle, the two-parameter VBGM provides the best fit. Although it was found the bestfit growth models for elasmobranchs were associated with their life history traits, exceptions were also noted. Therefore, a multimodel approach incorporating the best-fit model selected for each group in this study was recommended in growth estimation for elasmobranchs.

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

AUTHOR CONTRIBUTIONS
K-ML and C-BW conceived and designed the experiments and performed the experiments. C-BW, W-PT, K-YS, and K-ML analysed the data. S-JJ, W-PT, and K-YS contributed reagents, materials, and analysis tools. K-ML and C-BW wrote the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This study was financially supported by the Ministry of Science and Technology, Taiwan under Contracts no. MOST 105-2313-B-019-005-MY3.