Shifts in the developmental rate of spadefoot toad larvae cause decreased complexity of post-metamorphic pigmentation patterns

Amphibian larvae are plastic organisms that can adjust their growth and developmental rates to local environmental conditions. The consequences of such developmental alterations have been studied in detail, both at the phenotypic and physiological levels. While largely unknown, it is of great importance to assess how developmental alterations affect the pigmentation pattern of the resulting metamorphs, because pigmentation is relevant for communication, mate choice, and camouflage and hence influences the overall fitness of the toads. Here we quantify the variation in several aspects of the pigmentation pattern of juvenile spadefoot toads experimentally induced to accelerate their larval development in response to decreased water level. It is known that induced developmental acceleration comes at the cost of reduced size at metamorphosis, higher metabolic rate, and increased oxidative stress. In this study, we show that spadefoot toads undergoing developmental acceleration metamorphosed with a less complex, more homogeneous, darker dorsal pattern consisting of continuous blotches, compared to the more contrasted pattern with segregated blotches and higher fractal dimension in normally developing individuals, and at a smaller size. We also observed a marked effect of population of origin in the complexity of the pigmentation pattern. Complexity of the post-metamorphic dorsal pigmentation could therefore be linked to pre-metamorphic larval growth and development.

Patterning has a strong genetic component 16,21 although the formation of complex patterns can result from developmental alterations 11 , possibly influenced by environmental stressors 22 . Indeed, a higher pattern complexity in the black bib of red-legged partridges has been shown to be an indication of higher individual quality 23 . Likewise, ornament symmetry in little bustards 24 and facial patterning of paper wasps 25 also reflect the developmental conditions experienced.
While much of the research on colour patterning in vertebrates has focused on birds and teleost fish 1 , amphibians have been an important system in this topic as well 2,[26][27][28] . A characteristic feature of many amphibians is their complex life cycle, which most typically consists of an aquatic larval stage and a terrestrial adult phase 29 . In fact, amphibian larvae and adults differ physiologically, anatomically, and ecologically 30 to the extent that they are at times conceptually treated as separated units 31,32 . Nonetheless, multiple studies show that changes in the developmental trajectory of larvae have direct effects on adults [33][34][35] . This is particularly relevant for pigmentation, as the adult patterning forms during late stages of the larval development and metamorphosis 36 . Hence, if development of pigmentation is altered during this critical period, the adult phenotype will be directly affected. This has been shown in zebrafish, which undergo similar complex life cycles, where environmental stress at the larval stage affects adult pigmentation 37 . Amphibian larvae can decouple growth and differentiation to adjust their development to local environmental conditions 35,38 , constituting an ideal system to test for consequences of plastic alterations of development on pigmentation. Understanding the evolution of colour patterning in adult frogs, regardless of its function in signalling or camouflage, is therefore contingent on understanding how larval development affects the correct formation of patterns.
Here we test whether the dorsal pigmentation pattern of juvenile spadefoot toads (Pelobates cultripes) is affected by developmental alterations during the larval phase. Adults are semi-fossorial and nocturnal, typically living in sandy substrate where they can bury themselves easily. Their dorsal colour patterning can be described as darkly mottled with uneven blotches on a light uniform background and likely serves camouflage by background matching. Its larvae can grow for extended periods of time, reaching large sizes at metamorphosis under benign conditions, but can accelerate development and precipitate metamorphosis if at risk of desiccation from pond drying 39 . Developmental acceleration is in part achieved through increased corticosterone levels and metabolic rate and although critical for evading drying ponds, it comes at the cost of smaller size at metamorphosis, higher oxidative stress, reduced immunocompetence and allometric changes in body shape [39][40][41] . To test whether postmetamorphic pigmentation patterning is impacted by environmentally induced acceleration of pre-metamorphic larval development in spadefoot toads, we reared individualized tadpoles from multiple egg clutches and from four populations in standardized conditions, and experimentally induced accelerated development in half of them by reducing their water level. Upon metamorphosis, we compared body mass, snout-to-vent length, and dorsal patterning between accelerated and non-accelerated individuals. Patterning was quantified using digital imaging to estimate parameters describing different aspects of the texture and complexity of pigmentation.

Results
Effects of water level on growth and development. Experimentally decreasing water level led tadpoles to accelerate their development by 21% on average (F 1,150 = 24.8, p < 0.001), reducing the age when Gosner stage 46 was reached from 171 ± 55 days to 132 ± 44 days. Larvae subjected to reduced water levels metamorphosed at lower body mass and shorter body length: body mass was decreased from 1.84 ± 0.97 g to 1.24 ± 0.31 g, representing 34% decrease (F 1,150 = 86.36, p < 0.001), and snout-to-vent length was decreased from 26.62 ± 12.13 mm to 23.39 ± 4.79 mm, representing a 14% decrease (F 1,150 = 101.4, p < 0.001).
Principal component analysis of pigmentation pattern characteristics. Eleven variables have been calculated to describe 3 features: frequency distribution of grey scales (mean, variance, skew, and kurtosis), homogeneity (angular second moment, inverse difference moment, entropy, contrast, and correlation), and complexity (lacunarity and fractal dimension). We observed substantial variation in pigmentation patterns in the metamorphosed toads emerging from our experiment. A principal component analysis showed that most of the calculated pattern descriptors loaded considerably on the first principal component, with entropy, mean, variance, and contrast of the grayscale distribution of the pattern showing similar magnitudes but opposite directions to kurtosis, correlation and inverse difference moment (S1 in Supplementary material). Fractal dimension and lacunarity loaded heavily on the second principal component while loading in opposite directions, angular second moment loaded on the third principal component, and skewness on the fourth component. The numerous variables showing similar loadings on PC1 indicated high collinearity among variables, which we further visualized using the ggpair function of the R package GGally 42 (S2 in Supplementary material). Principal components 1, 2, and 3 together explained over 89% of the variance in the pigmentation patterns across toads (S1 in Supplementary material). The PC scores largely overlapped across populations and across treatments (Fig. 1).
Effects of developmental acceleration, population and sibship on pigmentation pattern. We used Bayesian hierarchical linear models to test for the effects of experimentally decreasing water level on multivariate pigmentation pattern, specifically on the mean of the grayscale distribution, angular second moment, and fractal dimension chosen based on their high loading on different principal components and low collinearity among each other ( Fig. 2; see "Material and methods" section for more detail). Our model comparisons based on measures of the prediction accuracy by expected pointwise predictive density (ELPD) showed that the model including sibship nested within population as random (group) effect was the better fit than the one excluding it: the predictive ability was improved by adding the random effect according to both the loo (ΔELPD: fit1 = 0. www.nature.com/scientificreports/ Therefore, it could be inferred that sibship and population effect explained a considerable fraction of the variance of each of these variables and we included population and sibship (clutch) as random effects in the subsequent models. Since aspects of the pigmentation pattern seemed to be attributable to sibship, we estimated broad sense heritability (H 2 ) of these three variables (Table 1) using the sommer 43 package in R. The Bayesian generalized coefficients of the first model including only water treatment as fixed effect (fit1) showed that the experimental treatment (i.e. water level) explained a considerable fraction of the variance of the variables, while more so in mean and angular second moment than in fractal dimension (mean: R 2 = 0.29, fractal dimensions: R 2 = 0.17, angular second moment: R 2 = 0.23). Such proportion of variance explained increased when larval period was added as an additional fixed effect (fit2), and again, mean and angular second moment were better explained than fractal dimension (mean: R 2 = 0.36, fractal dimensions: R 2 = 0.20, angular second moment: R 2 = 0.31). Consistently, including body mass at metamorphosis as additional fixed effect improved the model fit (fit3), and mean and angular second moment were better explained than fractal dimension (mean: R 2 = 0.37, fractal dimension: R 2 = 0.19, angular second moment: R 2 = 0.31). In all models, 100% of the observations had a lower pareto k value than 0.7, indicating optimal effective sample size.
The conditional effects of this best fitting model ( Fig. 3) showed that decreased water levels led to lower mean of the grayscale distribution, lower fractal dimension and higher angular second moment. Similarly, the duration of the larval period was positively associated with mean grayscale and fractal dimension and negatively with angular second moment.

Discussion
The dorsal pigmentation pattern of P. cultripes metamorphs varied in texture and complexity across populations and sibships, but also in response to altered developmental rate. The among-sibship and among-population variation is intriguing because it suggests a heritable component. Our estimated broad-sense heritability for brightness, homogeneity and complexity aspects of the pigmentation patterns ranged between 0.15-0.30, although our full sib design prevented us from disentangling additive genetic from non-additive and maternal components 44 and should only be considered a general indication of a hereditary component of the pattern.
As expected, spadefoot toad tadpoles accelerated their development in response to reduced water level, resulting in tadpoles metamorphosing earlier but at a smaller size. Interestingly, this developmental acceleration of the larvae also caused alterations in different aspects of the dorsal pigmentation patterns of the emerging toads. Toads that had accelerated their larval development (i.e. experienced stressful conditions) showed a less contrasted and more homogenous pattern (less random in its intensity, as interpreted based on low fractal dimensions and high angular second moment values; Figs. 2, 3) compared to control toads (i.e. benign conditions), with the blotches being darker and more highly connected amongst each other. Toads induced to accelerate their development also showed a lower complexity of their pigmentation pattern and lower fractal dimension, or "gappiness" 6 .
Developmental acceleration in response to pond drying is known to take a considerable physiological toll on amphibian larvae 41,45 . In order to achieve such acceleration, spadefoot toads increase the corticosterone and www.nature.com/scientificreports/ thyroid hormone secretion and their metabolic rate 39,46 . The formation of the adult pigmentation pattern in the epidermis of amphibians is induced during metamorphosis by the melanocyte-stimulating hormone (MSH), in interaction with thyroid hormone 47 , to the extent that hypophysectomy causes bleaching of tadpole skin 48 . For instance, in oriental fire-bellied toad tadpoles (Bombina orientalis), experimental inhibition of thyroid hormone led to a marked reduction in the size of melanophores on their dorsal skin 26 . We thus hypothesize that enhanced pituitary activity in response to perceived risk of predation may trigger over-melanization of tadpole skin as When raised under environmentally stressful (i.e. low water level treatment) conditions, toads metamorphose with an overall darker dorsal pattern (i.e. lower mean grey value), with decreased heterogeneity (i.e. lower fractal dimension) and with a more homogeneous dorsal pattern (i.e. higher angular second moment). www.nature.com/scientificreports/ they approach metamorphosis. Sexual hormones are also known to affect the development of pigmentation pattern in amphibians, and disruption of these hormones cause profound changes to the pigmentation of postmetamorphic amphibians [49][50][51] , although it is still unclear how induced developmental acceleration may affect the secretion of these hormones. Pigments such as carotenoids are known to prevent oxidative stress 52 and enhance the coloration of metamorphosing frogs 53 , and their availability is reduced in birds when raised under stressful conditions 52,54 . In amphibians, however, a direct link between pigmentation and oxidative stress is not yet clearly established as not enough studies have explored amphibian pigmentation in the context of environmental stress experienced during their larval period. Rather, past studies have focused on pigmentation as aposematism 55,56 or camouflage 57 . Interestingly, deformed individuals originating from stressful conditions showed larger pigment spots despite being smaller in size in comparison to individuals originating from benign conditions in northern leopard frogs (Rana pipiens) 58 . Newts have been shown to be capable of background matching 59,60 , and newts that expressed more pigments have been shown to have a higher metabolic rate, suggesting a metabolic cost entailed in pigment expression 60 , but it is still unknown whether such increased metabolism has carry-over consequences.
While we find that P. cultripes metamorphs originating from stressful conditions (i.e. accelerated development) show darker and more continuous blotches, at this stage it is hard to tease apart whether this was to compensate for smaller size resulting from limited time of growth or rather a side effect of developmental acceleration and insufficient time for growth. Although variation in coloration can often play a role in thermoregulation, UV protection, predator avoidance, or sexual signaling 27,61,62 we do not yet know whether variation in the pigmentation pattern of this burrowing, nocturnal species with aquatic acoustic communication is of any adaptive value. Adaptive coloration in amphibians remains a vast field to be explored, as for instance, subterranean caecilians have evolved bright coloration 63 . Pelobates cultripes is mostly an amphibian of nocturnal habits (although males may stay by the clutch during daytime on the first day post-oviposition) and its dorsal pigmentation likely serves as camouflage rather than communication 64 . However, whether the observed variation in pigmentation and complexity across populations is an evolved response to different levels of predation pressure and whether developmental alterations of such patterns increase their risk of predation remain open questions. For now, our finding that tadpoles induced to accelerate their development emerged with a more homogeneous and darker pigmentation pattern (Figs. 2, 3) suggests two hypotheses. First, the development of lighter, well-contrasted, complex dorsal patterns in normally developing individuals may be disrupted due to physiological stress, indicating deterioration of body condition that carries over to the juvenile phase. Second, the darker and more homogenous pigmentation of accelerating individuals may simply reflect heterokairic effects (i.e. environmentally induced shifts in the timing of developmental events 65,66 ), wherein skin development did not proceed at the same pace as www.nature.com/scientificreports/ that of other organs in accelerating individuals. Individuals subject to developmental acceleration may not have had enough time and/or resources to allocate for the optimal expression of pigmentation, which is what might have led to their darker and less complex patterns. Amphibians show remarkable variation in pigmentation patterns, both within and among species. Here we show that patterns can be altered due to plastic alterations of development. Future studies would have to determine whether the texture and complexity of the pattern are heritable, the relative importance of parental effects, whether it is static or changes ontogenetically, and whether it shows geographic variation or phylogenetic signals.

Material and methods
Study animals and preparation of dorsal pigmentation photos. We collected portions of five egg clutches of four P. cultripes populations from southern Spain (Huelva province) and reared them in climatic chambers of Estación Biológica de Doñana, CSIC (Seville). Upon reaching the free-feeding stage, 14 tadpoles per sibship (clutch) were individualized in 4 L round plastic containers. They were kept in identical conditions (20 °C, carbon filtered dechlorinated tap water) and fed a mixture of ground up rabbit chow (alfalfa) and spinach. At Gosner stage 35 (i.e. stage at which digits are identifiable on the hindlimb bud 65 ), we dropped the water volume to 650 mL in half the containers to simulate risk of pond drying. Gosner stage 35 is the optimal developmental stage at which spadefoot toad larvae are able to perform the most acceleration in comparison to other developmental stages 67 . The experimental procedures were approved by CSIC's IACUC committee (authorization #560-2017), and all experiments were performed in accordance with relevant guidelines and regulations.
Image acquisition and preprocessing. We recorded the duration of the larval period of each toad, marked by the completion of metamorphosis (full tail resorption, Gosner stage 46) 68 at which stage we also measured snout-to-vent length (to the nearest 0.1 mm) and body mass (to the nearest 0.1 mg). A dorsal photograph of each toad placed over laminated grid paper was taken with a Nikon D50 camera. Since photos were taken from slightly different distances, images were downscaled to the same resolution (13.2 px/mm) by bicubic interpolation. We focused our analysis of pigmentation patterns on the torso and head (except eyeballs), excluding the limbs to avoid noise from variation in posture. Regions of Interest (ROI) on the toad's dorsal surface were manually outlined.
The resulting ROI images were transformed to grayscale 8-bit images (256 intensity values). Intensity values of images were obtained from the image histogram, in which the reflectance of each pixel ranged from 0 to 255. Uneven illumination was corrected using a pseudo-flat field correction tool (blurring radius: 40 pixels). Contrast was enhanced by a histrogram stretching to a full 8-bit range of gray (0.1% saturated pixels). The unsharp masking filter tool of imageJ (weight value of 0.5 and blur radius of 1.0) was applied to improve sharpness. Due to artificial light sources, wet dorsal skin of toads showed specular reflection, which may potentially have resulted in errors in image analyses. We therefore applied threshold criteria to determine areas affected by specular reflection, and removed them from the ROIs (intensity values > 190).
Extraction of textural features. Image texture can be defined as the spatial arrangement of intensities 69,70 .
Skin texture was characterized by measuring the likelihood of observing an intensity value at randomly selected pixels in the images (i.e. first-order statistics) and the spatial relationships between them (i.e. second-order statistics) 69,70 . From the frequency distribution of gray scales in the images, the following first-order statistics were calculated: mean ( µ ), variance ( σ 2 ), skewness (S), and kurtosis (K). Second-order textural features were extracted using the Gray Level Co-occurrence Matrix (GLCM) [70][71][72][73] . GLCM is a matrix representation of an image in which the gray tone of each pixel is quantized to a set (G) of N θ levels. Information on texture is specified by the matrix of relative frequencies P (i,j|d,θ) in which two neighboring pixels separated by distance d and angle θ occur on an image, one with gray tone i and another with j 71 . GLCMs of irregular ROIs were calculated from a distance between pixel (d) = 1 at four directions (θ) = 0°, 45°, 90°, and 135°, and then averaged for each image. In order to characterize the textural features of the dorsal images of the toads, we obtained the following five GLCM measurements: ii) Contrast (CON): measure of local intensity variation between a continuous set of pixels.
iii) Inverse difference moment (IDM): weighted by the inverse of the contrast ( i = j ) which is influenced by homogeneity, high values indicate homogeneous textures: Spatial heterogeneity (fractal dimension and lacunarity) analysis. Fractal dimension of grayscale images 76 was measured on ROIs applying differential box-counting methods using the plug-in FracLac in ImageJ 77 A series of boxes (samples) of decreasing scale (ε) were projected over a gray-scale image with an intensity I (i,j) scaled to a range of 0 to 255 for each pixel (i,j). For each box, the difference between the maximum and minimum values of pixel intensity (δI i,j,ε ) was calculated and the slope of the log-log regression between box size and the sum of all the intensity differences (Iε = ∑ [1 + δI i,j,ε ]) was used to calculate the gray fractal dimension ( D B ) . Fractal dimension constitutes a measure of spatial complexity [78][79][80] : Lacunarity (λ) is a measure of spatial heterogeneity that complements fractal dimension; whereas D B measures how much space is filled, lacunarity describes the spatial size of gaps and measures how space is filled 81 Thus, a fractal with large and heterogeneous gaps has high lacunarity values and vice versa: where σ stands for the standard deviation, μ for the mean difference in intensity for pixels per box at the scale ε and orientation g.
Statistical analyses. Variables were normalized using the bestNormalizer 82 package in R (version 3.6.2).
Fractal dimension (D B ) and entropy required no transformation, whereas mean (μ), variance (σ 2 ), contrast, and lacunarity (λ) were log-transformed. A Box Cox transformation was applied to angular second moment (ASM) and correlation (COR), and a square root to inverse difference moment (IDM). We then standardized each variable to a mean of 0 and standard deviation of 1. We then performed a principal component analysis (PCA, R version 3.6.2.) on the transformed variables, using factor loadings to assess the relative contribution of each variable to the variation observed in each of the principal components. Some of these variables were complementary to each other and showed high collinearity. Thus, multiple variables provided similar information regarding the frequency distribution of grey scales (mean, variance, skewness, kurtosis) whereas others described homogeneity in the pattern (ASM, IDM, entropy, contrast, correlation) and yet others spatial complexity (lacunarity, fractal dimension). We visualized the pairwise relationship among these variables using the ggpair function of the R package GGally 42 , and selected three variables of interest with low collinearity among them and loading on different principal components: mean, angular second moment, and fractal dimension. We then fitted multivariate Bayesian hierarchical linear models on these three variables using the R package brms 83 . We began with comparing model performance including and excluding sibship and population as random (group) effects, and included those in further analyses as models performed better with random effects. We compared model performance of fitting either on i) only water treatment, ii) water treatment and larval period, or iii) water treatment and body mass upon metamorphosis as fixed (population) effects, including sibship nested within population as random (group) effect for each model. All three response variables in the models were fitted within a gaussian distribution with default prior settings, with 152 observations (i.e. metamorphs) and 4 group levels (i.e. populations) for each model. Two sampling chains ran for 2000 iterations with a warm-up period of 1000 iterations for each model, hence yielding 2000 samples for each parameter. As diagnostic checks, chain convergence, autocorrelations and posterior predictive distributions were visually inspected and effective sample sizes calculated. We also checked for low pareto k heritage (< 0.7) for each observation within each model. To explore how much variation in the response variables were explained by our models, we employed a Bayesian generalization of the R 2 coefficient. Finally, the best model fit (i.e. the model with the highest predictive accuracy based on the expected log pointwise predictive density [ELPD]) was evaluated across models by estimating the out-of-sample predictive fit using the leave-one-out cross-validation (loo) and widely applicable information criterion (waic) criteria.