Evolutionary and ecological processes influencing chemical defense variation in an aposematic and mimetic Heliconius butterfly

Study sites and field collection of butterflies and host plants

Heliconius erato demophoon and its typical host plant P. biflora were collected in Panama along a steep rainfall gradient (the Smithsonian Tropical Research Institute Physical Monitoring Program; https://biogeodb.stri.si.edu/physical_monitoring/) from three populations (Wet, Intermediate and Dry study sites; Fig. 1A, File S4), and H. erato lativitta from High and Low altitude sites on the Eastern slope of the Andes in Ecuador (Fig. 1A, File S4). Research and collecting permits were granted by the Ministerio de Ambiente, Republica de Panama (Permit SE/AP-21-16) and the Ministerio del Ambiente, Ecuador (under the Contrato Marco MAE-DNB-CM-2017-0058). Butterflies were caught using a hand net and transported live to the laboratory. All individuals caught from the Panamanian study areas (n = 92, of which 34, 28 and 30 individuals were caught from Dry, Intermediate and Wet sites, respectively) were separated by sex, weighed, and the body excluding wings and one half of the thorax was preserved in 1 ml 100% MeOH (methanolic preservation stabilizes cyanogenic glucosides (Gleadow et al., 2011), and was used with both butterfly and plant samples) on the day of capture (males) or after being allowed to oviposit in an insectary (females, see File S4). At the Ecuadorian sites, fertilized females were brought to a common garden environment and kept in separate cages. Eggs were collected daily, and larvae raised in individual containers all through development, placed in randomly assorted positions in common garden insectary conditions. The thorax of 21 F1 individuals from each altitude (total n = 42) were sampled in 1 ml 100% MeOH. We do not expect the sampled body part to influence estimates of cyanogenic glucoside compound concentrations, as these compounds are located in the hemolymph, and are thus expected to be evenly distributed throughout the body (de castro et al., 2020).

Passiflora biflora leaf samples (one to four leaves pooled per plant individual, see File S4) were collected from the Panama Dry (n = 20), Intermediate (n = 11) and Wet (n = 20) study areas (Fig. 1A), and preserved fresh in 1ml 100% MeOH. In addition, 50 cm cuttings, all from different P. biflora individuals, were collected from the Dry (n = 12) and Wet (n = 12) study areas for greenhouse cultivation.

Host plant greenhouse cultivation and treatments

Passiflora biflora were cultivated in greenhouses (75% relative humidity; 8-20 h: 25 °C, light; 20-8 h: 20 °C, dark) in a soil mixture (50% compost, 20% coir, 15% perlite, 15% sand or gravel) (see File S4 for details on cultivation methods and watering treatments). Standard (std)-treated P. biflora were used as oviposition plants for the H. erato parental generation and larval diet for the F1 families not included in feeding treatments (see below). The quality of host plants in H. erato feeding treatments was manipulated with watering treatments, motivated by the importance of water availability for many cyanogenic plants (Gleadow & Woodrow, 2002; Hayden & Parker, 2002; Liang, 2003; Vandegeer et al., 2013). The P. biflora cuttings collected from the Dry and Wet study areas in Panama were cultivated in a full-factorial combination of origins and water treatments, resulting in four host plant types. Plants in the dry treatment were kept at 75% relative humidity and watered two to three times weekly with about 1/3 of the water volume compared to the wet treatment (constant water availability, relative humidity up to 100%), such that the surface soil was allowed to dry and the leaves slightly drooped between watering events. The dry treatment thus aimed to mimic conditions of drought stress.

Butterfly rearing, breeding design and feeding treatments

Laboratory H. erato populations descendent of the wild Dry population were established in greenhouse conditions (see Supplementary File 4 for details on rearing methods). Larvae were reared on std-type P. biflora. All butterflies were individually marked, and were allowed to fly, mate, and oviposit freely in 2 ×  2 × 2 m mesh cages with potted std-type P. biflora and a standardized diet of sugar solution supplemented with amino acids provided ad libitum. The butterfly cages were observed throughout the day for mating pairs, and mated females were moved to individual 2 ×  2 × 2 m mesh cages and allowed to oviposit on std-type P. biflora. Eggs of 14 females (eight mated with a known father, average no. of full-sibs/family = 8, see File S2) were collected, and hatched larvae were reared in full-sib groups fed on std-type P. biflora until pupation.

The effect of larval dietary cyanogen content and quality on adult butterfly biosynthesized toxicity was studied with feeding treatments. We expected larval dietary quality to affect butterfly toxicity due to potential variation in the acquisition of nitrogenous resources for cyanogen biosynthesis (Cardoso & Gilbert, 2013), possible detoxification costs (Lindstedt et al., 2010; Reudler et al., 2015) and the general importance of larval diet for fitness traits in holometabolous insects (e.g., Boggs & Freeman, 2005). The offspring of six mated pairs (both parental identities known, average no. of full-sibs/family = 29, File S2) were included in the feeding treatments, where eggs from each mother were divided equally into four groups and transferred onto each of the four P. biflora treatment types. Individuals of the F1 generation were collected within three days following emergence (unmated and unfed), weighed, and preserved in 100% MeOH (body excluding wings and one-half of thorax).

Butterfly toxicity analyses with ¹H-NMR

The concentrations of the two cyanogenic compounds biosynthesized by Heliconius larvae and adults, linamarin and lotaustralin, were analyzed from the butterfly samples using nuclear magnetic resonance (¹H-NMR) with Bruker Avance III HD NMR spectrometer (Bruker BioSpin, Germany) (see File S4 for details on sample extraction and ¹H-NMR methods). The data were processed and analyzed using Bruker TopSpin software (versions 3.2.6; 2014 and 3.5.6; 2016). The peaks of cyanogen compounds were identified based on comparison of sample ¹H spectra with the ¹H spectra of authentic reference samples of linamarin and lotaustralin. Note that we analyzed only the biosynthesized cyanogenic compounds (not sequestered cyanogenic compounds which differ in chemical structure from the two biosynthesized compounds linamarin and lotaustralin, nor other chemicals potentially influencing unpalatability in Heliconius; see (Jiggins, 2017)), and Heliconius “cyanogenic toxicity” or “toxicity” refer here to biosynthesized cyanogenic toxicity. Our study species has been considered generalist or oligophagous on Passiflora, largely acquiring its cyanogenic toxicity through de novo biosynthesis of cyanogens (Engler-Chaouat & Gilbert, 2007; Merrill et al., 2013), and the cyanogen profile of its main host plant P. biflora, also used in our experimental setup, consists of complex compounds which Heliconius larvae are unable to sequester (Engler-Chaouat & Gilbert, 2007; de castro et al., 2019).

All statistical analyses of the study were undertaken using R 3.4.4 (R Core Team, 2018). The values of cyanogen concentrations in the wild-collected samples were normalized with a square-root-transformation before applying ANOVA to test for population differences. Distributions of cyanogen traits were also explored by inspection of truncated weighted density distributions (accounting for the skewness of data towards near-zero values) using the “sm” R package (Bowman & Azzalini, 2014). For the analyses using common-garden data, the data were normally distributed and untransformed values were used.

Host plant toxicity and quality analyses

The cyanogenic content of the Passiflora biflora samples was analyzed using liquid chromatography-mass spectrometry (LC-MS/MS) with Agilent 1100 Series LC (Agilent Technologies, Germany) hyphenated to a Bruker HCT-Ultra ion trap mass spectrometer (Bruker Daltonics, Germany) following the procedure of de castro et al. (2019) (see File S4 for details on sample extraction, spectrometry and analyses). Mass spectral data were analyzed with the native data analysis software (Compass DataAnalyses, Bruker Daltonics). The cyanogen compounds of P. biflora were detected and quantified following de castro et al. (2019).

Plant quality traits assumed to be proxies of nutritional value, including nitrogen/carbon ratio (a proxy of protein/carbohydrate-ratio), were measured from the greenhouse-cultivated P. biflora treatment types (average number of individual leaves sampled per P. biflora treatment type: n = 63) with Dualex^® Scientific+ leaf clip meter (Cerovic et al., 2012; File S4). We tested for differences in cyanogen levels and plant quality traits between the treatment groups in an ANOVA including plant origin, treatment and their interaction.

Estimation of variance components, heritability, evolvability and maternal effects in Heliconius biosynthesized cyanogenic toxicity

The concentrations of the two tested biosynthesized cyanogens, linamarin and lotaustralin were strongly correlated with each other (r =0.897). Cyanogen concentration was not correlated with adult body mass in either sex (F_1,193 = 0.52, P = 0.470 for total cyanogen concentration; mass data available for n = 196 individuals). Adult mass also did not differ between the feeding treatments (F_4,186 = 0.88, P = 0.480; model including feeding treatment, sex and their interaction), and was not included as a variable in further analyses.

Variance components of Heliconius erato biosynthesized cyanogenic toxicity were analyzed based on the two-generation full-sib/half-sib breeding design (Lynch & Walsh, 1998a) by fitting linear mixed models applying the “animal model” framework (Kruuk, 2004; Wilson et al., 2010), in which the full pedigree (relatedness matrix) is incorporated as an individual-level random effect. The models were fitted using ASReml-R Version 4 (Butler et al., 2017), an R package that fits linear mixed models using Residual Maximum Likelihood (REML) estimation. The two-generation pedigree consisted of 20 broods and their parents, and included mother-offspring, father-offspring, full-sibling and paternal half-sibling relationships (n = 322 individuals, File S2). Of the 20 broods, six were included in the feeding treatments receiving four different P. biflora types as larval diet (File S2). Individuals of the parental generation were assumed to be unrelated. With the materials available, we maximized our power to detect genetic variation by generating as many family groups as possible. Because the number of families was limited, we also analyzed a large number of offspring per family, which allowed precise estimates of family means. Because our breeding design resulted mostly in full-sib families, we have limited power to distinguish additive-genetic from other genetic variance components (dominance, epistasis, maternal/paternal effects). We thus consider our genetic-variance estimate as a ‘broad-sense’ genetic variance, although the estimate approaches the narrow-sense estimate more than a simple full-sib design due to the parent–offspring and half-sib relations in the breeding design (Lynch & Walsh, 1998a).

Animal models were fitted separately for the three biosynthesized cyanogenic toxicity traits (concentration of linamarin, lotaustralin and total cyanogens). The primary models (referred to as “Animal Model1”) included the individual breeding value as a random effect, and the fixed effects of sex and feeding treatment (variance explained by fixed effects computed using the method of de Villemereuil et al. (2018); see also File S4). “Animal Model2” included additionally the random effect of maternal identity. We computed broad-sense heritabilities (H²) as the proportion of the total phenotypic variance (V_P) explained by the (broad-sense) genetic variance component (V_G). When calculating heritabilities we included the variance due to sex in the total phenotypic variance, but we did not include the variance due to the treatment because this component of the variance is unlikely to represent natural variation (de Villemereuil et al., 2018). Heritabilities are useful measures of the resemblance between offspring and parents, but are poor standardized measures of evolutionary potential (Hansen, Pélabon & Houle, 2011). As a measure of evolutionary potential, we therefore computed mean-scaled evolvabilities (e_μ) as e_μ = V_G/µ², where µ is the trait mean (of common garden brood data trait values), which gives the expected percent change in the trait mean per generation under a unit strength of selection (Hansen & Houle, 2008; Hansen, Pélabon & Houle, 2011).

Given that our design included mostly full-sib groups, we obtained an additional estimate of the component of toxicity variance explained by brood by estimating the intra-class correlation coefficient (ICC), i.e., repeatability, with a linear mixed model approach. For this, we implemented the “rptR” package (Stoffel, Nakagawa & Schielzeth, 2017), with brood identity as the grouping factor, a Gaussian distribution and 1000 parametric bootstraps to quantify uncertainty.

To further explore the role of maternal contributions to the inheritance of cyanogenic toxicity, we estimated maternal-effect coefficients following Falconer’s dilution model, which assumes that offspring phenotype is a direct, linear function of maternal phenotype, beyond the contribution of maternally transmitted genes (Lande & Price, 1989; Lynch & Walsh, 1998b; Walsh & Lynch, 2018b). Thus, the interpretation of this ‘maternal effect’ is related to specific measured traits (here, cyanogenic toxicity), and is different from the typical formulation of ‘maternal environmental effects’ (influence of the maternal environment on offspring traits). We estimated maternal effect coefficients (m) as the difference in slope between a mother-offspring (mo) and a father-offspring (fo) regression, i.e., m = β _mo − β _fo (Lande & Price, 1989; Lynch & Walsh, 1998b). We fitted the parent–offspring regressions as linear mixed models with the lme4 R package (Bates et al., 2014), including family as a random factor and maternal or paternal toxicity phenotype as explanatory variables of offspring toxicity phenotype.

Wide natural variation in toxicity within populations and between regions, but not along environmental gradients

The distribution of biosynthesized cyanogen levels, as estimated from the concentrations of biosynthesized compounds linamarin and lotaustralin, was similar across natural populations of H. erato sampled along an environmental rainfall gradient in Panama (Fig. 1C). A large proportion of individuals had low toxicity values (e.g., 40% of individuals had total cyanogen concentrations falling within the lower 25% of toxicity values), while a small proportion of individuals had a considerably higher toxicity level (e.g., 5% of individuals had total cyanogen concentrations within the highest 25% of toxicity values). Total cyanogen concentration and concentration of linamarin did not differ among the three Panamanian populations (Fig. 1B, Table 1, File S1). The concentration of lotaustralin differed among populations (F_2,89 = 7.97, P < 0.001), with the Intermediate population having higher concentrations than the two populations at the opposite extremes along the rainfall gradient (File S1).

The total cyanogen levels of the host plant Passiflora biflora differed among the natural populations along the rainfall gradient in Panama (F_2,48 = 3.65, P = 0.034; Fig. 1B, Table 1). The Dry population had substantially lower cyanogen levels compared with both the Intermediate and Wet populations (Fig. 1B, Table 1). In contrast, levels of cyanogen did not vary substantially across populations of H. erato (Fig. 1B, Table 1, File S1).

The populations along the altitudinal gradient in Ecuador did not differ in biosynthesized cyanogenic toxicity from each other, but differed markedly from their Panamanian conspecifics (Fig. 1B, Table 1). Most Ecuadorian samples lacked both linamarin and lotaustralin (File S1). Linamarin was detected in only 29% of samples, most of which originated from the high-altitude site (75% of those with any signal), and the average biosynthesized cyanogen level in these samples was very low (0.063%; Fig. 1B, Table 1, File S1), which is below the minimum level observed in the Panamanian populations.

Biosynthesized toxicity does not vary with diet but varies among broods

The level of biosynthesized cyanogenic toxicity varied two-fold among the 20 common-garden reared broods (range of brood averages =0.48%–0.88% of dry mass, Fig. 2A). Phenotypic variance components of cyanogenic toxicity estimated within an “animal model” framework (Kruuk, 2004; Wilson et al., 2010), which included the full pedigree as an individual-level random effect and fixed effects of sex and feeding treatment (Model1), revealed a statistically significant genetic variance component (V_G) for total biosynthesized cyanogen concentration level (Fig. 1D, Table 2, File S1), corresponding to a broad-sense heritability (H²) of 0.115 (H² (linamarin) =0.139, H² (lotaustralin) =0.067; Table 2, File S1). We confirmed the presence of detectable among-brood variance in a simpler brood-level analysis, in which brood identity explained 12.5% of the variation in biosynthesized cyanogen concentration across common-garden reared broods (repeatability R = 0.125, SE =0.059, CI = [0.015, 0.245], P = 6.02e−05). Among-brood variance could originate from several sources including (1) additive genetic variance, (2) maternal environmental or genetic effects, and (3) shared environment, and below we explore these options further.

First, we consider the estimated broad-sense genetic variance component and estimate evolutionary potential of toxicity, assuming the genetic component consists of largely additive genetic variance. The broad-sense evolvability of total cyanogen concentration was e_μ = 1.55%, SE ± 0.89 (e_μ(linamarin) =1.89%, e_μ(lotaustralin) =1.04%; Table 2, File S1), meaning that the trait mean is expected to change by 1.55% per generation per unit strength of selection (Hansen & Houle, 2008; Hansen, Pélabon & Houle, 2011).

Our baseline estimates of repeatability, as well as broad-sense heritability and evolvability (Model 1) thus showed that there is moderate but detectable variation among broods in biosynthesized toxicity. We then estimated maternal and paternal contributions to this variation. When including maternal identity as a second random effect in the animal model (Model2), this maternal variance component explained most of the variance explained by the genetic variance component in Model1, whereas V_G became negligible, illustrating that the two random terms explain essentially the same variance component (i.e., the among-brood variance; Table 2, File S1). Parent-offspring regressions revealed that offspring cyanogen concentration was positively associated with the cyanogen concentration of the mother (slope =0.181 ± 0.071, t =2.539, P = 0.012 for total cyanogen concentration; Fig. 2C), while the association between fathers and offspring showed a negative, but statistically non-significant trend (slope = −0.122 ± 0.083, t =  − 1.469, P = 0.143 for total cyanogen concentration). In turn, the midparent–offspring regression was nearly flat (slope =0.046 ± 0.107, t =0.435, P = 0.664 for total cyanogen concentration; Fig. 2C). These regression slopes translate into a maternal-effect coefficient for total cyanogen concentration of m =0.261, CI = [−0.099, 0.614] (m (linamarin) =0.288, CI = [−0.094, 0.687]; m (lotaustralin) =0.193, CI = [−0.135, 0.520]), which implies the existence of positive maternal effects in biosynthesized cyanogenic toxicity. However, uncertainty in these maternal effects remains notable, as illustrated by the wide confidence intervals.

Finally, we investigated other sources, particularly environmental origins of variance. Larvae reared on different host plant types (Fig. 1A) did not differ in biosynthesized cyanogenic toxicity as adults (Fig. 2B), and the host plant type explained a negligible proportion of variance in butterfly toxicity level (V_trm/V_P = 0.011; Fig. 1D, Table 2, File S1). This is despite the strong effect of the watering treatments on the cyanogenic toxicity and quality indicators of the host plants (e.g., nitrogen-carbon balance; see File S3), and also indicates that among-brood variance is unlikely to be explained by common environment effects related to dietary quality. Butterfly sex was also a negligible source of variation in toxicity (V_sex/V_P = 0.011; Fig. 1D, Table 2, File S1). Despite the common-garden conditions (outside of dietary treatments) that all individuals experienced throughout their life cycle, the residual (unexplained) variance component of cyanogen concentration remained large (V_res/V_P = 0.874 for total cyanogens, Fig. 1D, Table 2, File S1). This large amount of unexplained variance is also demonstrated by a difference in toxicity level between the parental and offspring generation (F_1,320 = 27.42, P < 0.001, effect size (Hedges’ g) =0.962 for total cyanogens), despite the two generations sharing common garden conditions.

Implications of evolutionary potential in chemical defenses for aposematic and mimetic systems

Studies of mimicry often focus on color pattern, i.e., the warning signal (Benson, 1972; Mallet & Barton, 1989; Kapan, 2001; Merrill et al., 2015; Jiggins, 2017). However, it is important to note that contrasting selective agents could be acting on the defense trait associated with the warning signal, both traits equally required for aposematism and mimicry. The selection pressures acting on chemical defense traits can be expected to be complex and dynamic (Speed et al., 2012; Jones, Speed & Mappes, 2016; Briolat et al., 2019). The role of predators as the major selective agent is clearly a key component of possible adaptive evolutionary processes herein (Ruxton et al., 2018; Briolat et al., 2019). Even high evolutionary potential will only lead to evolutionary changes if the defense trait variation influences predation rates (Speed et al., 2012). In general, empirical data on the impact of variation in defenses and its signaling on predators is lacking, but there is evidence for variation among predators in traits including perceptual abilities, individual motivation and tolerance to chemical defenses (Halpin et al., 2017; Briolat et al., 2019; Hämäläinen et al., 2020). Predators may also respond not only to average defense levels, but to the variation in them (Skelhorn & Rowe, 2005; Barnett, Bateson & Rowe, 2014). Furthermore, prey toxicity level does not necessarily correlate with actual unpalatability, the ultimate cue for predator choices (Winters et al., 2018). Predator-mediated selection is therefore expected to be dynamic and difficult to estimate. In addition, bottom-up selection related to co-evolutionary relationships with host plants, as well as selection on energy optimization may also play a role (Turner, 1984; Ruxton et al., 2018).

Species may only respond adaptively to such selection pressures given heritable genetic variation in these traits. Evolutionary potential in chemical defenses could allow species to respond to the expected dynamic top-down and bottom-up selection, resulting in diverse patterns of adaptive evolution in chemical defense traits. This could have important consequences for the evolutionary trajectories of aposematic and mimetic species, and ultimately play a role in the evolution of warning signal diversity. This is because evolutionary changes in chemical defense levels could alter mimetic relationships among groups of species across geographic space, including shifts between Müllerian, Batesian and quasi-Batesian relationships, which could in turn influence advergent vs. convergent selection and contribute to the observed diversity in mimicry systems (Mallet & Joron, 1999; Speed, 1999; Briolat et al., 2019). With regards to Heliconius, the combination of a locally dynamic selective environment, flexible genetic mechanisms for generating warning signals (Wallbank et al., 2016; Van Belleghem et al., 2017), and evolutionary potential in toxicity could explain the rapid origin of inter- and intraspecific variation. In this light, the ability of Heliconius to biosynthesize defense compounds could allow for increased genetic control of chemical defenses compared to most other chemically defended Lepidoptera, which rely solely on host plant compound sequestration (Bowers, 1992; Nishida, 2002; Zvereva & Kozlov, 2016). Because Heliconius can sequester cyanogenic compounds of only some Passiflora species (Engler, Spencer & Gilbert, 2000; de castro et al., 2019), selection would be targeted toward de novo synthesized defenses particularly in environments lacking such host species. Furthermore, previous investigations both confirmed that the cyanogenic glucoside compounds in Heliconius deter vertebrate predators (Cardoso, 2020) and that predator avoidance learning rate increases when comparing low to moderate levels of cyanogenic compound concentrations (Chouteau et al., 2019). Thus, the level of cyanogen defenses is likely to be a direct target of selection by predators in some circumstances.

By altering the nature of mimetic relationships or selection for increased defense levels on a new warning signal form, evolvability of defense traits could drive the evolution and spread of new color patterns and the reproductive isolation that may follow, potentially even leading to the formation of new species (Mallet, McMillan & Jiggins, 1998; Jiggins et al., 2001). A phylogenetic signal in the cyanogenic glucoside chemical profile in the Heliconiini tribe supports this idea (de castro et al., 2019; Sculfort et al., 2020) and underscores the potential for adaptive evolution in chemical defense traits to influence evolutionary trajectories in mimicry systems. However, more empirical data is called for in future studies to estimate levels of evolutionary potential and patterns of selective responses in defensive traits in aposematic and mimetic species.

Evolutionary processes shaping chemical defense variation: the role of automimicry

A further complication in chemical defense evolution is the interplay between individual and public good (Jones, Speed & Mappes, 2016). In this respect, our data demonstrating a broad-sense genetic component in the defense trait coupled with patterns of natural variation of cyanogen levels in Heliconius erato populations are consistent with balancing selection related to automimicry, in which automimic “cheaters” exploit the protection given by their warning coloration with investing less in chemical defenses themselves. Automimicry may occur if producing and maintaining chemical defenses is costly (Brower, Pough & Meck, 1970; Speed et al., 2012), leading to selection reducing individual-level toxicity (Bowers, 1992; Zvereva & Kozlov, 2016; Ruxton et al., 2018), and if there are also individual benefits of investing in chemical defenses (e.g., Ruxton & Speed, 2006; Svennungsen & Holen, 2007). Previous studies indicate that de novo synthesis of defenses may generally be more costly than the sequestration of plant compounds (Zvereva & Kozlov, 2016), which may make automimicry more likely in species with biosynthesized defenses. The occurrence of automimicry could lead to an increase in population average palatability, which could theoretically dilute the protection conveyed by the warning signal (Gibson, 1984; Speed et al., 2012; Jones, Speed & Mappes, 2016), as in Batesian and quasi-Batesian mimicry, and influence the evolutionary trajectories of the species sharing the warning signal. Thus, empirical data on within-population defense variation and automimicry can give important insights into the diversity of mimetic systems (Jones, Speed & Mappes, 2016).

Our results show that all three butterfly populations along the studied rainfall gradient exhibited similar skewed distributions of biosynthesized cyanogens. The majority of individuals in these populations had relatively low, although not zero, levels of toxicity (24% of individuals had less than 0.5% concentration of biosynthesized cyanogenic compounds from an overall range of 0.1–2.9%), whereas a smaller proportion had substantially higher levels, and the pattern was surprisingly similar in all populations. Such a pattern could be explained by age distributions (Heliconius have maximum adult lifespans in excess of six months and continuous reproduction throughout life; Gilbert, 1972; Dunlap-Pianka, Boggs & Gilbert, 1977), as previous studies suggest that Heliconius accumulate cyanogenic toxins with age (Nahrstedt & Davis, 1985; de castro et al., 2020), or by individual differences in the balance between sequestration and biosynthesis of cyanogenic compounds (Engler-Chaouat & Gilbert, 2007; de castro et al., 2019; de castro et al., 2020; de castro et al., 2021). More generally, chemical defense polymorphisms could also be maintained by patterns of predator selection targeting defense variation itself (Skelhorn & Rowe, 2005; Barnett, Bateson & Rowe, 2014). However, as many individuals in the study populations had near zero levels of biosynthesized toxins, the pattern could also be indicative of automimicry. In Heliconius, the costs of chemical defense remain poorly known, but some cases of automimicry have previously been reported (Arias et al., 2016). Here, toxicity and mass were not associated, which could imply that potential costs of cyanogen biosynthesis are not related to or mediated by a trade-off with body size. In terms of benefits of chemical protection, requirements for automimicry are likely to be met, as empirical data on Heliconius imply direct benefits (surviving predator attacks) for chemically defended individuals (Boyden, 1976; Pinheiro, 1996; Pinheiro & Campos, 2019). Furthermore, an increase in protection level is found especially when comparing low to moderate cyanogenic toxicity levels (Chouteau et al., 2019), applicable to the observed cyanogen levels of H. erato (Arias et al., 2016; de castro et al., 2019). However, future studies on natural defense variation and automimicry in Heliconius should measure the entire range of cyanogens (both sequestered and biosynthesized) and other chemical defense compounds potentially interacting to produce the unpalatability level perceived by relevant predators (Jiggins, 2017; Winters et al., 2018). In more general terms, future research in this framework should investigate the costs of defensive compound biosynthesis (see also Zvereva & Kozlov, 2016) and its role in the evolution of automimicry, as well as automimicry-originating defense variation as a potential source introducing diversity into mimetic systems.

The roles of phenotypic plasticity and ecological sources of defense variation in wild populations

It has been suggested that the often wide variation in defense toxins could be explained by the non-adaptiveness of toxicity profiles, in which case the variation would stem from chance effects and plastic responses to the environment (Speed et al., 2012). Empirical data on chemical defense variation in wild populations e.g., along environmental gradients could help to evaluate ecological vs. evolutionary origins of chemical defense variation, and the expected consequences for diversity in aposematic and mimetic systems (Briolat et al., 2019). Although the large amount of unexplained variance in our quantitative genetics model may suggest the presence of unknown environmental sources of variation in biosynthesized toxicity, we found very little evidence for a role of known ecological factors in explaining variation in toxicity profiles along our studied environmental gradients, despite the steepness of these gradients. For example, the chemical profiles were very similar in the Panamanian Wet and Dry populations, which experience a five-fold difference in rainfall. Similarity among populations is unlikely to be caused by large-scale dispersal along the gradient, because Heliconius erato has a restricted home range and a dispersal distance of around 1–2 km per generation (Turner, 1971). For the same reason, the collection locations of adults are also likely to reflect the locations in which larvae fed. Furthermore, despite water availability inducing substantial differences in the cyanogen content and quality indicators of P. biflora, both our field and experimental data indicate that variation within the host plant species does not contribute to variation in H. erato biosynthesized toxicity. This is in contrast to our expectations of dietary quality as a major source of plastic chemical defense variation.

However, there were striking differences in the chemical profiles between the two distinct warning color forms of H. erato that we sampled from geographically distant populations in Panama and Ecuador. Compared to their Panamanian conspecifics sampled along the gradient of rainfall at sea level, the Ecuadorian populations inhabit a very different and geographically distant biotope along the eastern slope of the Andes (Montejo-Kovacevich et al., 2020). These environments differ in larval host plant and nectar/pollen availability, the community of predators and of co-mimics, which are far more diverse in Ecuador. Lower levels of cyanogen biosynthesis in Ecuadorian populations of H. erato may suggest that these populations have specialized on a Passiflora host which provides more cyanogens for sequestration (Engler-Chaouat & Gilbert, 2007; de castro et al., 2019; de castro et al., 2020; de castro et al., 2021). Alternatively, the more diverse mimetic ring in Ecuador could also provide greater protection by more toxic model species, relaxing selection on defenses in co-mimics (Turner, 1987; Speed, 1993; Speed, 1999; Mallet, 1999). Differences in the collection methods between populations in Panama and Ecuador are unlikely to explain the extremely low levels of biosynthesis in the Ecuadorian individuals, as biosynthesis levels in common-garden reared individuals are generally shown to be similar with the levels of corresponding wild individuals in the current and previous studies on Heliconius (Engler-Chaouat & Gilbert, 2007; de castro et al., 2020).

In this light, our data may provide an example of how ecological variation could lead to differential selective environments between conspecific populations, thereby leading to local adaptation of the acquisition mechanisms and levels of defensive compounds. Such adaptive radiations would be facilitated by evolutionary potential in defense traits. More data on different aposematic and mimetic systems will be needed to confirm whether ecological variation may have a greater influence on chemical defense variation through introducing spatially and temporally varying selection pressures, rather than by introducing “ecological noise” through plastic responses in defense traits.