Understanding metrics of stress in the context of invasion history: the case of the brown treesnake (Boiga irregularis).

Invasive species can exert rapid depletion of resources after introduction and, in turn, affect their own population density. Additionally, management actions can have direct and indirect effects on demography. Physiological variables can predict demographic change but are often restricted to snapshots-in-time and delayed confirmation of changes in population density reduces their utility. To evaluate the relationships between physiology and demography, we assessed metrics of individual and demographic stress (baseline and 1-h corticosterone (CORT), body condition and bacterial killing ability) in the invasive snake Boiga irregularis on Guam collected in intervals of 10-15 years. We also assessed potential discrepancies between different methods of measuring hormones [radioimmunoassay (RIA) versus enzyme immunoassay (EIA)]. The magnitude of difference between RIA and EIA was negligible and did not change gross interpretation of our results. We found that body condition was higher in recent samples (2003 and 2018) versus older (1992-93) samples. We found corresponding differences in baseline CORT, with higher baseline CORT in older, poorer body condition samples. Hormonal response to acute stress was higher in 2018 relative to 2003. We also found a weak relationship between circulating CORT and bacterial killing ability among 2018 samples, but the biological significance of the relationship is not clear. In an effort to develop hypotheses for future investigation of the links between physiology and demography in this and other systems, we discuss how the changes in CORT and body condition may reflect changes in population dynamics, resource availability or management pressure. Ultimately, we advocate for the synchronization of physiology and management studies to advance the field of applied conservation physiology.


Introduction
Due to biogeographic isolation, perturbation of island ecosystems can have rapid, cascading, destabilizing effects (Nilsson and Grelsson, 1995). A particularly salient type of ecosystem perturbation occurs when invasive species are introduced to islands (Fritts and Rodda, 1998;Spatz et al., 2017). Beyond direct impacts on native species, invasive species also affect the environments where they are introduced, influencing resource availability and potentially their own population dynamics (O'Dowd et al., 2003). These effects may be magnified on islands such that after establishment, invaders must respond to resource and population changes more rapidly than in mainland invasions, where population expansion and carrying capacity may not be geographically restricted to the same degree (Ricciardi, 2007). Additionally, the added pressure of removal efforts from managers may trigger changes in demography, similar to responses to commercial harvest (de Buffrénil and Hémery, 2002;Hamilton et al., 2007). Invasion success in such dynamic environments may be influenced by plasticity in traits influencing reproduction, growth and dispersal (Richards et al., 2006;Parker et al., 2013). Due to their capacity to rapidly respond to dynamic ecosystem changes, successful island invaders are an interesting case study in understanding ecological resilience of animal populations (Sax et al., 2007). To facilitate the successful efforts made to both conserve native species and manage invasive species, it is integral to understand the resilience (or lack thereof) of invasive species to ecosystem change.
Ecological resilience is difficult to measure, as long-term monitoring of animal population responses to resourceavailability, management or environmental change is costly and time-consuming (Wintle et al., 2010). Individual physiological and morphological responses can be useful to highlight demographic trends in populations (Young et al., 2006). In some cases, physiological metrics may reveal changes before population-level effects are apparent, e.g. predicting increased likelihood of mortality in salmonid fish before an increase in death rate is observed (Ham and Pearsons, 2000;Cooke et al., 2012). Unfortunately, evidence for direct links between individual-level physiology and population change are rare (Cooke et al., 2012;Bergman et al., 2019), perhaps owing to the rarity of long-term studies due to uncertainty in long-term funding stability, caution in interpreting variables obtained by different methods (Szeto et al., 2011) and/or the requirement of many funding agencies and publication outlets that studies be 'novel' to advance scientific knowledge. Although problematic, a 'snapshots-in-time' approach is often the only post hoc option for evaluating how physiological data inform later population change in wildlife conservation or invasive species management (McCormick and Romero, 2017). Because demographic signals often lag behind physiological metrics (McCormick and Romero, 2017) it is important to re-assess published physiological data in the context of later observed demographic change to better understand their utility.
Commonly employed metrics for assessing individual and demographic stress, i.e. perturbations to individual animal health and subsequent population trends, include condition indices (Stevenson and Woods, 2006), glucocorticoid hormone levels (Bonier et al., 2009) and indices of immune function (Downs and Stewart, 2014). Condition indices serve as measures of an individual's ability to convert resources into mass relative to conspecifics. When measured across a segment of the population or across time, condition indices can be interpreted as indicators of resource availability (Johnson, 2005;Smith and Iverson, 2016) or population trends (Jennings et al., 2006). Glucocorticoid hormone levels are another commonly employed metric. Glucocorticoid hormones have a wide range of functions across species, as they can influence glucose and protein metabolism, reproduction and offspring fitness (reviewed in Wingfield et al., 1998;Wingfield and Sapolsky, 2003;Landys et al., 2006;Love et al., 2013;Romero and Wingfield, 2015). Glucocorticoids are often equated with 'stress' because circulating levels of these hormones can be altered during food shortage, after aggressive encounters and after acute confinement (Vera et al., 2017;MacDougall-Shackleton et al., 2019). We strive to avoid confounding glucocorticoids with stress in this manuscript. Understanding of the direct effects of glucocorticoid hormones across vertebrate species is still a work in progress (Jessop et al., 2013;Vera et al., 2017;Romero and Gormally, 2019), yet these hormones have some demonstrated utility as an indicator of population status in wild animals (Kitaysky et al., 2007;Escribano-Avila et al., 2013;Sorenson et al., 2017). Another implicated effect of fluctuating glucocorticoids is their influence on immune function (Martin, 2009), which itself is gaining utility as a metric of individual and population health (Demas et al., 2011;Downs and Stewart, 2014). Immune function is important for immediate defense and survival against pathogens. Investment in immune defenses may be altered when other energy demands take priority (i.e. reproduction, dispersal, self-maintenance during starvation; reviewed in Downs and Stewart, 2014). Immune investment may also serve as an important indicator of ongoing disease processes in a population (Hawley and Altizer, 2011). Condition, glucocorticoid levels or immune function are not necessarily expressed independently within individuals, as these factors may respond to similar ecological pressures and physiological pathways. Interpretation of multiple synergistic physiological variables as a 'stress syndrome' may thus better inform changes in wild animal health or population trends rather than relying on a single metric (Madliger and Love, 2015;Sandfoss et al., 2020). We thus chose to investigate changes in synergistic physiological variables over time in an infamous invasive predator.
The brown treesnake (Boiga irregularis, hereafter, BTS) on the island of Guam is a well-documented example of an introduced vertebrate influencing its own resource dynamics via extirpation of native avian prey (Savidge, 1987  . BTS remained largely unstudied until the mid-1980s . BTS research largely began after they were implicated as the cause of plummeting native bird populations (Savidge, 1987), and focus turned to the now island-wide snakes as invasive predators (Engbring and Fritts, 1988). Since many studies have been conducted on the snake with the goal of preventing the spread of BTS to other Pacific islands and the control, reduction and potential eradication of BTS on Guam (reviewed in Engeman et al., 2018) to facilitate the conservation and re-introduction of near-extinct avian species such as the flightless Guam rail (Gallirallus owstoni). BTS control has been difficult, complicated by the introduction of other nonnative prey species to Guam that may have subsidized the BTS population and allowed for continued growth after the loss of native birds (Savidge, 1988;Fritts and Rodda, 1998;Rodda et al., 1999;Christy et al., 2007a;Campbell et al., 2012;Olson et al., 2012;Wostl et al., 2016;Fig. 1). Management pressure has recently increased on portions of the BTS population (Engeman et al., 2018;Siers et al., 2018;Fig. 1). Potentially, as a result of prey availability, management or other factors, BTS population estimates have fluctuated over time (Savidge, 1991;Fig. 1). BTS on Guam are perhaps the longest continuously studied vertebrate invasion, including more recent studies on stress and reproductive physiology (e.g. Mathies et al., 2001;Moore et al., 2005;Waye and Mason, 2008;Aldridge et al., 2010;Mathies et al., 2010). The nature of this wellstudied system affords opportunities to integrate historical and recent data, specifically to assess physiological responses over time.
By combining previously and recently collected data, we investigate changes in metrics commonly associated with demographic stress over time to assess the response or resilience of island invaders as they react to a changing environment. Specifically, we explore the changes in BTS body condition and corticosterone (CORT; a glucocorticoid hormone) at-capture and following acute confinement stress (Moore et al., 2005;Waye and Mason, 2008) to recently collected data (reported herein). We hypothesize that body condition and CORT are negatively correlated and discuss how these variables and their relationships may be altered by changes in population density, prey availability and management pressure. We also investigate the hypothesis that CORT influences immune function by assessing the relationship between CORT and plasma bactericidal capacity. We discuss our results in the context of the dynamic environment on Guam and offer hypotheses for future research avenues using this infamous island invasion as a system for continued eco-physiological investigation and conservation efforts.

Snake sampling
We collected adult, reproductively active BTS (greater than 90-cm snout-vent length, hereafter SVL; Savidge et al., 2007) in 1992 and 1993 (Moore et al., 2005), 2003 (Waye and Mason, 2008) and, in 2018, by visual searching with headlamps for 3 h immediately following sunset. In all years, sampling was concentrated in the northern portion of the island in primarily limestone forest habitats and adjacent to rural roads. Immediately upon capture, we bled snakes (hereafter referred to as baseline) from the caudal vessels by heparinized syringe (2003,2018) or via collection of blood after decapitation . The specific timing of handling before blood draw was not recorded for 1992-93 and 2003 but is estimated to be within 3-5 min of capture (Moore et al., 2005;Waye and Mason, 2008). In 2018, we used stopwatches from disturbance of snake to blood collection, resulting in baseline samples collected in a mean of 4 min and 46 s (± 1 min and 56 s SD) after capture. In 2003 and 2018, snakes were subjected to acute confinement in a cloth bag for 1 h, after which we obtained a second blood sample (hereafter, 1 h). In 2018, we centrifuged blood to separate plasma and then froze samples in liquid nitrogen vapour in the field following blood draws to prevent complement degradation for our immune assay. In 1992-93 and 2003, we separated plasma within 12 h of sample collection, then stored samples frozen at −20 • C until shipment from Guam and at −70 to −80 • C until analysis (Moore et al., 2005;Waye and Mason, 2008). We measured the SVL of snakes with cloth measuring tape and mass using a spring scale. In 2018, snakes were euthanized via isoflurane inhalation followed by intracardiac injection of potassium chloride. All procedures were conducted according to approved animal care and use protocols (Moore et al., 2005;Waye and Mason, 2008; University of Florida IACUC 201709774).

CORT assay
We assessed total blood CORT via radioimmunoassay (RIA) for 1992,1993, 2003. RIA methods for plasma from 1992-3 and 2003 are described in Moore et al. (2005) and Waye and Mason (2008), respectively. Hormone measurements can vary both by laboratory and among methods (e.g. Nizeyi et al., 2011;Szeto et al., 2011;Yadav et al., 2013;Fanson et al., 2016). Thus, we assessed plasma from 2018 using two methods [RIA and enzyme immunoassay (EIA)] from identical aliquots to evaluate the potential issues with multistudy comparisons of BTS hormones. We assessed plasma from 2018 via RIA using the same reagents and techniques as described in Moore et al. (2005) at Virginia Tech. Intra-assay coefficient of variation for 2018 RIA samples was 5.5%. The within-year interassay variation for RIA was <18%.  Note that estimates of snake density should be interpreted with caution, as studies may not be directly comparable due to different methods. Introduction dates are estimates based on the literature. References: (i) Savidge, 1987;(ii) Austin et al., 2011;(iii) Fritts and Rodda, 1998;(iv) Christy et al., 2007a; Savidge, 1988;(vi) (vii) Engeman 1998;(viii) Wiewel et al., 2009;(ix) Engeman et al., 2018;(x) Wostl et al., 2016;(xi) Olson et al., 2012;(xii) (xiii) Siers et al., 2018. for parallelism and quantitative recovery. We diluted plasma samples 1:100 with assay buffer and plated them in duplicate according to the protocol provided by the kit manufacturer. We determined the optical density of each well at 450 nm with a plate reader (Epoch, BioTek); results are reported in ng/mL. Sample pairs (baseline and 1 h) from all individuals were assigned haphazardly to each plate. Average recovery was 99.14%, and the kit limit of detection is 16.9 pg/mL. The mean intra-assay coefficient of variation was 3.93%.

Bacterial killing assay
We employed a functional immune assay that evaluates the ability of complement and other anti-microbial components of the blood plasma to kill or inhibit growth of bacteria (modified from French et al., 2012). We performed all bacterial killing assays within 30 days of collection. Briefly, we diluted a working solution of Escherichia coli (ATCC no. 8739) to a 10 3 colony forming units solution in sterile phosphate buffered saline (1 M PBS; Lonza VWR catalog no. 12001-67). We plated samples in triplicate in a sterile 96-well plate. In each well, we added 5 μL of E. coli solution to 18 μL of diluted plasma and incubated for 30 min at 37 • C. After this initial incubation period, we added 125-mL sterile tryptic soy broth (TSB) and scanned plates for initial absorbance at 300 nm using a microplate reader. We then incubated plates for 12 h at 37 • C, after which we scanned them again to quantify growth of E. coli. Each plate had six positive control wells containing E. coli, PBS and TSB and six negative control wells containing PBS and TSB matching the volume and concentrations of the wells with samples. Rather than plating a single, 'optimal' plasma dilution (e.g. the dilution where an average of all samples kills 50% of bacteria; French et al., 2012), we used a dilution series to determine the dilution at 50% killing capacity of BTS plasma. We initially diluted each plasma sample to 1:16 with PBS (where most individuals showed 100% killing in pooled samples), then serially diluted by half until 1:128 (where most individuals showed 0% killing in pooled samples) and conducted the assay as above. To calculate dilution at 50% bacterial killing capacity of BTS plasma, first, the difference in absorbances for each sample set are calculated (Equation 1), followed by calculation of percent bacterial killing ability, accounting for potential variation in TSB absorbance across plates (Equation 2). We assigned samples with calculated percent killing above 100 or below 0 as 100 or 0 for further analyses, respectively. For each individual sample dilution series, we generated a 4-parameter logistic regression curve using the nlme package in R (Pinheiro et al., 2020;version 3.4.0 Patched), setting dilution as x values and the percent killing ability at each dilution as the y values, then extracting the x-value from each curve where y = 50% killing. This x value (or dilution at 50% bacterial killing ability) was used in all subsequent bacterial killing ability analyses. A sample with a lower x value indicates that a lower concentration of plasma was necessary to achieve 50% killing of E. coli; thus, a 1:32 dilution represents better bacterial

Conservation Physiology • Volume 9 2021
Research article killing ability than a 1:16 dilution.
where abs 12 represents the absorbance in a single well at 12 h and abs 0 at 0 h and n represents the number of replicates for each sample set within a plate (n = 3 for each individual sample dilution and n = 6 for positive control or negative control, respectively).

Percent Bacterial Killing Ability
where pos, neg and samp represent the results of Equation 1 applied to the set of positive control wells, negative control wells and individual sample dilution wells for each plate, respectively. Subtraction of neg is applied as a correction for variation in TSB absorbance across plates.

Body condition
We calculated body condition index by extracting the residuals from a cubic regression of natural log-transformed mass and natural log-transformed SVL from a larger morphometric dataset comprised 526 BTS collected on Guam from 1991-93, 2003 and 2018. Fat content was unavailable from 1992-93 snakes and specimens are unavailable for years prior to 2018; thus, we were unable to use more sophisticated methods of calculating body condition (Falk et al., 2017) to conduct this long-term comparison.

Statistical analyses
Raw data on snakes from 1992-93 and 2003 were acquired from I.T. Moore and H. Waye, respectively. Data for baseline CORT and body condition were available for all years; to account for potential body condition variation among seasons, we restricted data to samples collected at the end of the dry season (March through May: 1992 n = 18; 1993 n = 14; 2003 n = 15; 2018 n = 23). Data for 1 h CORT was available for 2003 (baseline n = 46; acute stress n = 14) and 2018 (baseline n = 23; 1 h n = 24); analyses restricted to samples collected in spring versus year-round did not change the results. We report year-round analyses as this increased the number of 1 h samples from 2003. We conducted bacterial killing assays on 2018 samples (baseline n = 21; 1 h n = 22). We used linear models in R (R Core Team, 2019) to model body condition (Model 1) and baseline CORT across years (Model 3; Table 1). We used linear mixed effects models using the lme function in package lme4 in R (Bates et al., 2015) to assess CORT methods (RIA versus EIA, Model 2), acute CORT response (Model 4) and bacterial killing ability (Model 5), with animal ID as a random effect to account for repeated sampling. To allow interpretation of baseline and 1 h CORT levels in Models 2, 4 and 5, we included sampling timepoint as a factor (Table 1). We did not include time-tobleed as a covariate in models because (i) these data were only available for 2018 samples and (ii) we did not observe a relationship between time-to-bleed and baseline CORT for 2018 samples (linear, quadratic, cubic all P > 0.2). To assess the differences between groups, we used Tukey's post hoc tests. We transformed CORT values via natural log-transformation after adding 1 (i.e. log1p transformation) and used angular transformation on 50% bacterial killing dilution. We checked data to meet assumptions of normality and homoscedasticity for respective analyses. We ran identical Models 3-5 for each 2018 CORT method (EIA and RIA) and report results of both to compare differences in interpretation related to different locations or hormone measurement methods.
The EIA model showed similar results except that baseline CORT between 2003 and 2018 are now statistically distinguishable (2.57 ng/mL greater in 2018; P = 0.0002). EIA models estimated a smaller difference between 2018 and both 1992 and 1993 CORT than RIA models (average 2.59-2.6 ng/mL less). Standard errors for RIA and EIA model estimates were comparable (i.e. within 0.05 ng/ml of one another).
The EIA model showed similar results, with slight differences in mean differences calculated between groups. EIA  CORT in 2018, 1.52-ng/mL greater difference between each year's baselines and 0.12-ng/mL greater difference between each year's 1-h acute confinement samples. Standard errors for RIA and EIA model estimates were comparable (i.e. within 0.05 ng/ml of one another).
The EIA model showed similar results, slightly differing in the magnitude of the estimated relationship of CORT and bacterial killing ability, with an additional increase of 0.02 ng/mL for every unit decrease in killing dilution.
The y-axis is inverted to illustrate that samples with a lower concentration of plasma necessary to achieve 50% killing of E. coli represent better bacterial killing ability. Dashed lines connect baseline and 1 h stressed samples for each individual. Points without lines indicate individuals for which the corresponding sample did not have enough plasma for this assay. Although sampling timepoint was not a significant predictor, many lines show a positive trend for lower dilution samples from 1 h stressed individuals.

Figure 5:
Bacterial killing ability (the arcsine transformed dilution of plasma required to kill 50% of bacteria) was positively related to log total plasma CORT levels in B. irregularis (measured by RIA) sampled in 2018, such that higher CORT was associated with better bacterial killing ability. Blood samples were drawn immediately after capture (baseline, dark shading) and after snakes were subjected to acute confinement in a cloth bag for 1 h (light shading). Letters indicate significantly different groups.

Discussion
To maximize transparency and to avoid conflation of hypothesis generation with conclusion, we have opted to discuss our results in two contexts. First, we briefly discuss our results within the context of data collected herein (i.e. how our metrics were related to one another). Finally, we discuss our results in a broader ecological context to generate hypotheses for continued investigation of the intersections between physiology and demography in invasive species, using BTS as a focal species. Our work represents a long-term effort compiled from multiple shorter studies in a natural setting; thus, we were unable to directly control population density, resource availability or management pressure post hoc. This is the reality of many conservation and management studies: scientists and managers lack a crystal ball to predict data needed to interpret the past and the scope to direct the collection of future data with consistent methodologies over time. Nonetheless, both must make use of the imperfect data that do exist. Thus, hypothesis generation is an important and useful exercise in furthering the field of conservation physiology and especially in clarifying the utility of nebulous but commonly used metrics such as CORT.

Findings at face value
There were negligible differences between EIA and RIA measurement methods. Some studies assessing differences between the two methods have observed much higher levels via EIA compared with RIA (Nizeyi et al., 2011;Szeto et al., 2011;Yadav et al., 2013). The slightly higher levels in baseline CORT measured by EIA in this study are likely due to increased sensitivity of EIA relative to RIA, rather than differences in laboratory locations (Fanson et al., 2016). Indeed, recent studies assessing the two methods have found differences only at low concentrations (Burraco et al., 2015;Wheeler et al., 2017). This slight difference did not change overall interpretation of trends in our data regarding the 2018 samples. Given this result, long-term monitoring of BTS plasma CORT is possible using either RIA or EIA and may be possible for other species as well.
We observed a negative association between baseline CORT and body condition in the later years. This relationship has been reported many times before in colubrid snakes (Moore et al., 2000;Waye and Mason, 2008;Lind et al., 2018; but see Dayger et al., 2013) and other reptiles (reviewed in Moore and Jessop, 2003). Our data challenge the general usefulness of body condition as a correlate of baseline CORT levels in small-bodied colubrid snakes. The differences between years (i.e. a positive relationship in 1992), highlights that context is important to interpret these relationships, and a negative relationship between baseline CORT and body condition in snakes is not a hard-and-fast rule (discussed in Sandfoss et al., 2020).
We observed increased baseline CORT in 1992-93 relative to 2003 and 2018 across all sexes. We speculate on proximate causes of this observation below. The response of CORT to acute confinement differed between 2003 and 2018. A 2000 study assessing CORT in free-ranging BTS on Guam versus 2-h confinement more closely matches CORT results from 2018; however, the snakes subjected to acute confinement in the 2000 study had spent one night in a trap, so CORT levels could be confounded by increased CORT while in the trap (Mathies et al., 2001). Thus, whether lower 1 h CORT in 2003 represents a dampened acute response or higher 1 h CORT in 2018 represents a more sensitive response is difficult to interpret. Although body condition was not different between 2003 and 2018, previous work demonstrates a link between dampened 1 h CORT and low body condition in a viperid snake (Sandfoss et al., 2020). The observed difference in response to acute confinement between years may be partially explained by other interacting processes (discussed below).
A weak positive correlation between plasma bactericidal capacity and CORT levels was observed, but it is unknown whether the magnitude of the correlation is biologically significant. Interestingly, the 1 h and baseline samples were not statistically different, although CORT was higher at 1 h. This may be an artefact of small sample size, as the direction of the relationship between CORT and bactericidal capacity within individuals was relatively consistent (Fig. 5). Alternatively, while CORT can affect plasma bactericidal ability in as little as 10 min after handling in house sparrows (Gao et al., 2017), CORT can also take longer to ultimately exert downstream effects (Buttgereit and Scheffold, 2002). If this is the case in our study, the CORT level at 1 h may not necessarily reflect its effect on immune function in the same sample. In other words, the immune function at 1 h may be influenced by baseline levels of CORT and acute confinement would affect immune function at a future timepoint, which we did not examine. Regardless, the weak positive correlation of bactericidal capacity and CORT is interesting; it lends evidence towards CORT (baseline and 1 h) as potentially preparative for immunity as proposed by Sapolsky et al. (2000), rather than suppressive, as is assumed in many studies (reviewed in Vera et al., 2017). The relationship between CORT and immune function is context-dependent, may be species-specific (McCormick et al., 2015) or sometimes may be a misinterpretation. For example, CORT and immunity may be correlated but the variables may not influence each other; CORT and immune function may simply be responding at the same time (i.e. upon capture) to different pathways that are activated by capture stress (reviewed in Vera et al., 2017). There is no previous immune data for BTS, so we are unable to speculate if this relationship is characteristic of the species, if this is a plastic or adaptive response to current or past events or if the CORT-immune correlation is biologically relevant; this may be resolved with continued study.

Findings in an ecological context
Regarding CORT and body condition, our findings lead to three hypotheses to explain the patterns, which are not mutually exclusive and which require further experimental testing. These processes affect one another and carefully designed or timed experiments (e.g. in conjunction with future and ongoing management efforts) may elucidate these relationships. We discuss the potential contributions of each process to our observed data given what has been reported for other vertebrate populations.

Density of snakes
In Guam, the density of BTS has decreased from surveys just prior to sampling in 1992-93 , to the most recent survey in 2016  Fig. 1). Decreasing density estimates in animal populations can result from human intervention (e.g. culling as management effort; Gordon et al., 2004), decreased resource availability (Salamolard et al., 2000), disease processes (reviewed in de Castro and Bolker, 2005) and/or variable population estimation methods (reviewed in Freckleton et al., 2006). We will discuss potential influences of management and resource availability on our results separately. While a discussion of the comparability of data from different population density estimation methods is beyond the scope of this article, it is important to note that there is some margin of error in density estimate comparisons. Because there is no density estimate available for the middle of our sampling period, we restrict our discussion to general differences observed between the initial (1991-92) and two latter periods (2003 and  Changes in density can affect interspecific interactions (e.g. Wal et al., 2014), which may influence transmission of infectious disease (reviewed in Lloyd-Smith et al., 2005), and competition for mates (Jirotkul, 1999). There is currently little evidence to suggest the influence of disease processes on population density for BTS on Guam. Exotic cestode parasites have been documented in BTS from Southern Guam, where infected snakes had higher body condition; parasites are not fatal and are likely transmitted by consumption of an intermediate host so may not be related to density (Holldorf et al., 2015). We did not observe the parasites in snakes in our study that were collected outside the prevalence area (Holldorf et al., 2015); thus, it is not relevant to body condition or CORT differences in our study. To our knowledge, fatal disease has only been documented in BTS that were captive for at least 1 year (Nichols et al., 1999). Snakes included in our study were not noted as diseased, excepting one individual collected in 2018 excluded from analysis above. This animal had skin lesions, but histological analysis was determined by UF Veterinary Diagnostic Laboratories to be related to ecdysis issues rather than viral, bacterial or fungal infection (Case no. A18-0343-R). Overall, disease processes related to population density in BTS are unlikely to explain the pattern in our data, but other interspecific interactions are implicated.
Individuals in high density populations may interact with conspecifics more often, which may lead to increased aggression (Knell, 2009). Aldridge et al. (2010) suggested that elevated CORT levels in male snakes observed in 1992-93 as discussed in Moore et al. (2005) may be a result of increased combat at higher densities. While male BTS compete for mates via combat in captivity (Greene and Mason, 2000), and presumably in the wild, it has not yet been documented in the native or introduced range. In some reptiles, losers in combat show increased CORT levels 1-h post fight (Schuett et al., 1996;Schuett and Grober, 2000) or maintain increased CORT levels up to 30 days (Greenberg et al., 1984), while in one lizard species, CORT is increased in winners for an unknown duration (Baird et al., 2014). It is unknown how long CORT is altered after combat in snakes or how BTS respond hormonally to combat. If 'loser' males exhibit increased CORT in BTS, this is a possible partial explanation of our observed higher CORT at higher density. There appears to be few receptive or gravid females at any given time (Rodda et al., 1999;Savidge et al., 2007), even in the native range (Whittier and Limpus, 1996;Trembath and Fearn, 2008). It is possible that at higher densities, more males compete for the same female, leading to more 'losers'. In addition to combat, increased conspecific availability in high-density populations may increase stimulation of mate-searching behaviour (Greene et al., 2001;Mathies et al., 2013). Some male snakes exhibit lower body condition in breeding season because breeding precedes feeding (O'Donnell et al., 2004). The authors are unaware of this phenomenon in BTS but concede the possibility of lower body condition in male BTS pursuing females. While there is some evidence for density-dependent interactions influencing our findings for males, increased male combat and mate-searching do not explain why we observed similar CORT and body condition results in female snakes.
Both receptive and possibly non-receptive female BTS may be harassed by male BTS attempting to breed (Mathies et al., 2013). At high densities, increased interactions with males appears to cause physiological stress via increased lactate in female garter snakes (Shine et al., 2003), but we are unaware of similar studies quantifying CORT in female snakes. Even when harassment by males is severe (e.g. in female garter snakes emerging from a communal den site) females are only delayed from feeding by a few days (Shine et al., 2000). We thus doubt that female body condition is affected by densitydependent male mating attempts, assuming food availability is constant. In female snakes, anorexia is typically associated with gravidity (Gregory et al., 1999), which, as noted above, is rarely observed and unlikely to affect our body condition or CORT results. As we did not observe differences between the sexes in body condition or CORT, it is difficult to reconcile the possibility of density-dependent increased interactions as wholly explanatory of our data. However, until now, we have ignored another key density-dependent effect. Increased density certainly puts pressure on limited food resources, which may aid in explaining our observations.

Food availability
Our understanding of fluctuating prey availability on Guam is based on a combination of surveys that studied prey species presence and BTS diets (see Fig. 1 timeline). By the time of our study, endothermic prey was likely consistently at low levels. Local avifauna were mostly extirpated by 1986 due to depredation by BTS (Savidge, 1987). At the time of our study, native birds were not a major component of BTS diets, although some consumption of domestic chicken chicks and eggs along with introduced avian species is documented (Savidge, 1988;Siers et al., 2017a). Introduced small mammals, consumed by large BTS (Savidge, 1988), also decreased relative to estimates before our study period (Wiewel et al., 2009). The species composition and population levels of ectothermic prey have been in flux. Ectothermic preys are common in the diet of BTS (Savidge, 1988;Siers, 2015) and comprise a majority of the diet in small BTS, which prefer geckos to neonate rodents (Lardner et al., 2009). About the time BTS were introduced, a skink (Carlia ailanpalai) also established, becoming common around 1968 (Austin et al., 2011), and a house gecko (Hemidactylus frenatus) was already established . In the early '90s, many native lizard species were noted as scarce (pers comm in Savidge, 1991;Campbell et al., 2012). Carlia ailanpalai and H. frenatus increased by the late '90s, apparently in response to a decline in nonnative shrews (Suncus murinus; Fritts and Rodda, 1998), but perhaps also to a synchronous decline in BTS population Campbell et al., 2012). It is difficult to provide concrete estimates on how lizard populations have since changed, due to unresolvable confounds in the only long-term study conducted (Rodda et al., 2015). In addition to lizards, BTS consume frogs in their native range (Shine, 1991) and on Guam (Christy et al., 2007a;Cook, 2012;Mathies et al., 2012;Siers, 2015); although they were uncommon in the diet overall by 2012 (Siers, 2015). Five frog species were introduced to Guam in the 2000s (reviewed in Christy et al., 2007a,b), two of which were established and common by 2012 (Olson et al., 2012;Wostl et al., 2016). To our knowledge, there have been no published diet contents of BTS since Siers' (2015) investigation and no BTS prey preference experiments using frogs. Cannibalism is also documented in BTS (Engeman et al., 1996) but may be uncommon. Without concrete estimates of ectothermic prey populations, it is difficult to interpret our data in terms of overall food availability. However, knowledge of BTS's cosmopolitan diet combined with other studies relating food availability, body condition and CORT responses lead us to speculate on the effects of relative prey abundance.
Based on the information chronicled above, our samples were collected when endothermic prey were likely consistent across years, while ectothermic prey were relatively scarce in 1992-93, increased by 2003 and further increased in 2018 with introduced frogs. Fluctuations in availability of ectothermic prey are likely to influence the juvenile size class (Savidge, 1988(Savidge, , 1991. Although we sampled adult snakes, their previous experience as juveniles may inform their survival, reproduction and physiology as adults (Marcil-Ferland et al., 2013;Holden et al., 2019). Interestingly, Rodda et al. (1999) found that lizard (but not frog) abundances corresponded to BTS abundances of all size classes. It is important to note that Rodda et al. (1999) occurred before multiple frog introductions (Fig. 1). Our data show lower body condition in 1992-93, which coincides with scarcity of ectothermic prey. Food availability was previously implicated as a proximate cause of decreased body condition (Moore et al., 2005;Waye and Mason, 2008). Lack of food, i.e. starvation, can lead to low body condition as animals metabolize fat, carbohydrates and muscle to sustain life (McCue, 2007(McCue, , 2010. Movement in search of prey can accelerate these processes (Higginson and Ruxton, 2015). Experimental removal of rodents on open plots led to greater BTS movements and increased BTS activity outside of areas where rodents were removed, indicating BTS alter food-searching behaviour when prey is scarce (Christy et al., 2017, but see Rodda et al., 2008. As more nonnative lizard prey became available by 2003, BTS body condition increased and was indistinguishable in 2018 when presumably additional anuran prey was available. Some posit that larger, endothermic prey were not available to sustain larger snakes, leading to abundant smaller BTS that fed on abundant lizards (Rodda et al., 1999). Indeed, in areas where endothermic prey is available (i.e. urban habitat), BTS tend to be larger and in better condition (Savidge, 1991;Siers et al., 2017a). However, our sampling was not conducted in urban areas, and we observed body condition differences across years so endothermic prey is less likely to explain better body conditions in 2003 and 2018. Several other studies investigat-ing resource availability in snakes find lower body condition with restricted prey availability (Beaupre, 2008;Sandfoss et al., 2018). The relationship between food availability and CORT is more complex.
Glucocorticoids such as CORT are first and foremost metabolic hormones that mediate glucose availability and protein catabolism when energy demands change (Jacob and Oommen, 1992;Jimeno et al., 2018;MacDougall-Shackleton et al., 2019); these actions may be especially apparent during starvation (Dallman et al., 1993;Romero et al., 2010). There is a body of evidence in avian and non-avian reptiles that baseline CORT is elevated when resources are scarce (Sapolsky et al., 2000;Romero et al., 2010;Dickens and Romero, 2013;Sorenson et al., 2017), coinciding with our observation that scarce resources in 1992-93 relate to increased baseline CORT. In other colubrid snakes, baseline CORT is increased in populations with fewer food resources (Palacios et al., 2012) or after food restriction (Holden et al., 2019). CORT reactivity is also affected by resource availability (Jessop et al., 2013). We observed an increased magnitude of 1 h CORT in 2018 relative to 2003, potentially corresponding to greater availability of frogs in 2018. It is possible that mild toxicity of amphibian prey may lead to increased baseline  and 1 h CORT (Mohammadi et al., 2017) in 2018; however, we have no evidence of this type of response in BTS, and given the scarcity of anuran prey reported in gut contents (Siers, 2015), it is less likely that many individuals in 2018 had recently consumed frogs. Thus, resource availability may better explain lower 1 h CORT in 2003. Dampened acute CORT corresponds with reduced resources in many avian and non-avian reptiles (see Dunlap and Wingfield, 1995;Kitaysky et al., 2007;Romero, 2001;Romero and Wikelski, 2001;Romero et al., 2010), including other snakes (Sandfoss et al., 2020, but see Neuman-Lee et al., 2015. The apparent correlations between CORT and resource availability in BTS deserve experimental attention.

Management pressure
Management of BTS has increased in effort and diversity of techniques since the initial sampling period in our study (reviewed in Engeman et al., 2018). Techniques include but are not limited to trapping, visual searches, barriers and baited oral toxicants. Much of the intense removal efforts are implemented in experimentally closed populations (i.e. no emigration to surrounding environment, snakes in our study unaffected) or in targeted areas such as military bases, harbors and airports (e.g. Engeman et al., 1998;Siers et al., 2018;Nafus et al., 2020). Although less severe, multiple agencies have increased targeted removals on Guam since 2004 (reviewed in Engeman et al., 2018). Removing snakes from the system could have obvious direct implications for population density and demography. In other harvested species, removals affect life history strategy of the survivors, which commonly lead to maturation at smaller sizes or earlier ages (reviewed in Kuparinen and Festa-Bianchet, 2017)

Conservation Physiology • Volume 9 2021
Research article bait methods miss small snakes Boyarski et al., 2008), which may similarly produce selective pressure for breeding at smaller sizes (Siers et al., 2017b). Our data may better reflect 'indirect' effects of management efforts, that is, the physiological effects on individuals that are not captured or avoid capture (Hollins et al., 2018). Even in game species, indirect effects of management and harvesting are only very recently being considered (Kuparinen and Festa-Bianchet, 2017;Hollins et al., 2018). With that in mind, we consider the possibility of indirect effects on management of BTS in our data and advocate for its further study in vertebrate invasions.
We postulate that 'missed' snakes may feel some effects of management. After heavy removals in one area, other individuals may immigrate to that area, as evidenced by high gene flow across Guam in BTS (Kierepka et al., 2019). On the other hand, nearby management may affect snakes that have emigrated to surrounding, less-targeted areas. Given this, snakes not in immediate management areas may have some prior experience or exposure to management techniques. Managers may be thought of as predators in the ecological context of invasion management. Increased presence or threat of predators can cause transgenerational behavioural and physiological changes in prey species (Sheriff et al., 2010;Sheriff, 2015). In snowshoe hares, faecal CORT was highest when threat of predation was also the highest (Sheriff et al., 2011). We observe the opposite in BTS-lower baseline plasma CORT in years corresponding to increased management activities. It is likely that BTS consider humans a predator, as defensive behaviour is increased in sites where snakes had previous experience with humans (Spencer et al., 2015); however, human presence is not constant, even during management control activities. The hormonal responses to human presence in BTS are also unknown. Other reptiles have contrasting CORT responses to human activity or presence (reviewed in Injaian et al., 2020), so it is difficult to interpret our results in this context. Clearly, this is an avenue for further experimental research. Body condition may also be affected by predation. If snakes pursue more active avoidance behaviours or are hesitant to feed, they may decrease in body condition, as seen in snowshoe hares (Sheriff et al., 2011) or heavily harvested ungulates (Proaktor et al., 2007). We also observe the opposite in BTS: improved body condition correlated with increased management efforts and exposure to humans (i.e. living in urban areas). Exposure to novel objects in the environment can affect CORT responses in other animals (Dinces et al., 2014;Baugh et al., 2017). Traps or toxic bait deployment materials could be considered novel items to BTS. Traps may additionally influence snakes that are not captured as they may convey information about captured, stressed conspecifics via pheromones or other scent (musk) if not thoroughly sanitized between captures. Unfortunately, snake responses to novelty are poorly understood (Holding et al., 2014;Heiken et al., 2016). In rats, exposure to novelty early in life led to greater magnitude of acute CORT responses as adults (Dinces et al., 2014). Interestingly, in largemouth bass, individuals with lower acute CORT response were more vulnerable to capture (Louison et al., 2017). Another important consideration is how increased human development on Guam influences snake hormonal responses and responses to novelty. It is possible that our observed greater acute CORT response in adult snakes in 2018 relative to 2003 corresponds to increased exposure of 2018 snakes to novel objects as neonates. Further study is needed to seriously evaluate these and other potential indirect effects of management and human activity in BTS and other species.

Putting interactive processes in context
The increases in ectothermic prey from 1990s to mid-2000s coincide with decreases in snake density. These decreases in density may be reactive responses to previous food scarcity (2003 from 1992-93), the result of increased management pressure removing more snakes (from 2003 to 2018) or a combination of the two. Changes in CORT, body condition and the relationship between them over the years may thus be due to increased management pressure, which decreases effective population size, which increases effective food availability that is also subsidized by potential increases in invasive prey, but the interacting processes are not quite so simple. For example, increased management and the resulting increased food availability or decreased density may have opposing effects on CORT and body condition. To disentangle the effects of each, targeted sampling will be necessary. For example, sampling in areas in Southern Guam with longer term frog abundance versus more recent, smaller frog populations may elucidate the effect of prey availability on body condition and CORT expression. Likewise, sampling closed populations or areas with intense, regular management pressure versus wild areas not frequented by management agencies may reveal effects of management pressure on these morphological and physiological metrics. Finally, sampling of high-density, closed populations versus low-density, closed populations with equal management pressure could show how BTS respond physiologically to population density. Immune metrics can be applied to these questions to grow our understanding of immune investment in wild ectotherms. All the above scenarios may be conducted in conjunction with, or adjacent to, ongoing management efforts to help advance how physiological tools are interpreted in conservation ecology.

Conclusion
We showed that different assay methods can be used to evaluate long-term responses in BTS. Additionally, our data add to growing literature supporting the use of CORT and body condition for use in interpreting demographic health. Although intertwined, we find potential relationships between changes in population density, food availability, management pressure and physiological and morphological metrics. For BTS, a species that is regularly sampled for control and research purposes, at minimum the regular collection of morphology information and blood samples before, during and after direct Overall, the field of applied conservation physiology may be improved by the synchronous study of physiology with management actions for species of conservation concern such as the BTS.

Funding
This work was supported by multiple grants and fellowships over multiple years. The 1992-93 data collection was supported by the US Fish and Wildlife Service (USDI 14-16-0009-1577), the Young Investigator Award (IBN-9357245

Data Availability
Data generated during this study are available as a USGS data release (Claunch and Reed, 2021).