Variation in size and shape of toxin glands among cane toads from native-range and invasive populations

If optimal investment in anti-predator defences depends on predation risk, invading new regions (and thus, encountering different predators) may favour shifts in that investment. Cane toads offer an ideal system to test this prediction: expensive anti-predator toxins are stored mainly in parotoid glands whose dimensions are easy to measure, and toad invasions have changed the suites of predators they encounter. Although plasticity may influence parotoid morphology, comparisons between parents and progeny revealed that gland dimensions were highly heritable. That heritability supports the plausibility of an evolved basis to variation in gland dimensions. Measurements of 3779 adult toads show that females have larger glands than males, invasive populations have larger glands than in the native-range, and that parotoid sexual size dimorphism varies strongly among invaded areas. Geographic variation in parotoid morphology may be driven by predation risk to both adult toads and offspring (provisioned with toxins by their mother), with toxins allocated to eggs exacerbating the risk of cannibalism but reducing the risk of interspecific predation. Investment into chemical defences has evolved rapidly during the cane toad’s international diaspora, consistent with the hypothesis that organisms flexibly adjust resource allocation to anti-predator tactics in response to novel challenges.

history may affect optimal levels of investment into chemical defences. For example, a toxic toad in its native range is likely to be confronted with predators that either are unaffected by its toxic arsenal (small genetic changes to predator physiology confer many thousandfold resistance 20 ) or have evolved to exclude large toads from the diet 21,22 . Thus, increased investment into toxin production may confer little benefit to toad fitness in the native range. In contrast, greater investment into chemical defences may enhance survival in the invaded range, because toad-naïve predators cannot tolerate bufonid toxins 19 . In keeping with that hypothesis, interspecific comparisons within the Bufonidae show that possession of potent chemical defences is a significant correlate of range expansion 23 .
In the case of the cane toad, the invaded range in Australia contains many predators that are susceptible to the toxin, and willing to attack adult toads (especially at the invasion front 19 ), favouring increased allocation of resources to toxin production. Additionally, cane toads in many parts of Australia are highly dispersive, stimulating rapid evolution of morphological traits (such as relative limb lengths) associated with higher mobility 24 . In contrast, cane toads in Hawai'i are relatively sedentary 25 and the colonized area contains few species that are predators of adult toads 26 , favouring reduced investment. However, toxins do not protect only adult (terrestrialphase) toads; reproducing female toads allocate toxins to eggs to protect aquatic stages also 27,28 . As a result, the magnitude of toxin stores in reproductive female toads may comprise two components: that needed for defence of the adult toad, and that needed to provision her offspring.
Toads can change the magnitude of their investment into defensive chemicals via two mechanisms: adaptation and phenotypic plasticity. Under the former mechanism, the size of parotoid glands is expected to show significant heritability; gland sizes vary geographically because individuals with specific sizes of glands survive, and pass that trait on to their progeny. The alternative mechanism, phenotypic plasticity, requires an individual toad to modify its investment into toxin production based on circumstances that it encounters earlier in life. For example, a larval cane toad that is exposed to cues predicting high predation risk develops larger parotoid glands after metamorphosis 29 . Both canalised adaptation and plasticity can increase individual fitness, and can co-occur. Phenotypic plasticity in toxin investment is best demonstrated with experimental studies [30][31][32][33] , whereas the likelihood of longer-term selection on parotoid morphology can be clarified by quantifying heritability of this trait.
To test these ideas, we measured the sizes of parotoid glands in field-caught cane toads from the species' native range in South America, as well as invasive populations on three Hawai'ian Islands and in four regions of Australia. Because climates diverge between windward and leeward shores within each island, and the two sides are separated by inhospitable habitat (making it more difficult for toads to move from one side of an island to another), we treated the "wet" (windward) and "dry" (leeward) sides of each Hawai'ian island as a separate sampling unit. We also bred and raised toads from three Australian populations under standard ("common garden") conditions, to measure heritability of parotoid morphology.

Results
The relationship between ln parotoid area and ln SVL was linear, with SVL explaining 75.4% of the variance in ln parotoid area (t = 106.52, 1 df, P < 0.0001). Our measure of parotoid shape (width divided by length) was significantly (P < 0.0001) correlated with our measure of overall relative parotoid size (residual scores from the linear regression of ln parotoid area to ln SVL) but explained very little variance in the latter trait (r 2 = 0.006). Thus, we treat relative parotoid size and parotoid shape as separate dependent variables in the analyses below.
Variation in morphology of the parotoid macroglands among field-collected adult toads. In seven of the 11 regions that we sampled, female toads had larger parotoid glands than did males at the same body length. At three other sites dimorphism was minimal, and at one site (Maui wet-side) the dimorphism was reversed (larger parotoids in males than females: Fig. 2a). The glands of females were rounder in shape (greater width relative to length) than were those of males, at all sites (Fig. 2b). However, the magnitude of sexual dimorphism in parotoid morphology differed significantly among populations (interaction region*sex from ANOVA, for size F 10,3669 = 4.70, P < 0.001; for shape F 10,3681 = 3.06, P < 0.001). Thus, we examined geographic patterns separately in the two sexes, as well as in the combined dataset.
Parotoid size. Relative to body length, toads from the native range had smaller parotoids than did conspecifics from any of the invasive populations except those on the dry side of the large island of Hawai'i ( Fig. 2a). Parotoids were much larger in other sites within Hawai'i (especially in samples from the wet regions) and in the coolest-climate population in Australia (NSW: see Fig. 2a). Sexual dimorphism in parotoid size was low in the native range (2.3% difference in scores), high in Hawai'i (dry-side 9.2% divergence, wet-side 6.1%) and NSW (9.2%), variable in O'ahu (dry-side < 0.1%, wet-side 7.2%) and low in Maui (dry-side 2.0%, wet-side 3.7%) and in all tropical populations within Australia (0.2 to 2.2% in QLD, NT, WA: Fig. 2a). Because of the significant interaction between sex and region (above), we examined data separately by sex. For females, geographic variation in relative size of the parotoids was significant (F 10,37.63 = 2.29, P < 0.04), but with no significant posthoc Tukey tests (all comparisons P > 0.05). In males, relative parotoid size also differed significantly among regions (F 10,39.7 = 6.98, P < 0.0001; posthoc tests show that wet-side Maui and dry-side O'ahu toads differed significantly from Hawai'i dry-side conspecifics).
Parotoid shape. Broadly, sexual dimorphism in shape of the parotoids was low in the native range (6% difference in scores), variable in Hawai'i (3.3 to 8.3%), high in Queensland (7.4%), and decreased over the course of the Australian invasion, such that invasion-front populations exhibited similar sexual dimorphism in this trait as did native-range toads (e.g., 5.9% in NT, 5.2% in WA toads: Fig. 2b). Looking only at female toads (because of the significant sex*region interaction, see above), we found significant geographic variation in shape of the parotoids www.nature.com/scientificreports/ (F 10,41.24 = 5.53, P < 0.0001; posthoc tests show that parotoids were wider relative to length in toads from French Guiana than in most other groups (QLD, Hawai'i dry-side, O'ahu and Maui both sides), and Maui dry-side toads had significantly smaller parotoids than did those from the native range, NT, WA, NSW or Hawai'ian wet-side conspecifics. Shape variation was also significant in male toads (F 10,44.58 = 9.83, P < 0.0001), with parotoids significantly more rounded in French Guianan animals than in those from QLD, both sides of O'ahu and dry sides of Hawai'i and O'ahu. Toads from the latter three areas also had significantly smaller parotoids than did toads from WA and the NT. Across both sexes and in both wild-caught and captive-raised toads, parotoid glands were more rounded in individuals whose limbs were relatively short compared to the body, based on linear regression comparing our parotoid shape index to residual scores from the general linear regression of ln limb length vs. ln SVL. Although limb length explained relatively little variation in parotoid shape in some groups, the relationship was statistically significant for all comparisons (wild-caught females arm length N = 1804, r 2 = 0.01, P < 0.0001, leg length N = 1907, r 2 = 0.04, P < 0.0001; wild-caught males arm length N = 1907, r 2 = 0.04, P < 0.0001, leg length N = 1907, r 2 = 0.05, P < 0.0001; captive-raised females arm length N = 66, r 2 = 0.20, P < 0.0002, leg length N = 65, r 2 = 0.29, P < 0.0001; captive-raised males arm length N = 86, r 2 = 0.12, P < 0.0015, leg length N = 86, r 2 = 0.10, P < 0.0025).
Effect of toxin expulsion on dimensions of the parotoids. ANCOVA (with SVL as covariate to remove body size effects) indicated that when measured five days post-manipulation, parotoid size and shape were indistinguishable between toads from which we had manually expelled toxin from the glands compared to sham-manipulated controls (both F 1,32 < 1.42, both P > 0.24).

Discussion
Cane toads provide an excellent system in which to study the factors influencing investment into anti-predator defences: the toad's toxins are expensive to produce 16 , the magnitude of toxin stores is easily quantified 17 , and recent range expansions have modified the toad's rate of dispersal as well as its exposure to predation in both aquatic and terrestrial phases of the life-history 19 . Variation in dimensions of the parotoid macroglands among individuals is not attributable to recent expulsion of toxins (i.e., as measured over 5 days) and exhibits significant heritability (present study), consistent with the hypothesis that geographic variation in gland morphology reflects adaptive responses to novel challenges. Our sampling design does not allow us to determine whether the heritability of parotoid dimensions is underpinned by genetic versus epigenetic changes.
The most notable changes in shape of the parotoid glands are the evolution of a more elongated parotoid in invasive populations (relative to the situation in the native range) and a decrease in the degree of sexual dimorphism in shape as toads colonized the Australian tropics (Fig. 2b). The former shift might be due either to founder effects (genetic drift) or to novel selective forces. The latter shift mirrors evolutionary changes in relative head size (and the decline in sexual dimorphism in that trait) that has taken place in Australia, perhaps as an adaptation to increased rates of dispersal 34 . Consistent with that hypothesis, length of the limbs relative to the body (a trait that affects dispersal speed, and has evolved rapidly in invasive populations of cane toads [35][36][37] ) was linked to parotoid shape. A large rounded parotoid gland extends well down behind the shoulder, and thus may interfere with mobility of the forelimbs (see Fig. 1a). Under this interpretation, the evolution of greater mobility in invasive populations of cane toads was accompanied by elongation of the parotoid glands as part of a suite of morphological features related to enhanced mobility (see 36,37 for other traits). Such a shift might have evolved either through adaptation (a fitness benefit to higher rates of dispersal 38 ) or through spatial sorting (accumulation of dispersal-enhancing genes at an expanding range edge 39 ).
Geographic variation in relative size of the parotoids exhibits more complex patterns (Fig. 2a). Parotoids were larger in the invaded range than in the native range, and varied in size even among invaded sites within both Hawai'i and Australia. These results support Phillips and Shine's 40 demonstration of geographic variation in parotoid size within Australian toads, but do not reveal a strong cline associated with invasion history (contra the previous analysis).
Variation in parotoid size may be driven at least partly by phenotypic plasticity; exposure to predation cues during larval life stimulated recently metamorphosed cane toads in Australia to develop larger parotoid glands 29 . However, subsequent work in which we have raised toads to maturity revealed no significant difference in parotoid sizes of adult toads as a function of exposure to predation cues during larval life (Sharma et al., unpublished data), suggesting that impacts of larval experience on parotoid dimensions do not persist through to adulthood. High repeatability of parotoid measures in our study of captive-raised toads indicate that the morphology (size and shape) of the glands were consistent across an individual's lifetime. Additionally, we found significant heritability of parotoid sizes and shapes among offspring raised in standardised conditions (above), a result that would not be expected if gland morphology was sensitive to disruption by environmental factors (although such effects likely were minimised in our study, because offspring were raised under standardised conditions). Our estimates of heritability (26 to 32%) are higher than those we have calculated for other morphological traits of these toads 37 , but are similar to published estimates of heritability of morphological traits more generally 41 .
Even if variation in parotoid size is driven by heritable factors, it might not be adaptive. Biological invasions often involve successive founder effects, increasing the likelihood that gene frequencies will be affected by nonadaptive processes such as genetic drift 42  www.nature.com/scientificreports/ by a maximum of 150 individuals, and the Australian populations by a maximum of 101 animals 43 . However, existing populations both in Hawai'i and Australia contain thousands or millions of individuals, reducing the impact of random effects on gene frequencies. Thus, the consistent pattern for smaller and more rounded parotoids in the native range than in the invaded range (see Fig. 2a) suggests that geographic variation in parotoid size likely is driven at least partly by deterministic processes (i.e., adaptation). If variation in levels of investment into parotoid macroglands is indeed driven by natural selection, what selective forces are likely to have been important? The evidence that these glands function to deter predator attacks is conclusive (see above), but the parotoids might have other functions as well. For example, the glands contain hydrophilous glycosaminoglycans that could provide a store of water during long dry periods 30,44 . We doubt the importance of hydroregulation, however, because the two regions where parotoids were largest (Maui wet-side and New South Wales; Fig. 2a) experience moister climates than most of the other regions from which toads were sampled (see 45 for climatic data). Thus, the most likely selective force driving geographic variation in size of the parotoid glands involves defence against predators.
We first consider the case of males, because their investment into toxins should be driven only by vulnerability to predators of adult toads (whereas investment by females reflects the need for additional allocation to the offspring). The overall pattern for males is straightforward, with low investment in the native range and on the dry sides of Hawai'ian islands, and higher everywhere else (Fig. 2a). Increased parotoid size in Australian populations is consistent with the putatively higher predation pressure in these areas (see above), but available information about predation on cane toads in Hawai'i is too meagre to suggest an explanation for the higher investment into parotoids on wet-side Maui populations.
A female toad's investment into toxin stores should incorporate two components: the amount needed for her own defence against predation, and the toxins that she needs to allocate to her eggs 28 . However, the optimal allocation of toxins to eggs is complicated by geographic variation in the incidence of cannibalism, which is frequent in Australian cane toads 46,47 but is infrequent in the native range (DeVore et al., unpublished data). Cannibalistic conspecific tadpoles locate eggs by detecting toxins exuded during late stages of egg development 46 ; and thus, a higher investment of toxins into the eggs may render a female's offspring more vulnerable to cannibalistic attack by conspecifics but less vulnerable to predation by other taxa such as fish and turtles 48,49 .
The disparity between parotoid sizes in female versus male cane toads offers an approximate index of the allocation of toxins to eggs. That disparity is low in the native range (where local predators have evolved to tolerate bufonid toxins, reducing the effectiveness of toxins as a defence) and in Queensland (where larvae are highly cannibalistic, and most clutches are laid in waterbodies that already contain such larvae). However, females from Maui and the dry side of O'ahu also have relatively small parotoids (Fig. 2a), for no obvious reason. In contrast, the sexual disparity in favour of females is maximal in New South Wales, where cannibalism may be infrequent at the invasion front (due to low densities of toads) and embryonic development is slow (due to low water temperatures), increasing the duration of time for which eggs are vulnerable to other types of predators. Field studies are needed to identify selective forces on optimal allocation of toxins to eggs, by quantifying rates of predation (including cannibalism) of toad eggs across the cane toads' geographic range.
Throughout the manuscript, we have considered cane toad parotoid gland morphology as a proxy for toxin production, and hence anti-predator defence capabilities. While it is likely that toads with larger glands produce more toxin by volume (and invest more energy in doing so), it is possible that divergent selective forces are operating on toxicity or chemical composition of the macrogland secretions as well. For example, in environments www.nature.com/scientificreports/ where toads frequently encounter naïve predators (i.e., the invaded range), they may require greater volumes of toxin if they exude often as a deterrent in response to harassment. Conversely, in habitats where predators are resistant to toad toxins 20 , as is expected with co-evolved predators in the native-range, selection may act to increase toxicity or alter the composition of parotoid exudate, leading to reciprocal adaptation between predator and prey. Previous work has shown that local adaptation to predation regime or climatic factors can influence both the chemical composition of amphibian toxins, and their potency 50,51 . Therefore, comparisons of the parotoid secretions of cane toads throughout their native and invasive ranges would be an illuminating avenue for future research, and may provide a more nuanced picture of the evolution of anti-predator investments during colonization of new habitats. In summary, invasive species provide excellent models for studying anti-predator adaptations; the invader arrives with phenotypes that have evolved elsewhere, and may experience different selective forces in colonized habitats than in their native range. The fitness benefits of investment into anti-predator traits can vary greatly in novel environments, generating selection to fine-tune levels of investment depending on cost-benefit trade-offs in local environments. By comparing introduced populations to those within the native range, we can study phenotypic adaptations to geographically variable predation regimes, and explore rapid evolution in response to novel selection pressures. Research on invasive species also allows us to explore the impacts of range spread www.nature.com/scientificreports/ on phenotypic traits, with investment into anti-predator traits likely trading off with dispersal rate as well as with protection against predators. Overall, strong geographic variation in size, shape and sexual dimorphism of the parotoid macroglands in cane toads supports the hypothesis that investment into defence against predation is fashioned by a complex interplay among selective forces, acting across multiple life-history stages, and that shifts in selective forces can rapidly change investment optima, and thus modify heritable variation in allocation of resources towards anti-predator adaptations.

Materials and methods
Study system. Many bufonids ("true toads") possess paired postorbital macroglands known as parotoids 18,52 that produce cardiotonic steroids such as bufogenins and bufotoxins 50,[53][54][55] . These toxins have cardio-acceleratory and vasopressor effects on vertebrates 52,56 , such that ingestion can be lethal for predators [57][58][59] . In cane toads (Rhinella marina), this toxin is contained both within the parotoids and a series of smaller glands that are distributed on the dorsal surface and limbs 17 . The cane toad uses the toxin as a passive defence, exuding it when stressed (Fig. 1); predators also may be poisoned by biting into the glands. The cane toad was introduced from French Guiana to Puerto Rico in the early 1900s, and from there to Hawai'i in 1932, and from there to Australia in 1935 (see 60 for genetic evidence of this translocation route). Since these introductions, the cane toad has become widely distributed across Hawai'i 25 and Australia 43,61,62 . Invertebrate predators of eggs and larvae are common in both countries, but predation on adult toads is higher in Australia than in Hawai'i 19,25 . We treated the "wet" and "dry" sides of the Hawai'ian islands as separate sampling units, because the climatic divergence between windward and leeward shores within an island is greater than inter-island differences in climatic regimes 25,63 ; sample sizes Hawai'i dry-side N = 286, wet-side 244; Maui dry-side N = 95, wet-side 66; O'ahu dry-side N = 126, wet-side 200). Toads were captured by hand, and we used Vernier callipers (± 0.1 mm) to measure the SVL of the toad, and the length (PL) and width (PW) of its right parotoid macrogland. We also measured limb lengths with a plastic ruler (see 36,37 for details). Sex was determined by secondary sexual characteristics 64 .

Effect of toxin expulsion on dimensions of the parotoids. If expulsion of toxin alters the size or
shape of the parotoids, measurements of these glands on field-caught toads might be affected by local rates of predator attack. Such an effect could introduce significant variation into data for any given population, weakening our ability to detect broader patterns. As part of another study, we manually expressed toxins from the parotoid macroglands of 16 toads, and sham-manipulated another 20 animals to serve as controls. The two groups were matched for overall body sizes. Following toxin removal treatment, toads were radio-tracked for 5 days and then recaptured and their glands and body size measured (see 16

for details).
Assessment of heritability of parotoid morphology. Australian toads captured from the wild were brought to our field station 66 km southeast of Darwin, NT (12°37′S 131°18′E) where they were bred. We raised the resultant offspring under standard conditions, to reduce the influence of environmental factors on morphology. The adult toads were collected from three sites in north-eastern Queensland (Townsville, 19°26′S 146°82′E, Innisfail, 17°52′S 146°03′E, and Tully, 17°93′S 145°92′E) and four sites in northern Western Australia (El Questro, 16°00′S 127°98′E, Purnululu, 17°53′S 128°40′E, Wyndham, 15°46′S 128°10′E, and Oombulgurri, 15°18′S 127°84′E). Toads had been present in all of the QLD sites for > 70 years, but had only recently colonized the WA sites (< 2 years 65 ). We induced spawning by injection of leuprorelin acetate (Lucrin; Abbott Australasia; using 1 mL of Lucrin diluted 1:20 with saline) and raised the progeny in captivity using the protocols described by 66 . After post-metamorphic toads attained body lengths > 20 mm, we toe-clipped them for identification (a procedure that causes minimal stress to the animals 67 ) and kept them in outdoor enclosures in groups of 30 (with mixed parental origins). We raised "common garden" toads from 31 egg clutches (16 QLD, 15 WA) totalling 489 offspring (287 QLD, 202 WA). The offspring were measured three times (at approximately 8, 14, and 17 months of age); 196 animals (132 from QLD, 64 from WA) had reached adult size (> 90 mm SVL) by the end of this study. Our heritability analyses are based on data from 61 parents and 317 offspring > 60 mm SVL, each measured one to three times.

Calculations and statistical analyses. Morphology of the parotoid glands.
To obtain an estimate of overall size of the glands, we calculated an approximate area by treating parotoids as rectangular in shape (length*width). That method slightly overestimates area for rounded shapes but the error is small, and is similar among all samples (calculating the area as a circle, based on mean radius, yields a score that is almost perfectly correlated with the rectangular approximation: r 2 = 0.997). We then ln-transformed values of parotoid area and body length to achieve normality, and regressed ln parotoid area against ln body length to obtain residual scores (deviations from the average parotoid size expected at any given body length) that we used as indices of relative parotoid size (i.e., a high score represents a larger-than-usual parotoid for an animal of that body length). To quantify shape of the parotoids, we divided the width measurement by the length measurement; thus, high scores indicate more rounded shapes. We did not correct for allometry in this measure because our index of parotoid shape was not strongly related to body length in either sex (r 2 = 0.006 for females, r 2 < 0.0001 for males). www.nature.com/scientificreports/ For analyses of geographic variation, we recognized 11 major regions (although we retained information on collection sites within each region, and used this as a random factor in our statistical analyses to avoid pseudoreplication): (1) all French Guiana sites; (2) Wet (windward) sides of each of three Hawai'ian islands (Hawai'i, O'ahu, Maui); (3) Dry (leeward, rain shadow) sides of each of the same three islands; and four states within Australia: (4) Queensland, (5) New South Wales, (6) Northern Territory, and (7) Western Australia. Table 1 shows sample sizes and associated information for all 54 collection sites. We compared parotoid dimensions between sexes and among regions using a two-factor ANOVA with sex and region as factors, and population as a random variable (to avoid treating toads from the same collection site as statistically independent). All other analyses were performed using JMP 14 software (SAS Institute, Cary, NC). We assessed residuals from all analyses to detect violations of assumptions.
Common-garden offspring. We measured the animals in the same way as for field-collected specimens, and calculated the same variables as above to describe parotoid size and shape. We excluded measurements of individuals < 60 mm SVL (because of imprecision in measuring such small animals). To estimate heritability and repeatability of size and shape of the parotoids, we used ASREML software (VSN International Ltd., Hemel Hempstead, UK) to run an animal model 68 incorporating individual ID and family ID as random effects.

Ethics statement.
All procedures in the current study were approved by the University of Sydney Animal www.nature.com/scientificreports/