Somatic Growth and Maturity for Four Species of River Cooter Including Pseudemys concinna suwanniensis, P. nelsoni, P. peninsularis, and P. texana

Simple Summary River Cooters commonly occur in many freshwater ecosystems, but little is known about their growth and maturity. We used nearly 20 years of turtle tracking to estimate the growth, length at maturity, age at maturity, and sexual dimorphism of four River Cooters in five freshwater springs, four in Florida and one in Texas. We found that River Cooter juveniles grow exceptionally fast, doubling or tripling in size in the first year, while adults grow very slowly after maturing. Female River Cooters grow larger than males and mature later for the four taxa we studied but much more so for Suwannee River Cooters. Environmental differences among springs affected the growth rate of turtles more than the maximum size. Our efforts elucidate more about the life history of River Cooters and define baseline information to track potential changes to growth from natural and human disturbances in the future. Abstract Pseudemys is a genus of commonly occurring freshwater turtles with limited growth information across their long lifespans. We used 11,361 mark-recapture events to estimate the somatic growth rates of P. nelsoni, P. peninsularis, P. concinna suwanniensis, and P. texana from freshwater springs and developed a Bayesian growth model to estimate the species-specific, site-specific, and individual effects on growth. We corroborated evidence for fast juvenile growth and slower adult growth in Pseudemys but found uncommonly fast growth rates, with turtles doubling or tripling in size in the first year. P. texana males had the smallest average maximum size (L∞, 243 mm), while P. c. suwanniensis females had the largest (423 mm). Environmental conditions at springs had significant effects on k, the growth coefficient, but not L∞. We derived, using a ratio of length at maturity to L∞ (71.7% and 87%, males and females), that females matured 1.15–1.57 times older than males except for P. c. suwanniensis, which matured three times older. Given the local abundance declines in many Pseudemys from anthropogenic impacts, this study provides important baseline life history information for Pseudemys species for use in ongoing conservation efforts and presents a novel hierarchical modeling approach using a long-term mark-recapture dataset.


Introduction
The genus Pseudemys (River Cooter) is a group of freshwater turtles (family Emydidae) distributed throughout the south and eastern United States and northern Mexico [1]. The genus has long been considered a taxonomic "quagmire" and has been the subject of numerous revisions [2]. Currently, there are nine recognized species and subspecies with most species, when found, being locally abundant [3][4][5]. Despite this, little is known about the growth and, subsequently, age at maturity for most of the Pseudemys species [3][4][5][6]. With observed and perceived declines in the local abundance of many Pseudemys species due to anthropogenic threats [7,8], this critical information on the species' growth and maturity is imperative for understanding their resilience to anthropogenic threats [9,10]. Here, we focus on providing growth rates for three Pseudemys species and one subspecies occurring in Florida and Texas: Florida Red-bellied Cooter (P. nelsoni), Suwannee River Cooter (P. concinna suwanniensis), Peninsula Cooter (P. peninsularis), and the Texas River Cooter (P. texana).
The Pseudemys genus can be found in a variety of habitats, including large rivers, smaller tributaries, lakes, manmade water features, and freshwater spring habitats [3][4][5]11,12]. In Florida, both the Florida Red-Bellied Cooter and Peninsula Cooter prefer habitats with slow water flow with abundant vegetation [3,5], with Florida Red-Bellied Cooter preferring faster water flows than Peninsula Cooter [12]. The Texas River Cooter and Suwannee River Cooter are more adapted to higher flow systems [4,13,14]. All four species/subspecies occur in freshwater springs and spring run systems [12,[14][15][16][17], where higher flows occur in the spring runs, and slower flow areas occur in the spring boils, lagoons, and associated sloughs.
Freshwater springs are excellent for comparing growth across the four Pseudemys species/subspecies given that all of these taxa occur in high densities in spring systems [12,14,16,18]. The springs of Florida and Texas are largely limestone-based, which dampens spring water temperature variability resulting in a near-constant temperature reflecting the average annual air temperature. This constant temperature and water flow allows for a year-round growing season for aquatic plants, which provides a constant food source for species of turtles that are strongly herbivorous, like Pseudemys [3][4][5]. However, freshwater springs are under increased threat by local to global anthropogenic stressors, including habitat alteration, abundant recreational use, groundwater and recharge zone depletion, pollution, invasive species, and climate change [12,19,20]. Additionally, the water quality and quantity of many springs have decreased drastically over the past several decades, resulting in environmental contrast across spring systems despite their thermal stability [19,20]. Given that these systems are bastions of high local abundance for Pseudemys, it is important to establish baseline growth rates in spring systems to monitor potential changes from continued anthropogenic stressors.
A well-documented Pseudemys trend has been fast juvenile growth with slower growth (typically < 5 mm/year) following the onset of maturity [3][4][5][6]16]. Only P. texana has a formal growth model using length-at-age data, where age was estimated from scute growth annuli [6]. The lack of growth modeling is directly related to the ambiguity in the growth annuli being deposited on an annual schedule [21] and a loss of the growth annuli beyond a few years [4]. As a result, there is high uncertainty in key demographic parameters such as age at maturity and longevity. Generally, females have been noted to grow larger and mature later than males, with Suwannee River Cooter females documented to grow 150% larger [4].
With limited length-at-age data, mark-recapture-based models can be a suitable replacement for estimating growth [22,23]. These models are often based on the von Bertalanffy growth model that represents individuals growing towards an asymptotic size [24] and is typically used for taxa exhibiting indeterminate growth. This model assumes all individuals grow identically [22], but this assumption has been relaxed with many modifications [25], most assuming individual variation occurs independently in the asymptotic size, L ∞ (see [26] for a review). Another limitation that has been noted is that the von Bertalanffy  [27,28]. In teleosts, this limitation has spurred a series of biphasic growth models generally aimed at addressing changes in the growth rate due to maturation [28]. In testudines, Armstrong and Brooks [23,29] modified the von Bertalanffy to change the Brody growth coefficient, k, at a fixed age at maturation. However, the age at maturation is largely unknown for Pseudemys species resulting in a critical assumption and easily resulting in bias if the age at maturity is misspecified or if variation exists between groups (e.g., sex, location, subspecies).
Here, we make use of long-term, ongoing, mark-recapture efforts at multiple sites of Florida Red-Bellied Cooters, Suwannee River Cooters, Peninsula Cooters, and Texas River Cooters from five freshwater springs systems, four in Florida and one in Texas, to model somatic growth for these taxa. We compare our estimates to the existing sparse literature on Pseudemys growth, sexual dimorphism, and the first-year growth increment. We also produce estimates of age at maturity to compare to existing estimates based on limited growth annuli studies. To do such, we make use of life-history invariants that relate the time necessary to reach some proportion of L ∞ [30][31][32]. These approaches have been used to derive life history parameters from von Bertalanffy growth extensively in teleosts [33] but recently in data-deficient elasmobranchs [34,35]. For three of the four species, our estimation of growth is the first for the species and represents a major contribution to understanding the life history of these taxa.

Sampling Sites
Wekiwa Springs in Orange County, Florida ( Figure 1) is a second magnitude spring that opens up to a large slow-flowing pool (149 × 98 × 6.5 m, length x width x height) then narrows (6-12 m wide) becoming the 1.3 km long Wekiwa Springs Run [36]. Invasive aquatic weeds, principally hydrilla, have been an issue in the past, while algal bloom issues persist [12,16,37]. Volusia Blue Spring in Volusia County, Florida ( Figure 1) is a first magnitude spring with a 0.72 km long spring run (30 to 46 m wide) which has been impacted anthropogenically over time through pollution, invasive species, recreational use, and algal growth. The spring acts as the largest manatee (Trichechus manatus latirostris) winter refuge on the St. Johns River, and this has left the spring largely denuded of native vegetation [17]. Manatee Springs in Levy County, Florida ( Figure 1) is a first-magnitude spring with a public swim area dominated by green algae that becomes a hardwood swamp spring run before eventually emptying into the Suwannee River. The study site also includes several sinkhole ponds, which are connected to the main spring. Fanning Springs (~0.70 ha), located in Levy County, Florida ( Figure 1) has historically been considered a first-magnitude spring, but is now a second-magnitude spring. Fanning Springs is the only site in the current study with a presence of motorized boats. The spring area is largely devoid of vegetation, while the spring run is a hardwood swamp. Comal Springs system feeds directly into Landa Lake, located in New Braunfels, Comal County, Texas ( Figure 1). American eelgrass dominates large portions of the lake. Depending on water levels, the dense vegetation and algal blooms make parts of the lake become almost impossible to survey [18]. Additional site covariates, including size, discharge, temperature, nitrate levels, vegetation density, algal bloom frequency, and annual visitors, are provided in Table 1.

Sampling Protocol
At each of the five study sites, researchers (a mix of professional biologists, students, and interested citizen scientists) conducted multi-day annual or semi-annual snorkel surveys between 1999 and 2020 led by the Turtle Survival Alliance North American Freshwater Turtle Research Group [38]. For each sampling period, a variable number of snorkelers (typically between 5 and 30, depending on the study site) hand-captured turtles from ca. 0800-1600/1900 h, depending on the time of year and weather conditions. All turtles observed were captured by hand while snorkeling and placed into canoes and brought back to a central location for processing. All turtles were tagged and measured (details below) and then released. For each turtle, we recorded maximum straight-line measurements of carapace length (CL), plastron length, carapace width, and shell height to the nearest 1 mm. We determined the sex of turtles based on secondary sexual characteristics, notably tail length and girth and placement of the cloaca in relationship to the edge of the carapace, according to Ernst and Lovich [39]. We noted any unique features or physical anomalies, such as damage, scars, or coloration for each turtle, which helped confirm individual identity. All turtles were weighed either using hanging Pesola spring

Sampling Protocol
At each of the five study sites, researchers (a mix of professional biologists, students, and interested citizen scientists) conducted multi-day annual or semi-annual snorkel surveys between 1999 and 2020 led by the Turtle Survival Alliance North American Freshwater Turtle Research Group [38]. For each sampling period, a variable number of snorkelers (typically between 5 and 30, depending on the study site) hand-captured turtles from ca. 0800-1600/1900 h, depending on the time of year and weather conditions. All turtles observed were captured by hand while snorkeling and placed into canoes and brought back to a central location for processing. All turtles were tagged and measured (details below) and then released. For each turtle, we recorded maximum straight-line measurements of carapace length (CL), plastron length, carapace width, and shell height to the nearest 1 mm. We determined the sex of turtles based on secondary sexual characteristics, notably tail length and girth and placement of the cloaca in relationship to the edge of the carapace, according to Ernst and Lovich [39]. We noted any unique features or physical anomalies, such as damage, scars, or coloration for each turtle, which helped confirm individual identity. All turtles were weighed either using hanging Pesola spring scales (Pesola AG, Baar, Switzerland) or Ohaus top-loading digital scales (Ohaus Corporation, Parsippany, NJ, USA), depending on turtle size. All turtles were marked using a variation of the technique described by Cagle [40]. In 2009, passive integrated transponder (PIT) tags were used for turtles with CL greater than 70 mm and injected under the right bridge of the turtle into the inguinal cavity [41].

Modeling Growth
The lack of accurate aging structures late into life and the lack of maturation information severely limit the application of simpler modeling approaches to Pseudemys. This limitation is due to the nature of mark-recapture growth information, which incorporates individual variation. Without accounting for this variability, growth model parameters have increased uncertainty, and the estimation of population-level growth variation is inflated. Below, we present a modification of the von Bertalanffy growth model to parse the population-level and the individual-level effects on somatic growth. In doing so, we account for the individual variation by using a hierarchical multivariate approach to draw a L ∞ and k for each individual. This approach retains strong correlations between L ∞ and k [34,42] that previous individual variability approaches focused solely on L ∞ have disregarded. Additionally, we pool the mark-recapture datasets across species and sites, which facilitates sharing environmental effects on L ∞ and k between similar species monitored across a range of environmental conditions. These density-independent factors are expected to impact somatic growth by changing metabolic rates [43,44]. Pooling of the data facilitates species-specific effects to set the baseline site-specific growth while allowing shared environmental effects on L ∞ and k, which is not an unreasonable assumption for this set of species/subspecies given the similarity in habitats they occupy. Without pooling species across sites into one model, it would not be possible to separate out site-specific effects from environmentally driven effects.

Mark-Recapture von Bertalanffy Growth Model
We wished to understand the sex-specific growth of Pseudemys spp. across our study sites while accounting for various site covariates that might impact growth. We solely used captures that resulted in zero or positive growth and of individuals that eventually had a sex determined at some point in their capture history. We modified the Fabens [22] implementation of the von Bertalanffy growth model [24] commonly applied to testudines growth modeling [23,29] (Equation (1)).
whereL t is the expected length at capture in time t for individual i, L ∞ is the asymptotic length (also the average maximum length), L t−1 is the length in the past capture event, k is the Brody growth coefficient (the proportion of L ∞ − L t obtained per unit time), and ∆t is the amount of time between t and t − 1 in fractions of a year. We assumed a normal likelihood as has been done in previous testudines growth studies [23,29] (Equation (2)).
where L t,i is the observed length at capture in time t for individual i and σ is the variability in growth around the expected length. In order to estimate separate growth parameters across sexes and species/subspecies, we estimated sex-and species-specific growth parameters, L ∞ , k, and σ (Equation (3)). We also wished to incorporate individual variability in growth and did so by incorporating a random effect for individual Equations (3)- (5).
Biology 2023, 12, 965 where s denotes species, j denotes sex, i denotes individual, and β i is the random effect for L ∞ and k. Unfortunately, the nonlinearity of the von Bertalanffy growth model results in growth parameters that are highly correlated [34,42]. We chose to directly estimate this correlation as part of the growth model and structured the individual growth parameter random effects as a set drawn from a multivariate normal distribution (Equation (6)).
where the mean of the individual random effects is zero for each transformed growth parameter, and Σ is the covariance matrix of the multivariate normal distribution. For the sake of easing the fitting of this correlation, we decomposed the covariance matrix into the Cholesky factorization (Equations (7) and (8)).
where τ is the scale of the individual random effects, Ω is the correlation matrix between log(L ∞ ) and logit(k), and W Ω W Ω is the Cholesky factorization of the correlation matrix (it is worth noting that W Ω is typically written as L Ω , which we replaced to reduce confusion with L t ) [45].

Incorporating Site Covariates
We also wished to incorporate site covariates into the von Bertalanffy growth model. We did this by modifying Equations (4) and (5) to include the part-worth effects of the site covariates (Equations (9) and (10)).
where Γ is the loading matrix that loads site characteristics, X g , at site g on the respective transformed growth parameter. The site covariates we included were site size (hectares), flow (liters per day), temperature ( • C), nitrate (mg/L), a categorical variable for vegetation density (low, moderate, high), a binary variable for bloom frequency (low, high), and the number of visitors per year (Table 1). We simplified this model such that the individual growth parameters incorporated the mean of the transformed growth parameters for each sex, the part-worth site effects, and the individual random effects (Equation (11)).
In this formulation, β i,{L ∞ ,k} vary independently, but the diag(τ)W Ω diag(τ) modifies this independence based on the correlation between the transformed L ∞ and k. This is equivalent to hierarchical prior on β i,{L ∞ ,k} that occurs in the multivariate formulation (Equation (6)) but allows vectorization speeding up model fitting. We assumed that the correlation between L ∞ and k was shared across species inducing a hierarchical prior on the individual random effects across Pseudemys spp. We made this assumption for a practical reason-some species were more heavily sampled across sites than others-and under the broader assumption that closely related species should share similar growth rates [46].

Predicted Length at Age and Derived Age at Maturity
Within the model, we desired to estimate the expected length at age and had to assume a size at hatch (birth), L 0 . We collated studies with hatchling sizes to obtain L 0 for our study's species, including the mean and standard deviation in L 0 . Several studies were lacking standard errors, and we estimated these standard errors given the mean, range, and assumed quantile of the range (coming from a shared significance level, α). To do such, we built an optimization function that minimized the difference between the probabilities from a normal probability distribution function of the known range and the assumed quantile. We first used this optimization function to estimate a typical significance level, α, for standard errors of hatchling size using a known mean, range, and standard error from the Lindeman [6] study on P. texana. This significance level was then used to predict the standard errors for P. concinna and P. nelsoni from Jackson [47] and Heinrich et al. [48], respectively. The Aresco [49] data on P. floridana was used to inform P. peninsularis ( Table 2). The mean, µ L 0 , and standard error, σ L 0 , of the size at hatch for each species was used to inform the prior (and posterior) distribution on L 0 (Equation (12)).
and this L 0 was used in the original formulation of the von Bertalanffy growth model [24] to predict the length at age (Equation (13)).
where Equation (13) is analogous to Equation (3) but uses t instead of ∆t where t is equal to a vector of ages to predict length at age at, ranging from age zero to 50 and an interval of a tenth of a year. Table 2. The mean (µ), standard deviation (σ), and range of the size at birth (L 0 ) for each Pseudemys species and the source. Standard deviations noted with asterisks were estimated given the species' mean and range, and assuming an α equal to that would result in the range for P. texana given its mean and standard deviation. We also wished to predict length (L mat ), age at maturity (t mat ), and age at capture t cap . Life-history invariants [30,32,50] have been used to estimate maturity parameters from the von Bertalanffy growth model by assuming that maturity occurs at some proportion of L ∞ [34,35]. Gibbons et al. [51] estimated this proportion for males and females of Trachemys scripta for two populations, and we used the average sex-specific proportions of L ∞ , L mat L ∞ , for Pseudemys (Equation (14)).

Species
We did not find a similar dataset for any Pseudemys species, and with the genus Trachemys sister to Pseudemys [52], it stands to reason that the two genera might mature at similar proportions of L ∞ . With this length at maturity or the length at capture, we can then derive the expected age at the respective length (Equation (15)).
where t x is the age at maturity (t mat ) or the age at capture t cap , and L x is the corresponding length. By predicting length at age and estimating a derived L mat , t mat , and t cap within the model, we ensure that all error in the model estimates of the growth parameters propagates to these predictions/derivations [34,35].

Bayesian Implementation
We fit the mark-capture von Bertalanffy growth model [22] using STAN and the cmdstanr package [53]. We specified priors on L ∞ , k, and σ by estimating Equation (1) for each species and sex using maximum likelihood estimation to estimate the location and setting a coefficient of variation of 50% to determine the respective scale following Caltabellotta et al. [54] and Rolim et al. [34]. We assumed that β i,{L ∞ ,k} and Γ in Equations (6)-(11) followed standard normal priors. We assumed W Ω had a prior following the Lewandowski, Kurowicka, and Joe (LKJ) [45] correlation matrix Cholesky factor with a shape of 2 and set a prior on τ following [55]. We fit the model in STAN using four chains with 5000 warmup and 1000 sampling iterations per chain using the NUTS sampler. We assessed chain convergence using the Gelman-Rubin statistic [56].  Figure 1A). Comal Springs had the highest average discovery rate of new individuals (7.2 individuals/trip), with similar rates for Wekiwa Springs, Fanning Springs, and Manatee Springs (3.05-4.09 individuals/trip) and the lowest rates for Blue Springs (0.81-1.95 individuals/trip). We also sampled turtles above the noted maximum size for several species: 330 mm versus 300 mm in Jackson [47] for P. nelsoni males, 342 mm versus 320 mm in Thomas and Jansen [5] for P. peninsularis males, 349 mm versus 330 mm in Jackson [4] for P. c. suwanniensis males, and 344 mm and 294 mm versus 241 mm and 161 mm in Lindeman [6] for P. texana females and males, respectively.

Growth
The mark-recapture von Bertalanffy growth model of species-and site-specific growth rates converged with the Gelman-Rubin statistic below 1.1 for all parameters after 6000 iterations per chain. The median correlation between log-transformed L ∞ and logistic-transformed k individual random effects was −0.31 (−0.4−0.23, 90% credible interval).

Species-Specific Growth Rates
Among the four species modeled, P. nelsoni and P. texana were the most similar in median female asymptotic size, with L ∞ at 311 and 313 mm, respectively, while P. penin-sularis and P. c. suwanniensis females grew to larger sizes, 363 and 423 mm, respectively (Table 3, Figure 2). The asymptotic size of males did not follow this trend as P. peninsularis and P. c. suwanniensis L ∞ were very similar (309 and 311 mm, respectively) while P. nelsoni males had L ∞ of 274 mm and P. texana males, the smallest of the four species, had L ∞ of 243 mm (Table 3, Figure 2). The Brody growth coefficient, k, for P. c. suwanniensis males was almost double females and the only species where k M > k F (Table 3, Figure 2). All standard deviations in growth, σ, were similar between males and females (Table 3) except for P. c. suwanniensis whose females had the highest variability ( Figure 2). The predicted growth increment in the first year, ∆L a=1 , for females was highest in P. peninsularis (56.9 mm) and similar for P. texana (55.2 mm), somewhat lower for P. c. suwanniensis (49.5 mm) and lowest in P. nelsoni (41 mm). In contrast, the male predicted growth increment in the first year was far greater for P. c. suwanniensis (61 mm) and similar among the other three species (31.3-36 mm). Table 3. von Bertalanffy growth parameters (L ∞ , k, σ, t 0 ) for the four Pseudemys taxa in this study as well as the first-year growth increment (∆L a=1 ) and the age at maturity (t mat ). The median and 90% credible interval (in parentheses) are reported for each parameter. Note that t mat was derived from the life-history invariants and was not directly informed by maturity measurements. All standard deviations in growth, , were similar between males and females (Table 3) except for P. c. suwanniensis whose females had the highest variability ( Figure 2). The predicted growth increment in the first year, Δ , for females was highest in P. peninsularis (56.9 mm) and similar for P. texana (55.2 mm), somewhat lower for P. c. suwanniensis (49.5 mm) and lowest in P. nelsoni (41 mm). In contrast, the male predicted growth increment in the first year was far greater for P. c. suwanniensis (61 mm) and similar among the other three species (31.3-36 mm).

Species-Specific Derived Age at Maturity and Longevity
From Gibbons et al. [51], the ratio of L mat to L ∞ ranged from 63.5% to 93%, with a mean ratio of 71.7% for males and 87.2% for females. The result of this ratio is that females mature larger and later than males for all species (Figure 2). Across the study taxa, P. c. suwanniensis had the biggest difference between female and male length at maturity (∆L mat , 145 mm SCL), followed by P. texana (99 mm SCL), then P. peninsularis (94 mm SCL), and then P. nelsoni (74 mm SCL) ( Table 3). For P. peninsularis, the difference between female and male median age at maturity, t mat , was small (∆t mat , 1.15 years) but was greater for P. nelsoni and P. texana (3.88 and 2.96 years, respectively) and greatest for P. c. suwanniensis (9.79 years). For all species except P. c. suwanniensis, this results in females maturing 1.13-1.57 times older than males while, for P. c. suwanniensis, females matured 3.25 times older. The oldest median female age at maturity was P. c. suwanniensis, 14.12 years, while the youngest was P. texana, 8.14 years. The oldest median male age at maturity was P. peninsularis, 9.07 years, while the youngest was P. c. suwanniensis, 4.36 years (Table 3).

Site-Specific Effects on Growth
None of the site covariates had significant effects, but, when aggregated to the site level, the totality of the environmental conditions had significant effects on k but not L ∞ (Figure 3). Among the sites, Fanning Springs had the greatest impacts on L ∞ and k, with negative impacts on L ∞ and significant positive impacts on k. The remaining sites had negative impacts on k, with all but Manatee Springs having significant effects. Comal Springs and Blue Spring had positive median effects on L ∞ , while the remaining sites, Wekiwa Springs and Manatee Springs, had slightly negative median effects (Figure 3). The result of these effects on L ∞ and k on P. nelsoni and P. peninsularis was that turtles at Blue Spring had slightly larger asymptotic sizes than Wekiwa Springs and lower k values leading to slightly lower length at maturity in Wekiwa Springs ( Figure 4A-D). For P. c. suwanniensis, turtles at Manatee Springs had slightly larger asymptotic sizes than Fanning Springs and lower k values leading to slightly lower length at maturity in Fanning Springs ( Figure 4E,F). The result is that relative to the species-specific growth rate, the environmentally informed L ∞ increased slightly on average, and k decreased on average ( Figure 4G). (Table 3).

Site-Specific Effects on Growth
None of the site covariates had significant effects, but, when aggregated to level, the totality of the environmental conditions had significant effects on bu (Figure 3). Among the sites, Fanning Springs had the greatest impacts on and negative impacts on and significant positive impacts on . The remaining s negative impacts on , with all but Manatee Springs having significant effects Springs and Blue Spring had positive median effects on , while the remaini Wekiwa Springs and Manatee Springs, had slightly negative median effects (Fi The result of these effects on and on P. nelsoni and P. peninsularis was tha at Blue Spring had slightly larger asymptotic sizes than Wekiwa Springs and lowe ues leading to slightly lower length at maturity in Wekiwa Springs ( Figure 4A-D c. suwanniensis, turtles at Manatee Springs had slightly larger asymptotic sizes th ning Springs and lower values leading to slightly lower length at maturity in Springs ( Figure 4E,F). The result is that relative to the species-specific growth environmentally informed increased slightly on average, and decreased age ( Figure 4G).

Sexual Dimorphism
Across the suite of , , and parameters, P. c. suwanniensis and P. peninsularis had the least overlap between female and male posteriors (2% and 6% on average, respectively), with P. nelsoni and P. texana having more overlap (24% and 21% on average, respectively) (Figure 2). In the first year of life (Δ ), P. peninsularis was the most sexually dimorphic (20% posterior overlap) in SCL, while P. c. suwanniensis was the least sexually dimorphic (68%). However, P. c. suwanniensis was the most sexually dimorphic (4% posterior overlap) in straight carapace length at their maximum size ( ). These results held between species across sites, but the environmental effects on growth resulted in some within-species site-specific differences in sexual dimorphism. For P. nelsoni and P. peninsularis, sexual dimorphism in was greater at Wekiwa Springs (0% posterior overlap for both species) than Blue Spring (47% and 36%, respectively) but similar for and Δ . For P. c. suwanniensis, Fanning Springs exhibited more sexual dimorphism in (0% posterior overlap) than Manatee Springs (10%) but no differences in and Δ .

Discussion
Somatic growth is a key life history parameter that correlates with a variety of physiological processes, including maturity, fecundity, and mortality [31,57]. As such, somatic growth is often used as a proxy or surrogate to inform these, at times, more difficult-to- . Predicted length at age for female (purple) and male (green) Pseudemys species sampled at five spring systems. For each species and site, the median (solid line) and 90% credible interval of the prediction (shaded region) are shown from zero to fifty years and are overlaid on top of the back-calculated time series of length at age for each individual's capture history. On the right side of each plot is the posterior distribution of the species and site-specific asymptotic length (L ∞ , dark shaded region) and length at maturity (L mat , light shaded region). Note that length at maturity is derived from life history invariants (see Section 2.3.3).

Sexual Dimorphism
Across the suite of L ∞ , k, and σ parameters, P. c. suwanniensis and P. peninsularis had the least overlap between female and male posteriors (2% and 6% on average, respectively), with P. nelsoni and P. texana having more overlap (24% and 21% on average, respectively) ( Figure 2). In the first year of life (∆L a=1 ), P. peninsularis was the most sexually dimorphic (20% posterior overlap) in SCL, while P. c. suwanniensis was the least sexually dimorphic (68%). However, P. c. suwanniensis was the most sexually dimorphic (4% posterior overlap) in straight carapace length at their maximum size (L ∞ ). These results held between species across sites, but the environmental effects on growth resulted in some withinspecies site-specific differences in sexual dimorphism. For P. nelsoni and P. peninsularis, sexual dimorphism in L ∞ was greater at Wekiwa Springs (0% posterior overlap for both species) than Blue Spring (47% and 36%, respectively) but similar for k and ∆L a=1 . For P. c. suwanniensis, Fanning Springs exhibited more sexual dimorphism in L ∞ (0% posterior overlap) than Manatee Springs (10%) but no differences in k and ∆L a=1 .

Discussion
Somatic growth is a key life history parameter that correlates with a variety of physiological processes, including maturity, fecundity, and mortality [31,57]. As such, somatic growth is often used as a proxy or surrogate to inform these, at times, more difficult-to-measure processes resulting in somatic growth frequently being critical for understanding a species response to disturbance, natural or anthropogenic [9,10]. Testudines, and more specifically, freshwater turtles, are a group of species whose general life history characteristics, slow-growing, long-lived, late maturing, and low fecundity, put them at risk from human impacts [29,58]. The habitats these species live in are also frequently disturbed, degraded, or destroyed by human activity [59], elevating the need for measuring a species compensatory ability to such impacts. Here, we present somatic growth relationships from a large, multi-site mark-recapture survey of three Pseudemys species (P. nelsoni, P. peninsularis, and P. texana) and one subspecies (P. concinna suwanniensis). We corroborate the broad finding in Pseudemys of rapid juvenile growth and very slow adult growth that begins with the onset of maturity but find exceptionally high juvenile growth rates. The hierarchical Bayesian von Bertalanffy growth model we developed was essential for estimating speciesspecific growth rates from individual capture histories while accounting for environmental effects on asymptotic length (L ∞ ) and growth coefficients (k) across sites. This method also allowed us to share the uncertainty in the growth parameters with the estimates of length at maturity and age at maturity we derived using life history invariants. Together, the somatic growth rates and maturity information represent a substantial information gain on the life history of these species.
The somatic growth patterns we estimated broadly corroborate many of the growth findings for past studies on Emydidae species [6,11,16]. For P. nelsoni, P. c. suwanniensis, and P. peninsularis, our somatic growth relationships are the first estimated for the species/subspecies, while for P. texana they are the first from mark-recapture surveys. Lindeman [6] estimated a von Bertalanffy growth model for P. texana using growth annuli on 26 turtles that were less than six years old, and estimated k M = 0.191, k F = 0.129 and L ∞,M = 169, L ∞,F = 255 mm SCL, where M and F indicate males and females. We estimated higher k and L ∞ values for P. texana based on 1769 captures of 914 females and 855 males with k M = 0.201, k F =0.232 and L ∞,M = 243, L ∞,F = 313 mm SCL. There is a chance that the differences in k between the studies are from site-specific differences, as Lindeman [6] sampled the South Llano River, whereas our mark-capture study occurred in the Comal Spring-Landa Lake system. A co-occurring explanation is the Lindeman [6] study lacked older turtles, with their largest female being 274 mm SCL and ours being 344 mm SCL. Poor sampling of larger, older individuals results in lower k as L ∞ is overestimated as there is little data to inform the asymptotic size [42]. For P. concinna concinna, Dreslik [60] estimated k M = 0.136, k F = 0.087 and L ∞,M = 239, L ∞,F = 327 mm SCL from 15 male and 34 female turtles in Round Pond, Illinois. Similar to the findings for P. texana, our estimated P. c. suwanniensis growth parameters were considerably higher (k M = 0.258, k F = 0.138 and L ∞,M = 311, L ∞,F = 423 mm SCL).
As such, it is perhaps coincident that the first-year growth rates we estimate are exceptional, with females of P. peninsularis and P. texana as well as males and females of P. c. suwanniensis increasing their length more than 49 mm SCL in the first year from hatching size. These rates are high compared to previously published juvenile rates, with Bancroft et al. (1983) [62] estimating 13 mm yr −1 for P. peninsularis. For female P. texana and male P. c. suwanniensis, our growth rate estimates result in a one-year-old turtle growing to over 100 mm SCL and already close to a third of its asymptotic size. P. c. suwanniensis standouts with the highest growth coefficient among the study species with one-year-old males predicted to grow almost 12 mm SCL larger than females in the first year, the only species in our study to do so. Males in the other Pseudemys species/subspecies were typically 58-76% of female size in the first year of life. For P. nelsoni and P. texana, males and females had similar growth curves with smaller differences between the sex-specific k and L ∞ , while P. peninsularis males exhibited the slowest growth among the males and females of all the species/subspecies in our study. These elevated growth rates are potentially due to the spring systems we sampled. Springs provide a thermally stable environment reducing the effects of seasonality on primary production and providing turtles with a consistent food supply year-round. Huestis and Meylan [61] found elevated average growth rates (35 mm yr −1 ) for P. c. suwanniensis in Rainbow Run, another spring system in Florida.
We found strong evidence for sexual dimorphism among the species/subspecies, with females larger than males, similar to past studies. Sexual dimorphism was the most extreme in P. c. suwanniensis, with the L ∞ of females 1.36 times greater than that of males, which is similar but lower than the Jackson [4] estimate of 1.5 times. Interestingly, the smallest species among the set, P. texana, had the next highest sexual dimorphism with a ratio of L ∞,F /L ∞,M of 1.29. Both P. nelsoni and P. peninsularis had similar ratios (1.14 and 1.17, respectively) but were consistent with past studies which found females were larger than males [5,63]. This results in median sexual dimorphism index (SDI; [64]) of 0.135, 0.174, 0.36, and 0.288 for P. nelsoni, P. peninsularis, P. c. suwanniensis, and P. texana, respectively, based on L ∞ . These differ from previous studies: for P. peninsularis, this is far lower than Aresco's [49] SDI for P. floridana of 0.50, for P. c. suwanniensis, our SDI is slightly lower than that of Jackson and Walker's [65] SDI of 0.45, for P. texana, this is very different from the Lindeman [6] estimate of 0.88, though Lindeman notes their SDI is likely skewed by sampling. In all, our study quantified long-observed similarities between the Pseudemys species and subspecies we sampled while corroborating some of the observed differences and providing evidence for some additional differences that have not been widely noted.
By using life-history invariants of somatic growth, we were able to derive estimates of length-and age-at-maturity for the Pseudemys species/subspecies in our study. Given that these quantities are not estimated from direct measurements of maturity, we chose to generate these quantities within the model to pass any uncertainty in the growth parameters onto the uncertainty of the length-and age-at-maturity estimates. To facilitate this derivation, it was necessary for us to generate a percent of L ∞ that length at maturity occurs, L mat :L ∞ . We calculated this ratio as 71.7% and 87.2% of L ∞ for males and females, respectively, from Trachemys scripta from measurements of maturity and estimates of L ∞ made by Gibbons et al. [51]. These are higher than the ratio commonly used in teleosts, 66%, but smaller than has been used for sea turtles, 97.5% [32,66]. It should be noted that this ratio is a fundamental assumption of our derived length at maturity and age at maturity estimates, and we shared this ratio across species. Future studies should focus on generating this life-history invariant ratio of length at maturity to asymptotic length for a variety of Pseudemys and other turtle species.
For P. nelsoni, Bancroft et al.
[7] estimated males mature at three years old and females mature at seven to eight years old. We estimated a far higher male age at maturity (8.26 years) but only a slightly higher estimate for females (12.2 years). However, our lengths at maturity, 197 mm SCL for males and 271 mm SCL for females, are consistent with the reported length at maturity for P. nelsoni between 187-231 mm SCL for males and 275-290 mm CL for females [47]. In P. c. suwanniensis, Bancroft et al. [7] estimated males mature at five years and females at 16 years, while Jackson and Walker [65] estimated females mature at 10-13 years. Our estimates were close to Bancroft et al. [7] for males at 4.36 years and between the two studies for females at 14.12 years. Additionally, our lengths at maturity, 223 mm SCL for males and 369 mm for females are consistent with past report values of 216 mm and 360 mm SCL [4]. For P. peninsularis, Congdon and Gibbons [67] estimated males mature at three to four years and females mature at five to seven years, while we estimated far higher age at maturity, with males at 9.07 years and females at 10.25 years. Coincidentally, we estimate higher lengths at maturity, 222 mm SCL for males and 316 mm SCL for females, than previously reported 130-150 mm SCL for males and 240-300 mm SCL for females [62,67]. While Lindeman [6] does not report age at maturity for P. texana, the reported length of mature individuals were 87-178 mm SCL for males and 234-265 mm SCL for females, with our reported lengths at the high end for males (174 mm SCL) and above the maximum for females (273 mm SCL).
The unexpected late age at maturity of males for P. nelsoni and P. peninsularis likely results from separate phenomena. In P. nelsoni, the agreement between our estimated lengths at maturity and those published from observed maturation indicates that the ages estimated in Bancroft et al. [7] are underestimated. In P. peninsularis, it is likely that this comes down to the split of P. floridana where Congdon and Gibbons [67] sampled P. floridana in South Carolina, whereas we sampled P. peninsularis (syn. P. floridana peninsularis) in peninsular Florida [1]. Another potential source of discrepancy is that the L mat to L ∞ ratio was derived from the Gibbons et al. [51] study of T. scripta and, without a study of the variability in life-history invariants across Deirochelyinae or Emydidae, there is some reason to interpret our L mat and t mat with circumspection. However, the close agreement between many of the lengths at maturity across our study species suggests that the Gibbons' et al. [51] invariants are likely appropriate for Pseudemys. Life-history invariants have been useful for understanding the variation of life-history strategies across taxa [31,32], and future studies should be pursued in Testudines to understand how these invariants may differ from other indeterminate growing clades.
The joint von Bertalanffy growth model made some key assumptions. Principally, we assumed the environmental effects on L ∞ and k were shared among species. Ideally, we could relax this assumption with further overlap between species and sites (or with additional sites) but with at most any species occurring at two sites, we would have suffered from separability issues in the current model. We concluded this pooling of effects was appropriate as no one species at a site massively outnumbered another in terms of captures, and we still resolved significant site-specific effects on k. Interestingly, we did not estimate any significant effects by any given environmental covariate on L ∞ or k but solely when in aggregate at a site ( Figure 3). We attribute this to the variability in individual growth as well as the low overlap between species and sites.
We also made a structural assumption that individual variation occurs in L ∞ and k, which differs from past approaches that assume individual variation only occurs in L ∞ (see [26] for a review). We made this assumption for biological, empirical, and structural reasons. Growth is the net result of anabolic (energy gain) and catabolic (energy loss) processes, and the von Bertalanffy model [24,68] formulates k parameter to depend on the catabolic rate while the L ∞ depends upon the ratio of anabolism and catabolism. However, there is some debate on how von Bertalanffy [24,68] defines catabolism and whether catabolism includes processes related to routine behaviors such as obtaining food [43]. Empirically, the spring systems studied here are remarkably thermally stable [36], reducing any direct temperature-related effects on k as have been shown in other studies [44]. We were additionally unwilling to discern whether a specific environmental covariate would solely affect anabolic or catabolic processes as none directly measured density-dependent effects or food availability. Structurally, L ∞ and k are highly correlated parameters demonstrated by our negative correlation of −0.31 even in transformed space. Allowing for individual or environmental variation in one parameter while fixing the other seemed unlikely to garner robust growth curves. Failing to estimate one parameter properly has been shown to have severe, negative biases on the corresponding parameter [42]. Thus, we chose to let the data decide where environmental covariates played a role, in L ∞ or k.

Conclusions
Given these assumptions, the novel estimates of Pseudemys species generated with the hierarchical Bayesian somatic growth model provide critical information for the ongoing management and conservation of these similar species [2,52]. Somatic growth rates and age at maturity in wild populations are major data gaps for many freshwater turtle species preventing a robust understanding of the ability to compensate for anthropogenic impacts [3][4][5]58]. Thus, increasing the understanding of how freshwater turtles grow can aid in management of the species at local, state, regional, and range-wide levels [9]. Changes in growth rates at a population level can act as an early warning sign for changes in environmental conditions and or impacts on habitat, such as water quality, food availability, and climate change. Our study serves as a baseline reference for the species, and at these impacted spring sites and provides a benchmark to track the population into the future.
These benchmarks are essential for measuring restoration or management intervention successes in the highly impacted spring systems these species inhabit [20].