Temporal stability of δ2H in insect tissues: Implications for isotope-based geographic assignments

Hydrogen isotope geolocation of insects is based on the assumption that the chitin in the wings of adult migratory insects preserves the hydrogen isotope composition (δ2H) of the larval stages without influence of adult diet. Here, we test this assumption by conducting laboratory feeding experiments for monarch butterflies (Danaus plexippus) including: (1) a starvation treatment where adults were not fed and (2) an enriched treatment where adults were fed a diet isotopically enriched in deuterium (~ +78‰) compared to the larval diet. The δ2H values of adult wings were measured at different time steps along the 24-day experiment. We also investigated intra-wing differences in δ2H values caused by wing pigmentation, absence of wing scales, and presence of major wing veins. We conclude that, although the magnitude of the changes in δ2H values are small (~6‰), wing δ2H values vary based on adult diet and insect age, particularly early after eclosion (i.e., 1–4 days). We found that wing shade, wing pigmentation, and the presence of wing scales do not alter wing δ2H values. However, wing samples containing veins had systematically higher δ2H values (~9‰), suggesting that adult diet influences the hemolymph that circulates in the wing veins. We hypothesise that there is a stronger influence of adult diet on the isotope signal of wings during early adult life relative to later life because of increased metabolic and physiologic activity in young insect wings. We argue that the influence of the isotopic contribution of adult diet is generally small and is likely minimal if the wings are carefully sampled to avoid veins. However, we also demonstrated that wings are not inert tissues, and that adult feeding contributes to some of the intra-population δ2H variance. We conclude that δ2H geolocation using insect wings remains valid, but that adult feeding, butterfly age and wing vein sampling generate an inherent uncertainty limiting the precision of geolocation.


Introduction
Investigating, understanding and predicting migration patterns has critical implications for species conservation, biosecurity and ecology (Meretsky et al., 2011;Lohmann, 2018;Satterfield et al., 2020). One effective tool to track migratory animals is through endogenous biochemical markers such as stable isotopes (Hobson et al., 2010). Biological and physical processes drive predictable variations of isotope compositions across the landscape (West et al., 2010), which can be used to generate isotopic baseline maps, referred to as isoscapes (e.g., Bowen and Revenaugh, 2003;West et al., 2010;Bataille et al., 2018Bataille et al., , 2020. The principle behind stable isotope geolocation is that metabolically inert tissues (e.g., enamel and keratin) will preserve the isotopic signature of the location where the tissue developed through the isotopic composition from local diet (Hobson and Wassenaar, 1996;Holder, 2012). To track the mobility and determine the natal origin of migratory individuals, the isotope composition of the metabolically inert tissue is compared to isoscapes of the study area to generate probabilistic maps of potential origin (Wunder, 2010).
Hydrogen isotopes have been used to track animal migrations for decades (e.g., Hobson and Wassenaar, 1996;Cryan et al., 2004;Britzke et al., 2012;Holder, 2012;Popa-Lisseanu et al., 2012;Flockhart et al., 2013Flockhart et al., , 2017Talavera et al., 2018;Hobson et al., 2020;Reich et al., 2021). The hydrogen isotope composition (δ 2 H) is provided in a delta notation: ( ) with R corresponding to the ratio of heavy isotopes to light isotopes of the sample (R sample ) or standard (R standard ) (Coplen, 2011). Hydrogen isotopes fractionate in precipitation and vary across the landscape following climate patterns as hydrogen fractionates due to the rainout effect and temperature (Bowen and Revenaugh, 2003). Hydrogen isotopes tend to become more depleted in heavy isotopes in cooler (higher latitudes), more inland and higher altitude conditions (Bowen and Revenaugh, 2003;Bowen et al., 2005), creating distinctive and predictable δ 2 H precipitation patterns. These unique δ 2 H precipitation patterns are inherited by local plants and animals in a predictable manner through diet, allowing tissue-specific δ 2 H isoscapes to be calibrated using the relationship between known-origin tissue δ 2 H values and δ 2 H in precipitation (Yapp and Epstein, 1982;Bowen, 2008;Hobson, 2019). Consequently, in migratory individuals, the δ 2 H values of inert tissues that developed at a specific location can be compared to tissuespecific δ 2 H isoscapes to estimate an individual's potential geographic origin (Bowen, 2008;Hobson, 2019;Hobson et al., 2010). Studies using the δ 2 H values of inert tissues including enamel or keratinous tissues (e.g., hair, feather and nail) have been used to track the mobility of many vertebrate species (e.g., Britzke et al., 2012;Popa-Lisseanu et al., 2012;Nordell et al., 2016;Voigt and Lehnert, 2019;Hu et al., 2020;Fauberteau et al., 2021;Kruszynski et al., 2021). However, estimating the natal origin of migratory invertebrates using δ 2 H values requires the existence of inert tissues in which the δ 2 H value of the earlier stages of life (i.e., larval stage) is preserved as the animal migrates (Hobson et al., 2010). In insects, neither enamel nor keratin is present (McKittrick et al., 2012). When applying isotope geolocation to track migratory insects, it has long been assumed that the δ 2 H values of chitinous insect wing tissues reflect the location of natal origin (i.e., larval stage; e.g., Hobson, 1999;Flockhart et al., 2013;Stefanescu et al., 2016;Talavera et al., 2018;Vander Zanden et al., 2018;Hobson et al., 2019). However, the assumption that insect wing tissues preserve the isotope signal of the larval stage has not been properly verified (but see Holder, 2012;Hungate et al., 2016;Morra et al., 2021). This assumption is based on the idea that insect wings are analogous to bird feathers and mammal hair, which are metabolically inert after tissue formation and are mainly composed of keratin, a structural scleroprotein (Hobson, 1999;Holder, 2012). Hair and feathers are inert (i.e., "dead") tissues that preserve the δ 2 H of the location of tissue formation, allowing isotope geolocation principles to be applied (Ehleringer et al., 2008;Hu et al., 2020). As hair grows continuously and feathers moult regularly, isotope geolocation principles can be applied to the hair and feathers to track adult migration and movement (e.g., Cryan et al., 2004;Fauberteau et al., 2021). In the case of insect wing tissues, it is assumed that isotope geolocation principles can also be applied because the tissues developed at the larval/pupa stage and were subsequently inert during adult life (Wassenaar and Hobson, 1998;Hobson, 1999). However, chitin, a structural polysaccharide (Ravi Kumar, 2000;Resh and Carde, 2009), fundamentally differs from keratin.
Although insect wings do not continuously grow or repair themselves following formation, these tissues are known to be metabolically active and do not qualify as inert tissues (Pass, 2018;Tsai et al., 2020). Insect wings are irrigated with hemolymph, the equivalent of blood in insects, which circulates in the insect for thermoregulation (Rawlins, 1980;Tsai et al., 2020) and for replacing water loss from evaporation (Wasserthal, 1983). Insect wings also require oxygen, water, nutrients, hormones, and, waste removal due to the presence of active "living" tissues, such as the epidermis and pheromone glands (Pass, 2018). The cuticle that makes up the wing membrane remains active during early adult life and undergoes sclerotization following eclosion. In the cuticle, hydrogen (in the form of OH) is used for cross-linking cuticular proteins during sclerotization (Wright, 1987;Hopkins and Kramer, 1992;Andersen, 2010). Due to the continued development of wing structures in the adult phase (i.e., tissue turn-over and tissue activity; Wright, 1987;Hopkins and Kramer, 1992;Andersen, 2010;Pass, 2018), food intake (Wyatt, 1961), and effects from metamorphosis (Morra et al., 2021), it is likely that the isotopic composition of the wing will change through time. As hydrogen isotopes, other isotope systems and trace elements are increasingly used in insect ecology studies, it is critical to investigate the potential influence of adult diet, insect age and, metabolic processes on insect wing isotope values.
In this study, we propose to revisit the assumption that the δ 2 H values of wing tissues of migratory insects are not altered by adult diet and/or insect age. We use monarch butterflies (Danaus plexippus) as our study species for a series of laboratory feeding experiments. Monarch butterflies are an ideal species for this study as: (1) they are easily sexed, (2) they have large wings allowing the same individual to be sampled multiple times to test intra-wing differences, (3) they can live up to 8 months allowing changes in isotope composition to be tested across several weeks of adult feeding, and (4) caterpillars develop on a single host plant (Asclepias spp.). Additionally, hydrogen isotope geolocation has been used for decades to investigate the multi-generational annual round trip of the species between Mexico, the United States, and Canada (e.g., Hobson, 1999;Flockhart et al., 2013;Vander Zanden et al., 2018;Hobson et al., 2019;Reich et al., 2021). The assumption behind these studies is that the wings of monarch butterflies are inert after eclosion and preserve the δ 2 H value of the larval stage throughout adult life . Adult monarchs often travel across long distances and feed throughout their adult life on flower nectar found along their migratory paths (Brown and Chippendale, 1974). As the monarch butterfly migrates, the isotopic composition of the nectar along its path likely diverges from that of the site of larval development (Yapp and Epstein, 1982;Bowen, 2008;Hobson, 2019). We hypothesise that wing δ 2 H may change through time and if some of the hydrogen ingested during adult feeding is incorporated into the wing tissue, it could influence the δ 2 H value of the wing.
The secondary focus of this study is to further test if wing δ 2 H values vary based on the wing sampling location by testing if there are intrawing differences in δ 2 H values in wings of migratory insects (i.e., monarch butterflies). As migratory insect wings are used for δ 2 H geolocation, the wing sampling procedure and sampling location on the wing can vary based on wing wear, wing shade, wing pigmentation, insect species and wing size. To know if wing δ 2 H can reliably be used  . In this study, we test if δ 2 H values vary based on wing pigmentation (orange vs. black), wing shade (dark orange vs. light orange), the presence of major veins and the absence of wing scales in the wing sample. As insect wing tissues are not metabolically inert (Pass, 2018;Tsai et al., 2020), we predict that wing δ 2 H values of migratory insects, specifically butterflies, will vary with the sampling location on the wing.

Feeding experiment
To test if the wings of migratory insects preserve the δ 2 H values of the larval stage, feeding experiments for monarch butterflies were designed. In summary, monarch caterpillars were fed a diet with a distinct δ 2 H value, then upon emerging from the chrysalis as adults (i.e., eclosion) they are either provided with no food or provided with an enriched deuterium diet. The approach used in this study is summarised in Figure 1.

Larval stage
In August 2021, 67 eggs were collected from laboratory-raised monarch butterflies from the Department of Biology at the University of Ottawa. All of the 67 eggs were collected from the same female monarch butterfly. The eggs were harvested and transferred to fresh milkweed leaves. At the larval stage, caterpillars were fed every day with wild common milkweed (Asclepias syriaca) collected daily from the same urban green area in Ottawa, Ontario (45°24′48″ N, 75°39′53″ W) between August and September 2021. The caterpillars were contained individually in mason jars with a screen netting and frass was removed daily. The caterpillars were kept at room temperature and exposed to natural daylight cues. A humidifier was placed in the room to mimic natural humidity to promote successful ecdysis. We expect that any surficial absorption of hydrogen by the milkweed leaves from the humidifier will have little to no influence on the δ 2 H of the caterpillars. Holder (2012) showed that the δ 2 H in insect wings, specifically Lepidoptera, are largely derived from plant or leaf solid and only minimally influenced by plant or leaf water. Additionally, the hydrogen provided by the leaf solid and plant water will far outweigh contributions of hydrogen absorbed on the leaf surface from the ambient environment, thus minimising any contribution of hydrogen from water vapour provided by the humidifier. One day after the caterpillars entered the pupal phase, the chrysalises were tied at the base and strung to ensure that the chrysalis would not fall or be squished by the side of the mason jar. As the mason jar may not be sufficient for the butterflies to grip, a thin piece of fabric was placed in the mason jar. By tying the chrysalises and placing fabric in the mason jar, the butterflies will have ample space to hang and pump their wings following eclosion. The larval and pupal stages took place from the end of August to mid-October 2021. The monarch butterflies enclosed from the chrysalises throughout September and October 2021. Of the 67 eggs, 56 monarchs reached the adult stage.

Adult stage
To test if an isotopically distinct adult diet from that of the larval stage alters the δ 2 H values in the wings of monarch butterflies, monarchs Flow chart summarising the methodology used in this study into three main stages with sub-steps depicting the procedure used to rear monarch butterflies in the feeding experiment, prepare wing samples, measure δ 2 H and perform statistical analyses.
Frontiers in Ecology and Evolution 04 frontiersin.org were assigned to one of two adult dietary treatments. In the first adult dietary treatment, the "enriched treatment, " adult monarch butterflies were fed a solution of 2:1 deuterium-enriched water (δ 2 H = +78‰): New Zealand honey. The deuterium-enriched water was obtained by boiling sequentially several litres of water to distil the heavy isotopes and regularly measured δ 2 H values using a Los Gatos Cavity Ringdown Spectroscopy instrument at the University of Ottawa. The δ 2 H value of the adult diet water (δ 2 H = +78‰) was highly enriched in deuterium relative to estimates for the larval diet water (δ 2 H = −62 ± 2.1‰; Bowen et al., 2005). We estimated the larval diet water based on the δ 2 H value for the Ottawa region of the global growing season precipitation hydrogen isoscape (Bowen et al., 2005). We assumed that plant tissues preserve the δ 2 H value of local precipitation with minimal change (Hobson et al., 2010;Holder, 2012). In the second adult dietary treatment, the "starvation treatment, " adult monarch butterflies were not provided with food or water. This treatment was used to test the influence of physiology on the wing δ 2 H values with no influence from the adult diet. To test if wing δ 2 H values are influenced by the butterfly's age following eclosion, henceforth referred to as butterfly age (days), the monarch butterflies were assigned to a specific sampling date (0-24 days following eclosion). In total, six monarch butterflies (three females, three males) were sampled on the day of their eclosion (i.e., day 0), which is the reference δ 2 H value for the larval stage without alteration from butterfly age and/or adult diet. The monarch butterflies in the starvation treatment were sampled at 2, 4 and 6 days following eclosion, whilst the monarch butterflies in the enriched treatment were sampled at 4, 6, 12, 18 and 24 days following eclosion. The design was unbalanced to accommodate the expectation that monarchs in the starvation treatment would show increasing rates of natural mortality with time. In total, 31 monarch butterflies (14 females, 17 males) were assigned to the enriched treatment and 19 monarch butterflies (10 females, 9 males) were assigned to the starvation treatment. Upon eclosion, the monarch butterflies that were not assigned to a day 0 sampling date were carefully transferred to netting cages. Each small and medium-sized cage contained one monarch butterfly, whereas the large cages contained two monarch butterflies that were separated by a piece of fabric. The adults were kept at room temperature and were not exposed to a humidifier or natural daylight, so there was a strict lights schedule paired with lights on a timer to mimic natural daylight cues (~12:12). The monarch butterflies in the enriched treatment were encouraged to feed on the 2:1 deuterium-enriched water: New Zealand honey mixture daily for a minimum of 2 minutes. Since the butterflies do not always recognise the artificial food source, they were hand-fed to ensure they would ingest the deuterium-enriched diet. To feed the monarch butterflies, the deuterium-enriched water and honey solution was put into small faux-flower tubes at the centre of a butterfly feeder. The monarch butterfly was gently captured and placed on the feeder, where a thin paintbrush was used to expose the proboscis and hold it in the deuterium-enriched water and honey mixture. If the butterfly began to feed on its own, the paintbrush was removed. No monarch butterflies were fed on the day of eclosion as the wings and proboscis were too delicate to be handled.
On the assigned sampling date, sampled monarch butterflies were not fed, irrespective of treatment, and were placed in an envelope. The envelope and butterfly were then weighed to obtain the butterfly wet weight. The monarch butterflies were killed by freezing and were stored at −18°C until stable isotope analysis. The monarch butterflies sampled on day 0 (i.e., day of eclosion), were killed by freezing 2-4 hours after emergence from the chrysalis.

Sample preparation
Before the δ 2 H values of the wings were measured, the wings underwent a sampling preparation procedure ( Figure 1) and a location on the wing was selected to undergo stable isotope analysis ( Figure 2). The wings were assigned a wing wear score (1-5) which represents the physical condition of the wing, 1 being wings in perfect condition and 5 being wings in which the majority of scales were removed and there is significant tearing of the wing (Flockhart et al., 2013). The wings were then removed from the thorax. The wings and bodies were desiccated at 65°C and the bodies were weighed (i.e., butterfly dry weight), including only the head, thorax, and abdomen to account for differences in physical damage (e.g., loss of scales and torn wings).
Next, the wings were cleaned using the 2:1 chloroform:methanol solution to remove lipids, dust, and contaminates from the wings that could have impacted the δ 2 H values of the wing samples (Paritte and Kelly, 2009;Hobson et al., 2017). The wings were washed three times for 60, 30, and 15 minutes, respectively, with agitation and were flipped halfway through each wash. The wings were air dried in a class-100 fume hood before the wing sample is taken.
A standard wing sample was taken in the lower distal portion of the wing that is dark orange, veinless, and scaled (location S, Figure 2). To test if there are intra-wing differences in δ 2 H values based on wing pigmentation, wing shade, presence of major veins, or the absence of scales in the wing sample, for each butterfly, a second sample of the wing was taken to be compared to the standard wing sample (Figure 2). The butterflies were accordingly assigned a secondary sample (i.e., black wing sample, light orange wing sample, wing sample containing vein and wing sample without scales) to ensure that age (i.e., 0-24 days), sex and dietary treatment (i.e., starvation or enriched) were distributed for all of the secondary samples. To test for differences in pigmentation, a black Annotated photo depicting the location of samples on the right hindwing. S is the location of the standard samples (dark orange, no vein and scaled), A is the location of black samples, B is the location of samples with vein, C is the location of light orange samples, and D is the location of samples with no scales (scales were manually removed). Results from samples from location S were used for linear regression modelling and the Tukey HSD test for the δ 2 H of wing samples analysis. Results from samples from locations S, A, B, C, and D were used for the paired t-tests for the analysis of intra-wing differences.
Frontiers in Ecology and Evolution 05 frontiersin.org pigmented wing sample was taken (location A, Figure 2). A sample was taken from a light orange section of the wing (location C, Figure 2), in contrast with the dark orange shade of the standard sample, to look for differences due to wing shade. To explore the influence of wing veins on intra-wing δ 2 H values, samples were taken across a prominent wing vein (location B, Figure 2). Finally, scales were manually removed with a paintbrush (location D, Figure 2) to detect any influence of wing wear scores on δ 2 H values.

Statistical analysis
We subdivided the dataset into four categories: (1) starvation, which were butterflies in the starvation treatment aged 2-6 days (n = 19), (2) early enriched, which were butterflies in the enriched treatment aged 2-6 days (n = 14), (3) late enriched, which were butterflies in the enriched treatment older than 6 days (n = 15) and (4) day 0, which can be considered the control treatment as they were not aged or exposed to the adult diet (n = 5). These four categories subdivide the dataset into comparable groups based on adult dietary treatment and time and are henceforth referred to as treatment-time. To test if the deuterium-enriched diet and/or insect age altered the δ 2 H values of wing tissues of monarch butterflies, a multivariate linear regression and Tukey HSD test were applied to compare the wing δ 2 H values of the monarch butterflies in the four treatment-time categories.
To test if there are intra-wing differences in δ 2 H values based on wing pigmentation (orange vs. black), wing shade (dark orange vs. light orange), the presence of major veins, and the absence of wing scales in the wing sample, paired t-tests were applied to compare the δ 2 H values of these groups. All of the data collected and used for statistical analysis can be found in Supplementary Table 1.

Multivariate linear regression
We ran a multivariate linear regression model to test the correlation between categorical, continuous and discrete predictor variables (i.e., adult dietary treatment and time categories, sex, dry butterfly weight (g), percent hydrogen of the wing sample and date of eclosion) and wing δ 2 H. For the multivariate linear regression samples, EL20, EL22 and EL23 were removed as outliers. EL20 is a monarch butterfly that was sampled on the day of eclosion and yielded a δ 2 H value of −115‰, which is an outlier from the rest of the dataset as seen in Supplementary Figure 1. EL22 and EL23 were lost during isotope analysis due to autosampler malfunction. The linear regression model was performed using the ols linear model function of the statsmodels. formula.api package in Python3 (Seabold and Perktold, 2010).
To test whether wing δ 2 H values (from sampling location S; Figure 2) are correlated to adult dietary treatment and time categories, sex, dry butterfly weight (g), percent hydrogen of the wing sample, and date of eclosion, a multivariate linear regression model was applied. Model assumptions of normality and homoscedasticity were visually confirmed before running the multivariate linear regression model (Supplementary Figure 2).
As there is a sample size of 53, five independent variables were selected for the linear regression model (abiding by the recommended 10 samples:1 variable ratio; Peduzzi et al., 1996). Sex and dry butterfly weight were selected as independent variables to test if the sex or size of the butterfly were correlated to wing δ 2 H values. Dry butterfly weight and wet butterfly weight have a direct linear correlation (Supplementary Figure 3), so either can be used as an independent variable in the model. The independent variables all had a normal or close to a normal distribution, but a log transformation was applied to dry butterfly weight (Supplementary Figure 4). The percent hydrogen of the wing sample (Supplementary Figure 5) was selected as an independent variable to test if the proportion of hydrogen in an individual wing sample is correlated with δ 2 H values, as the quantity of hydrogen in a given sample may differ based on adult dietary treatment. As individuals enclosed across 3 weeks and were fed milkweed collected during their larval stages, we tested whether the date of eclosion affected wing δ 2 H values (Supplementary Figure 6). The date of eclosion variable was first converted to a DateTime data type using the datetime function in Python3 (Van Rossum and Drake, 2009). Weight of wing sample and wing wear score were excluded as independent variables as they showed no univariate correlations to wing δ 2 H values, with an r 2 of 0.019 and 0.000, respectively (Supplementary Figures 7, 8).

Tukey HSD test
A post hoc multiple comparison Tukey HSD test was used to identify differences between the four treatment-time categories: (1) starvation, (2) early enriched, (3) late enriched and (4) day 0. Since the δ 2 H values of wing samples had a close to normal distribution, as seen in Supplementary Figure 2, a Tukey HSD test could reliably be applied to this dataset. The outliers removed for the multivariate linear regression model were also removed for the Tukey HSD test. The Tukey HSD test was performed using the stats pairwise Tukey function from the statsmodels. Regression.linear_model.OLSResults.t_test_pairwise package in Python3 (Seabold and Perktold, 2010).

Paired t-tests
Paired t-tests were applied to test the intra-wing differences in δ 2 H values based on different wing sample locations. We compared the δ 2 H values of the standard sample (i.e., dark orange, vein-less, scaled; Frontiers in Ecology and Evolution 06 frontiersin.org location S, Figure 2), to the δ 2 H of either a black wing sample (n = 12, location A, Figure 2), a light orange wing sample (n = 13, location C, Figure 2), a wing sample containing vein (n = 13, location B, Figure 2), or a wing sample with no scales (n = 12, location D, Figure 2), on the same individual. In addition to the previously removed outliers, EL2, EL19 and, EL54 were removed for the paired t-tests. EL2 had a sample weight of 0 mg likely due to an error during sample capsule packing. EL19 was removed as an outlier as it had a much higher wet butterfly weight (0.64 mg) compared to that of the others. EL54 was lost during isotope analysis due to autosampler malfunction. The paired t-tests were performed using the stats t-test function of the scipy package in Python3 (Virtanen et al., 2020).

Trends in wing δ 2 H based on adult diet and insect age
The trends for wing δ 2 H values from the standard location (location S, Figure 2) are observed in Figure 3. The mean δ 2 H value of wings from monarch butterflies sampled on the day of eclosion (day 0) was −104 ± 1.7‰ (n = 5; Figure 3). This is the reference δ 2 H value of the larval stages with no influence from adult diet, butterfly age and/ or physiology. The δ 2 H values of the wings of monarch butterflies from the enriched treatment were higher than those of the starvation treatment ( Figure 3). The mean wing δ 2 H value from the enriched treatment was −100 ± 4.8‰ (n = 29) and for the starvation treatment was −104 ± 4.0‰ (n = 19). The δ 2 H values of the wing samples appeared to vary based on the butterfly's age for both adult dietary treatments (Figure 3). The most notable changes to δ 2 H values occurred during the first 4 days following eclosion (Figure 3). δ 2 H values in both adult dietary treatments peaked at days 2 and 4, before leveling off towards the expected δ 2 H value from the larval stage (i.e., −104‰; Figure 3). The linear regression model explained approximately half of the variance in wing δ 2 H values (Adj. r 2 = 0.51; Table 1). The correlation of wing δ 2 H values with sex, dry butterfly weight, and percent hydrogen were not statistically significant (i.e., p-values of 0.958, 0.267 and 0.168, respectively; Table 1). Day of eclosion was significantly associated with wing δ 2 H values, with an average increase of ~0.4‰ per day (β = 0.37 ± 0.09, p < 0.001; Table 1 and Supplementary Figure 6). The treatment-time categorical variable was also significantly associated with wing δ 2 H values (Table 2; Figure 3). The comparisons of treatment-time groups using the post-hoc Tukey HSD test found significant differences between the early enriched and day 0 categories (p = 0.004), and the early enriched and starvation categories (p = 0.03) ( Table 2). The early enriched category was on average more enriched than the day 0 category by ~6‰ and more enriched than the starvation category

Intra-wing differences in δ 2 H values
Tests of intra-wing δ 2 H variability found that the δ 2 H values between the paired variables of wing shade (dark orange vs. light orange), wing pigmentation (orange vs. black) and scaled vs. non-scaled wing samples were not significantly different (t-test; p-values of 0.276, 0.065 and 0.950, respectively; Table 3), with mean wing δ 2 H differences of −1.1‰, −2.1‰, and, 0.06‰, respectively, (Table 3; Figure 5). The mean δ 2 H of wing samples with veins compared to wing samples without veins on the same individual were significantly different (t-test; p < 0.001; Table 3), with samples with veins having δ 2 H values ~9‰ more enriched than those without veins (Table 3; Figure 5).

Identifying non-significant controls on δ 2 H values in butterfly wings
As indicated by the multivariate linear regression model, sex, percent hydrogen of wing sample and dry weight of the butterfly do not have a significant influence on δ 2 H values in wings (Table 1). This experiment demonstrates that if metabolism influences the fractionation of hydrogen isotopes in wings, this influence is not related to the adult body mass and/or sex of the individual, which supports the findings of Hobson et al. (2017). When using isotopes to trace the migration of butterfly species, the sample set often includes individuals of different sizes. For example, monarch butterflies have different phenotypes depending on the summer or overwintering generations (Altizer and Davis, 2010;Flockhart et al., 2017). Similarly, sex is often challenging to identify for some butterfly species (e.g., Papilio victorinus and Papilio glaucus; Glassberg, 2017). The lack of influence of morphology and sex validates the use of δ 2 H geolocation tools for a sample set composed of randomly collected individuals (e.g., Flockhart et al., 2013).
Wing shade (dark orange vs. light orange), wing pigmentation (orange vs. black) and the absence of scales show an insignificant and negligible isotopic effect on wing δ 2 H values (Table 3; Figure 5). Our findings differ from those of Hobson et al. (2017), who demonstrated that wing pigmentation (i.e., black vs. orange) and the absence of wing scales influenced wing δ 2 H values, though this effect was small with little variation between wing quadrants. Hobson et al. (2017) also demonstrated that proper sampling procedures and wing washing (i.e., 2:1 chloroform:methanol solution) minimised the intra-wing variations of δ 2 H values . We demonstrate that worn wings have approximately the same isotopic composition as fresh wings ( Figure 5). This is important because migratory insects travel great distances, and their wings tend to be worn with significant scale loss (Flockhart et al., 2017). Similarly, sampling different portions of the wings with distinct pigmentation and shades will have only a small impact on the measured δ 2 H values as long as veins are avoided ( Figure 5). The combined findings of this study and Hobson et al. (2017) further support the applicability of hydrogen isotope geolocation of migratory insects' wings.  Frontiers in Ecology and Evolution 08 frontiersin.org

Identifying significant controls on δ 2 H values in butterfly wings
On day 0 (day of eclosion), most individuals have a similar wing δ 2 H value of −104 ± 1.7‰ (n = 5; Figure 3). If we calculate the expected δ 2 H values in wings using an existing calibration equation (i.e., δ δ . ) from known-origin monarch butterfly individuals from Hobson et al. (2019) and local growing season precipitation (GSP) in Ottawa (−62 ± 2.1‰; Bowen et al., 2005), we find an expected δ 2 H value of −101 ± 3.5‰ for local monarch butterfly wings. This value is within the uncertainty of the measured δ 2 H values, confirming the major and predictable influence of larval diet on the wing δ 2 H value. However, we also notice that at a single location with a well-controlled laboratory setting the intra-population δ 2 H variance at day 0 is only slightly larger (~2.27‰, n = 5) than the analytical uncertainty (>2.00‰; Figure 3). Interestingly, we also notice that the δ 2 H variance is larger on most days (e.g., starvation treatment days 4 and 6, enriched treatment days 4, 6, 12 and, 24) than on day 0, underlining the probable role of metabolic processes in fractionating the isotopes differently in each adult (Figure 3).
We found that wing δ 2 H values increase with later eclosion dates (Table 1; Supplementary Figure 6). As the milkweed leaves were collected daily throughout August and September 2021, changes in milkweed senescence and local weather events (i.e., thunderstorms, changes in temperature and/or humidity) may have altered milkweed δ 2 H values throughout the feeding experiment, potentially altering larval diet δ 2 H values and driving the observed changes in wing δ 2 H values associated with the date of eclosion. Considering that all milkweed leaves were collected within a month and on the same 30 m 2 patch of land, this is a surprising result but underlines the potentially rapid changes in δ 2 H values in soil water influencing larval diet and δ 2 H values of known origin individuals. Alternatively, but less likely, the effect of eclosion dates on δ 2 H values might be a result of the conditions in which the larvae were reared which were not perfectly controlled, such as the δ 2 H value of ambient water vapour (e.g., from the humidifier). To avoid such an effect, we recommend that future studies rear the caterpillars and butterflies in growth chambers where temperature and humidity can be tightly controlled.
As the early enriched and starvation categories occur over the same time span (6 days or less), the significant difference of ~4‰ in wing δ 2 H values between these categories provides evidence that adult diet influences wing δ 2 H values and thus wings are not inert tissues (Tables 1,  2; Figures 3, 4). Our work investigating the δ 2 H values based on wing sampling location further supports the hypothesis that adult diet   Figure 5). We suggest that wing veins have greater variation in δ 2 H values, and greater influence from the adult diet because they are made of metabolically active endocuticle tissues. Endocuticle tissues undergo active physiological processes (Pass, 2018;Tsai et al., 2020) and contain a large volume of metabolically active hemolymph (Wyatt, 1961;Emmel, 2012;Pass, 2018) which provides constant hydration to the wing veins (Wasserthal, 1983); therefore, it is expected that the δ 2 H value of the wing veins may vary from that of the rest of the wing, as observed in this study. This second experiment confirms that hemolymph and veins have a significant portion of newly metabolised hydrogen influenced by the adult diet and/or metabolic processes. We recommend that, as much as possible, veins should be avoided when using insect tissues for geolocation. This might not be currently feasible for small species where entire wings are needed for δ 2 H analysis or for other isotope analyses that require larger wing samples (e.g., Sr wing analysis; Reich et al., 2021). However, with advances in mass spectrometry, the required sampling mass is rapidly decreasing and avoiding veins should become a more common practise. We suggest that the pronounced δ 2 H change during early adult life (i.e., 2-6 days post-eclosion) is strengthened by early adult wing formation processes (Figure 3). The early enriched and day 0 categories had the greatest significant mean δ 2 H difference (~6‰) relative to any of the other compared categories, but the early enriched category was only enriched over the starvation category by ~4‰ (Table 2; Figure 3). We suggest that the remaining difference between the day 0 and starvation categories (~2‰), although not significant, is likely due to an additional influence from non-dietary related processes, which also further enriched the wing δ 2 H values in the early enriched treatment. Following eclosion, the wings are wet and fragile, and an increased hemolymph volume is pumped through the wing veins to facilitate hardening (Andersen, 2010;Emmel, 2012), whereas hardened mature wings do not have as great a volume of hemolymph circulating through the veins (Wasserthal, 1983;Emmel, 2012). We argue that during early adult wing formation, physiological processes (e.g., sclerotization, wing hardening) also fractionate hydrogen isotopes creating further isotope variations in the wing in early adult life. This is demonstrated as wing δ 2 H values appeared to vary the most between 2 and 4 days following eclosion before converging towards the day 0 δ 2 H mean (Figure 3). The hemolymph circulating in the wings is a fluid partially metabolised from the recent adult diet (Wyatt, 1961). We suggest that the decreased hemolymph volume in hardened mature wings (Wasserthal, 1983;Emmel, 2012) limits the influence of adult diet and/or physiology on wing δ 2 H values. Thus, the δ 2 H values of wings in the early enriched treatment increasingly diverge from their initial values towards those of the adult diet due to the high volume of deuterium-enriched hemolymph filling the wings (Wyatt, 1961;Wasserthal, 1983;Emmel, 2012). Box and whisker plots of intra-wing differences in δ 2 H (‰), comparing the standard wing samples (from location S) to (A) wing pigmentation (orange vs. black), n = 12; (B) wing shade (dark orange vs. light orange), n = 13; (C*) vein or no vein wing samples, n = 13; and (D) scaled or non-scaled wing samples, n = 12. There is a mean difference of ~9‰ between samples containing veins and samples with no veins. Panel C* has an asterix to indicate statistical significance between the compared intra-wing differences in δ 2 H (‰; i.e., vein or no vein wing samples). Black diamond points indicate outliers.
Frontiers in Ecology and Evolution 10 frontiersin.org However, as the individual ages, the δ 2 H value of the wings trends back towards the initial values (Table 2; Figure 3). Fortunately, most sampling of migratory insects occurs at later stages of adult life after the individual migrated or dispersed. Alternatively, as observed in previous studies (e.g., McCue and Pollock, 2008;Doi et al., 2017) starvation conditions can lead to an isotopic starvation effect, where in our study this could account for the difference in δ 2 H values between the day 0 and starvation categories. However, the starvation treatment was short (fewer than 6 days) and, prior to euthanizing the butterflies by freezing, all individuals were active and in overall good health. Therefore, we do not expect that the effects of starvation are driving our isotopic results. Nevertheless, as insect wings are not inert, as shown in this study and by Pass (2018), Tsai et al. (2020), and Morra et al. (2021), when an individual travels across large distances with varying isotopic values, adult diet and physiology could contribute a small amount to the isotope variability.

Isotope element cycling and isotope-based geographic assignments for migratory insects
We conclude that there is cycling of hydrogen isotopes in the wings of monarch butterflies which is influenced by the butterfly's age and adult diet. Our results support those of a previous feeding experiment performed using moths (Helicoverpa. armigera; Holder, 2012). We calculated that, for monarch butterflies, up to 8% of hydrogen isotopes in wing tissues are sourced from the adult diet (3.0 ± 3.2%; Supplementary Figure 9). We also found that early adult wing formation alters hydrogen isotopes by up to 2% (0.7 ± 0.9%; Supplementary Figure 10). These proportions may differ from those found by Holder (2012) for moths due to physiological differences between the insect species, and/or differences between experiments. However, the results from these studies support our hypothesis that wing tissues of migratory insects are not metabolically inert and that wing δ 2 H will be altered by the isotopic composition of the adult diet and/or the insect's age.
We found that the influence of adult diet and physiology is small (average ~6‰ for an isotopic difference of ~140‰ between larval and adult stage diets) relative to the uncertainty generally associated with hydrogen-based geographic assignments of migratory insects. For example, the spatially explicit map of isoscape uncertainty in Reich et al. (2021) ranges from 6 to 19‰. This high uncertainty is due, in part, to the high intra-site δ 2 H variance of known-origin individuals used to calibrate the tissue isoscape (Ma et al., 2020). Adult diet and physiology likely contribute to the high intra-site δ 2 H variance observed in previous studies (e.g., Hobson et al., 1999Hobson et al., , 2019Morra et al., 2021). Whilst this high intra-site variance does not preclude using δ 2 H geolocation, it sets inherent limits of this geolocation tool to provide specific geographic assignments for migratory insects. Other possible sources of intra-site isotope variance such as variations in the δ 2 H of plants and water at one given site, should be further explored to develop a better understanding and prediction of this uncertainty (e.g., Magozzi et al., 2020). To minimise the uncertainty identified in this study, researchers should be careful to avoid sampling wing portions with veins, and they should also consider sampling individuals of similar ages (using for instance similar wing worn indices). Additionally, appropriate uncertainty propagation should be performed when relating the δ 2 H of wing chitin and precipitation isoscapes to account for the inherent uncertainty of this geolocation tool (Ma et al., 2020). Using simple transfer equations relating δ 2 H of wings with that of precipitation without accounting for the full uncertainty propagation (Ma et al., 2020) will yield overly optimistic geographic assignments.

Limitations and future studies
This study used a sample size of 53 monarch butterflies collected from the same mating pair. Future studies could test the isotope effect in a feeding experiment on individuals selected from different mating pairs to evaluate if natural population diversity would incorporate additional isotope uncertainty. As the date of eclosion had a significant effect on our study, individuals with identical dates of eclosion should be used in future studies to better control for potential isotopic differences in larval diets. Additionally, the monarch butterflies raised in this study could have been raised in growth chambers to minimise physiological and environmentally driven isotope fractionation. Studies testing the role of adult diet and physiology could also be conducted for other insect species of interest, given the diversity of natural-history strategies, including diets, in this large taxonomic group. The mechanisms of δ 2 H changes over time could also be explored by expanding on recent work by Morra et al. (2021) by looking at how compound-specific (e.g., amino acids, lipids) δ 2 H values of different life stages and tissues (e.g., caterpillars, frass, exuviae and wings) change over time. We also recommend testing other isotope systems [e.g., Sr (e.g., Reich et al. this issue), C, O, S etc.] in wings to further ensure that wings preserve the isotopic signal of the location of larval development over time. As there were only two adult dietary treatments, we propose further studies with multiple adult dietary treatments of varying isotopic compositions to have a deeper understanding of the role of adult diet in wing tissues of migratory insects, especially under natural conditions.

Summary
Through a feeding experiment with two adult dietary treatments, we tested whether the δ 2 H values of wing tissues of migratory insects (i.e., monarch butterflies) were altered by adult diet and/or insect age. We found that the δ 2 H value of the wings did not vary based on sex and body mass but showed significant effects from the adult diet, age, and the date of eclosion. We also showed that the intra-wing δ 2 H values did not vary based on condition or pigmentation/shade of the wing but varied based on the presence of veins. We conclude that the effect of adult diet and age remains small and that the δ 2 H values of wings can reasonably be used as an endogenous biochemical marker to track migratory insect species. We predict that for wild monarch butterflies where the difference between larval and adult diet is likely much smaller than in this study, the influence of adult diet will be minimal and potentially negligible (<4‰). However, our results suggest that wings cannot be assumed to be "dead" or "inert" tissues, and should no longer be referred to as such in isotope literature. We also caution researchers that adult diet and physiological processes contribute to a fraction of the large intra-site isotope variance observed in this isotopic system. Interestingly, the changes in δ 2 H values based on eclosion dates underline the rapid isotopic variation of the larval diet and its impact on δ 2 H values of known-origin individuals. In summary, we recommend that researchers carefully select wing samples from individuals of approximately the same age, collect samples from the same location and time period, avoid sampling veins, and take samples from the same location on the wing to minimise δ 2 H variation and the differential isotopic contribution from adult diet and insect age. We also recommend that researchers fully propagate the δ 2 H uncertainty in their isoscapes as proposed by Ma et al. (2020). Whilst δ 2 H geolocation has clear limitations and uncertainties, it remains a powerful tool for tracking the mobility and the natal origin of monarch butterflies and other migratory insect species in biological, ecological and conservation studies.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary material; further inquiries can be directed to the corresponding authors.
Author contributions CB, PH, and MR designed the experiment inspired by PH, 2012 Ph.D. Thesis. EL and MR reared monarch butterflies and carried out the feeding experiment. Sample preparation for stable isotope analysis and statistical analyses was performed by EL. EL led the writing of the manuscript. All authors contributed to the article and approved the submitted version.