Age estimation and body size of the Parsley Frog, Pelodytes caucasicus Boulenger, 1896 from Lake Borçka Karagöl, Turkey

. In this study, we described age structure, body size, body mass and the relationships among these parameters for a population of P. caucasicus from Lake Borçka Karagöl, Artvin, Turkey. The mean SVL with standard error was 45.87 mm ± 0.55 (range: 39.98-50.28 mm) and the mean weight with standard error 8.81 g ± 0.39 (range: 6.10-11.47 g) in females whereas 48.16 mm ± 0.45 (range: 43.64-54.78 mm) and 11.32 ± 0.25 (range: 9.56-14.80 g) in males, respectively. We found a significant male-biased difference reflecting sexual dimorphism and statistically significant positive relationships between these variables. According to the results, the age ranged between 2-5 years in females and 2-6 years in males. The mean age distributions significantly differed between the sexes (females: 3.28 years ± 0.19; males: 3.94 years ± 0.20). The mean ages and maximal ages were found identical to the previously reported results from Turkey, but the mean ages were higher than in Georgian populations. Von Bertalanffy growth models demonstrated similar curves, and the growth rate was faster up to 3 years in both sexes. To conclude, this study was the first to determine age structure and growth patterns in Borçka Karagöl population and weight data for P. caucasicus was presented for the first time in the literature.


INTRODUCTION
The determination of individual age is essential for studies like demographic and life history, including developmental biology and population dynamics of a species. This provides the researchers basic ecological data related to population structure such as sexual maturity time, the age structure, growth rate, and life span. Moreover, this is a parameter essential to infer the life history traits of a species and to compare it to other species (de Buffrénil et al., 2021;Ma et al., 2022). In this context, skeletochronology is an effective and reliable method for estimating the growth and age of many amphibian species as the growth of amphibians is not independent of environmental conditions (Guarino et al., 2019;Üzüm et al., 2020). In addition, the method can be also applied to animals such as mammals (Nacarino-Meneses et al., 2016), lizards (Beşer et al., 2019), turtles (Guarino et al., 2020), and even fossils (de Buffrénil et al., 2021). Skeletochronology calculates the lines involving the formation of calcium carbonate, called "annual rings" or "growth markers" in bone tissues (Castanet, 1994;Ma et al., 2022). For this, the diaphyseal region of their long bones, which show weaker vascularity, provides the best result for calculating the age of individuals (Castanet et al., 1993). Rozenblut and Ogielska (2005) showed that lines of arrested growth (LAGs) are most complete in the middle part of the phalangeal diaphysis in European water frogs and thus pointed out that the middle part of the long bone is optimal for age studies. Moreover, the skeletochronology helps to calculate the lifespan of the population in amphibians, explain the sexual size dimorphism, and reveal the differences between the sexes in terms of age and size. Also, the knowledge based on the skeletochronological studies tends to show population dynamics (Peng et al., 2021).
Pelodytes caucasicus Boulenger, 1896, the Caucasian parsley frog, is a native species of the Caucasus fauna. The species is distributed throughout northwest Azerbaijan, Georgia (southwest and South Ossetia), Russia (Krasnodar district), and Turkey (Blacksea region) (Zazanashvili et al., 2012;Litvinchuk and Kidov, 2018;Çiçek et al., 2019). The species is considered as near threatened because of natural and anthropogenic pressures (Ananjeva et al., 2009;Iskanderov, 2009;Kaya et al. 2009), so it is recommended that public campaigns should be conducted to raise the awareness for this endemic relict species in the border of Georgia and Turkey . From this aspect, it is important to reveal the population dynamics of P. caucasicus in a new population based on skeletochronology. Although the age structure of P. caucasicus was reported by two studies in Georgia (Gokhelasvili and Tarkhnisvili, 1994;Chubinishvili et al., 1995), there is only single study in the border of Turkey (Erişmiş et al., 2009). Given the importance of skeletochronological information to better understand the population dynamics, ecological and evolutional processes, we aimed to present the age structure of the Borçka Karagöl population for the first time and compare the results with the previous studies. We also provided weight data for the first time in this species, as well as assessed the relationship between this trait and SVL for age and growth rate.

Fieldwork
We sampled 50 specimens (18 ♀♀, 32 ♂♂) during the breeding season (2013) from Lake Borçka Karagöl territory, in the vicinity of Artvin, Turkey (Fig. 1). The collected frogs were anesthetized with MS222. Snoutvent length (SVL) was measured at the nearest 0.01 mm using a digital calliper and the weight was noted using the nearest 0.01 g an electronic balance for each sample. To count growth lines, the longest toe of the left hindlimb was clipped and preserved in 70% ethanol. After the sampling procedure, all frogs were released at the site of capture. Sex determination was possible due to the presence of secondary sexual characters, such as nuptial pads on the forearm, and presence of vocal sacs in males.
The standard skeletochronology procedure (Castanet and Smirina, 1990) was applied to calculate the number of the lines of arrested growth (LAGs). After removing soft tissues, the phalanges were washed in tap water, then decalcified in 5% nitric acid for 2 hours. Lastly, the samples were washed in distilled water overnight to ensure the removal of nitric acid. The 18 μm cross-sections were obtained from the diaphyseal part of each phalanx using a cryostat (Thermo Shandon Cryotome, Germany) and stained with Ehrlich Haematoxylin approximately for 15 minutes. The stained cross-sections were treated with glycerol to control under a light microscope. The images of selected sections were acquired using Olympus BX51 microscope with an integrated camera to estimate the precise age of each specimen. The number of LAGs was counted by two researchers (S. Gül and N. Özdemir). The distance between two adjacent LAGs is a well-known indicator demonstrating individual growth in a given year (Kleinenberg and Smirina, 1969;Kurnaz et al. 2018). Following previous studies (Özdemir et al., 2012; Üzüm et al., 2014), the endosteal resorption was evaluated by comparing the diameters of eroded marrow cavities with the diameters of non-eroded marrow cavities in sections from the youngest specimens.

Statistical analyses
To summarize data and present basic features of the measurements we calculated descriptive statistics using the psych package (Revelle, 2019). The normality assumption of the variables was checked using the Kolmogorov-Smirnov test. Since the variables followed a normal distribution (P > 0.05), we used parametric tests in downstream analyses. The homogeneity of variances was compared using Levene's Test in the car (Fox and Weisberg, 2019) package. We utilized the student t-test to compare the mean differences of the variables between males and females. Thereafter, Pearson's product-moment correlation test was used to estimate the relationships among SVL, weight, and age.
The relationship of SVL and weight was analysed using a linear regression model. Additionally, we used an ANCOVA test to explore the patterns of SVL and weight between sexes, with age as the covariate. Thereafter, posthoc tests were applied using the emmeans package (Lenth, 2021) under Bonferroni correction (estimated marginal means: aka least-square means or adjusted means).
We estimated growth curve models under the typical Von Bertalanffy's equation modified by Beverton and Holt (1957): Lt=L ∞ {1-exp[-k(t-t 0 )]} where Lt is the expected or average length at the time (or age) t, L ∞ is the asymptotic average length, k is the so-called Brody growth rate coefficient and t 0 is a modelling artifact that is said to represent the time or age when the average length was zero. To estimate parameter values and run the analyses FSA (Ogle et al., 2021), FSAdata (Ogle, 2019), FSAsim (Ogle, 2020) and nlstools (Baty et al., 2015) packages were used following the guide "fishR Vignette" prepared by Derek Ogle (2013). To visualize growth curve, we also added a hypothetical individual to the dataset by the reference of Erişmiş et al. (2009) under the presented parameters: SVL 0 at metamorphosis is fixed to mean 20.15 and t 0 (age at metamorphosis) is 0.3 year.
Statistical analyses were carried out using the stats package. Data visualization was performed using ggplot2 (Wickham, 2016), ggpubr (Kassambara, 2020), and ggally (Schloerke et al., 2021) packages. All analyses were run in R Programming Language (R Core Team, 2020). RESULTS We successfully aged 50 individuals using skeletochronology technique. Ages ranged between 2-5 years in females, and 2-6 years in males. Descriptive statistics of both sexes and age groups are presented in Table 1.
Growth curves estimated by Von Bertalanffy's model fit adequately described the relationship between age and SVL, and the curves indicated similar shapes in both sexes (Fig. 4). The final models were found statistically significant for all parameters (P < 0.01). According to the constructed models, the estimated asymptotic SVL was not higher than maximum recorded SVL values (Males: 54.78 mm; Females: 50.28 mm). The growth parameters were presented in Table 2. DISCUSSION The sexual size dimorphism (SSD) is an observable characteristic in animals, and it is known that the direction of size dimorphism in most of the amphibians is female biased (Kupfer, 2007). Shine (1979) noted that females are larger than males in 61% of urodeles and 90% of anurans. The female biased SSD is generally explained by the fecundity advantage hypothesis (Shine, 1989;Andersson, 1994) when the selection is supporting large females due to the larger energy storage capacity for reproduction and more offspring production. However, a male biased SSD is relevant to the dominance in contests of strength, the extent of endurance, mate choice and higher sperm competition success in animal kingdom (Hudson and Fu, 2013). The male-biased SSD corresponding to 10% amphibian species is especially associated with territoriality behaviours and male-male competitions (Nali et al., 2014).
The family Pelodytidae is including five different species from the single genus Pelodytes distributing in Western Europe especially in the Iberian Peninsula, Caucasia, and north-eastern Turkey (Amphibiaweb, 2022). The previous studies have reported the female biased sexual size dimorphism in Pelodtytes species inhabiting in Iberia (Talavera, 1990;Escoriza, 2017). ing their morphological measurements, all the species (P. punctatus, P. hespericus, P. ibericus and P. atlanticus) were characterized by larger average body length in females. However, our results revealed that unlike Iberian species, P. caucasicus males have a larger and heavier body than females. Our results are also consistent with the findings presented in previous studies. For example, Erişmiş et al. (2009) have comprehensively assessed the age structure and growth in Pelodytes caucasicus from Uzungöl, Turkey and. The mean SVL of adult males with SD were reported as 47.16 ± 2.87 mm (n = 44; range 41.48-52.58 mm), and significantly smaller in females (45.79 ± 2.29 mm; n = 31; range 40.28-50.62 mm). Arıkan et al. (2007) noted the range of SVL in sexually mature males between 45.06-52.08 mm, and 46.70-49.62 mm in mature females in Uzungöl population. Yildirimhan et al. (2009) also recorded the mean SVL 50.6 ± 3 mm (range: 43-57 mm) in the population of Çaykara, Trabzon including 47 males, 7 females. Erişmiş et al. (2009) emphasized that the bigger size of older males in their study may deviate the results linked to SSD, but the size differences could also be derived from the biotic or abiotic selective pressures in poikilotherm animals. They validated their findings based on the former studies conducted in climatically different regions in Caucasia (Gokhelashvili and Tarkhnishvili, 1994;Chubinishvili et al., 1995) and they suggested male biased SSD is the species characteristic of P. caucasicus. The common SSD pattern in anura is generally female biased (Monnet and Cherry, 2002). Recently, Pincheria-Donoso et al. (2021) have investigated the SSD of amphibian species in global scale. As a result, they observed 90.8% female biased SSD 7.5% male-biased SSD Fig. 3. The matrix is demonstrating the correlation coefficients, scatterplots, and density plots of binary variables. Corr values indicate the correlation coefficients (r). The significance level of correlation coefficients was represented with asterisk (*, P < 0.05; **, P < 0.01; ***, P < 0.001). and 1.7% monomorphic in anurans. However, Han and Fu (2013) determined male biased SDD in 19 different anuran families distributing in six distinct continents and along tropical and temperate habitats. Moreover, constrained habitats and microhabitats can affect the female size and fecundity, and it can be resulted with more similar body size of both sexes as observed in Hylid frogs (Silva et al., 2020). From this aspect, it can be suggested that the directional difference of SSD between the Iberian species and Caucasian P. caucasicus may cause due to different ecological preferences in response to their habitats. The Iberian Peninsula which has complex orographic structure is surrounded by Atlantic Ocean and Mediterranean Sea. The Mediterranean coast and surrounding areas show dry and warm/hot summers, and mild and wet winters. The Atlantic coasts are characterized by oceanic type of climate, milder but humid winters and cooler/wetter summers, without large temperature variations. Lastly, the inner areas which is represented by continental climate type have hot and dry summers,  and cold and humid winters with large-scale precipitation, but also semi-arid areas extremely low precipitation and very hot temperatures especially in summer season (Carvalho et al., 2021) which are corresponding to the distribution of Iberian Pelodytes taxa (see Diaz Rodrigez et al., 2017). However, Caucasus Ecoregion has quietly high mean annual rainfall in the southwestern part over 2000 mm in the coastal area of the Black Sea. The mean annual temperature in the South Caucasus part of the Black Sea coast around 15 centigrade degree (Zazanashvili, 2009). Besides, the distribution of P. caucasicus species is restricted to the humid subtropics in Caucasia and Turkey with Colchic vegetation type (Tarkhnishvili, 1996;Iskanderov, 2009;Beşir and Gül, 2019). Gül (2014) also noted that the species prefers wet and warm microhabitats in Turkey. Additionally, the amount of precipitation is known as one of the most important factors constructing the distribution of P. caucasicus (Lukina and Koneva, 1996;Litvinchuk and Kidov, 2018). Pincheria-Donoso et al. (2021) proposed that the underlying impact of geographical variation in climatic pressures can shape largescale patterns of SSD in synergy with natural and sexual selection such as intensification of fecundity selection to shorten breeding season in anurans. They also implied that the different selection forms can shaped by macroecological factors climate, geographical gradients and temperature seasonality which are triggering the evolution of life-history traits associated with fecundity. Therefore, we suggested that the ecological preferences of P. caucasicus are potential reason causing male-biased SSD comparing to the representatives of Iberian Peninsula. Radojićić et al. (2002) also revealed similar SSD pattern in the Bombina species. They noted that Bombina bombina showed male-biased SSD while the larger body size of females was observed in B. variegata subspecies, and the differences were associated with reproductive behaviours and possible ecological differences. This pattern can also be observed in Pelobates syriacus (Bülbül et al., 2020 and references therein), Rana arvalis populations (Glandt and Jehle, 2008) and Rana nigrovittata which has similar ecological needs with P. caucasicus such as streams in shaded forest environments (Khonsue et al., 2000). On the other hand, the male-male competition was reported in Pelodytes species (Pargana, 2003;Marquez et al., 2004 and unreported mating ball of P. caucasicus observed by S. Gül). Therefore, it should be also taken into consideration that sexual selection which maximize the fitness in reproductive traits may be an alternative force to describe the SSD pattern in our study as reported in different anuran species e.g., Paa spinosa (Gen-Yu et al., 2010) and Hypsiboas atlanticus (Camurugi and Juncá, 2013). SSD can be driven by life-history traits, so the accuracy of the age determination is critical to estimate these characteristics. The post-metamorphic terrestrial growth is a major part of total growth (over 90%) in amphibians (Werner, 1986) and the variation in terrestrial growth rates and age at maturity is playing an important role in the intersexual adult size variability (Marangoni et al., 2012). Furthermore, the age determination of amphibians yields crucial information on the demographic features such as size at sexual maturity growth, longevity, and growth rate (Duellman and Trueb, 1994;Otero et al., 2017;Baraquet et al., 2021). In this study, we found that the age is ranged between 2-6 years. The mean age is 3.28 years in females and 3.94 years in males. Our constructed Von Bertalanffy's growth curves adequately fitted age/SVL relationship. The models demonstrated similar curve shape for both sexes, but the growth coefficient was higher in females.
Previously, Gokhelasvili and Tarkhnisvili (1994) studied the age distribution of the reproductive population of P. caucasicus in Borjomi canyon during two consecutive years (1992)(1993). According to their results, the mean ages were 2.96 and 2.74 for males, 2.74 and 3 for females, respectively. They also reported the dominancy of young specimens in the reproductive portion of the populations and a very-high annual mortality rate index (0.78-0.83). The males reached maximum of 6 years, and 4 years in females, which is in accordance with what Gokhelasvili and Tarkhnisvili (1994) reported. In the following year, Chubinishvili et al. (1995) studied certain aspects of the population ecology of P. caucasicus. They reported sexual maturity between 2-3 years. Contrary to their findings, the mean age is found approximately one year above the one measured in the Borçka Karagöl population. Erişmiş et al. (2009) reported the mean age of males and females 3.61 ± 0.9 years and 3.03 ± 0.7 years, respectively and the age structure difference is statistically significant (P < 0.05). The oldest male was 5 years while the female was 4 years. The sexual maturity is reached at two years of age. We also obtained approximate values from the mean SVL and age (see table 1). The differences observed in SVL and age between sexes were also statistically significant. Given the identical results for Uzungöl and Borçka Karagöl populations, we can assume that there is a separation between Georgian and Turkish populations related to age structure. Erişmiş et al (2009) also calculated a survival rate at 0.78 in males and 0.76 in females, and these figures can be taken as a reference to explain the high mortality rate noted by Gokhelasvili and Tarkhnisvili (1994). Based on the literature, the maximum age reported in the genus Pelodytes is 10 years for females and 8 years for males in Pelodytes punctatus species (Esteban et al., 2004). Additionally, the mean age (Burgos: 1.71 years ± 1.41 years, Valencia: 3.82 years ± 1.22 years) and mean SVL (Burgos: 34.36 mm ± 2.22 mm, Valencia: 36.41 mm ± 2.50 mm) of males in the study of Esteban et al. (2002) were lower than our mean values. The authors noted that males were aged between 1-6 years old, and the age structure was highly skewed corresponding 50% of the sample of males being 1 year old. This situation can cause uncertainties when age/SVL relationships in the genus Pelodtyes. In this study, we also did not age any individual in 1 year old from Lake Karagöl population. However, it is known that the low average age and high proportion of younger individuals might be associated with rapid decrease of local population of the species such as Rana porosa (Misawa et al., 2002) and Esteban et al. (2002) pointed that this pattern is more relevant to conservation status due to anthropogenic pressures in the Burgos population.
According to the growth curves estimation of Erişmiş et al. (2009), asymptotic SVL and growth coefficient were very similar between the sexes (Males: SVL max : 53.42 mm ± 1.01 mm; K = 0.42 ± 0.03; Females: SVL max = 52.04 mm ± 0.75 mm; K = 0.38 ± 0.01). In both sexes, the estimated asymptotic SVL was slightly higher than the maximum SVL record and the age-specific growth rate under the 3 years reported faster than higher ages. On the contrary, the estimated SVL max values were lower, but the estimated K values were higher in our constructed growth model. From this aspect, we think that the more bell-shaped curve of our models is likely due to lack of the data from one-year individuals, and this may affect parameter fitting. The second reason is likely the preferred formula equation differences. The common point of the models is a remarkable decrease of growth rate after 3 years following sexual maturity in both sexes.
Data on body mass is limited in amphibians, but Santini et al. (2018) emphasized that weight data can be used to assess ecological and evolutionary processes such as dispersal distance, reproduction, population abundance and energy intake and it was also highly correlated with SVL in various anuran families (Bufonidae, Hylidae, Myobatrachidae, Ranidae). Moreover, they said that rapid body mass rise is associated with SVL rise for terrestrial and semi-aquatic species because allometric relationships between length and mass vary in amphibian species based relevant to the different habitat preferences. Lastly, they suggested the geometric similarity hypothesis (Hill, 1950) is fitting these assumptions better because if two organisms are geometrically similar their linear dimensions can be made equal by multiplying those of one of them by a constant. In this study, we contributed to the literature by first recorded weight data of P. caucasicus species. In the genus Pelodytes the available weight data was presented by Esteban et al. (2002) for male P. punctatus specimens (Burgos: 2.36 g ± 0.28 g, Valencia: 3.71 g ± 0.58 g). Contrary to the mean weight in our male specimens, Spanish populations have lower values in parallel with mean SVL and age differences. The linear regression model and ANCOVA results also indicated the allometric relationship between SVL and weight, and the weight differed both between sexes and among age groups. Habitat preference of P. caucasicus are terrestrial and freshwater systems , thus our findings are suggesting the geometric similarity hypothesis and habitat preference effect proposed by Santini et al. (2018) to describe the relationship between SVL and weight. In addition, Xiong et al. (2020) noted that age is affecting the intersexual difference on the weight for the Shangcheng Stout Salamander Pachyhynobius shangchengensis. From this aspect, the synchronized rise of SVL with aging may contribute to the intersexual difference associated with weight in P. caucasicus species. Camurugi and Juncá (2013) said that malebiased sexual size dimorphism corresponding to length and weight is unusual pattern and among frogs and even if it is generally described by male-male competition, the reasons are not fully understood yet. They also described the male-biased SSD as synapomorphyc characteristic in Hypsiboas punctatus species group. Monnet and Cherry (2002) indicated that Bufo achalensis, Rana cascadae, R. nigrovittata species showed male-biased SSD because males were older sex as we determined that the mean age of P. caucasicus males. This pattern is causing due to delayed maturity which is determinant factor when larger body is contributing to mating success and males are not accelerating the growth to breed earlier than females. Tarkhnishvili (1993) said that the lower fecundity, smaller eggs, the volume of clutch and smaller body are the differences in spawning mood associated with female body size, and they are causing to rapid maturation for the species including Pelodytes caucasicus as a strategy to shorten the period between generations and to increase the number of adult individuals. Therefore, we assume that the male biases SSD observed in weight data may also occur in response to the life-history and reproductive traits.
To sum up, we determined age structure and growth patterns in Borçka Karagöl population. Future studies can comprehensively investigate the relationship between these variables and ecological conditions, life-history traits, and reproductive characteristics.