The impact of multiple climatic and geographic factors on the chemical defences of Asian toads (Bufo gargarizans Cantor)

Chemical defences are widespread in nature, yet we know little about whether and how climatic and geographic factors affect their evolution. In this study, we investigated the natural variation in the concentration and composition of the main bufogenin toxin in adult Asian toads (Bufo gargarizans Cantor) captured in twenty-two regions. Moreover, we explored the relative importance of eight climatic factors (average temperature, maximum temperature, minimum temperature, average relative humidity, 20–20 time precipitation, maximum continuous precipitation, maximum ground temperature, and minimum ground temperature) in regulating toxin production. We found that compared to toads captured from central and southwestern China, toads from eastern China secreted higher concentrations of cinobufagin (CBG) and resibufogenin (RBG) but lower concentrations of telocinobufagin (TBG) and cinobufotalin (CFL). All 8 climatic variables had significant effects on bufogenin production (ri>0.5), while the plastic response of bufogenin toxin to various climate factors was highly variable. The most important climatic driver of total bufogenin production was precipitation: the bufogenin concentration increased with increasing precipitation. This study indicated that the evolution of phenotypic plasticity in chemical defences may depend at least partly on the geographic variation of defensive toxins and their climatic context.

50 different bufagenin compounds have been identified in Asian toads to date 24,25 . Among them, gamabufotalin (GB), telocinobufagin (TBG), bufotalin (BFL), cinobufotalin (CFL), bufalin (BL), cinobufagin (CBG) and resibufogenin (RBG) (Fig. 1) are the dominant components, representing more than 80% of the total amount of bufagenins 26 . These compounds are considered major contributors to the role of toxins in defence 27 . Despite similar chemical structures, the different bufagenins may exhibit functional differences in the potency and selectivity of their toxicity to different types of predators 20,28 . For example, many arthropod Na + /K + -ATPases lack the binding site for certain bufagenins, whereas certain vertebrate Na + /K + -ATPases have high affinity for the same bufagenins 29,30 . Therefore, the concentrations and composition of such highly toxic components in the toxin blend are more relevant to the deterrent efficiency of the chemical defence than the amount of toxin 20 .
In this study, we aimed to explore the natural variation in chemical defences and its relationships with assorted climatic and geographic factors in adult Asian toads. We described the variation in chemical defences with three variables: the concentrations of single bufogenins, the total bufogenin concentration (the sum of the concentrations of the main seven bufogenins) and the constituent ratio of bufogenin compounds (relative proportions of various bufogenins). In addition, we investigated the extent of variation in these defensive toxins from different geographical regions and explored the relative importance of the influence of various climatic factors on toxin variation. Our study confirms that climatic factors may substantially regulate bufogenin plasticity and affect the particular geospatial distribution of these defensive chemicals. To our knowledge, this study is the first to determine the variation in chemical defences among Asian toads in different geographical regions and its relationship to climatic factors.

Results
Constituent ratios and correlation of the bufogenins. The total bufogenin concentrations at the twenty-two test sites ranged from 8.43% to 19.34% (Table 1). A percentage stacked column chart of the various component concentrations is shown in Fig. 2. Clearly, the constituent ratios of the seven bufogenins from the twenty-two sites were significantly different. The highest constituent ratios of GB, TBG, BFL, CFL, BL and CBG were more than 4.0, 56.8, 5.6, 13.7, 5.4 and 31.6 times the lowest ones, respectively. However, the RBG concentration was too low to detect in Zhangjiajie, Hunan Province ( Table 1).
The concentrations of various bufogenins were related to each other (Table 2). Generally, with increasing CBG concentration, the concentrations of both BL and RBG significantly increased, but the TBG and CFL concentrations decreased.
Geospatial distribution characteristics of the concentrations and proportions of bufogenins. The variation in bufogenin concentrations was significantly related to geographic longitude (p < 0.01) but not to latitude (p > 0.05). Specifically, the concentrations of TBG and CFL were negatively correlated with longitude (r = −0.787, p < 0.01 and r = −0.836, p < 0.01, respectively), while conversely, the concentrations of BL, CBG and RBG were positively correlated with longitude (r = 0.708, p < 0.01; r = 0.685, p < 0.01 and r = 0.851, p < 0.01, respectively). The total bufogenin concentration was positively correlated with longitude (r = 0.560, p < 0.01).
The composition and proportions of defensive chemicals can markedly affect deterrent effects 20 . This study investigated the geospatial distribution characteristics of the GB/CBG, TBG/CBG, BFL/CBG, CFL/CBG, BL/ CBG and RBG/CBG ratios, as shown in Fig. 3. With increasing longitude, the line chart of all the ratio values     www.nature.com/scientificreports www.nature.com/scientificreports/ except the RBG/CBG presented an increasing trend at the southwestern China sites, and the line chart of all the ratio values displayed cross-fluctuation trends at the central and eastern China sites. However, except for RBG/ CBG, the fluctuation trends of all ratios were consistent at the central China sites, and only the fluctuation trends of BFL/CBG and BL/CBG were identical at the eastern China sites. Moreover, TBG/CBG and CFL/CBG were always significantly higher than GB/CBG, BL/CBG and RBG/CBG at the southwestern and central China sites. In addition, at the twelve eastern China sites, except for PX, SR and TX, all six ratio values from different sites were less than 1, which indicated that CBG was the richest component among the bufogenins in the eastern regions.
Geographic clustering of toad toxins based on chemometrics. The principal component analysis (PCA) score plot of the toad toxin samples (Fig. 4) showed a clear separation among the samples from different geographical origins, with the first two components explaining 75.75% of the total variance. Based on their chemical profiles, the toad toxin samples were classified into two main clusters (I and II) encompassing the twenty-two geographical origins: 1) C2-C6, SW1-SW4 and E8 were cluster I, and 2) C1, E1-E7 and E9-E12 were cluster II. The toxin samples C1-C6 and SW1-SW4 were obtained from central and southwestern China, and samples E1-E12 were obtained from eastern China. From a geographical perspective, C1 was obtained from Anyang, which is located in north of the Huaihe River, in the northernmost part of Henan Province; E8 was collected from Pingxiang, located in the western part of Jiangxi near Hunan Province. The concentrations of bufogenins from these two locations (Anyang and Pingxiang) were similar to those in the adjacent regions but different from those in other regions of the provinces.
The classification results based on the hierarchical clustering heatmap (HCH) were consistent with the PCA results (Fig. 5). The distributions of the seven bufogenins in twenty-two geographical origins was clearly displayed by the heatmap. The distribution characteristics of bufogenins in different clusters were distinctly different. The concentrations of TBG and CFL were significantly higher in cluster I than in cluster II, whereas the concentrations of CBG and RBG were higher in cluster II than in cluster I. In six out of ten samples in cluster I, CFL constituted the largest proportion, ranging between 26.04% and 40.41%. However, CBG was the predominant active constituent in cluster II except in E5 and E9, ranging from 25.05% to 41.75%.
The results of orthogonal partial least squares discriminant analysis (OPLS-DA) were consistent with those of PCA and HCH. The OPLS-DA score plot (see Fig. 6) could also be divided into two zones according to the twenty-two sites with R2Y = 0.871 and Q2 = 0.831, demonstrating the good classification and prediction ability of the model.  www.nature.com/scientificreports www.nature.com/scientificreports/ The effects of climatic factors on bufogenins. All 8 climatic variables had significant effects on bufogenin concentration (r i > 0.5; Table 3). Among them, PRE_Time and PRE_Max showed the most significant effect on bufogenin accumulation, followed by RHU_Avg. The results from January 2018 (r i -Jan.) were essentially in agreement with those from February 2018 (r i -Feb.).
Based on the coefficient plots of the OPLS model (Fig. 7A,B), the extent and direction of the influence of the eight climatic factors on the seven bufogenins were significantly different. The concentrations of GB and BFL were  www.nature.com/scientificreports www.nature.com/scientificreports/ minimally affected by the eight climatic factors (|r| < 0.2). In the other five bufogenins, because of the correlation between the components, the directions of the influence of the climatic factors on TBG and CFL were the same. The directions of influence on BL, CBG and RBG were also the same, but the intensity of the influence was not completely consistent in our study. Moreover, these two groups of toxic ingredients were affected by climatic factors in opposite directions. Specifically, PRE_Time and PRE_Max were both negatively associated with TBG and CFL but positively correlated with BL, CBG and RBG, and the effect on total bufogenin concentration was consistent with the effects on the latter three components. That is, the total bufogenin concentration increased with increasing precipitation. The main climatic factors influencing the concentration were PRE_Time and PRE_Max, which had variable importance in the projection (VIP) values greater than 1 (Fig. 7C,D). These results were consistent with those of the GRA.

Discussion
We found that the concentrations and proportions of the main defensive chemicals in adult Asian toads from different geographical regions were highly variable. Compared to toads from central and southwestern China, the toads from eastern China had higher concentrations of CBG and RBG but lower concentrations of TBG and CFL. This variation in bufogenin concentrations showed a specific spatial structure related to the origins of the bufogenins, and the degree of variation increased significantly with increasing geographical distance, suggesting that toads with different bufogenin profiles are nonrandomly distributed among geographical regions.
Interestingly, despite the geographical variation in chemical defence, the concentrations of the various compounds in the toxin blend were related to each other. Therefore, we hypothesized that positive or negative correlation of various toxin components result from constitutive defences in toads, whereas inducible plasticity determines the dominant compounds of the defensive toxin blend. In terms of why the dominant compounds in the toxin blend from different habitats were significantly different, we argue that environmental factors in different regions may regulate the synthesis of these compounds. Many studies have demonstrated that the biosynthesis of the bufogenin component starts with common precursor molecules 20,31,32 , and the negative correlation between  Table 3. The grey relational grades r i between climatic factors and the total bufogenin concentration. www.nature.com/scientificreports www.nature.com/scientificreports/ the production levels of certain toxins may be the result of competition for these limited precursors. Additionally, the synthesis of these compounds may involve different biosynthesis pathways, which may be selectively accessed in different environmental conditions. Some studies have indicated that bufogenins may also be transformed by certain microorganisms 33,34 , and the correlations between different chemical components may result from biotransformation of the toad skin microbiome in different habits.
In addition to elucidating the variation in chemical defences, the relationships between defensive chemicals and environmental factors are also relevant for understanding the consequences of this variation. As the mixed defensive toxins were obtained from different adult toads at each geographical site, this regular pattern of chemical defence phenotypic plasticity induced by geographical factors was not related to the gender, age, or size of the individual toads. Although predators and competitors in different regions may have certain effects on the variation in defensive toxin production 8,12,34 , the density and type of these species in different regions vary greatly 7,34,35 . Accordingly, we conclude that this regular variation in geographical environment was not mainly induced by the presence of these species.
The variation may be primarily attributable to the climatic factors in different geographic locations. We found that the relationships between climate variables and bufogenin concentration, despite various climatic drivers, exert different influences on toxin production. Our results also showed that the total bufogenin concentration was positively correlated with precipitation, demonstrating that toxin production exhibits climate-induced plasticity.
Although the biosynthesis pathway of these defensive compounds is unknown 32 , previous studies found that the amount of biosynthesis-related energy absorbed by toads mainly depends on weather conditions as modified by habitat structure 36 . The habitats chosen by toads will determine their body temperature and hydration state, which can have direct physiological, behavioural and functional consequences [37][38][39] . As an amphibian, toads rely particularly on humid environments. In addition, toads obtain water primarily by osmotic absorption across their integument under a favourable osmotic gradient because moist skin is highly permeable 40 . The highly vascularized abdominal skin contains a seat patch or pelvic patch 41,42 , which is specialized for rapid water absorption to restore plasma volume and osmolality 40 . Therefore, the precipitation and humidity of the external environment will directly affect the water balance in toads, which in turn will have a certain impact on the metabolism. Previous studies found that toad tadpoles tended to contain higher concentrations of defensive toxin in ponds that were less likely to dry out 7 . As an alternative explanation for the positive correlation between total bufogenin concentration and precipitation, it is possible that toads retain the trait of tadpoles by producing more bufogenins in a humid environment.
Climate-induced chemical defence plasticity may have major ecological consequences. First, the defensive chemicals produced by animals can provide protection against natural enemies 11,15 , and the composition of the defensive toxin mixture may substantially affect its function or deterrence efficiency 20 . Such climate-induced plasticity may reflect a diversity of chemical defence strategies, thereby possibly enhancing anti-predator efficacy. For example, individuals with certain compositions of toxic chemicals may defend against enemies with corresponding toxin susceptibilities 28 . However, the survival probability of animals may decrease if toxin production is reduced by unfavourable climate conditions at a particular time. Based on our current results, extreme drought may lead to a reduction in the main toxin concentration, and prey species exposed to the environment may increase their predation risk. Recent studies have suggested that most species are facing a high risk of extinction 18,43,44 , and climate change is an important factor affecting the survival of various species [45][46][47] . Relatively little is known about the demographic mechanisms through which climatic variation affects population dynamics in species 48 . Our results show that climate might indirectly and partially reflect the species survival probability by regulating chemical defences. In addition, climate-induced the plasticity of chemical defence may impact the fitness of species. A previous experimental study showed that the synthesis of defensive chemicals imposes a substantial cost related to energetic trade-offs between compound production and reproduction 31 . Animals that invest more in toxin production do so at the expense of reproductive investment 31 . Further testing is needed to investigate whether climate-induced high levels of toxin production will reduce the reproductive success of both males and females.
In conclusion, our results are the first to document that toads phenotypically adjust their chemical defences in response to variable climate and geographical conditions. The observation that toads produced a high concentration of the defensive toxin in regions with considerable precipitation indicated that bufogenin production may be related to water in the environment. These results, combined with those of previous studies 5,7,10 , suggest the existence of multiple approaches to inducible chemical defences. Studies on the environmental factors regulating toxin synthesis and the costs and mechanisms of toxin production, clarifying the deterrent efficacy of toxin mixtures constituted by diverse compositions in facing different environmental risks will be favoured for further understanding of the evolution and ecology of chemical defences.

Materials and Methods
Study species and sites. The study species was Asian toads (Bufo gargarizans), the most common toad in The study was further approved by the Ethics Review Committee of Experimental Animal Welfare, Zhejiang University. All methods were carried out in accordance with relevant guidelines and regulations.
The parotid glands of the toads were scraped moderately with an aluminium scraper. The milk-like fresh secretions obtained were placed into clean glass dishes. After this procedure, the wounds of the animals were gently wiped with erythromycin ointment to reduce inflammation. Then, the animals were returned to the land. The

Preparation of standard solutions.
A stock solution of seven standards was prepared by dissolving accurately weighed standards in methanol. Working standard solutions for calibration curves were obtained by diluting the mixed standard solution with methanol to obtain a series of different concentrations of these analytes. All solutions were stored in a refrigerator at 4 °C.
Preparation of sample solutions. 0.2 g fresh secretions sample obtained from the parotoid glands of toad was accurately weighted and extracted with 20 mL methanol through heating reflux for one hour. The extraction solutions were cooled at room temperature and then the lost solvent were replenished with methanol. All of the samples solutions were centrifuged at 4,000 rpm for 5 min and then the supernatant was filtered through a 0.45 μm filter membrane before HPLC analysis.
Chromatographic analysis. Methanol extracts of samples were prepared as described above. The bufogenins concentration were determined by the method established previously 26    www.nature.com/scientificreports www.nature.com/scientificreports/ climatic factors were obtained from the China Meteorological Administration. Given that the biosynthesis of defensive chemicals takes time 31 , we selected daily climate data for the first two months before the toxin samples collection (January and February 2018) and calculated the monthly average values of the eight climatic factors at the twenty-two geographical sites (Table S1). Data analysis. Climate system can usually be considered as grey system due to their high complexity and unpredictable change. It's considerably difficult to explore the combined influence of multiple climatic factors on toxic components. In view of this, this paper utilized grey relational analysis (GRA) and multivariate analysis (MVA). GRA, as a powerful tool dealing with problems of small samples and imprecise information, has been widely applied in many fields 49,50 . But the correlation coefficient in the GRA can only indicate the influence degree of the compared series to the reference series, and cannot reflect the influence direction (positive or negative) 51 . Therefore, the influence direction of various climatic factors on bufogenin accumulation were further explored by MVA. MVA can re-establish a new set of features by extracting highly correlated projections of data representations from input and output spaces. MVA includes PCA, OPLS-DA, and OPLS [52][53][54] . PCA is an unsupervised method and is very favourable for preliminary classification 54 . OPLS is a popular supervised method and is being increasingly adopted as an alternative to partial least squares regression (PLS) due to better generalization 55 . The multiple climatic factors and variety of bufogenin concentrations, as high-dimensional variables, can be analysed by these methods.

Grey relational analysis.
To evaluate the importance of various climatic factors on the concentration of bufogenins, independent input variables, i.e.TEM_Avg, TEM_Max, TEM_Min, RHU_Avg, PRE_Time, PRE_ Max, GST_Max and GST_Min, were chosen as the compared series and defined as: i i The output variable (the total bufogenins concentration) was set as the reference series and expressed as: where m is the total number of indexes to be considered, and n is the total number of samples (m = 8, n = 22 in this experiment). Due to the different dimensions of the original data, it is necessary to standardize the original data. In this experiment, raw data of multiple climatic factors and the total concentration of seven bufogenins were normalized by SPSS 20.0 software for Z-Score standardization to obtain a series of dimensionless normalized data.
The grey relational coefficients reflects the relationship between climatic factors and the total bufogenins concentration. The calculation formula of grey relational coefficients is: are the standardization value of X 0 (k) and X i (k), respectively. ρ (0 < ρ ≤ 1) is a distinguishing coefficient to adjust the range of the comparison environment, normally used value is ρ = 0.5.
Finally, the value of the grey correlation degree were obtained by calculating the mean values of all the grey relational coefficients.
i k n i 1 Above steps were carried out by Matlab R2018a.

Statistical analyses.
All experiments were performed in triplicate, and the values are expressed as the mean ± standard deviation. The raw data on multiple climatic factors and the bufogenin concentrations were normalized by Z-score standardization to obtain a series of dimensionless normalized data. The obtained data were analysed by one-way ANOVA of Duncan's multiple comparison test and Spearman's two-tailed correlation analysis using SPSS 20.0 software. The level of significance for all hypothesis tests (p) was 0.05. The bufogenin concentrations were used as the variables, and the toxin samples were the observation IDs. The obtained data were imported into SPSS 20.0 software for PCA to reveal the intrinsic relations in the data set and to diagnose any possible outlier. In the HCH analysis, the data were processed based on the average linkages using HemI (Heatmap Illustrator, version 1.0). Then, OPLS-DA was performed to evaluate the robustness and accuracy of the classification results using SIMCA-P version 14.0. Finally, the relationships between climatic factors and bufogenin concentration were evaluated by the regression coefficient and VIP of the OPLS. The regression coefficient can reflect the positive or negative effect of the X variables (eight climatic factors) on the Y variables (the concentrations of single bufogenins and the total bufogenin concentration). In general, VIP values of the X variables that are greater than 1 are considered decisive indicators with significant influence on the Y variables 56 .