Modelling the growth of the brown frog (Rana dybowskii)

Ethics statement

All frogs used in this study were handled in strict accordance with Northeast Agricultural University (NEAU) Institutional Animal Care and Use Committee (IACUC) protocols (IACUC#09-012) and tissues-of-opportunity waivers were approved by NEAU. Tissues-of-opportunity are defined as samples collected (1) during the course of another project with an approved IACUC protocol from another institution; (2) during normal veterinary care provided by appropriately permitted facilities; or (3) from free-ranging animals at appropriately permitted facilities. No endangered or protected species were harmed during this study.

Breeding environment

The experiments were carried out from April to October of 2010 at a farm located in Jiankou Town (46°51′54″N, 130°17′32″E; 80 m above sea level) owned by the Folin company in Jiamusi City, Heilongjiang Province, China. The area experiences a temperate continental monsoon climate with a long winter and a short summer, a frost-free season that typically lasts for 130 d, average yearly precipitation of approximately 510 mm, and an average temperature of 2.8 °C.

The open-air breeding pens were enclosed by a sealed fence, a sprinkler was installed, and some low vegetation was planted. This simulated environment designed for rearing brown frogs resembled an actual forest habitat to ensure the necessary conditions for growth. Age 1 and age 2 brown frogs were reared separately in different pens, but all pens were maintained under the same conditions, including the air temperature, photoperiod, and air moisture as well as the construction materials and physical properties of the pen. The temperature and humidity were recorded every 2 h by temperature–humidity metres. The air humidity in the pens varied slightly but was maintained at approximately 60∼80% by spraying water into the air, and the ground-level humidity of the pens was controlled within a range of 25∼35%. During the trial, the average diurnal, maximum, and minimum temperatures recorded inside the facility were 18.73 ± 3.73 °C, 24.36 ± 3.40 °C, and 13.10 ± 2.29 °C, respectively. The temperature of the pens rarely exceeded 35 °C due to ventilation and spraying, thus preventing stress from rapid changes in the frogs’ body temperature.

The age 1 and age 2 frogs were reared in three different breeding pens, each of which was approximately 80 m² in size. The density of the age 1 brown frogs in the pens was 8/m², and that of the age 2 brown frogs was 4/m². During the experiments, the frogs and their behaviour were periodically examined. When the average temperature was above 10 °C and the frogs were observed to be moving, we provided food (Tenebrio molitor). The frogs were fed once a day at 09:00, and the feeding amount was approximately 4% of their mean body weight. Insects in the natural environment were also available for consumption. The age 2 frogs used in this experiment had been cultured the previous year (2009) using breeding methods consistent with this experiment (2010). In the fall of 2009, the groups of age 2 frogs were placed in separate wintering ponds to avoid errors in discriminating between the different ages. The initial body weights for the age 1 and age 2 brown frogs were 3.47 ± 0.74 and 0.48 ± 0.05 g, respectively.

Sample collection

The snout–vent lengths (SVLs) and body weights of 40 brown frogs were measured each week. Callipers with a precision of approximately 0.02 mm were used to measure the body sizes (SVLs) of the specimens, and their body weights were measured using electronic scales (accurate to 0.01 g). Age 1 brown frogs included frogs in the developmental stages from metamorphosis to out-hibernation in the second year, and age 2 brown frogs included frogs in the stages from out-hibernation in the second year to out-hibernation in the third year. Additionally, the one-year-old brown frogs completed their metamorphosis by 5th June and began eating between 10th June and 15th June. Brown frogs cannot eat after 1st October due to low temperatures, which terminates growth, so they must enter hibernation one month later. For age 2 brown frogs, the hibernation period ended on 15th May, and they began to feed from 25th May to 1st June before entering hibernation at the same time as the one-year-old brown frogs.

Males and females were mix-reared in the first year because the brown frogs did not exhibit sexual dimorphism, but late in the second year, males, and females were reared separately.

Growth functions

Frog growth was studied and mathematically modelled using various nonlinear, linear, and polynomial functions. The nonlinear mathematical models included the logistic, Gompertz, von Bertalanffy, Brody, Janoschek, and Richards models.

These models contain several common parameters and can associate any biological meaning to each of them.

W is the measurement value (g or cm); t is the number of experimental weeks; A is the body weight or length at maturity, B is an integration constant related to the initial weight of the animal, and C is the maturation rate. The value m is the parameter that gives shape to the curve by indicating where the inflection point occurs.

Logistic Model. The following equation describes the Logistic (Fekedulegn, Mac Siurtain & Colbert, 1999) growth model: (1) $$W=\frac{A}{\left(1+B{e}^{\left(-Ct\right)}\right)}$$

Gompertz Model. The following equation describes the Gompertz (Wellock, Emmans & Kyriazakis, 2004) growth model: (2) $$W=A{e}^{-b\exp \left(-Ct\right)}$$

Von Bertalanffy Model. The following equation describes the Von Bertalanffy (Sarasola-Puente et al., 2011) growth model: (3) $$W=A{\left(1-B{e}^{-Ct}\right)}^3$$

Brody Model. The following equation describes the Brody (Hossein-Zadeh, 2015) growth model: (4) $$W=A\left(1-B{e}^{\left(-Ct\right)}\right)$$

Janoschek Model. The following equation describes the Janoschek (Gille & Salomon, 1995) growth model: (5) $$W=A-\left(A-W_0\right)e^{\left(-C{t}^m\right)}$$

Richards Model. The following equation describes the Richards (Fekedulegn, Mac Siurtain & Colbert, 1999) growth model: (6) $$W=\frac{A}{{\left(1+B{e}^{\left(-Ct\right)}\right)}^{\frac{1}{m}}}$$

To estimate body weight and body size at a certain age, first-, second-, third-, and fourth-order polynomial functions were fitted to the brown frog body weight and body size data.

Polynomial Model. The following equation describes the Polynomial (Hadeler, 1974) growth model: (7) $$W=d_0+{\displaystyle \sum _{i=1}^r{d}_i}\times t^i$$

The value of r is the second to the fourth order of fit; d₀ is intercept; d_i is the regression coefficient.

Prediction consistency and quality

The models were assayed for prediction consistency and quality using an adjusted coefficient of determination (R²_adj). Furthermore, the residual standard deviation, or the RMSE, DW, AIC, and BIC were used.

R²_adj was calculated as follows: (8) $$R^2{}_{\text{adj}}=1-\left[\frac{n-1}{n-p}\right]\left(1-R^2\right)$$where R² is the coefficient of determination (R² = 1−(SSE/TSS)); the total sum of squares is denoted as SST; and the residual sum of squares is denoted by SSE. Additionally, n and p are the number of observations or data points and the number of parameters in the equation, respectively. The coefficient of determination is denoted by R² and is a statistical measure of how well a regression line can approximate real data points.

The importance of testing a regression line in this study to resolve polynomial models renders the determination coefficient a useful statistic. The R² value is a measure of the relationship between the general deviation around the average trend expressed by the growth curve model, and its value lies between 0 and 1. If the regression line of a model has a coefficient of determination of 1, it is considered a perfect approximation of the real data. Therefore, the predictions of a model are considered satisfactory if the coefficient of determination approaches unity.

The root mean square deviation is a particularly generalized type of standard deviation that may occur among sets of observations and represents the sample standard variation of differences between the values predicted by a model and the observed values (Wikipedia, 2017). Individual differences between predicted values and observed values are referred to as residuals. The RMSE was calculated as follows: (9) $$\text{RMSE\hspace{0.17em}}=\sqrt{\frac{\text{SSE}}{n-p-1}}$$where SSE is the residual sum of squares; n is the total number of observations (data points); and p is the number of parameters in the equation. The SSE is the squared sum of residuals, which is the squared sum of the individual differences between values estimated from the model and the observed values. The SSE determines the level of discrepancy between data and a model and is therefore a vital criterion for model selection.

The RMSE value is a major criterion for analysing the predictability and adequacy of the growth curve model of a desired function. Consequently, the growth curve of a function with the smallest error in its RMSE is considered the most adequate model.

The DW is used in regression analyses to detect relationships between values separated from one another by a particular time lag; these relationships are generally referred to as autocorrelations. The DW was employed as discussed above, and the results are discussed below.

The DW lies between the values of 0 and 4; a value close to 2 suggests no autocorrelation; a value trending towards 0 implies a positive autocorrelation; and a value trending towards 4 implies a negative autocorrelation (Hossein-Zadeh, 2015). The DW was calculated using the following formula: (10) $$\text{DW\hspace{0.17em}}=\frac{{{\displaystyle {\sum }_t^n\left(e_t-e_{t-1}\right)}}^2}{{\displaystyle {\sum }_{t=1}^n{e}_t{}^2}}$$where e_t is the residual at time t, and e_t−1 is the residual at time t−1.

Furthermore, AIC was calculated using the following expression (Hossein-Zadeh, 2015): (11) $$\text{AIC}=\text{\hspace{0.17em}nln}\left(\text{SSE/n}\right)+2p$$

Akaike’s information criterion is a good statistic for comparing models or functions of varying complexity due to its ability to modify the SSE for a different range of parametric inputs in a particular model. A smaller AIC value implies greater consistency between the data and the model. The basic use of AIC is in model selection, which has its foundation in information theory and is carried out by implicit comparison of several models. Notably, if all the contrasting models fit the data poorly, AIC selection is affected but not disrupted during usage.

Meanwhile, the BIC (Lindberg, Schmidt & Walker, 2015), which is also a model selection criterion, uses the sum of the most likely outcome of data fitting and the selected model but penalizes the logarithm (Log) of the most likely outcome with a term associated with the complexity of the model, as shown in the following equation: (12) $$\text{BIC}=\text{\hspace{0.17em}nln}\left(\text{SSE/n}\right)+\text{\hspace{0.17em}pln}\left(n\right)$$

Similar to AIC values, smaller BIC values from the above expression suggest better consistency between the data and the model.

The model input variables (parameters) were estimated, and statistical analyses were conducted using R software version 2.14.2 (R Development Core Team, 2010). The nonlinear and linear models were implemented by employing the nonlinear system (NLS) and linear modelling (LM) procedures. The LM procedure involved the least-squares method, and the NLS procedure involved the Gauss–Newton algorithm.

Body weight/size ratio

The body weight/size ratio (K_wl) was calculated as: (13) $$K_{\text{wl}}=W*L^{-1}$$where W is the body weight (g), and L is the body length (cm).

Fulton’s condition factor

Fulton’s condition factor (K) was calculated as follows: (14) $$K=100W*L^{-3}$$where K is fatness (g/cm³); W is the weight (g); and L is the length (cm). Generally, larger fatness values indicate that the frog is in good shape (Jin et al., 2015; Mozsár et al., 2015).

The K_wl was calculated for each frog to assess its relative condition. The numerical value of the ratio decreases as the condition of the individual improves (Heikura, 1977).

Growth curves

The body weights and body sizes of the age 1 animals from weeks 1 to 18 are illustrated in Figs. 1A and 2A, which show S-shaped or sigmoidal trends. Additionally, body weights and sizes were assumed to be dependent on age (weeks) according to several growth models for all frogs. The estimated growth curves were typically sigmoidal; in other words, the growth curves were shaped in the form of the letter S (sigmoidal).

The observed body sizes and shapes of the age 2 frogs from week 1 to week 20 are shown in Figs. 1B and 2B, with an obvious sigmoidal trend in body weights with increasing age. Additionally, the estimated body size (cm) as a function of age (weeks) was examined with several growth models for all frogs, and the original estimated growth curve was sigmoidal.

Goodness of fit statistics for the growth models

Low AIC, BIC, RMSE, and RMSE values and a high R²_adj value indicate good fit of the model to frog growth, which occurred for the same model.

The results showed that the third- and fourth-order polynomial models provided the most consistent predictions of body weight for age 1 and age 2 brown frogs, respectively, and the Gompertz and third-order polynomial models yielded similarly adequate results for the body size of age 1 brown frogs. Additionally, the Janoschek model generated a similarly adequate result for the body size of two-year-old brown frogs due to decreases in the AIC, BIC, RMSE, and R²_adj values and an increase in the DW value relative to those of the other models (Tables 1 and 2).

Among all the nonlinear functions, the Richards model also provided the most satisfactory results for the body weight of one-year-old frogs. Meanwhile, the most satisfactory model for the body weight of age 2 frogs was the logistic model, and both the Gompertz and the Janoschek models provided the most consistent predictions of body size for age 1 and 2 frogs, respectively, due to decreases in the AIC, BIC, and RMSE values and increases in the DW and R²_adj values compared to those of the other models (Tables 1 and 2).

The DW of all models was between 0.22 and 0.54. The linear model had the smallest DW values among all considered models of body weight and body size, and the fourth-order polynomial function presented the largest sets of values with the highest DW values among all considered models of the body weight and body size of age 2 frogs (Tables 1 and 2).

The R² and R²_adj values of all models showed slight discrepancies between the models of body weight and those of body size. The R²_adj value ranged from 0.762 to 0.822, and the R²_adj values ranged from 0.777 to 0.909. However, the R² and R²_adj values of body size were higher than those of the models of body weight (Tables 3–6). The Richards model, which is a modification of the sigmoid logistic that allows greater flexibility in S-shaped curves, yielded the greatest R² and R²_adj values, and the linear model yielded the smallest R² and R²_adj values for all frogs.

Initial body weights and body sizes

For the linear model, d₀ (intercept) is interpreted as an estimate of initial body weights or body sizes, although parameter d₁, which is defined as an average body weight or body size, changed over the course of a week. Table 4 shows negative d₁ values for the body weight of age 1 frogs, indicating imprecise estimates of this parameter.

Asymptotic weight and maturation rate

The parameter A is an estimate of asymptotic weight, which can be described as the mature or adult weight. The Brody and Janoschek models yielded the highest and lowest estimates of A, respectively, for the body weights of all frogs. The Janoschek growth model is similar to the Richards model in terms of the level of flexibility (Tables 3 and 5).

Another parameter of interest is the speed of growth to the asymptotic weight, which is represented by the parameter C and is called the maturation rate. In this research, age 2 brown frogs generally showed greater values for this parameter than the age 1 frogs (Tables 3 and 5).

K_wl value and K value

Figure 3 shows that the K_wl values of both the age 1 and 2 frogs increased with growth, and no major changes were observed during the first four weeks or from weeks 13 to 18. We can also see that the growth of the frogs was concentrated during weeks 4–14. The K_wl value of the one-year-old frogs increased from 0.40 to 1.08, and the K_wl value of the age 2 frogs increased from 1.08 to 3.18.

Figure 4 shows that the K value of the one-year-old frogs exhibits uncertain variation, and it decreased from 23.81 to 9.45 in the first four weeks, indicating development of growth characteristics after metamorphosis. The K values increased before declining. The K value of the age 2 frogs decreased from 10.54 to 9.10 within six weeks and then increased until week 15 before decreasing again.

Growth of the brown frogs

The growth of the brown frog, similar to all animals, depends on the genotype-environment interaction (Mansano et al., 2017a). Animal growth is associated with genetics, biotype, race, weight, age, and body state. However, culture management and nutrition supply are also the key factors in frog culture (Mansano et al., 2017b; Olvera-Novoa, Ontiveros-Escutia & Flores-Nava, 2007). Modelling animal growth curves is necessary for optimizing the management and efficiency of animal production (Köhn, Sharifi & Simianer, 2007), and growth modelling provides many benefits for frog production. Many similar growth modelling studies of pigs (Köhn, Sharifi & Simianer, 2007; Wellock, Emmans & Kyriazakis, 2004), sheep (Hossein-Zadeh, 2015), poultry (Rizzi, Contiero & Cassandro, 2013), fish (Baer et al., 2011), crabs (Durán, Palmer & Pastor, 2013), and bullfrogs (Mansano et al., 2012, 2017a; Pereira et al., 2014) are available. Pereira et al. (2014) tested Gompertz and logistic models for captive bullfrogs during the fattening phase (the terrestrial phase). In the aquatic phase, different models were applied for evaluation and simulation of bullfrog tadpole growth (Mansano et al., 2012, 2016). Monitoring the growth and development of brown frogs is essential to better understand the general relationship between production and management on frog farms (Pereira et al., 2014) as well as the relationships between growth and age or growth and feeding over time.

Recent frog studies frequently used the basic unit of year (López et al., 2017; Sarasola-Puente et al., 2011), and research based on days-old or weeks-old animals are rare (Sarasola-Puente et al., 2011; Tessa et al., 2017) even though this scale is very important to studies of froglet growth. However, skeletochronological methods have many deficiencies in terms of studying the growth of amphibians based on days or weeks and many animals must be sacrificed and the work is too difficult and time-intensive.

The artificial breeding technology used for frogs in this study is specific to R. dybowskii, which was domesticated in China (Bing, 2012). Similar to many fully terrestrial frogs in the temperate zone, R. dybowskii lives on land in summer but hibernates in water in the winter. In this study, we built suitable habitats for R. dybowskii and investigated its growth in breeding areas. The time points for metamorphosis and out-hibernation and the proximity of the feeding initiation and termination (for hibernation) dates were determined. The ages of the frogs were known from management of the frog farm. The above time points are pertinent to the intrinsic growth rate, which provides a better understanding of growth over weeks or the explicit relationship between growth and feeding.

Studying frog growth and development is essential in the first two years of development (Iturra-Cid, Ortiz & Ibargüengoytía, 2010) from metamorphosis to sexual maturity because frogs can reach more than 50% of their maximum body weight during this period (Yang et al., 2011). A weight gain of approximately 5–10 times the initial weight was recorded each year. The weights and sizes of the age 1 frogs reflect the growth and development in the current year and affect the next breeding cycle through overwintering survival and slaughter weight for the next year (Hedeen, 1972).

The growth curve for frogs is sigmoidal, or S-shaped, from out-hibernation to hibernation, showing a slow rate of increase in the initial stage before accelerating at the metaphase and then slowing again at a later stage (Hedeen, 1972; Martof, 1956). From the results, we can see that body weight increased throughout the lifetimes of the animals, with only slight decreases in later periods.

The most rapid rate of growth in body weight and size occurs in the first two years (Iturra-Cid, Ortiz & Ibargüengoytía, 2010), and the intrinsic growth rates are higher in the first year than those at other times (Semlitsch, Scott & Pechmann, 1988). In the first year before hibernation, a weight gain of approximately 7.02-times the initial weight was recorded, and the weight increased during the second year to approximately 5.41-times that of the first year.

As shown in Figs. 1 and 2, growth in terms of body length was faster than that of body weight in the initial stage. In Figs. 3 and 4, the relative fatness index of the brown frogs is low; the growth rates for body weight and body size were unsynchronized, and the growth of body weight was slower than that of body size.

In brown frog breeding, growth requires a reasonable food supply and balanced nutrition, but the frogs often showed dietary deficiencies from the 4th to the 10th week after initiating feeding (Densmore & Green, 2007). Therefore, vitamins A and D must be supplied during this period (Ogilvy, Preziosi & Fidgett, 2012). Furthermore, growth in terms of body size began prior to the increase in body weight, and growth in terms of body weight started after body size stopped increasing (Von Bertalanffy, 1957). The final weight was restricted by body volume and was shown to be self-inhibiting.

Environmental factors, such as temperature (time constraints), obviously restrict growth, as occurs with the gradual decline in temperature before hibernation. The brown frog can experience hypothermia, which directly affects its metabolism (Ziegler, Arim & Bozinovic, 2016), followed by a decline in feed intake until feeding and growth stop. Frogs are then faced with six months of hibernation, and they sustain themselves by utilizing energy stored in the liver and fat.

Making sense of the parameters

Knowledge regarding how parameters, such as the growth rate, differ among all developmental stages in frog breeding is important for formulating culture management strategies. Management must be appropriate for populations at different stages.

In this study, the logistic and Gompertz functions were used because of their ability to simply describe vital factors. Some of these simple descriptions include modelling the growth rate as a basic function of size, continuous growth, asymptotes, parameters with biological points, points of inflection, and sigmoidal trends. According to Richards (1959), the four-parameter flexible inflection points of the Richards function that are suitable for modelling animal growth were developed by modifying the logistic and Gompertz functions (Richards, 1959). In contrast, the Janoschek function (Janoschek, 1957) can mainly be used to describe the postnatal growth of individuals, but this function is also flexible in terms of its points of inflection, all of which are desirable properties in nonlinear growth models (Wellock, Emmans & Kyriazakis, 2004). The fourth-order polynomial model was selected due to its similarity to nonlinear growth models. In this study of brown frogs, the Richards function was found to be most suitable for modelling growth according to the low AIC values, and this conclusion was also determined in a study of the growth of ducks and other poultry (Knižetova et al., 1991) after evaluating other studies that used the Richards, logistic, and Gompertz functions to model poultry growth. The R²_adj values for the Richards growth functions were slightly higher than those of the remaining nonlinear models, which all had similar R²_adj values; the Richards growth function also showed the lowest AIC values. The R²_adj value measures linearity and is consequently a more suitable and appropriate factor for linear models.

According to Tompić et al. (2011), the calculated AIC and R²_adj values suggest that the fourth-order polynomial model was more suitable. The estimated A and C parameter values are interpreted as mature body weight and the instantaneous relative growth rate, respectively.

A linear model, unlike nonlinear models, lacks parameters with biological meaning (Brown, Fitzhugh & Cartwright, 1976), but nonlinear models with such parameters can be linearized and evaluated by linear regression (Bates & Hunter, 1985). The results of the current study suggest that the logistic model is adequate for explaining brown frog growth in terms of both body size and body weight, but further research is required to fully understand the growth of brown frogs in China. Ideally, these studies would involve large frog data sets, each with a similar number of body weight records collected and regulated across intervals that are sufficiently similar to permit sufficient sample sizes.

Obtaining frog descendants with rapid growth and sexual maturity rates (Dmitriew, 2011) requires understanding the growth rate, implying that a process exists that alters the curvilinear growth of frogs. There are two strategies for changing growth curves: altering the parameters of a curve while leaving its basic configuration untouched, and changing the basic shape of the curve. The first strategy is the best approach for mature body weight, and the latter strategy is more desirable for the growth rate alone. The second strategy is known to apply to frog culture; the growth in the early stage is exponential, and the frog grows more rapidly if the process is strongly implemented.

The fitness of the captive environment

R. dybowskii is a rare species among breeding frogs. It usually lives on land in a forest environment during its terrestrial stage, which differs from other frogs that exhibit a more typical, amphibious lifecycle.

For R. dybowskii, an amphibian that lives in terrestrial habitats, influencing the micro-climate factors of the captive environment was easy, which is unlike cultivation of other frogs or snakes and turtles that have complex lives both in water and on land. The rearing cycle of R. dybowskii was up to two- to three-years long, and the frog has the characteristics of ‘three adaptions and nine inadaptations’; the three adaptions are gloomy, dampness, and clean, and the nine inadaptations are high light, rain, dry, wind, high temperature, predators, dirty, noise, and off-flavour. These features reflect the high risk of frog culture due to complicated breeding technology and high morality. High-density breeding is the main factor affecting the occurrence of disease and high morality. The present study was conducted at a site within the natural distribution of R. dybowskii, and the density of the captive frogs was relatively low, which can largely prevent the occurrence of disease and ensure normal growth of the brown frog.

Many factors affect the growth of frogs, including not only a captive environment but also the quality and quantity of food and health status (Mansano et al., 2017a). Culturing frogs in an artificial environment is difficult even if some aspects of the habitat, such as temperature and humidity, satisfy the needs of the animals. For example, the average diurnal temperature recorded inside the facility was 18.73 ± 3.73 °C. The temperature of the pens rarely exceeded 35 °C during the trial. If other necessary parameters for optimal culturing are not satisfied, such as density and food composition, growth and health status may change. Typically, animals with free access to food mature faster and grow to larger sizes, and when resource levels decline, animals tend to grow more slowly, reach maturity later, and grow to smaller sizes (Day & Rowe, 2002; Dmitriew, 2011; Lind, Persbo & Johansson, 2008). The breeding system used in this experiment included an appropriate culturing environment, sufficient quality and quantity of food, and proper management, fulfilling the basic requirements for frog growth. The K, K_wl, and fatness of the cultured frogs were greater than those of wild frogs or frogs grown in a much more natural environment (Wang et al., 1999), suggesting that frogs can live in a cultured environment and still obtain a not bad nutrition status.