Expanding the Isotopic Toolbox to Track Monarch Butterfly (Danaus plexippus) Origins and Migration: On the Utility of Stable Oxygen Isotope (δ¹⁸O) Measurements

Introduction

The monarch butterfly (Danaus plexippus) is an iconic migratory insect that navigates over thousands of kilometers between natal sites in the eastern USA and Canada and overwintering sites in the Oyamel fir (Abies religiosa) forests of central Mexico (Urquhart and Urquhart, 1976; Ries et al., 2015; but see Vander Zanden et al., 2018). This outstanding migration involves multiple generations whereby those individuals returning to well-established, long-term overwintering sites do so without any previous experience of their locations. Despite being a high-profile conservation issue among Mexico, Canada and the USA, monarch butterflies have declined considerably over the last two decades (Semmens et al., 2016). Potential causes for the decline are many but likely factors include loss of milkweed (Asclepia spp.) on the breeding grounds, declining extent and quality of overwintering sites, climate change and the myriad of challenges faced during the migratory passage (Flockhart et al., 2015; Thogmartin et al., 2017). Key to understanding declines of this and other migratory species is the establishment of migratory connections that link breeding, stopover and wintering sites (Hobson and Wassenaar, 2019). However, this is challenging for small insects. Previous attempts have used an impressive mark recapture effort involving tagging Monarchs on breeding and stopover sites with labels affixed to wings (Monarch Watch; https://www.monarchwatch.org/tagmig/index.htm), an approach requiring recovery en route and at the few publicly accessible wintering sites in Mexico. However, during the mid-1990s Hobson et al. (1999) and Wassenaar and Hobson (1998) developed an approach using ratios of naturally occurring stable isotopes in Monarch wings to infer natal origins. Those studies marked a turning point in the conservation of migratory monarchs since it readily identified the US Midwest as the center of production of Monarchs recruited into the overwintering population in Mexico. Since then, the approach has been used by other researchers interested in Monarch patterns of spring recolonization in North America (Miller et al., 2012; Flockhart et al., 2013), the effects of natal origin on parasite loads (Altizer et al., 2015). Other applications have involved the role of wing coloration in flight distance (Hanley et al., 2013) and general conservation concerns related to where most individuals are being produced (Flockhart et al., 2017).

The stable isotope approach to tracking migrant animals is based on the fact that naturally occurring stable isotopes of several elements in nature can provide information on provenance and/or habitat. Spatial patterns of stable isotopes (i.e., “isoscapes”) are ultimately transferred to animal tissues through local foodwebs and so spatial information related to the period of tissue growth can be inferred providing we know the nature of such isoscapes and how stable isotope values change or discriminate from the source to the tissue of interest (Hobson and Wassenaar, 2019). The early studies on stable isotope patterns in eastern monarchs indicated that both stable–carbon (¹³C/¹²C expressed as δ¹³C) and stable-hydrogen (²H/¹H, δ²H) isotope ratio measurements of wings could provide information on monarch natal origins in Canada and the USA. However, δ²H measurements provided a more powerful means of inferring origins than δ¹³C measurements. Hydrogen in foodwebs is ultimately derived from precipitation. The amount-weighted average precipitation δ²H (expressed as δ²H_p) varies across continents according to well-established principles (Clark and Fritz, 1997) and these isoscape patterns are transferred to plants and higher trophic-level organisms (Bowen and West, 2019). Metabolically inert tissues such as animal keratins and chitins are related to local growing season averaged δ²H_p and can be later measured to infer natal origins once the relationship between δ²H_p and tissue δ²H is established (Hobson, 2019). This principle has been used effectively to infer origins of numerous taxa including insects, birds and mammals on several continents (Hobson et al., 2018a; Hobson and Wassenaar, 2019).

Water contains both hydrogen and oxygen but measurement of the stable isotope ratios in oxygen (¹⁸O/¹⁶O; δ¹⁸O) have not been used extensively to study animal movement (but see Bryant and Froehlich, 1995; Kohn, 1996; Hobson et al., 2004, 2009a, 2012; Bowen et al., 2005; Ehleringer et al., 2008; Chesson et al., 2013; Hobson and Koehler, 2015; Pekarsky et al., 2015). This derives from the fact that oxygen has a compressed isotopic scale compared to hydrogen and that δ²H and δ¹⁸O values in environmental waters are highly correlated via the global meteoric water relationship (Craig, 1961). Also, previously, measurement of δ¹⁸O in organic materials has been challenging analytically (Hobson and Wassenaar, 2019). Oxygen in animal tissues can be derived from air, drinking water, and diet whereas hydrogen is derived from only diet and drinking water. However there is potential for δ¹⁸O when combined with δ²H measurements to provide additional information on provenance of animals because environmental waters can record the degree of evapotranspiration through the evaluation of deuterium excess (defined as δ²H-8^*δ¹⁸O; Dansgaard, 1964; Rozanski et al., 1993). In highly evaporative systems, we expect departures from the global meteoric water line whereby deuterium excess values are higher than the global average of 10‰, a phenomenon highly associated with low relative humidity. So, again, the relationship between δ²H and δ¹⁸O values in insect wing chitin (δ²H_w and δ¹⁸O_w) has the potential to contain important environmental and climate information.

The evaluation of both δ²H_w and δ¹⁸O_w values of monarch butterflies was first reported by Fourel et al. (1998) as a proof of principle for the measurement of these isotopes via continuous-flow isotope ratio mass spectrometry. That study showed a positive correlation between the two isotopes. However, since then, no studies have further investigated the relationship between these isotopes in monarch wings in particular, and only a few researchers have reported on δ²H and δ¹⁸O values among higher-order taxa. Hobson and Koehler (2015) reported that while feathers of American Redstart (Setophaga ruticilla) showed a good relationship between δ²H and mean growing-season precipitation δ²H (r² = 0.77), the same was not true for feather δ¹⁸O (r² = 0.32). However, for dragonflies, insects that form wings from nutrients derived from an aquatic larval stage, Hobson et al. (2012) found an excellent relationship between δ²H_w and δ¹⁸O_w values (r² = 0.92). The question remains, then, whether δ²H_w and δ¹⁸O_w values of terrestrial insects are also well-correlated and indeed reflect underlying relationships in environmental waters that could be used to increase the accuracy of assignment to origin compared to the use of δ²H measurements alone. Here, we report measurements of δ²H_w and δ¹⁸O_w values of monarch butterflies sampled at known origin natal sites in eastern North America. Our objective was to determine if a dual isotope approach could provide additional environmental information related to regional patterns in evaporative conditions that could ultimately contribute to more accurate delineation of natal origins of monarchs and other terrestrial insects.

Materials and Methods

Monarch Sample

Monarch wing material was obtained from samples reported in Wassenaar and Hobson (1998) and Hobson et al. (1999). Those were monarchs raised on milkweed exposed only to precipitation at 31 known sites throughout the range of the eastern monarch population and which formed the calibration relationship between mean growing season δ²H_p and δ²H_w reported by Hobson et al. (1999) (Table 1). Additionally, monarchs collected from roadkill mortality (n = 92) at a single site in northeast Mexico during the autumn migration of 2015–16 were used to investigate assignment to natal origins using both δ²H_w and δ¹⁸O_w values.

Table 1. Mean isotope values (± 1 SD ‰) of wild reared monarch butterflies from eastern North America and mean growing-season δ¹⁸O_p (GSO) and δ²H_p (GSD) at each location plus the 95% CI shown in parentheses from Bowen et al. (2005).

Stable Isotope Analyses

All monarch wing samples were soaked and rinsed in a 2:1 chloroform:methanol solution and air dried. Subsamples were cut from the same region of the hindwing to reduce inter-sample variance due to isotopic effects from pigmentation (Hobson et al., 2017) and weighed (0.35 ± 0.02 mg) into silver capsules. All samples were prepared for δ²H and δ¹⁸O analysis at the Stable Isotope Laboratory of Environment Canada, Saskatoon, Canada. Our approach involved the analysis of both δ²H and δ¹⁸O on the same analytical run (i.e., both H₂ and CO gases were analyzed from the same pyrolysis) from samples and standards. All measurements were performed on a HTC system (Thermo Finnigan, Bremen, Germany) equipped with a Costech Zero-Blank autosampler. The helium carrier gas rate was set to 120 ml/min through a 0.6 m ¼-inch 5-Å molecular sieve (80–100 mesh) GC column. The HTC reactor was operated at a temperature of 1,400°C and the GC column temperature was set to 90°C. After separation, the gases were introduced into a Delta V plus isotope-ratio mass spectrometer via a ConFlo IV interface (Thermo Finnigan, Bremen, Germany). The eluted N₂ was flushed to waste by withdrawing the CF capillary from the ConFlo interface. We used keratin reference standards CBS and KHS to calibrate sample δ²H (−197 and −54.1‰, respectively) and δ¹⁸O (+2.50 and +21.46‰, respectively; Qi and Coplen, 2011). These standard values were used vs. the newly reported values of Soto et al. (2017) simply to maintain consistency with our earlier published work for monarch assignments and we note that this will not affect the assignment per se (providing appropriate calibration algorithms are used). Based on replicate (n = 5) within-run measurements of keratin standards, sample measurement error was estimated at ±2‰ for δ²H and ±0.4‰ for δ¹⁸O. All H results are reported for non-exchangeable H and for both H and O in typical delta notation, in units of per mil (‰), and normalized on the Vienna Standard Mean Ocean Water – Standard Light Antarctic Precipitation (VSMOW-SLAP) standard scale.

Isotopic Assignments

We created δ¹⁸O_w and δ²H_w isoscapes based on derived transfer functions using wing stable isotope values of known origin monarchs found in this study. This was accomplished by calibrating amount-weighted growing-season average precipitation δ¹⁸O and δ²H (δ¹⁸O_p and δ²H_p) (Terzer et al., 2013; Hobson et al., 2014) surfaces into separate δ²H_w and δ¹⁸O_w isoscapes. We then depicted potential origins of a sample of roadkilled monarchs salvaged during their autumn migration through northeastern Mexico. That depiction was done using techniques described previously (Hobson et al., 2018b). Briefly, we used a likelihood-based assignment method (Hobson et al., 2009b; Wunder, 2010; Van Wilgenburg et al., 2012) to assign individual monarchs using δ²H_w and δ¹⁸O_w to each calibrated wing isoscape (see results), separately. We used the residual SD error of 9.33‰ for δ²H_w and 1.50‰ for δ¹⁸O_w from regressions in our assignments. We limited assignments to the current known geographic range for the eastern monarch population and used this as a spatial mask (i.e., clip) to limit our analysis. We estimated the likelihood that a cell (pixel) within the δ²H_w and δ¹⁸O_w isoscapes represented a potential origin for a sample by using a normal probability density function (pdf) based on the observed δ²H_w and δ¹⁸O_w, and thus depicted the likely origins of each monarch by assigning individuals to the δ²H_w and δ¹⁸O_w isoscapes, separately. We arbitrarily chose a 2:1 odds ratio to include only those pixels (coded 1) with at least a 67% probability of origin vs. all others (coded 0). This resulted in a binary map per assigned individual of presence and absence. We then summed the results of individual assignments by stacking the surfaces for a final depiction. We conducted geographic assignments to origin using functions within the R statistical computing environment (R Core Team., 2016) employing the “raster” (Hijmans, 2016) and “maptools” (Bivand and Lewin-Koh, 2016) packages. Thus, the final assignment surface depicted the number of individuals co-assigned at each pixel based on the odds criteria. We also conducted assignment to origin analysis for individual Monarch samples using a dual-isotope (δ²H_w, δ¹⁸O_w) approach applied with a multivariate normal probability density function (mvnpdf; see Hobson et al., 2014 for details). Similar to the univariate pdf approach, the mvnpdf method calculates the probability that a particular spatially referenced cell represents a potential origin in calibrated raster isoscape space.

Results

Monarch Calibration

We found a positive relationship between monarch δ²H_w and δ¹⁸O_w values for known-origin, outdoor-raised individuals (δ²H_w = 5.21 ^* δ¹⁸O_w – 211, r² = 0.42; Figure 1). The relationship between δ¹⁸O_w and δ¹⁸O_p was positive (δ¹⁸O_w = 0.54 ^* δ¹⁸O_p + 22.5; Figure 2) but relatively weak (r² = 0.30, p < 0.001). In contrast, the relationship between δ²H_w and δ²H_p was positive (δ²H_w = 0.78 ^* δ²H_p – 77.42; Figure 3) and stronger (r² = 0.61, p < 0.001).

Figure 1. Relationship between δ²H_w and δ¹⁸O_w for monarchs raised outdoors at 31 known sites (Table 1) throughout their range in 1996. Sample originally reported for δ²H_w and δ¹³C_w by Hobson et al. (1999) but these are new δ²H_w values measured simultaneously with δ¹⁸O_w.

Figure 2. Relationship between monarch δ¹⁸O_w and amount-weighted mean growing-season average δ¹⁸O_p (Bowen et al., 2005) for monarchs raised outdoors at 31 known sites (Table 1) throughout their range in 1996.

Figure 3. Relationship between monarch δ²H_w and amount-weighted mean growing-season average δ²H_p (Bowen et al., 2005) for monarchs raised outdoors at 31 known sites (Table 1) throughout their range in 1996.

Assignment Evaluation

Using the relationships established between δ¹⁸O_w and δ²H_w values and δ¹⁸O_p and δ²H_p values, respectively, we created predicted monarch wing isoscapes for each isotope (Figures 4A,B). These surfaces are superficially similar due to the correlation between isotopes. We then used these surfaces to assign origins of a sample of monarchs from roadkill mortality during their fall migration. That sample included monarchs from across their range and we contrasted depictions of origins of these individuals using only δ²H_w values, only δ¹⁸O_w values and both δ¹⁸O_w and δ²H_w values in a multivariate normal assignment. These depictions were similar in that they showed highest probability in the southwestern portion of the range but differences were also apparent (Figures 5A–C). The dual isotope approach constrained the assignment compared to the use of δ²H_w values alone.

Figure 4. Depiction of (A) the expected monarch δ¹⁸O_w isoscape based on the relationship derived between δ¹⁸O_w and amount-weighted mean growing season average δ¹⁸O_p (Bowen et al., 2005) shown in Figure 2 and (B) the expected monarch δ²H_w isoscape based on the relationship derived between δ²H_w and amount-weighted mean growing-season average δ²H_p (Bowen et al., 2005) shown in Figure 3. Legend is the number of individuals assigned to a pixel based on the odds ratio criterion used.

Figure 5. Depiction of origins of a subsample of roadkilled monarchs during fall migration near Monterrey, Mexico, using (A) only δ²H_w values, (B) only δ¹⁸O_w values, and (C) both δ²H_w and δ¹⁸O_w values in a multivariate normal likelihood approach. Legend is the number of individuals assigned to a pixel based on the odds ratio criterion used.

Deuterium Excess

We contrasted deuterium excess values calculated for individual monarch wings and compared these to the growing-season average deuterium excess calculated for each natal site. No correlation was found between these values (Figure 6, r² = 0.02). This suggests that little information associated with precipitation environmental deuterium excess was transferred to monarchs per se.

Figure 6. Relationship between measured deuterium excess in known-origin monarch wings (Table 1) and predicted growing season precipitation deuterium excess for these locations.

Discussion

Our results confirm a strong relationship between wing chitin δ²H values and precipitation δ²H values at known natal sites. This result underlines the strong assignment power of using Monarch wing δ²H values to infer origins and confirms that even though the earlier investigations of Monarch δ²H values to infer origins used different analytical techniques (i.e., offline zinc reduction and steam equilibration, Wassenaar and Hobson, 1998; Hobson et al., 1999), the power of that isotopic assignment holds. The different analytical approaches resulted in different calibrations between wing chitin δ²H (δ²H_w) and predicted amount-weighted growing season δ²H_p (earlier method: δ²H_w = 0.62 ^* δ²H_p – 79; r² = 0.69; recent method: δ²H_w = 0.78 ^* δ²H_p – 77.4; r² = 0.61) and we recommend that from now on authors use the recent calibration relationship reported here. The δ²H_w results contrasts with the much weaker relationship we found between δ¹⁸O_w and δ¹⁸O_p. Thus, while the measurement of δ¹⁸O_w possibly confers additional information on origin, the strong relationship between δ²H_w and δ¹⁸O_w values likely precludes the effective use of a dual isotope approach in most cases. This result has been seen previously in assignment exercises involving several taxa (reviewed by Hobson and Koehler, 2015). Our comparison of deuterium excess values calculated for monarch wings and those of the mean growing season average deuterium excess for natal sites confirms that little environmental information related to evaporative conditions was available from our sample by running δ¹⁸O analyses. Nonetheless, we recognize the potential utility of using both δ¹⁸O_w and δ²H_w values, especially for areas corresponding to potential origins involving high deuterium excess values that were not included in our study. It may be possible then, to identify migrant individuals raised in xeric vs. humid environments and we encourage further work in this area for migrant insects in general. Clearly, further dual isotope analyses of a large sample of monarchs from known natal origins that spans large isotopic gradients in North America will be useful in testing situations where this more extensive analytical approach might be justified.

It is interesting to speculate why tissue δ¹⁸O values do not seem to add much additional information on origins of Monarch Butterflies and other terrestrial taxa based on δ²H analyses alone. The typical response to this question is that oxygen enters metabolic pathways from more diverse origins than hydrogen (Hobson and Koehler, 2015). Fundamentally, oxygen is derived metabolically from diet, drinking water and air whereas hydrogen is derived only from diet and drinking water. Atmospheric oxygen is added to the body water pool during metabolism and respiration and has a relatively constant isotopic composition (Luz and Barkan, 2011). Therefore, the oxygen isotopic coupling between tissues and environmental waters is dependent on the relative contribution of respirative processes to the oxygen isotopic composition of body water. Additionally, with insects, atmospheric oxygen for respiration is obtained primarily by diffusion through the spiricules and tracheal pathways (Chown et al., 2006) which may further alter its isotopic composition. It is worth noting that while the hydrogen isotopic compositions of chitin are broadly similar to those of keratins and other tissues of birds or mammals (Ehleringer et al., 2008; Hobson and Koehler, 2015), the δ¹⁸O values of chitin are generally more positive. This, along with the relatively poor fit between oxygen isotopic compositions of wing chitin and precipitation, likely indicates that slightly different sources, processes or pathways are utilized for oxygen in terrestrial insects than in other taxa. This would not necessarily preclude the use of both isotopes in typical assignments to origin but such assignment models require isotopes to be relatively orthogonal. For oxygen isotopes, current analytical methods result in a much higher relative measurement error than for hydrogen isotopes (0.4 vs. 2‰ for hydrogen). To use both isotopes to determine an inferred deuterium excess of environmental waters requires the Meteoric Water Line multiplicative factor of 8, so that currently we can only determine the deuterium excess values to precisions of about 3.2‰. Precipitation varies in deuterium excess from ~0 to +10‰ which means that such poor precisions will require large numbers of samples depending on which populations are being examined. This is not to say that measurement of both isotopes in monarch tissue will not be useful, but that, in general, the use of wing chitin δ²H is the most preferred for single isotope assignments.

In addition to the earlier proof of concept papers on the use of both wing δ¹⁸O and δ²H values to investigate natal origins of monarch butterflies (Fourel et al., 1998), our investigation has underlined the fact that there is much future work required to move the field forward. For example, we stress again that monarchs deriving from natal sites in areas of high deuterium excess (e.g., southern states) might still produce wings which reflect these relationships and so to improve our statistical power more sampling in those areas would be desirable.

Data Availability

All datasets generated for this study are included in the manuscript and/or the Supplementary Files.

Ethics Statement

Only dead monarch butterflies were examined in this study. These were derived from collections, natural overwinter mortality and roadkills, so no animal care or ethics approval was necessary.

Author Contributions

KH and GK conceived of the idea and performed stable isotope measurements. KK performed isotope assignments and assisted with figures. All authors contributed to the writing of the manuscript.

Funding

This study was funded by Environment and Climate Change Canada operating funds to KH and GK and by an NSERC Discovery Grant to KH.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Blanca X. Mora Alvarez for assistance with preparing samples for stable isotope analyses. Two reviewers provided useful comments on an earlier draft of the paper.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2019.00224/full#supplementary-material

References