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

The measurement of naturally occurring stable hydrogen (δ2H) and carbon (δ13C) isotopes in wings of the eastern North American monarch butterflies (Danaus plexippus) have proven useful to infer natal origins of individuals overwintering in Mexico. This approach has provided a breakthrough for monarch conservation because it is the only viable means of inferring origins at continental scales. Recently, routine simultaneous analyses of tissue δ2H and δ18O of organic materials has emerged leading to questions of whether the dual measurement of these isotopes could be used to more accurately infer spatial origins even though the two isotopes are expected to be coupled due to the meteoric relationship. Such refinement would potentially increase the accuracy of isotopic assignment of wintering monarchs to natal origin. We measured a sample of 150 known natal-origin monarchs from throughout their eastern range simultaneously for both δ2H and δ18O wing values. Wing δ2H and δ18O values were correlated (r2=0.42). We found that wing δ2H values were more closely correlated with amount-weighted growing season average precipitation δ2H values predicted for natal sites (r2=0.61) compared to the relationship between wing δ18O values and amount-weighted growing season average precipitation δ18O values (r2=0.30). This suggests that monarch wing δ2H values will be generally more useful in natal assignments than δ18O values. Spatial information related to the use of deuterium excess in environmental waters was similarly found to be not useful when applied to monarch wings likely due to the considerable variance in wing δ18O values. Nonetheless, we recommend further testing of monarch wing δ2H and δ18O values from known natal sites with an emphasis on field data across a strong gradient in precipitation deuterium excess.

The measurement of naturally occurring stable hydrogen (δ 2 H) and carbon (δ 13 C) isotopes in wings of the eastern North American monarch butterflies (Danaus plexippus) have proven useful to infer natal origins of individuals overwintering in Mexico. This approach has provided a breakthrough for monarch conservation because it is the only viable means of inferring origins at continental scales. Recently, routine simultaneous analyses of tissue δ 2 H and δ 18 O of organic materials has emerged leading to questions of whether the dual measurement of these isotopes could be used to more accurately infer spatial origins even though the two isotopes are expected to be coupled due to the meteoric relationship. Such refinement would potentially increase the accuracy of isotopic assignment of wintering monarchs to natal origin. We measured a sample of 150 known natal-origin monarchs from throughout their eastern range simultaneously for both δ 2 H and δ 18 O wing values. Wing δ 2 H and δ 18 O values were correlated (r 2 = 0.42). We found that wing δ 2 H values were more closely correlated with amount-weighted growing season average precipitation δ 2 H values predicted for natal sites (r 2 = 0.61) compared to the relationship between wing δ 18 O values and amount-weighted growing season average precipitation δ 18 O values (r 2 = 0.30). This suggests that monarch wing δ 2 H values will be generally more useful in natal assignments than δ 18 O values. Spatial information related to the use of deuterium excess in environmental waters was similarly found to be not useful when applied to monarch wings likely due to the considerable variance in wing δ 18 O values. Nonetheless, we recommend further testing of monarch wing δ 2 H and δ 18 O values from known natal sites with an emphasis on field data across a strong gradient in precipitation deuterium excess.

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 wellestablished, long-term overwintering sites do so without any previous experience of their locations. Despite being a highprofile 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 ( 13 C/ 12 C expressed as δ 13 C) and stable-hydrogen ( 2 H/ 1 H, δ 2 H) isotope ratio measurements of wings could provide information on monarch natal origins in Canada and the USA. However, δ 2 H measurements provided a more powerful means of inferring origins than δ 13 C measurements. Hydrogen in foodwebs is ultimately derived from precipitation. The amount-weighted average precipitation δ 2 H (expressed as δ 2 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 δ 2 H p and can be later measured to infer natal origins once the relationship between δ 2 H p and tissue δ 2 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 ( 18 O/ 16 O; δ 18 O) have not been used extensively to study animal movement (but see Bryant and Froehlich, 1995;Kohn, 1996;Hobson et al., 2004Hobson et al., , 2009aHobson et al., , 2012Bowen 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 δ 2 H and δ 18 O values in environmental waters are highly correlated via the global meteoric water relationship (Craig, 1961). Also, previously, measurement of δ 18 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 δ 18 O when combined with δ 2 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 δ 2 H-8 * δ 18 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 δ 2 H and δ 18 O values in insect wing chitin (δ 2 H w and δ 18 O w ) has the potential to contain important environmental and climate information.
The evaluation of both δ 2 H w and δ 18 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 δ 2 H and δ 18 O values among higher-order taxa. Hobson and Koehler (2015) reported that while feathers of American Redstart (Setophaga ruticilla) showed a good relationship between δ 2 H and mean growing-season precipitation δ 2 H (r 2 = 0.77), the same was not true for feather δ 18 O (r 2 = 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 δ 2 H w and δ 18 O w values (r 2 = 0.92). The question remains, then, TABLE 1 | Mean isotope values (± 1 SD ‰) of wild reared monarch butterflies from eastern North America and mean growing-season δ 18 O p (GSO) and δ 2 H p (GSD) at each location plus the 95% CI shown in parentheses from Bowen et al. (2005).  All data shown in Table S1.
whether δ 2 H w and δ 18 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 δ 2 H measurements alone. Here, we report measurements of δ 2 H w and δ 18 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.

Monarch Sample
Monarch wing material was obtained from samples reported in 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 δ 2 H p and δ 2 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 δ 2 H w and δ 18 O w values.

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 intersample variance due to isotopic effects from pigmentation  and weighed (0.35 ± 0.02 mg) into silver capsules. All samples were prepared for δ 2 H and δ 18 O analysis at the Stable Isotope Laboratory of Environment Canada, Saskatoon, Canada. Our approach involved the analysis of both δ 2 H and δ 18 O on the same analytical run (i.e., both H 2 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 2 was flushed to waste by FIGURE 1 | Relationship between δ 2 H w and δ 18 O w for monarchs raised outdoors at 31 known sites ( Table 1) throughout their range in 1996. Sample originally reported for δ 2 H w and δ 13 C w by Hobson et al. (1999)

Isotopic Assignments
We created δ 18 O w and δ 2 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 δ 18 O and δ 2 H (δ 18 O p and δ 2 H p ) (Terzer et al., 2013; IAEA/WMO, 2015) surfaces into separate δ 2 H w and δ 18 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 δ 2 H w and δ 18 O w to each calibrated wing isoscape (see results), separately. We used the residual SD error of 9.33‰ for δ 2 H w and 1.50‰ for δ 18 O w from regressions in our assignments. We limited assignments to the current known geographic range for the eastern monarch FIGURE 3 | Relationship between monarch δ 2 H w and amount-weighted mean growing-season average δ 2 H p (Bowen et al., 2005) for monarchs raised outdoors at 31 known sites ( Table 1) throughout their range in 1996.
Frontiers in Ecology and Evolution | www.frontiersin.org 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 δ 2 H w and δ 18 O w isoscapes represented a potential origin for a sample by using a normal probability density function (pdf) based on the observed δ 2 H w and δ 18 O w , and thus depicted the likely origins of each monarch by assigning individuals to the δ 2 H w and δ 18 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 (δ 2 H w , δ 18 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.
FIGURE 4 | Depiction of (A) the expected monarch δ 18 O w isoscape based on the relationship derived between δ 18 O w and amount-weighted mean growing season average δ 18 O p (Bowen et al., 2005) shown in Figure 2 and (B) the expected monarch δ 2 H w isoscape based on the relationship derived between δ 2 H w and amount-weighted mean growing-season average δ 2 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. Frontiers in Ecology and Evolution | www.frontiersin.org

Assignment Evaluation
Using the relationships established between δ 18 O w and δ 2 H w values and δ 18 O p and δ 2 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 δ 2 H w values, only δ 18 O w values and both δ 18 O w and δ 2 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 δ 2 H w values alone.

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 2 = 0.02). This suggests that little information associated with FIGURE 6 | Relationship between measured deuterium excess in known-origin monarch wings (Table 1) and predicted growing season precipitation deuterium excess for these locations.
precipitation environmental deuterium excess was transferred to monarchs per se.

DISCUSSION
Our results confirm a strong relationship between wing chitin δ 2 H values and precipitation δ 2 H values at known natal sites. This result underlines the strong assignment power of using Monarch wing δ 2 H values to infer origins and confirms that even though the earlier investigations of Monarch δ 2 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 δ 2 H (δ 2 H w ) and predicted amount-weighted growing season δ 2 H p (earlier method: δ 2 H w = 0.62 * δ 2 H p -79; r 2 = 0.69; recent method: δ 2 H w = 0.78 * δ 2 H p -77.4; r 2 = 0.61) and we recommend that from now on authors use the recent calibration relationship reported here. The δ 2 H w results contrasts with the much weaker relationship we found between δ 18 O w and δ 18 O p . Thus, while the measurement of δ 18 O w possibly confers additional information on origin, the strong relationship between δ 2 H w and δ 18 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 δ 18 O analyses. Nonetheless, we recognize the potential utility of using both δ 18 O w and δ 2 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 δ 18 O values do not seem to add much additional information on origins of Monarch Butterflies and other terrestrial taxa based on δ 2 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 δ 18 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 δ 2 H is the most preferred for single isotope assignments.
In addition to the earlier proof of concept papers on the use of both wing δ 18 O and δ 2 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.