Global synthesis of sea turtle von Bertalanffy growth parameters through Bayesian hierarchical modeling

Knowledge of sea turtle demographic rates is central to modeling their population dynamics, but few studies have quantitatively synthesized existing data globally. Here, we used a Bayesian hierarchical model to conduct a meta-analysis of published von Bertalanffy growth curve parameters (growth coefficient, K; asymptotic length, L∞) for chelonid sea turtles. We identified 34 studies for 5 of 6 extant chelonids that met minimum selection criteria. We implemented a suite of models that included a multivariate normal likelihood on the log-transformed values of the 2 parameters to evaluate the influence of species, population (regional management unit, RMU), parameter estimation method (mark−recapture, skeletochronology, length-frequency analysis), latitude, and sampled body size range (all sizes, no large, no small, no large or small) on growth parameter estimates. According to information criteria, the best model included a random effect of species. The second best model also included latitude as a fixed effect, but RMU, parameter estimation method, latitude, and sampled body size ultimately did not strongly influence the means or variances of K and L∞ among studies. The apparent lack of RMU effect on parameter estimates within species may be an artifact of the small number of RMUs with published growth parameter estimates. The species-specific, and in some cases RMU-specific, posterior means and standard deviations of K and L∞ from this study would be appropriate priors for future studies of growth in chelonid sea turtles or for models of population dynamics. We highlight the need for expanded study and synthesis of sea turtle somatic growth rates.


INTRODUCTION
A central feature of protected species conservation and management is the use of demographic models to evaluate population status and trends and inform management decisions. Yet, demographic models for long-lived, late-maturing migratory marine species such as sea turtles are notoriously difficult to parameterize due to their natural history, where long-term inaccessibility can limit data collection (Lascelles et al. 2014). For example, while sea turtle nest protection programs have in some cases provided robust counts of nests and nesting adult females through time, inadequate knowledge of juvenile abundances and demographic or vital rates have impeded understanding of the causes of trends (National Research Council 2010). Quantifying the means and variances in sea turtle vital rates, especially growth rates and age at maturation, therefore remain high-priority research areas in sea turtle population ecology (Hamann et al. 2010, National Research Council 2010, Rees et al. 2016. Such information is fundamental to accurately understanding and predicting the dynamics of sea turtle populations in addition to prescribing management targets, quantifying sustainable biological re moval levels, and evaluating likely responses to management actions or environmental change scenarios (Crouse et al. 1987, Crowder et al. 1994, Casale & Heppell 2016, Piacenza et al. 2019.
As free-living sea turtles are nearly impossible to age, somatic growth curves (i.e. size-at-age relationships) have served as the principal mechanism by which age, life stage duration, and age at maturity are estimated and ultimately integrated into population models (e.g. Crouse et al. 1987, Heppell et al. 1996, Chaloupka 2002, Mazaris & Matsinos 2006, Casale & Heppell 2016. Sea turtles are characterized by rapid growth during the first few years of life followed by monotonically decreasing growth rates until maturity, at which point energy is redirected to reproduction and growth becomes negligible (Omeyer et al. 2017(Omeyer et al. , 2018. Although studies of in-water life stages are dwarfed by those of adult nesting females, somatic growth rates have been measured across a wide range of populations through capture−mark−recapture of turtles on their foraging grounds (e.g. Frazer & Ehrhart 1985, Chaloupka & Limpus 1997, Tanaka 2009), analyses of stranded turtle humeri and scleral ossicles (e.g. Zug et al. 1986, Snover & Hohn 2004, Avens et al. 2009), and analyses of length-frequency data (e.g. Bjorndal & Bolten 1995, Bjorndal et al. 2001, Casale et al. 2011b). However, a comprehensive, cross-population synthesis of existing sea turtle growth studies, particularly those that report somatic growth curves, has not been conducted for most species (but see Bjorndal et al. 2017, Omeyer et al. 2017. Given the variable recovery trajectories of sea turtle populations globally and poor understanding of underlying drivers (Wallace et al. 2011, Mazaris et al. 2017, Valdivia et al. 2019, there is a critical need to synthesize existing data to identify data gaps, aid parameterization of population models, illuminate mechanisms underpinning trends, and set research priorities. Somatic growth curves are most accurate and useful when they encompass all life stages, but individual growth studies tend to sample only a subset of a species' life cycle. Sea turtles exhibit wide geographic distributions and undergo multiple ontogenetic habitat shifts during their lifetime (Chaloupka & Musick 1997, Bolten 2003, which makes characterizing growth rates across all life stages within a single study challenging. For example, there is a dearth of growth rate data for sea turtle pelagic life stages, during which turtles occupy open-ocean habitats for multiple years (e.g. Reich et al. 2007, Avens et al. 2013. Although this life stage remains understudied, pelagic stage growth rates have in some cases been obtained through skeletochronological analyses of stranded turtle humeri (Avens et al. 2013, Ramirez et al. 2015, Turner Tomaszewicz et al. 2017, which retain multi-year growth records that can backfill this data gap, or through study of turtle growth rates near islands that are waypoints during developmental migrations (e.g. The Azores; Bolten et al. 1993, Bjorndal et al. 2000b, 2003. Data for large subadult turtles are similarly sparse due to their relatively low abundance but can be estimated through long-term mark−recapture programs, skeletochronological studies, or study of adult turtles (NMFS SEFSC 2001, Omeyer et al. 2018. As a result of these limitations, individual growth studies tend to present an incomplete picture of a population or species' growth dynamics. Integrating data across studies within a meta-analytic framework may help fill in data gaps of individual studies and provide more precise estimates of population-and speciesspecific growth parameters (Thorson et al. 2015).
Growth rates vary considerably in sea turtles due to a suite of biological and environmental factors that can cumulatively drive a wide variation in size and age at maturation among populations and species (Avens & Snover 2013, Omeyer et al. 2017. Major drivers of variation include temperature (Wallace & Jones 2008, Marn et al. 2017, habitat use (Hawkes et al. 2006, Hatase et al. 2010, population density (Bjorndal et al. 2000a, Balazs & Chaloupka 2004, prey availability and distribution (Balazs 1982, Diez & van Dam 2002, and diet quality and composition (Peckham et al. 2011, Ramirez et al. 2020b), among others. Temperature, which is in turn inversely related to latitude, is expected to be a particularly important driver of sea turtle somatic growth variation because they are ectothermic reptiles (Snover et al. 2015, Marn et al. 2017, Stubbs et al. 2020. For example, Bjorndal et al. (2013 identified distinct latitudinal changes in growth for loggerhead Caretta caretta and hawksbill Eretmochelys imbricata sea turtles across the western North Atlantic and Caribbean Sea. Integrating sea turtle somatic growth data within a metaanalytic framework may therefore provide important insight into population-and species-specific sources of growth variation, including potential latitudinal drivers of variation (Helser & Lai 2004, Gaertner et al. 2008. Moreover, such an analysis may elucidate the generality or specificity of sea turtle somatic growth curves within species, and thereby the suitability of extension of population-or species-specific parameters to populations lacking robust somatic growth data. Bayesian methods are often used as the framework for implementing meta-analysis as they allow for data integration within a flexible and consistent statistical framework and better accounting of uncertainty (Pilling et al. 2002, Gelman 2006, Helser et al. 2007. For example, within hierarchical meta-analysis, variations in parameters between populations and species are treated as random effects (Thorson et al. 2015), thereby allowing for sharing of information across groups through partial pooling. This tends to pull parameter estimates for 'data-poor' populations and species towards the group mean (pooled) and away from the observed data (unpooled; Helser et al. 2007). In contrast, parameter estimates for 'data-rich' populations or species are less influenced by the hierarchical structure. Thus, the method improves parameter estimates for less-studied groups, while only marginally influencing parameter estimates for more studied groups (Hoenig et al. 2016), making the best use of available data for all groups.
The primary goal of this study was to synthesize existing sea turtle growth data within a meta-analytic framework to develop species-and populationspecific parameter estimates, identify critical data gaps, and evaluate the potential influence of multiple sources of error (e.g. methodological vs. ecological). To this end, we first performed a global literature review of somatic growth studies across all sea turtle species, specifically identifying those that fit somatic growth functions to growth increment data. We then modeled the parameters of the von Bertalanffy growth curve (growth coefficient, K; asymptotic length, L ∞ ), which is the most commonly used growth function for sea turtles, by integrating data from these studies with a series of Bayesian hierarchical models to characterize means and variances of these parameters for each species and population. We specifically evaluated the potential influence of species, population, parameter estimation method (mark−recapture, skeletochronology, length-frequency analysis), study latitude, and body size range of turtles sampled in each study on growth parameter estimates by combining all studies into a single likelihood within each model. We believe outputs from this synthesis of sea turtle growth studies will provide important priors for growth parameters in demographic models, particularly for species and regions where somatic growth data may be lacking. When used in combination with additional demographic data, either existing or newly collected, they may aid in the study of sea turtle population dynamics and prioritization of research and conservation efforts for sea turtles globally.

Literature review
We first performed a structured literature search to identify and extract data from published sea turtle growth studies. The search was conducted in Web of Science, Google Scholar, the ProQuest Biological Science Database, and the Sea Turtle Online Bibliography (Archie Carr Center for Sea Turtle Research, University of Florida). The search terms 'growth rate,' 'maturity,' 'mark−recapture,' 'length-frequency analysis,' and 'skeletochronology' were used in conjunction with 'sea turtle' and the names of the 6 extant sea turtle genera (Caretta, Chelonia, Lepidochelys, Eretmochelys, Natator, Dermochelys). We then performed a secondary search for pertinent literature by reviewing the International Sea Turtle Society (ISTS) Conference Proceedings, Indian Ocean Turtle Newsletter, Marine Turtle Newsletter, Sea turtles of India website, and a suite of books and journals known to contain sea turtle publications (for a full list, see Table S1 in the Supplement at www. int-res. com/ articles/ suppl/ m657p191_ supp .pdf). Lastly, we performed an unstructured literature search of the reference lists of all papers and reports identified in the structured searches.
We initially identified 195 sea turtle growth studies that contained somatic growth rate data or growth models. To allow for comparisons across studies, this initial pool of studies was reduced to only those that reported growth models. Within these studies, growth was described with Gompertz (n = 11), logistic (n = 18), or von Bertalanffy growth models (n = 59), or with generalized additive models (GAMs; n = 33). Based on the results of this literature review, we selected the von Bertalanffy growth function as the focus of our meta-analysis given that it was the most commonly used functional form and one of the most flexible and easily applied approaches to describing length-at-age. In addition, GAM smoothing spline fits cannot be easily integrated across studies. The von Bertalanffy growth equation is: (1) where L t = average length at age t in straight carapace length (SCL), L ∞ = asymptotic length, t 0 = age when L t = 0, and K = Brody's growth coefficient, a measure of the rate at which animals grow when young. The von Bertalanffy curve is non-linear, assuming rapid growth at the beginning of life and a slow reduction in growth until sexual maturity, after which point growth rate asymptotes at a certain mean body length (L ∞ ). Studies subsequently needed to meet 4 primary criteria to be included in the meta-analysis. First, the publication needed to report L ∞ and K estimates, identify study locations, and provide sample sizes. Sample size was recorded as the number of data points used to estimate the von Bertalanffy growth curve or the number of sea turtles sampled; these could be different numbers in cases where age was estimated from multiple skeletal growth rings from the same individual or when there were multiple recaptures of the same individual. Study location was recorded as the latitude and longitude at the midpoint of the study area, which in many cases spanned hundreds of kilometers. Second, the publication needed to identify the types of growth data and methodology used to derive the parameter estimates (e.g. mark−recapture, skeletochronology, length-frequency analysis). Third, L ∞ and K needed to be estimated from data rather than given assumed values based on ad hoc relationships. Fourth, data must have come from wild sea turtles only (i.e. not captive individuals), or in the case of head-started sea turtles, derived from growth rates after the sea turtles were released from captivity.
The resulting list of studies was then reduced further to remove redundancies and anomalous data. When multiple studies used the same dataset, only the study that utilized the most comprehensive dataset was included in the meta-analysis. Region-or population-specific growth curves generated within single studies were kept separate (e.g. Avens et al. 2017). In one case , separate curves were provided for males and females. As we did not investigate the influence of sex on growth, we averaged the parameter estimates to generate a curve for both sexes combined; most studies did not fit separate growth curves for each sex. In 2 cases, i.e. Hart et al. (2013) for hawksbill sea turtles and Snover et al. (2007a) for loggerhead sea turtles, data on length at capture, length at recapture, and time at liberty were reported for individual turtles, which allowed us to estimate a von Bertalanffy growth curve using the Fabens method from these data (see Table S2 and Fig. S1). We used the Fabens method for simplicity and because it was used in many of the mark− recapture studies included in the meta-analysis.
However, methods that more explicitly model uncertainty in estimating ages may provide more accurate estimates of L ∞ and K in future analyses (Francis 1988, Zhang et al. 2009). Lastly, because we could not reasonably include temperature as a covariate in our analysis, we excluded data from Eguchi et al. (2012), who reported estimated growth coefficients (K = 0.20) 2 to 10 times higher than those reported in other green turtle studies (K = 0.02−0.09). These turtles occupied habitats in an unnatural state where heated effluent discharged from a local power plant acted as a thermal refuge during winter months, which likely extended active periods, enhanced local productivity, and promoted fast growth (Eguchi et al. 2012). This power plant ceased operation in 2010 and turtle growth rates are expected to slow as the system reverts to its natural state (Eguchi et al. 2012(Eguchi et al. , 2020. The full list of excluded studies can be found in Table S3. Most excluded studies were early analyses of datasets that were later expanded, fixed either L ∞ and K at a specified value, or did not report both L ∞ and K.
Leatherback Dermochelys coriacea and flatback Natator depressus sea turtles were excluded from our analysis. Our model assumes the species effect on each growth parameter is drawn from a common distribution. As leatherback sea turtles belong to a different taxonomic family and grow to much larger sizes than the other 6 sea turtle species, their species effect was presumed to be different than the other species. However, there were too few leatherback studies to estimate a family effect and they were therefore excluded from this study. Similarly, flatback sea turtles lack published von Bertalanffy growth parameter estimates and were necessarily excluded from the analysis.
Our final growth dataset comprised 37 separate growth curves from 34 studies (Table 1). At least 1 study was identified with published von Bertalanffy growth parameter estimates for 5 of the 7 extant species of sea turtle: loggerhead (n = 15), green (n = 9), Kemp's ridley (n = 9), hawksbill (n = 3), and olive ridley Lepidochelys olivacea (n = 1). Due to lack of published t 0 values in studies using mark−recapture or cohort tracking methods, and a few of the skeletochronology studies, we excluded t 0 from the metaanalysis. Some studies presented L ∞ in curved carapace length (CCL), while the majority used SCL. When necessary, we used species-and region-specific body size conversion equations to convert all CCL estimates to SCL (Table S4). Studies were categorized by species and regional management unit (RMU, Table  S5, Wallace et al. 2010), an ecologically relevant bio-geographical unit intermediate to nesting population and species. Wallace et al. (2010) identified a single RMU for Kemp's ridleys, but we used 2 RMUs in our analysis given known differences in somatic growth rates for Kemp's ridleys that inhabit the Gulf of Mex-ico (GoM) and the US Atlantic Coast (Avens et al. 2017, Ramirez et al. 2020a. Within the literature, 3 categories of methods were used to generate growth rates and size-at-age data and estimate the von Bertalanffy growth curve para - . Sample size is given as number of turtles, with the total sample size in parentheses if the method used multiple data points from the same individual. SCL: straight carapace length; L ∞ : asymptotic length; K: growth coefficient meters: mark−recapture, skeletochronology, and length-frequency analysis. Each approach introduces different sources of error that contribute to uncertainty in estimates and may bias results (e.g. turtle or bone measurement error, back-calculation error; Chasco et al. 2020). We therefore sought to account for these sources of uncertainty in our analyses.
(1) Mark− recapture methods used information on the length at first capture and length at recapture to calculate growth increments. Growth increments, along with information on the time at liberty, were used to fit a variation of the von Bertalanffy growth curve re-parameterized to fit growth increments, such as the Fabens method (e.g. Hart et al. 2013). When studies reported parameter estimates for different time at liberty intervals, only estimates from the longest time interval, generally ≥ 365 d, were used.
(2) Skeletochronological methods were based on counting and measuring humerus bone growth rings to backcalculate body length based on allometric relationships and estimated age (see Snover et al. 2007a). Growth curves were generated by fitting Eq.
(1) to the paired age and length data (e.g. Goshe et al. 2010). (3) Length-frequency methods, such as the MULTIFAN approach (Fournier et al. 1990), followed cohorts across years by tracking nodes in the length-frequency distribution in each year. Since this method combines measured lengths in each node with estimates of age in each node, the method can be used to estimate von Bertalanffy growth parameters (e.g. Bjorndal et al. 2000b). For all 3 methods, parameter estimates are sensitive to the size range of the sea turtles sampled. For example, a von Bertalanffy fit to data concentrated in early life stages may overestimate L ∞ , while the lack of early life stage data may result in a poorly estimated K value (Gaertner et al. 2008). Therefore, we classified each study into 1 of 4 categories of size ranges sampled: complete (minimum size sampled < 0.25 L ∞ and maximum size sampled > 0.75 L ∞ , 6 studies), no large (minimum < 0.25 and maximum < 0.75 L ∞ , 5 studies), no small (minimum > 0.25 and maximum > 0.75 L ∞ , 23 studies), and no large or small (minimum >0.25 and maximum <0.75 L ∞ , 5 studies), and included this variable as a categorical fixed effect in the models. Finally, because sea turtle growth rates are ex pected to broadly vary with temperature, we also included latitude (absolute value) as a numerical variable, recorded as the latitude at the midpoint of the study area reported in each paper. Specifically including temperature as a covariate in the models was not feasible given that most studies spanned multiple years and wide geographic areas, but here we considered latitude a proxy for mean temperature.

Hierarchical modeling
We used 9 linear mixed-effects models to estimate L ∞ and K for species and RMUs of turtles from the 37 growth curves. The response variables were the natural logs of the von Bertalanffy growth parameters L ∞ and K from each study because the 2 variables are correlated and have a non-linear relationship (Pilling et al. 2002). Species and RMUs were included as random effects so that we could separate the variance components between species and between RMUs within species. When an RMU effect was included, it was nested within species. The effect of parameter estimation method (skeletochronology, mark−recapture, length-frequency analysis), body size range sampled, and latitude were included as fixed effects. Model configurations are summarized in Table 2. All models had multivariate normal (MVN) likelihoods, defined as: ( 2) where θ i is the parameter vector of log(L ∞ ) and log(K ) for study i, μ i is the vector of mean values predicted for study i, and Σ is the error variance−covariance matrix. The models varied in how the means and variances were treated, as well as which fixed and random effects were included.
In Model A, the mean was: where individual studies are denoted by i = 1 to N and j = 1 (L ∞ ) or 2 (K ); μ is the expected mean log value of the parameter; α is the intercept; β is the fixed effect of parameter estimation method X i used in study i, where X i = 1 (mark−recapture), 2 (skeletochronology), and 3 (length-frequency analysis); and γ Z[i],j is the effect of RMU on the log mean, where Z i indicates the RMU associated with study i. Both α and β were given vague normal priors with mean 0 and standard deviation 1000. The RMU effect γ k for RMU k was drawn from a multivariate normal distribution associated with its species: where γ s [k] is the vector of species effects for the species s associated with RMU k, and Σ p is the variance− covariance matrix across RMUs, defined as: where σ 2 p,L is the variance among RMUs in log(L ∞ ), σ 2 p,K is the variance among RMUs in log(K ), and log ( ) The species effects γ s are similarly drawn from a multivariate normal distribution: with a mean of 0, and a variance−covariance matrix Σ s defined as: where σ 2 s,L , σ 2 s,K , and σ s,L,K are the variance in log(L ∞ ), variance in log(K ), and covariance among species. The variances parameters σ 2 p,L , σ 2 p,K , σ 2 s,L , and σ 2 s,K each had a vague inverse gamma prior, while the covariances σ s,L,K and σ p,L,K were defined as the correlation between the 2 variables, given vague uniform priors between −1 and 1, times the standard deviations of the 2 variables. The likelihood was multivariate normal (Eq. 2), with the mean vector defined by Eq. (3), and error variance−covariance matrix Σ given a vague prior on the inverse of the variance− covariance matrix using the Wishart distribution.
Model B was the same as Model A, except that instead of having an RMU effect nested within species, there was only a random effect on species, defined as in Eqs. (6) and (7). Model C was the same as Model A, except that the error variance− covariance matrix was allowed to be different for the parameter estimation method employed for each study, each given a vague prior on the inverse of the variance− covariance matrix using the Wishart distribution. Models D through I were based on either Model A or Model B, but varied which fixed effects were included (Table 2). For example, in Model H, the mean was: where β 2 is the fixed effect of the size range category X 2 and γ s was the random effect of species S i from study i on parameter j. In Model I, the mean was: where β 3 is the slope associated with the latitude X 3 from study i for parameter j.

Model evaluation and selection
The posterior distributions for each model were calculated using Markov chain Monte Carlo (MCMC) through the software JAGS run in R (version 3.5.3, R Core Team 2019) with R2jags (Plummer 2016). We ran each model with 2 chains of 110 000 iterations, discarding the first 10 000 as a burn-in and with a thin of 4. Each model was considered converged when the Gelman-Rubin diagnostic (R ) value was below 1.05 and effective sample size (n effective ) was above 300 (Lunn et al. 2013). Additional MCMC iterations were run as needed to achieve this standard. We then analyzed the residual plots of each model to ensure they adequately described the data. Model adequacy was also evaluated with the Pareto k diagnostic, where data points with values < 0.5 are considered to be well fitted (Vehtari et al. 2017). The best model was chosen using information criteria including the deviance information criterion (DIC; Spiegelhalter et al. 2014) calculated by JAGS, the widely applicable information criterion (WAIC; Watanabe 2013), and leave one out information criterion (LOOIC) calculated by the R package 'loo' (Vehtari et al. 2017).

Δ) between each model and the best model (where 0 is best) in the deviance information criterion (DIC), widely applicable information criterion (WAIC), and leave-one-out information criterion (LOOIC) and the LOOIC model weights (w).
Method: technique used to derive parameter estimates (skeletochronology, mark−recapture, or length-frequency analysis). RMU: regional management unit The differences in LOOIC between each model i and the best model (Δ LOOIC,i ) were used to compute the model weights w i as: which sum to 1 and indicate the relative support of the data for each model. The best model according to the information criteria was used to estimate the mean value of L ∞ and K for each RMU k as: To predict the mean values for each species, the same equation was used, except substituting the species effect γ s for the RMU effect.

RESULTS
The expected negative correlation between log(K ) and log(L ∞ ) was seen across the 37 von Bertalanffy growth curves we identified from the primary literature (Fig. 1, Table 1). However, individual studies varied from this relationship, especially for green sea turtles.
All 9 meta-analysis models converged adequately according to the Gelman-Rubin diagnostic and the effective number of parameters ( Table 2). The distribution of the residuals was quite similar between the models and appeared to be normally distributed except for some outliers with larger predicted values of both parameters (Figs. S2 & S3). The Pareto k diagnostic was < 0.5 for most data points in all models, but the simpler models performed better ( Table 2, Table S5). The DIC, WAIC, and LOOIC were consistent in finding that the best model was Model D (model weight = 0.71), which was the model with a random effect of species but no effect of RMU, parameter estimation method, sampled body size range, or latitude on either the mean or the error variance (Table 3, Fig. 2). This was the simplest model considered and had the best Pareto diagnostics (Table 2; Table S6). The only other models that were supported by the data according to the information criteria were the one with no fixed effects but including a random effect of RMU (Model E, Table 4 Marquez (1994), which summarized 26 yr of nesting female body size data from the species' primary nesting beach Table 3. Species-specific posterior means, standard deviations (SD), and 95% credible intervals of von Bertalanffy parameters (L ∞ : asymptotic length; K: growth coefficient) for the 5 species of sea turtles from the best model with species only (Model D). L max is the straight carapace length (cm) of the largest nesting female reported across the 15 most recent sea turtle studies for each species (see Table S7)  Fig. 1. Relationship between asymptotic length (L ∞ ) in cm straight carapace length, the Brody growth coefficient (K ), and sample size in each study (n = number of turtles). Abbreviations of regional management units (RMUs) as in Table 1. Data are given in Table 1 L ∞ with distance from the equator, but this effect was small. Model D estimated that the standard deviations among species for both variables and the covariance was roughly similar in magnitude to the error standard deviations and covariance. Model E estimated some variation among RMUs as well, especially for L ∞ (Table 5). From the posterior distributions of the parameters for species from Model D ( Table 3, Fig. 2) and RMUs from Model E (Table 4, Fig. 3), there was not much variation among the RMUs in our analysis for either parameter. For the 4 species that have multiple RMUs (green, loggerhead, hawksbill, Kemp's ridley), the posterior distributions of the RMU param-eters overlapped with each other with only slight apparent differences in the values of K for Kemp's ridley and L ∞ for green and loggerhead sea turtles. On the other hand, the species were quite different from each other for both parameters, as expected.
In general, the posterior distributions showed the expected pattern that well-studied species (green, loggerhead, Kemp's ridley) have more precise estimates of their parameters than the species that have only 1 to 3 studies (olive ridley and hawksbill). The exception to this pattern was the estimate of L ∞ for green sea turtles, which had a relatively large standard deviation ( Fig. 2; Fig. S5). This reflects the large variation in green sea turtle L ∞ values across studies (Fig. 1 Table 4. Population-specific posterior means, standard deviations (SD), and 95% credible intervals of von Bertalanffy parameters (L ∞ : asymptotic length; K: growth coefficient) for all 11 populations of sea turtles examined in this study from the model with populations nested within species (Model E). See Table 1 and Table S5 for regional management unit (RMU) definitions and additional information

DISCUSSION
Through Bayesian hierarchical modeling, we present a meta-analysis of von Bertalanffy growth parameters (L ∞ , K) for 5 chelonid sea turtle species (green, loggerhead, hawksbill, olive ridley, Kemp's ridley) and 11 sea turtle RMUs identified by Wallace et al. (2010). We observed considerable variation in parameter estimates among species, but lower variation as sociated with populations (among RMUs, across latitudes) and sources of uncertainty (parameter esti-mation method, size range sampled), such that the best model according to the information criteria only included a species effect. There was also a substantial error variance between studies even within the same species-specific RMU (Table 5), which contri buted to the wide credible intervals for many species. Because we found no strong differences in parameter estimates between RMUs within each species, the species-specific estimates herein are likely more robust than the RMUspecific estimates, and they are more precise for most species. However, these findings are based on limited data -only 22% (11/49) of RMUs for these species are represented in our dataset and only 10% (5/49) had multiple parameter estimates to synthesize. Therefore, RMU, latitudinal, and methodological effects may emerge as stronger drivers of sea turtle growth variation in future meta-analyses when additional data become available. The dearth of suitable data for this analysis highlights the need for expanded collection, analysis, and synthesis of sea turtle growth rate data globally.

Growth variation among sea turtles
Our meta-analysis identified clear differences in asymptotic length (L ∞ ) and growth (K ) among sea turtle species that largely align with known taxonomic variation. Of the species we investigated, green and loggerhead turtles are among the largest and longest-lived (maturity: > 30 yr, > 90 cm SCL), followed by hawksbills (maturity: >15 yr, > 75 cm SCL) and the ridleys (maturity: >10 yr and > 60 cm SCL; Avens & Snover 2013). Our results largely follow these patterns, with L ∞ decreasing and K increasing from the largest (green) to the smallest (ridleys) sea turtles. Moreover, mean L ∞ estimates are comparable to mean maximum lengths of nesting females (L max ) reported in the literature for all but olive ridleys (Table 3, Table S7). Therefore, our species-specific parameter estimates are ecologically realistic for green, loggerhead, hawksbill, and Kemp's ridley sea turtles. Until additional data become available, research applications for olive ridleys should utilize the von Bertalanffy growth curve reported by Petitet et al. (2015), or consider the sizeat-age smoothing spline growth function reported by Zug et al. (2006). The results of our model are most useful and robust for sea turtle populations in the western North Atlantic Ocean (green, loggerhead, and Kemp's ridley), where most studies included in our analysis were conducted. This large study count explains the higher precision in K estimates for green and loggerhead turtles, and L ∞ for Kemp's ridley and loggerhead turtles. The wider posterior density for K of Kemp's ridley and L ∞ for green turtles can be attributed to the differences between the estimates from the individual studies (Fig. 1,  Table 1). However, due to partial pooling, the information from those studies improves the precision of the estimates for other populations of the same species. Therefore, the species-and RMU-specific para meter estimates presented herein for green, loggerhead, and Kemp's ridleys may be useful priors for growth parameters in demographic models for populations of these species that lack published growth curves.
Although inferences were limited, our analysis identified some differences in species-specific parameter estimates among certain RMUs. For example, even with only 2 studies from the Mediterranean RMU, wider posterior density L ∞ and K for Mediterranean vs. Atlantic Northwest loggerheads suggested potential differences in estimates between these RMUs. Indeed, through skeletochronological and genetic analyses, Piovano et al. (2011) showed that loggerhead turtles of Mediterranean origin likely grow faster and mature at younger ages (and smaller sizes) than turtles of Atlantic origin, possibly due to unique physiological traits that enhance growth rates and lower time to maturity (i.e. higher assimilation efficiency) but reduce ultimate size (i.e. higher somatic maintenance costs, lower cumulative investment to maturation; Marn et al. 2019). Expanded analysis of loggerhead growth data from the Mediterranean, and other RMUs, may yield greater divergence in parameter estimates for these 2 RMUs and a wider range of parameter estimates across all loggerhead RMUs. Our analysis also identified apparent differences in K estimates between Kemp's ridleys that inhabit the Gulf of Mexico vs. the US Atlantic Coast as juveniles. These findings align with the growing body of evidence that Kemp's ridleys that inhabit the US Atlantic grow slower than those that inhabit the Gulf of Mexico (Avens et al. 2017, Ramirez et al. 2020a. Interestingly, despite having the most spatially comprehensive dataset and most variable studyspecific parameter estimates (Fig. 1), our model predicted that green turtle growth parameters vary little among RMUs. These findings may be an artifact of our limited sample size and the effect of partial pooling, which likely drew estimated values towards the species mean and away from the observed data for the 4 (of 5) green turtle RMUs with only a single study. Additionally, green turtle growth parameters may not be as variable as illustrated in Fig. 1. The 2 studies that deviated from the typical L ∞ and K relationship derived von Bertalanffy growth curves from incomplete size ranges (Boulon & Frazer 1990: 25.6− 62.3 cm SCL;Tanaka 2009: < 55 and > 85 cm SCL), which likely influenced their L ∞ and K estimates. Additionally, Goshe et al. (2010), who reported an anomalously high L ∞ estimate of 141.26 cm, found that both the logistic and Gompertz growth functions were better fits to their size-at-age data and yielded asymptotic length estimates closer to the other green turtle studies identified in our meta-analysis (logistic = 104.67 cm, Gompertz = 113.07 cm). This highlights that although the von Bertalanffy growth function  (Frazer & Ehrhart 1985, Bjorndal et al. 2000a, Snover et al. 2007b, different growth functions may prove better fits for certain species and populations. Of course, such variation among green turtle growth parameters may also be due to environmental factors such as temperature and diet (see Section 4.2). Ultimately, the question of whether growth parameters vary among populations cannot be adequately answered based on our analysis, and there is need for more data from most RMUs.

Methodological and ecological drivers of growth variation
Many sources of uncertainty contribute to sea turtle somatic growth variation, including those associated with measuring growth rates (i.e. observation error) and biotic and abiotic processes that lead to variability in growth among individuals. Within our analysis, we specifically evaluated the potential influence of parameter estimation method, body size range, RMU, and latitude on growth parameter estimates.
We found no difference in either means or variances of L ∞ and K estimates between studies that employed different parameter estimation methods, which may be reassuring for those who wish to integrate disparate growth datasets as inputs into population dynamics models. This result is primarily based on inferences from the 3 best-studied species, where there were studies that used 2 (Kemp's ridley) or 3 (green, loggerhead) parameter estimation methods. Importantly, these findings align with those of Chasco et al. (2020), who demonstrated that integrating multiple loggerhead sea turtle growth data types within a single modeling framework may help fill important data gaps in growth curves, reduce biases, and improve parameter estimates. We expect that these findings should apply to the whole Cheloniidae family, since biases specific to any of the methods used to study sea turtle growth rates should be similar across species for SCL measurements. We also found that size range sampled did not influence the results, which was somewhat surprising given the known biases in fitting the von Bertalanffy growth curve to samples missing either large or small animals (e.g. Ortiz de Zárate & Babcock 2016). This result may be an artifact of the fact that only 6 of the 37 included studies sampled a complete size range, so that, if there were biases, they were common across most studies.
The fact that the model including latitude was the second best model in our analysis indicates that temperature/latitudinal effects on growth may be important. Although sea turtle growth rates are expected to be strongly related to temperature/latitude (Wallace & Jones 2008), in practice, establishing direct links between temperature and growth has been challenging. For example, Bjorndal et al. (2013Bjorndal et al. ( , 2016Bjorndal et al. ( , 2017 identified latitudinal or climate-driven variation in green, loggerhead, and hawksbill sea turtle growth rates in the Caribbean, but similar effects were not observed for green sea turtles across Australia and Hawaii (Balazs & Chaloupka 2004 or Kemp's ridley sea turtles in the western North Atlantic (Avens et al. 2020, Ramirez et al. 2020a. Within our analysis, we observed weakly declining trends in K with increasing latitude for green and Kemp's ridley turtles (Fig. 4). However, limited spatial data within and among all sea turtle species contributed to high uncertainty in these relationships. Furthermore, marine ecosystems have experienced unprecedented, spatiotemporally variable change over the past century due to climate change, fisheries, and many anthropogenic stressors (Byrnes et al. 2007, Halpern et al. 2008, Rocha et al. 2015, Beyer et al. 2016, factors that may directly and indirectly influence sea turtle growth rates. We necessarily integrated growth curves based on data collected over multiple decades (publication dates: 1985− 2018). However, extensions of our model, when additional data become available, to include year effects or other environmental parameters will allow for greater analysis of correlations between growth rates and variable environmental stressors.
Diet variation could not be evaluated with our model but may also contribute to variation in L ∞ and K estimates within species, particularly for green sea turtles. Sea turtles display a wide range of foraging strategies both within and among species (Bjorndal 1997). Although loggerhead and Kemp's ridleys are generalist carnivores that forage primarily on crabs (but also other invertebrates and sometimes fish), individuals can have very narrow niches or unique foraging behaviors within populations (Hatase et al. 2002, McClellan & Read 2007, Vander Zanden et al. 2010, Pajuelo et al. 2016, Ramirez et al. 2020b. Moreover, there is a growing body of evidence that green sea turtles, classically considered herbivores, may be more omnivorous than previously thought. Although the consumption of seagrass and macroalgae remains widespread, consumption of animal protein has also been observed in benthic green turtles across the globe (Mauritania: Cardona et al. 2009; Piovano et al. 2020). While our modeling results suggested low variation in L ∞ and K estimates among RMUs, dietary plasticity could contribute to variation in green turtle L ∞ and K estimates among individual growth studies (Fig. 1).
More sea turtle growth studies are needed, especially syntheses of existing length-at-age or lengthfrequency data, in order to close data gaps and fully investigate sea turtle growth variation. Our analysis was conducted on a very narrow range of sea turtle growth studies (34/195) and RMUs (green: 5/17 RMUs; loggerhead: 2/10; hawksbill: 2/13; olive ridley: 1/8; Kemp's ridley: 1/1). However, there are considerable additional data that should be collated through region-wide collaborations among data holders (e.g. Bjorndal et al. 2017), such as those from the 61 studies that used Gompertz, logistic, or non-parametric growth functions or the 75 studies that presented annual growth rate data but not growth curves. Additionally, the Mediterranean (Greece, Turkey, Cyprus), Southwest Atlantic (Brazil), North Pacific (Japan, Hawaii), and South Pacific (Australia) likely have significant additional data that could be used to fit additional von Bertalanffy growth curves, as they include some of the longest-running sea turtle research and conservation programs in the world (Limpus et al. 2003, Margaritoulis 2005, Marcovaldi & Chaloupka 2007, Balazs et al. 2015. Fitting multiple new datasets together within a meta-analytic framework would allow for greater investigation and mitigation of biases caused by differences in study sample size and body size range (e.g. Chasco et al. 2020). We also encourage future analyses to evaluate and report fits of multiple model forms (e.g. Gompertz, logistic, nonparametric) to allow for better identification of the somatic growth function most appropriate for each sea turtle species and population and comparisons therein (e.g. Parham & Zug 1997, Goshe et al. 2010, Chasco et al. 2020).

Conclusion
Using a meta-analytic approach, we integrated existing von Bertalanffy growth data for 5 species of chelonid sea turtle to provide updated estimates of K and L ∞ through standard Bayesian hierarchical modeling methods. The estimated species-specific mean values of the growth parameters from this study (Table 3) would be appropriate priors for future growth studies and for growth parameters in population dynamics models, particularly for status assessments of populations that lack regional somatic growth curves. However, all available ecological information, including genetic, environmental, and ecological, should be used to assess the suitability of extending RMUspecific estimates to other populations. Data from the western North Atlantic Ocean dominate our dataset and thus provide the most robust estimates. Considerable effort is needed to expand growth studies to additional sea turtle populations through syntheses of existing data or initiation of new research projects. Our species-specific results will possibly be the most valuable to natural resource managers, sea turtle biologists, and modelers in regions lacking robust growth datasets, as they provide an important starting point for evaluating size-at-age relationships, constructing maturation schedules, and developing stage-structured population models.
Acknowledgements. This project was funded by the NOAA Educational Partnership Program through the Living Marine Resources Cooperative Science Center (NOAA Award No. NA16SEC4810007). Its contents are solely the responsibility of the award recipient and do not necessarily represent the official views of the US Department of Commerce, National Oceanic and Atmospheric Administration. We thank S. Heppell and 3 anonymous reviewers for their help in improving this manuscript.